Patients and samples

Preterm infants (<32 weeks gestation) from multiple births cared for in the Royal Victoria Infirmary (Newcastle upon Tyne, UK) had stool collected between April 2010 and February 2012. Stool samples (n = 173) were salvaged directly from the nappy from 12 sets of twins and 1 set of triplets (n = 27) from birth until discharge. EBM (n = 18) from 4 mothers was also analysed as EBM samples were only salvaged toward to latter stages of the study and often insufficient volume was acquired. EBM samples were obtained from the remaining liquid in the syringe after continual feed to the neonate through a catheter tube. All samples were frozen at -20°C within 24 hours. This unit has standardised feeding practices, antibiotic and antifungal use. Fluconazole is used as antifungal prophylaxis if <26 weeks/<1kg and in babies developing NEC. Antibiotics, specifically Penicillin and Gentamycin, were administered to all infants following birth and stopped at 48 hours unless proven infection. Clinical information was obtained from patient notes: NEC was confirmed independently by two experienced neonatal clinicians (JEB/NDE): sepsis was defined by positive blood culture, along with antibiotic treatment for a minimum of 5 days and clinical signs suggestive of infection.

DNA Extraction

Nucleic acid extraction of stool was carried out on 100 mg of sample using the PowerLyzer™ PowerSoil® DNA Isolation Kit (MoBio, CA, USA). Nucleic acid extraction of EBM was performed on 1.8 mL of EBM using the Powerfood™ Microbial DNA Isolation Kit (MoBio, CA, USA). Both kits were used in accordance with the manufacturer’s instructions.

PCR-DGGE

PCR amplification targeted the V3 region within the bacterial 16S rRNA gene and was performed with universal primers for total bacterial amplification; V3F-GC (5'-CGCCCGCCGCGCGCGGCGGGCGGGGCGGGGGCACGGGGGGCCTACGGGAGGCAGCAG-3') and V3R (5'-ATTACCGCGGCTGCTGG-3') as described by Muyzer et al. [11], according to the methods of Baxter and Cummings [12].

DGGE (denaturing gradient gel electrophoresis) analysis was performed using the D-Code DGGE system (Bio-Rad, CA, USA) and selected DGGE bands (n=17) underwent excision and sequencing as previously described [13] (Table S1). A ladder of amplicons derived from pure culture organisms was loaded on each DGGE gel for multiple gel alignment. The analysis of the DGGE banding patterns employed TotalLab Phoretix 1D and 1D Pro software (Nonlinear Dynamics, Newcastle upon Tyne, UK). Sequencing was carried out on an ABI 3730XL sequencer (Eurofins MWG Operon, London, UK) and BLASTn was used to determine the closest sequence homologies. DGGE sequence data was deposited to the National Center for Biotechnology Information GenBank under the accession numbers KC686617-KC686622. Analysis of sequence data was carried out at the genus level where possible. As the sequences generated were based on a ~200 bp fragment of the 16S rRNA gene, robust species level classification is not feasible.

454 Pyrosequencing

A subset of key samples (Stool = 36EBM = 5) was selected for metagenomic sequencing: primarily EBM and stool from the triplets (patients 145/147/148) and a discordant twin pair where one developed NEC (patients 139/140). Although other twin sets in the study were discordant for NEC and/or sepsis, patients 139/140 provided the most robust sampling pre, during, and post disease diagnosis. 454 pyrosequencing was carried out using the bifidobacteria-optimised primer set 357F (5'-CTCCTACGGGAGGCAGCAGAN-3') and 926Rb (5'-CCGTCAATTYMTTTRAGT-3') [14]. The DNA sequencing was performed using the 454 Genome Sequencer FLX Titanium platform (Roche, IN, USA) by the Research and Testing Laboratory (RTL, TX, USA) using previously described methods [15].

The raw sequencing reads were quality filtered in QIIME 1.6.0 [16] using the split-library.py script with the following criteria: 1) exact matches to bar code tags; 2) no ambiguous bases; 3) maximum of 5 mismatches to primer sequence; 4) read-lengths between 200–700 base pairs (bp); 5) average quality score of >25 in a sliding window of 50 bp. Remaining high quality sequences were clustered into operational taxonomic units (OTUs) at 97% similarity using UCLUST [17]. Representative sequences for each OTU were aligned using PyNAST [18] and taxonomic identities were assigned using RDP-classifier (version 2.2) [19] with 50% as confidence value threshold. Detection of potentially chimeric sequences was performed using ChimeraSlayer [20] and chimeric sequences were removed from downstream analysis prior to tree building using FastTree [21]. Sequences were deposited in MG-RAST under the accession numbers 4516545.3-4516585.3

Statistical analysis

The DGGE profiles obtained from both stool and EBM samples were analysed using multivariate partial least squares discriminant analysis (PLS-DA) using SIMCA 13.0 (Umetrics, Stockholm, Sweden) [22]. To check that data was adhering to multivariate normalities, Hotelling’s T² tolerance limits were calculated and set at 0.95. Shannon diversity indices (H'), a diversity index giving a measure of both species richness and the evenness of their abundance, was calculated using the PAST package [23].

The sequence reads generated using 454 pyrosequencing were analysed using weighted UniFrac [24] and visualised using principal coordinate analysis (PCoA). Rarefaction curves were produced to 5000 sequences for all samples (Figure S1). Martin’s P-test of significance [25] was applied to determine if significant differences were occurring in the profiles between and within patients. For all statistical analysis, a probability (P) value of <0.05 was deemed to be statistically significant.

Ethics Statement

Initial collection was part of routine service and all samples were collected during the course of normal treatment. Ethical approval was obtained from the County Durham and Tees Valley Research Ethics Committee to include molecular techniques in August 2010. For all infant stool samples, prospective parental informed consent was documented at the point of donation to have the samples stored for research purpose. For EBM, informed consent was documented at the point of donation to have the samples stored for research purposes from March 2011.

Patients

27 patients (12 twin pairs 1 triplet set) contributed to the study: five developed definite NEC and five developed sepsis (two of whom also had NEC). Demographic information is summarised in Table 1 and a detailed version based on single samples can be found in the supplementary information (Table S2). No infant received supplemental probiotics.

Stool profiles

PLS-DA of all samples from all patients based on DGGE data demonstrated twin pairs had comparable profiles which were distinct from unrelated individuals (Figure 1A). Samples which fell outside the ellipse (indicating Hotellings T² range, at 95% confidence) belonged to the triplets or were associated with late onset infection. Specifically, this comprised sample 139t collected one day post NEC diagnosis as well as samples 39c and 41b collected three days prior and on the day of sepsis diagnosis, respectively. Based on DGGE analysis of the cohort, a subset of longitudinal samples from the triplets (represented by red squares) and twin pair discordant for NEC (represented by orange circles) which showed significant variation in community development were selected for pyrosequencing. The pyrosequencing data was analysed at the genus level using weighted UniFrac. In accordance with the DGGE data, samples grouped with their related twin, showing high intra-sample similarities in the development of the gut microbiome (Figure 1B). This is reflected in the bar plots which show each set of twins and the triplets developed a distinct gut microbiome (Figure S2). Proteobacteria and Firmicutes dominated samples in both DGGE and pyrosequencing. From pyrosequencing, 5 OTUs were found in the core microbiome (present in over 85% of samples) in stools from the families Enterobacteriaceae, Enterococcaceae, and Staphylococcaceae.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 1. Community profiles of gut (stool) microbiome from preterm multiples.

Subjects are symbolised based on related multiples. a) Partial least squares discriminant analysis (PLS-DA) based on DGGE data of all stool samples. The ellipse indicates Hotellings T² range, at 95% confidence. Selected sample labels removed for clarity. b) Weighted UniFrac PCoA based on pyrosequencing data at genus level, generated from a subset of stool samples.

https://doi.org/10.1371/journal.pone.0073465.g001

Comparison of breast milk with respective stool samples

The DGGE profiles of all EBM samples were compared with the respective stool sample. EBM samples showed a relative lack of diversity making them cluster near the origin: despite this PLS-DA did reveal EBM samples clustered with the stool samples of the respective set of multiples (Figure 2A). There are cases where different EBM samples from the same mother cluster separately showing the EBM microbiome was not stable (e.g. BM148e / BM147b / BM145a / BM145i or BM139b / BM139f / BM139q).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 2. Comparison of breast milk with respective stool profiles.

Subjects are symbolised based on related multiples. a) Partial least squares discriminant analysis (PLS-DA) based on DGGE data of all breast milk (EBM) samples matched to respective stools. The ellipse indicates Hotellings T² range, at 95% confidence. b) Weighted UniFrac PCoA based on pyrosequencing data at genus level, generated from a subset of EBM and respective stool samples.

https://doi.org/10.1371/journal.pone.0073465.g002

EBM from the triplet set and twin set (patients 139/140) also underwent pyrosequencing. EBM samples from each mother showed high intra-sample similarity (Figure 2B). For twins 139/140 the EBM profiles (BM139/BM140) clustered closely with the respective stool (139s/140s). This clustering was less robust in the triplet set although triplet EBM (BM145/BM147/BM148) was still more similar to triplet stool (145a/147b/148b) than stool of other infants. The connection of EBM samples with stool can be visualised in the bar plot (Figure S2). Three OTUs were found in the core microbiome in EBM and, like in the stool core microbiome, were from families Enterobacteriaceae, Enterococcaceae, and Staphylococcaceae.

Bacterial community development in twin set 139/140

A reduction in diversity was observed in patient 139 at least 5 days before NEC was apparent clinically (Figure 3A), which was not shared by the twin. A less pronounced reduction in diversity did occur in the control twin (140) earlier in development (day 18) that coincided with antibiotic administration and diversity was quickly re-established when antibiotics ceased. A similar trend was observed when the H' of the pyrosequencing data was performed (data not shown).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 3. Development of gut microbiome in twin pair 139/140 mapped to life events.

P – Penicillin, G – Gentamicin, F- Fluctoxacillin, A - Amoxicillin, M – Metronidazole, ^s –Start of antibiotics, ^e – End of antibiotics, 72hr – full enteral feed (at least 150 ml/kg/day) sustained for 72 hours. a) Shannon Diversity indices (H') of twin pair based on DGGE data b) Turnover of the most prevalent bacterial OTUs throughout the first 36 days of life in twin 139 where antibiotics were prescribed for NEC. c) Turnover of most prevalent bacterial OTUs throughout the first 34 days of life in twin 140 where antibiotics were prescribed due to pyrexia (fever).

https://doi.org/10.1371/journal.pone.0073465.g003

Pyrosequencing of matched twin samples (full profiles in the bar plot; Figure S2) allowed the 12 most abundant OTUs to be incorporated into the analysis (Figure 3). These results are in accordance with the DGGE data and allow greater detail about specific contributors to the overall diversity to be observed. While community structure in the twins was initially comparable, twin 139 showed reducing diversity with an increasing abundance of Escherichia sp., before NEC diagnosis on day 28 (Figure 3B). Conversely, in the sibling there is an increase in diversity and Escherichia sp. was present in much lower abundance (Figure 3C). After antibiotic receipt (day 16 in twin 140 and day 28 in twin 139) both twins demonstrate reduced Escherichia sp. abundance and an increased abundance of other Enterobacteriaceae, rapidly reversing in twin 140 on antibiotic cessation. In twin 139 NEC and subsequent antibiotic treatment significantly (P=0.028) altered the bacterial community in comparison to its sibling, with a new notable bloom in Klebsiella sp. and a smaller increase in Parascardovia sp. (family – Bifidobacteriaceae).