Identical bacterial populations colonize premature infant gut, skin, and oral microbiomes and exhibit different in situ growth rates

Community profile

Two premature infants were recruited for this study with parental consent at Magee-Womens Hospital of the University of Pittsburgh Medical Center. Clinical information for these infants is summarized in Supplemental Table S1. Both infants were born via vaginal delivery to women with pregnancies complicated by chorioamnionitis. Each infant was treated immediately after birth with 7 d of initial antibiotic treatment (ampicillin and gentamycin). We collected 17 and 20 fecal samples during the first month of life for the two infants (referred to as #1 and #2). Additionally, two skin swabs and two tongue swabs were collected for each infant (see Fig. 1 for detailed information about timing of skin and oral swabs). For the skin and mouth, sample choice was based on the availability of sufficient DNA and, where more than two samples were available, to span the longest time period. In total, 183 gigabase-pairs of Illumina shotgun DNA sequencing were generated for 45 samples (Supplemental Table S2).

Historically, assembly-based metagenomics of the skin and mouth has been hampered by human DNA contamination (Liu et al. 2012; Tsai et al. 2016). While all but one skin or mouth sample consisted of >50% human DNA (range 34.2%–93.8%), the deep sequencing effort ensured that all samples had enough reads from microbial genomes for successful de novo genome reconstruction. Overall, an average of 98.7% (Infant 1) and 95.0% (Infant 2) of nonhuman reads could be assigned to assembled genomes (in some cases, by mapping to genomes reconstructed from another sample) (Supplemental Table S3).

Using ggKbase (Raveh-Sadka et al. 2015), we manually binned assembled DNA sequences to genome bins based on G+C content, coverage, and phylogenetic profile and subsequently used sequencing coverage patterns for binned scaffolds to verify the bins (Raveh-Sadka et al. 2015). For the bin refinement step, the clustering of fragments was analyzed using emergent self-organizing maps (ESOMs). Since identical genomes were assembled from different samples, the genomes were de-replicated based on similar average nucleotide identity (ANI), and the best genome for each strain was selected for downstream analyses. Across all three body sites for Infant 1, we recovered nine near-complete bacterial genomes and three partial genomes (all partial genomes were from the skin data sets). From Infant 2, we recovered 11 near-complete bacterial genomes and three partial genomes (all partial genomes were from the mouth). The community composition of all samples from all body sites and both infants are presented in Figure 1. It is interesting that the first two gut samples for Infant 2 that were collected during the antibiotic treatment course were dominated by Streptococcus agalactiae. The placenta showed heavy growth of S. agalactiae (group B streptococcus), yet, as noted above, the infant's blood cultures were negative.

The overlap in community composition between body sites is shown in Figure 2. Strains were considered identical if they had over 99.9% whole-genome ANI and were considered colonizers of a body site if they accounted for >1% of reads in any sample from the site (Supplemental Table S4). Despite obvious differences in habitat characteristics, we identified some identical strain populations in all three sampled body sites for both infants. Infant 1 body sites were heavily colonized by Citrobacter koseri, which comprised over 60% of all three communities. Six strains were colonists of more than one body site of Infant 2: Escherichia coli, Pseudomonas aeruginosa, Klebsiella pneumoniae, and Serratia marcescens colonized all three body sites, Enterococcus faecalis colonized the mouth and skin, and Staphylococcus epidermidis colonized the mouth and gut. Interestingly E. coli, which is traditionally thought of as a gut colonist (Tenaillon et al. 2010), accounted for the highest portion of Infant 2 reads at all three body sites. The infants were housed in different NICU rooms ∼3 mo apart, and no bacterial strains were shared between infants.

We also tested for the presence of organisms in multiple body sites at low abundance (<1% of the community) and, to the extent possible, evaluated the time periods in which specific strains appeared at these sites (for details, see Methods). For Infant 1, Enterococcus faecalis had appeared in the skin and gut by the time of collection of the first samples, persisted in the gut, and was present in both habitats at the later time point when samples from both body sites were collected (day of life [DOL] 23). A Haemophilus parainfluenzae population was not detected in the first collected samples from all body sites but was present in the gut on DOL 21 and had appeared in the mouth and skin 2 d later. For Infant 2, Staphylococcus epidermidis and E. coli were present in all body sites at most time points. Pseudomonas aeruginosa and Klebsiella pneumoniae were undetectable in gut samples collected during antibiotic administration (Fig. 1; Supplemental Table S3) but appeared in all three body sites the day after cessation of antibiotics (DOL 7). Staphylococcus sp. M0480 was present by the time the first skin (DOL 12) and mouth (DOL 14) samples were collected but was undetectable in the mouth at the second time point. E. faecalis was present in the first-sampled gut and mouth communities (and persisted there), was absent in the skin on DOL 12, but had appeared there by DOL 22. Serratia marcescens was a relatively late colonist, appearing first in the gut on DOL 19 and was present in all three body sites at later time points (Supplemental Table S3). In general, strains became detectable at all three body sites at around the same day of life, with the exception of E. faecalis in Infant 2, which persisted in the gut and mouth for over a week before being detected in the skin.

Growth rates are different across body sites

Recently, it was shown that accurate growth rates of microbial strains in their natural environment can be determined by measuring the ratio of the coverage of DNA at the origin and terminus of replication (Korem et al. 2015). However, this method requires complete closed circular genomes (which are rarely acquired from metagenomic studies) in order to locate the origin and terminus. We were able to circumvent that requirement by orienting each strain assembly (median number of contigs 67.5, range 12–1460) to a representative isolate genome in order to determine the order and orientation of the contigs. When available, we used multiple isolate genomes to confirm the best assembly, as some genomes in the RefSeq database were found to be incorrectly assembled around rRNA operons (Fig. 3A,B). The peak-to-trough ratio (origin to terminus coverage ratio) was then determined by mapping reads to the oriented assembly (see Methods).

Surprisingly, in all cases for which we could determine growth rates for the same strain from multiple body sites, growth was faster in the skin and mouth than in the gut (Fig. 3). When all growth rates were analyzed, microbes in the skin and mouth had significantly higher growth rates than microbes in the gut (P < 0.00001), whereas strains in the skin and mouth did not differ significantly in their growth rates (P = 0.12; Mann-Whitney U test). The strains exhibiting the fastest growth in the skin, mouth, and gut were P. aeruginosa, Streptococcus mitis, and Clostridium perfringens, respectively.

In general, growth rates measured for strains in the gut increased with increasing infant age, consistent with population recovery after early antibiotic treatment (Spearman rank correlation, R = 0.30, P = 0.0005, n = 135) (Fig. 4). In Infant 2, S. agalactiae, also known as group B Streptococcus, accounted for over 80% of the microbial community during antibiotic treatment in Infant 2 for a presumed group B Streptococcus infection. Despite being abundant, the organism exhibited a low and decreasing growth rate during the treatment period. Several microbial strains exhibited sharp changes in growth rates over the first month of life (Clostridium difficile, C. perfringens, P. aeruginosa). However, these events did not coincide with medical events indicated by the clinical metadata and were probably short lived, as they did not always lead to changes in relative abundance of these strains in the next-collected sample.

Microdiversity

We classified strains colonizing multiple body sites as the same based on >99.9% genome-wide ANI, but this analysis is insensitive to very small-scale differences that could be used to document subpopulation dynamics and constrain inoculum diversity. Thus, we performed high-resolution analyses of Citrobacter koseri, the strain with the highest coverage in all body sites of Infant 1. The 4.66-Mbp genome was initially assembled de novo into 27 scaffolds. By reference to an isolate genome, we confirmed potential joins supported by sequence overlaps and identified gap-filling reads so that these scaffolds could be reconstructed into a complete circular genome.

Single nucleotide polymorphisms were identified by mapping reads from each body site to the circularized C. koseri genome. No fixed mutations distinguished populations colonizing the skin, mouth, and gut. However, 47 polymorphic sites with a minimum frequency of 0.2 were identified along the genome. Of these, 21 occurred in intergenic regions and five within the coding region of yadA, a gene encoding adhesins with known pathogenic alleles (El Tahir and Skurnik 2001). Using hierarchical clustering, 10 cohorts of polymorphisms with similar variations in frequency over the sample series were identified, five of which had only one member (referred to as “singletons”) (Fig. 5). Each cohort is inferred to represent a strain subpopulation. Particularly interesting is cohort 1, which underwent a dramatic purge event around DOL 23, and singleton 3, which rapidly rose above detection level on DOL 17.

To evaluate differences in the populations across body sites, we focused on DOL 23, where gut, skin, and mouth samples are all available for Infant 1. Three variant positions had significantly different fractions of polymorphisms between the gut and both the skin and mouth (Fisher's exact test with Bonferroni correction, α = 0.01). These three positions make up the entirety of cohort 4 and singleton 3. Singleton 3 rises in abundance to comprise ∼40% of the gut population on DOL 17 but is only ever detected at ∼2% of the mouth and skin populations (Fig. 5B). There were no genomic positions with significantly different polymorphism levels between the mouth and skin on either DOL 6 or DOL 23. The detection of differences in subpopulation frequencies between body sites shows that there is limited gene flow between body sites and may indicate the start of in situ diversification.

CRISPR and phage

In addition to the bacterial genomes referenced above, our assembly-based metagenomic pipeline resulted in the recovery of 21 bacteriophage genomes and 18 plasmid genomes. An average of 2.1% of reads from all samples mapped to bacteriophage genomes, and the most significant bacteriophage bloom occurred in the gut of Infant 2 on DOL 7 (immediately following cessation of antibiotic administration) (Fig. 1). Three bacteriophage genomes alone accounted for ∼50% of the DNA sequenced during this bloom: K. pneumoniae phage A (7.9% of community; 12,700× coverage), K. pneumoniae phage B (39.1% of community; 19,400× coverage), and E. coli phage A (2.7% of community; 4100× coverage). A second bloom of K. pneumoniae phage B and E. coli phage A also occurred in the same environment on DOL 24. Interestingly, although the bacteriophages’ abundance in the second bloom was only a fraction of their abundance in the first bloom, both blooms dramatically shifted the microbial community composition (Supplemental Fig. S1). Overall, however, most phage, plasmid, and host population abundance patterns were highly correlated (partly due to integration). Consequently, bacterial distribution patterns across the three body sites generally predicted patterns for the associated phage and plasmids (Supplemental Fig. S2). For circularly recovered phage and plasmids (inferred to be nonintegrated), we also found GC skew and coverage patterns that could be indicative of DNA replication style (Supplemental Fig. S3). Unfortunately, no phage or plasmids with these coverage patterns were detected at multiple time-points, so the consistency of this pattern could not be evaluated. However, such analyses could possibly be used to elucidate plasmid and phage replication regulation and growth rate in future studies.

CRISPR-Cas loci confer bacterial phage resistance (Horvath and Barrangou 2010). Because CRISPR spacers are added uni-directionally, the loci provide a record of population history (Sun et al. 2016). We identified CRISPR arrays in 50% of bacterial genomes in Infant 1 and 43% of bacterial genomes in Infant 2, with a total of 111 and 182 unique spacer sequences, respectively (Supplemental Table S5). However, we found only four spacer targets in the samples from the infant from which the spacers were recovered, even using relaxed (1 mismatch allowed) search parameters. Thus, we broadened our search to include previously published gut metagenomes from the same NICU as this study (Raveh-Sadka et al. 2015, 2016) and the current NCBI database. This revealed 26 and 16 spacers/protospacer matches, respectively.

No changes in CRISPR spacer inventories over the study period of 28 d were found, but two coexisting E. faecalis populations were identified in Infant 2 based on distinct CRISPR spacer inventories. The variant A locus contained 13 spacers and the variant B locus contained 11 unique spacers, with one spacer appearing in the array twice (Fig. 6). The ratio of the two variants fluctuated over the study period, with variant B rising from 0% of the population on DOL 5 to over 60% during the DOL 12–24 period and declining to ∼5% by DOL 28. As the spacers that differentiate variants A and B had no protospacer matches, we sought other genomic features for which selection might explain the abundance changes. One polymorphism was significantly correlated with the variant A locus and another significantly correlated with the variant B locus (Pearson correlation with Bonferroni correction, α = 0.01) (Fig. 7). Both polymorphisms were single-nucleotide substitutions within the coding regions of separate hypothetical proteins. Bacteria with both CRISPR locus variants and the associated gene variants were detected in the skin on DOL 22 (the only sample from another body site with sufficient coverage for detection). CRISPR arrays reconstructed in this study were also compared to a number of previously isolated E. faecalis strains (Fig. 6; Palmer and Gilmore 2010). Surprisingly, we found remarkable similarity in spacer content (with no polymorphisms in the spacer sequences) for E. faecalis strains isolated as many as 82 yr apart. Additionally, we found that isolates from the 1970s still retain CRISPR that match the sequences of phage still found in the NICU today.

Patient recruitment and sample collection

Fecal samples from two preterm infants hospitalized in the NICU in Magee-Womens Hospital of UPMC (Pittsburgh, PA) were collected as available over the first month of life. Both infants were of low gestational age (<30 wk), and Infant 2 was of extremely low birth weight (<1000 g). See Supplemental Table S1 for additional clinical information.

Fecal samples were spontaneously expelled and collected from diapers or acquired directly using an established perineal stimulation procedure (Morowitz et al. 2011). Skin and oral samples were obtained by a member of the study team using a BD BBL Culture Swab EZ. Oral swabs were collected by rolling the swab head 5–10 times over the dorsal surface of the tongue. If intubated, the sample was collected by swabbing any exposed surface of the tongue. Skin swabs were collected by first dipping the swab into a 0.5 mL aliquot of a sterile solution of 0.15 M NaCl and 0.1% Tween 20. The swab head was then rolled 5–10 times over the left anterior upper chest wall. Stool samples were placed promptly into −20°C storage and transferred to a −80°C freezer for long-term storage as soon as possible. Swab samples were placed promptly in a −80°C freezer for storage.

DNA was extracted using either the MO BIO PowerSoil DNA Isolation kit (single tube extractions) or PowerSoil-htp 96-Well DNA Isolation kit. DNA extracted from stool using the single tube format followed the protocol as previously described (Raveh-Sadka et al. 2016). For DNA extracted from feces with the 96-well kit, fecal samples were added to individual wells of the bead plate and stored overnight at −80°C. The next day, the Bead Solution and Solution C1 were added, and the plates were incubated at 65°C for 10 min. The plates were shaken on a Retsch Oscillating Mill MM400 with 96-well plate adaptors for 10 min at speed 20. The plates were rotated 180° and shaken again for 10 min at speed 20. All remaining steps followed the manufacturer's centrifugation protocol. For swab samples, the swab head was cut off directly into the wells of the bead plate and stored overnight at −80°C. The next day, the Bead Solution and Solution C1 were added, and the plates were incubated at 65°C for 10 min. The plates were shaken on a Retsch Oscillating Mill MM400 with 96-well plate adaptors for 5 min at speed 20. The plates were rotated 180° and shaken again for 5 min at speed 20. The Solution C2 and C3 steps were combined (200 µL of each added) to improve DNA yield. All remaining steps followed the manufacturer's centrifugation protocol.

Metagenomic sequencing and assembly

Sample preparation and sequencing of skin and oral samples were performed at the University of Illinois at Urbana-Champaign sequencing facility, and sample preparation and sequencing of fecal samples were performed at the University of California at Berkeley Vincent J. Coates Genomics Sequencing Laboratory. Paired end reads of 160 bp with a combination of 1000- and 600-bp library insert sizes were sequenced using an Illumina HiSeq 2500 (Supplemental Table S2). Reads were trimmed with Sickle (https://github.com/najoshi/sickle). Reads that mapped to the human genome with Bowtie 2 (Langmead and Salzberg 2012) under default settings were discarded. An additional step of mapping with BBMap (Bushnell 2014) was performed on all projects with at least 10% of reads removed with Bowtie 2 mapping. See Supplemental Table S2 for depth of sequencing and levels of human contamination in each sample.

Reads were assembled using idba_ud (Peng et al. 2012) under default settings. Resulting scaffolds >1 kb in length were annotated using Prodigal (Hyatt et al. 2010) to predict open reading frames using default metagenomic settings. Annotated protein sequences were searched against KEGG (Kanehisa et al. 2014), UniReff100 (Suzek et al. 2007), and UniProt databases using USEARCH (Edgar 2010). All matches with bit scores greater than 60 were saved, and reciprocal best hits with a bit score greater than 300 were also cataloged. We identified rRNA sequences using Infernal (Nawrocki and Eddy 2013) by searching against databases from the SSU-Align package (Nawrocki 2009) and tRNAs using tRNAscan_SE (Lowe and Eddy 1997).

Genome binning was carried out using the online interface within ggKbase as described previously (Raveh-Sadka et al. 2015; http://ggkbase.berkeley.edu/). This method takes into account phylogenetic profile, GC content, and coverage information. Bins were refined based on differential coverage implemented using time-series emergent self-organizing maps as described previously (Sharon et al. 2013). The completeness of bacterial bins was evaluated based on the presence or absence of single-copy genes (Raes et al. 2007; Raveh-Sadka et al. 2015). Phage sequences were identified based on the presence of typical phage genes such as capsid, terminase, and tail-fiber and as distinct clusters in ESOMs; phage-host relationships were inferred based on phylogeny of annotated proteins and abundance patterns.

Genome recovery

All genome bins from all samples were pooled based on (1) the infant the genome was recovered from, and (2) the genome's identity as either of phage/plasmid origin or of bacterial origin. Genomes within each pool were then compared in a pairwise fashion based on ANIm (Richter and Rosselló-Móra 2009). For all clusters of genomes with high ANI values among members, a representative genome was chosen based on the highest total bin length, lowest scaffold fragmentation, and most complete complement of single-copy genes. Ambiguous genome clusters were visualized using Mauve alignments in Geneious (Kearse et al. 2012) to decide whether genomes could be included in the cluster.

Next, for each infant, ANIm was determined for all representative bacterial bins and all representative phage/plasmid bins together. The resulting ANI matrix was manually curated to resolve cases of overlap between the two genome sets. Most cases of overlap were of bacterial genomes containing prophage that were also represented in the phage list. These were resolved by removal of the prophage scaffold from the bacterial bin. The final genome list for each infant was verified by mapping reads from each project to the genomes to confirm strong coherence between the reads and the genomes, as well as to verify the completeness of the list based on total percentage of mapped project reads. Read-mapping data are available in Supplemental Table S3, and the final genome list is available in Supplemental Table S6.

Sample profiling

Reads from all samples were mapped to the corresponding infant's genome list generated above. SAMtools (Li et al. 2009) was used to convert mapping files (.sam) to mpilup format, and calculate_breadth.py and pileup_profile.py were used to determine the depth of coverage, breadth of coverage, and average nucleotide identify of each genome in each project. As SAMtools has an implicit coverage limit of 8,000×, coverage values from calculate_breadth.py were used. The results of both scripts were manually combined and are available in Supplemental Table S7.

A strain was considered a “colonist” of a body site if at least 1% of the reads from at least one sample from at least one body site mapped to the recovered genome. We chose to define colonization in this way to be consistent with previous studies of infant colonization (Ward et al. 2016) and because using a coverage-based threshold would have biased against samples with less sequencing reads (Supplemental Fig. S4). The same analysis was performed on all phage and plasmid sequences, using 99% breadth to define carriage (Supplemental Fig. S2). To identify when strains below the 1% threshold first appeared in body sites, we normalized to account for different sampling depths by using a read percentage cutoff which corresponds to 0.1× coverage of the most shallowly sampled data set (see Jupyter notebook, CallingColonisits, for details).

ANI calculation from metagenomic reads

Metagenomic reads from each sample were mapped to the genome list described previously, and nucleotide variants between reads and genomes were determined using pileup_profile.py. See the data availability section for full source code. Briefly, the script calculates ANI by masking regions of DNA near the ends of scaffolds, in conserved regions (tRNAs and rRNAs), or of insufficient coverage, locating all base pairs along the unmasked genome in which at least 80% of reads conflict with the reference genome, and calculating consensus ANI as (1−[# variant positions/unmasked genome length]). As shown in Supplemental Table S3, the number of SNPs found using this method was extremely low (average 34.2), with an average consensus ANI of 99.998%. We attempted to reduce the number of erroneously called SNPs and found that some appear to represent errors made during the process of metagenomic assembly (generation of the reference sequence) or unmasked regions of high sequence conservation (which recruit reads from other genomes) rather than real biological differences. Given this, and the extremely high reference ANI between strains on different body sites, we defined the strains as identical if they met the criterion of >99.9% consensus ANI.

Growth rate determination

To attain growth rates for the incomplete genomes recovered in this study, genome fragments were ordered and oriented to previously isolated reference genomes and the peak-to-trough coverage ratio was determined using bPTR.py based on the method previously described (Korem et al. 2015; Brown et al. 2016) (Supplemental Table S8). Circular reference genomes of the same species as draft genomes from this study were downloaded from NCBI GenBank. The expected form of the cumulative GC skew of genomes (Grigoriev 1998) was manually verified using the program gc_skew.py, and genomes with aberrant patterns were discarded. The ANI of each draft genome to all reference genomes was determined using the previously described ANIm method, and the reference genome with the highest ANI was chosen. Draft genome fragments were aligned to the reference genome using BLAST (Altschul et al. 1990), and any fragment with <20% alignment coverage was discarded. The remaining draft sequence fragments were then aligned to the reference genome using progressive Mauve (Rissman et al. 2009) (java -Xmx500m -cp Mauve.jar org.gel.mauve.contigs.ContigOrderer), resulting in an ordered and oriented draft “core genome.” All core genomes were verified by manual inspection of the cumulative GC skew and genome coverage plots generated by the script (Supplemental Figs. S5, S6). Projects with aberrant plots or coverage below 5× were excluded from analysis.

Microdiversity of C. koseri

To identify differences between the reads of specific data sets relative to reconstructed genomes (described above), VarScan (Koboldt et al. 2009) was run on all .pileup files using the pileup2cns command with the flag -min-coverage = 3. The frequency of each single nucleotide variant was tabulated for all samples using the script polymorpher2.py (Supplemental Table S9). Variants were then filtered based on a number of criteria and clustered based on changing relative frequency. Briefly, variants were required to have a minimum of 10× coverage in all samples, over 0.2 frequency in at least two samples, not be defined by polymorphisms present in conserved regions (rRNAs, tRNAs) (Supplemental Table S9), and pass an auto-correlation threshold. Clustering of variants was done using the Scipy hierarchical clustering package, with the cutoff threshold of 0.275 decided based on manual inspection of the resulting clusters. Full source code of analysis performed is available in the Jupyter notebook, CitroK_microdiversity (Supplemental Material).

CRISPR analysis

CRISPR arrays were identified in bacterial draft genomes using the program CRISPRFinder (Grissa et al. 2007). CRISPR spacer targets (protospacers) were identified by searching spacers for full-length BLAST hits in a sequence database, followed by filtering our results that also included full-length matches to CRISPR repeats (to remove instances of the CRISPR array itself). In addition to the assemblies of this study, we also searched for protospacers in previously published studies from the same NICU (Raveh-Sadka et al. 2015, 2016) and the NCBI nt database (accessed January 2016). DNA fragments with identifiable CRISPR arrays were excluded from the protospacer search. The alignments of E. faecalis CRISPR arrays were manually curated within Geneious (Kearse et al. 2012). The ratio of variant arrays was determined by comparing the number of reads that mapped (using Bowtie 2 default settings) to unique regions of both arrays. Mutations elsewhere in the genome that varied in frequency with CRISPR variants were identified based on a significant Pearson correlation. Source code is available in the Jupyter notebook, E. faecalis microdiversity (Supplemental Material).