Strain-resolved microbiome sequencing reveals mobile elements that drive bacterial competition on a clinical timescale

Abstract

Background: Populations of closely related microbial strains can be simultaneously present in bacterial communities such as the human gut microbiome. We recently developed a de novo genome assembly approach that uses read cloud sequencing to provide more complete microbial genome drafts, enabling precise differentiation and tracking of strain-level dynamics across metagenomic samples. In this case study, we present a proof-of-concept using read cloud sequencing to describe bacterial strain diversity in the gut microbiome of one hematopoietic cell transplantation patient over a 2-month time course and highlight temporal strain variation of gut microbes during therapy. The treatment was accompanied by diet changes and administration of multiple immunosuppressants and antimicrobials.

Methods: We conducted short-read and read cloud metagenomic sequencing of DNA extracted from four longitudinal stool samples collected during the course of treatment of one hematopoietic cell transplantation (HCT) patient. After applying read cloud metagenomic assembly to discover strain-level sequence variants in these complex microbiome samples, we performed metatranscriptomic analysis to investigate differential expression of antibiotic resistance genes. Finally, we validated predictions from the genomic and metatranscriptomic findings through in vitro antibiotic susceptibility testing and whole genome sequencing of isolates derived from the patient stool samples.

Results: During the 56-day longitudinal time course that was studied, the patient's microbiome was profoundly disrupted and eventually dominated by Bacteroides caccae. Comparative analysis of B. caccae genomes obtained using read cloud sequencing together with metagenomic RNA sequencing allowed us to identify differences in substrain populations over time. Based on this, we predicted that particular mobile element integrations likely resulted in increased antibiotic resistance, which we further supported using in vitro antibiotic susceptibility testing.

Conclusions: We find read cloud assembly to be useful in identifying key structural genomic strain variants within a metagenomic sample. These strains have fluctuating relative abundance over relatively short time periods in human microbiomes. We also find specific structural genomic variations that are associated with increased antibiotic resistance over the course of clinical treatment.

Keywords: Antibiotic resistance; DNA; Gut microbiome; Linked reads; Metagenomics; Read cloud assembly; Sequencing; Strain diversity; Structural variation.

Fig. 1

Patient condition, drug exposure, and intestinal microbiome composition during treatment. The study subject was admitted to Stanford Hospital with myelodysplastic syndrome and myelofibrosis and subsequently underwent hematopoietic cell transplantation (HCT, denoted by a red line). Stool samples were collected prior to HCT and over the following 5 weeks as the patient underwent chemotherapy, antibiotic treatment, and immunosuppression. Taxonomic classification of shotgun metagenomics sequencing reads (Illumina TruSeq Nano DNA) reveals pronounced dysbiosis emerging following HCT with gut domination by Bacteroides caccae, a commensal bacterium. Relative abundances of each species are determined using reads that are classified at the species level. The relative abundance of the species in the samples was determined after removing human and viral reads from the data. For this figure, the top 9 most abundant species in the shotgun metagenomics data are shown. Relative abundance data for all the species-level classification is shown in Additional file 2: Table S1. For a more detailed view of the species-level and genus-level classifications of all the samples, refer to Additional file 1: Fig. S1

Fig. 2

Read cloud and short-read genome drafts of Bacteroides caccae obtained from sequencing of stool samples. Circos plots of draft genome assemblies of B. caccae from the sequencing of the patient stool samples are ordered chronologically (A through D) from the outermost inward. Read cloud genome drafts assembled using Athena (blue) are more contiguous and complete than short-read drafts assembled using conventional assembly (gray). The read cloud B. caccae draft from timepoint C (third circos plot from the outside) was the most contiguous and served as the reference for all alignments. Contigs from each read cloud and short-read library are assigned a lighter color if they aligned to this reference, but did not belong to a draft classified as B. caccae. Short reads did not produce a draft annotated as B. caccae for timepoint A. The best read cloud and short-read drafts are both obtained through sequencing of timepoint C (read cloud: 414 kb N50, 5.5 Mb size, 99.3% complete, 1.6% contaminated; short read: 88 kb N50, 4.7 Mb size, 97.7% complete, 0.6% contaminated). The read cloud drafts include a total of 18 assembled integration sites of IS614 (red circles) and 25 assembled integration sites of a candidate insertion sequence (green circles) that are missing from all short-read drafts. Alignments of raw short reads to these sites indicated the presence of both strains without insertion sequence integration and strains with the insertion sequence integration. Estimated proportions of these strains for each site and timepoint are shown with different filled in areas of each circle, with an empty circle denoting predominance of ancestral strains lacking an IS at that location and shaded circles denoting predominance of strains with the IS at that location

Fig. 3

Co-occurrence of multiple B. caccae strains with differing IS614 integrations. a Alignments of short reads from timepoints B, C, and D to a representative IS integration site reveal domination of the strain without the IS (“ancestral strain”) in timepoints in B and C, and domination of the strain harboring the insertion in D. Short-read alignments from B and C show many reads spanning over both left and right junctions (red), indicating global alignment to the ancestral sequence, while short-read alignments from D show many reads supporting the IS integration (blue), indicating read pairs or single reads spanning the IS. This demonstrates that the IS is present at this locus at timepoint D but is undetectable at timepoints B and C. b Estimated relative abundances of B. caccae ancestral strains and strains with an IS integration for all 18 detected IS614 integration sites. Integrations upstream of annotated antibiotic resistance genes norM, thyA2, and per1 are shaded. Major shifts in abundances amongst the strains with and without integrations upstream of norM and thyA2 can be seen between timepoints B and C and timepoints C and D, respectively. Integration sites are sorted by ancestral strain fraction in timepoint C

Fig. 4

Metatranscriptomics support IS-mediated transcription within B. caccae. IS614 contains a putative outward-facing promoter. The relative contribution of the IS promoter to transcription was determined by comparing RNA sequencing read depths of genes upstream and downstream of it. a In timepoint B, which is dominated by ancestral strains without the promoter, RNA sequencing read coverage depth (relative transcript abundance) is relatively equal on both sides of the integration site. In timepoint C, which is dominated by strains with IS614 with its putative outwardly directed promoter, the transcription of the downstream gene norM is much higher than that of the upstream gene yidC. The relative transcript abundance of all neighboring genes increase in timepoint C relative to B, but this increase is 10-fold greater in genes immediately downstream of the introduced outward promoter. In later timepoints C and D, dominant strains harbor an introduced IS promoter positioned to upregulate norM. This is supported by read pairs spanning between the IS promoter and norM. This difference in coverage and domination by strains with this promoter both persist through timepoint D. Conversely, the earlier timepoint B is dominated by strains with no IS in this region. b PCR with primers flanking the above integration instance of IS614 yields amplicons without the insertion sequence in earlier timepoints A and B (400 bp), and with the insertion in later timepoints C and D (1.9 kb)

Fig. 5

Genomic alterations detected in B. caccae isolates support selection throughout treatment. The detected proportions of IS614 and large-scale genomic island integrations (columns) in 41 B. caccae isolates (rows) are shown (filled squares indicate presence). Hierarchical clustering of the 41 isolates by their IS614 integrations within each timepoint reveals four distinct subpopulations of B. caccae strains, three of which appear to shift in relative abundance between timepoints C and D. The IS614 integration upstream of norM is the only one that is absent from all isolates from timepoint A (before ciprofloxacin exposure), but appears in all isolates from timepoint C and D (after ciprofloxacin exposure). The IS614 integration upstream of per1 is present in all of the B. caccae isolates we obtained. Initial analysis of the isolate sequencing data was unable to detect IS614 integrations in front of per1 for two isolates, but manual inspection of the assembly graphs of these confirmed this integration to be present in these two as well (open squares in the per1 column). The IS614 integration upstream of thyA2 appears in a minority of strains in timepoint C (prior to trimethoprim exposure) and appears in the majority of strains in timepoint D (after trimethoprim exposure). The IS614 integration upstream of susC, which was detected in isolate sequence data, also appears to be under selection in timepoint C. Unlike the IS614 integrations, none of the large-scale genomic islands appear to be under selection between timepoints as they are broadly distributed across all subpopulations

Fig. 6

Antibiotic susceptibility testing of B. caccae isolates. Fifteen B. caccae isolates (3 isolates from timepoint A and 4 from each subpopulation present in timepoints C and D) were tested for their susceptibility to ciprofloxacin and trimethoprim, the two antibiotics that the patient was administered during treatment. The minimum inhibitory concentrations (MICs) of the two drugs against each isolate were determined as described in the “Methods”. a The presence (filled square) of different IS614 integrations within tested isolate strains that are upstream of annotated genes norM and thyA2, which are known to contribute to resistance to ciprofloxacin and trimethoprim. NorM is a known multidrug efflux pump that can confer resistance to ciprofloxacin. Upregulation of thyA2/dhfrIII has been shown to affect resistance to trimethoprim [61]. These IS integrations result in the potential introduction of likely bacterial promoters upstream of these genes likely leading to their upregulation increased expression and consequently as a result, increased antibiotic resistance. b MICs of ciprofloxacin against B. caccae strains. Strains from timepoints C and D with IS614 integrations upstream of norM and thyA2/dhfrIII predicted to increase resistance to ciprofloxacin (shaded) have a two- to fourfold increase in their MICs relative to strains from timepoint A. c MICs of trimethoprim against B. caccae strains. Overall, strains from all subpopulations showed high resistance to trimethoprim. Strains from the third subpopulation in timepoints C and D with IS614 integrations predicted to increase resistance to trimethoprim (shaded) show a twofold increase in resistance to trimethoprim compared to other subpopulations