High Prevalence and Genetic Diversity of Large phiCD211 (phiCDIF1296T)-Like Prophages in Clostridioides difficile

Most C. difficile isolates harbor one or more prophages in their chromosome, either as integrated prophages or as episomal circular DNAs (6–9, 12, 17, 28, 29). We sought to determine whether strain ATCC 9689 (DSM 1296^T, here referred to as CD211) carried functional prophages. To investigate this, we used a classical mitomycin C induction strategy. Early-log-phase cultures (optical density at 600 nm [OD₆₀₀] of 0.1) of strain CD211 were treated with various concentrations of mitomycin C, and we monitored the optical density over time. A typical lysis curve was observed at 1 μg/ml mitomycin C, suggesting the release of functional phage particles. We used this induction lysate in a spot test with multiple C. difficile isolates in search of a suitable replicating host. Unfortunately, despite screening over 50 different C. difficile isolates from various origins and ribotypes, none of them were susceptible to infection by the phage particles from the lysate.

We analyzed the crude induction lysate under the electron microscope and found only one type of phage particle, which we named phiCD211, in reference to the host. The phage is composed of an isometric capsid with a diameter of 76 ± 3 nm and a long noncontractile and flexible tail of 449 ± 24 nm in length (Fig. 1). The morphological characteristics thus are reminiscent of the Siphoviridae family of phages of the order Caudovirales (30).

Because no suitable host was available to propagate phiCD211, we could not amplify and purify the phage from single plaques. However, since we observed only one type of phage particle in the induction lysate, we purified the phage DNA directly from the cleared lysate and completed whole-genome sequencing using the Illumina HiSeq platform at the Pasteur Institute in Paris. The sequence was deposited in the ENA public database on 9 December 2014 under the accession no. LN681537, and its sequence was later used and cited in our paper published in September 2015 (19). Of note, the complete genome sequence of phiCDIF1296T, a 131,326-bp phage identified in C. difficile strain DSM 1296^T, which is the same strain as ATCC 9689, from which we induced phiCD211, was published in August 2015 (23, 31). Not surprisingly, phiCDI1296T is 100% identical to phiCD211 at the nucleotide level over the entire genome. Although phiCD211 and phiCDIF1296T are identical, the phiCD211 genome was released in public databases prior to phiCDIF1296T, and it was cited almost simultaneously with phiCDIF1296T (19). Moreover, two following publications made reference to phiCD211 (22, 26). Therefore, although phiCD211 is identical to phiCD1296T, we kept the original phiCD211 name in this work to avoid confusion with previous publications.

The double-stranded DNA genome of phiCD211 contains 131,326 bp and has a G+C content of 26.4%, which is slightly lower than the average G+C content of the host (∼29%) (32). Aside from phiCDIF1296T, which is 100% identical to phiCD211, nucleotide BLAST homology searches did not retrieve other sequences with significant coverage, suggesting that phiCD211 is a unique phage in its group. As proposed by Rashid et al. (22), phiCD211 would be the representing member of a new phiCD211virus phage genus. Of note, the phiCDIF1296T phage genome has been detected as an episome in the DSM 1296^T strain in the course of a whole-bacterial-genome sequencing project (31). Therefore, and despite the presence of integrase genes, the data suggest that this phage does not integrate into the chromosome of its host (23). We have not determined experimentally whether phiCD211 was integrated or not in CD211, but since this strain is identical to strain DSM 1296^T, phiCD211 is likely episomal, like phiCDIF1296T.

The phiCD211 genome encodes 190 putative open reading frames (ORFs), ranging in size from 34 amino acids (aa) to 2,176 aa. The few additional ORFs compared to the annotation of phiCDIF1296T (181 ORFs) (23) likely are due to differences in annotation algorithms and parameters. Of the 190 ORFs, 54 (28%) contained conserved domains according to InterProScan analyses, and a putative function could be assigned to 53 (28%) of them. Using PROKKA (24) and the PHASTER curated phage database (25), we obtained additional hits corresponding to proteins from known phages, even if distantly related to phiCD211. Overall, with this approach we could infer a function or find a related phage protein for nearly 40% of the predicted proteins. Detailed genome annotation is provided in Table S1 in the supplemental material, and a genomic map is shown in Fig. 2.

Annotation of the phiCDIF1296T genome has been detailed in the paper of Wittmann et al. (23), and here we will present only the features that were not previously described or for which further improvement can be provided. We identified 13 ORFs coding for structural proteins, and it is worth noting that ORFs 3, 5, 6, and 7 were related to the lactococcal phage 949, a large 115-kb Siphoviridae phage (33). Other ORFs were also related to large phages, like the structural protein ORF13 and the DNA polymerase III α subunit (ORF86), which were related to the 185-kb Clostridium botulinum phage c-st (34). Likewise, the hypothetical proteins ORF1, ORF161, and ORF189 were related to another large C. botulinum phage, D-1873 (35). Together, several gene products from phiCD211 matched to proteins found in other large phages, suggesting that structural features have been conserved among these large phages during evolution.

Among other interesting characteristics of the phiCD211 genome, 17 ORFs code for putative transcriptional regulators, suggesting that the regulation of the phage replication cycle is complex, and/or that the prophage could affect multiple bacterial functions. The phage also encodes a tRNA for serine (positions 128727 to 128817). Another striking feature of phiCD211 is the presence of 6 transposase genes, 3 of which (orf56, orf152, and orf153) are adjacent to 2 of the 3 integrase genes identified in the genome. The large ∼55-kb genome fragment delimited by the transposase genes orf56 and orf153 is flanked by regions containing several transcriptional regulators, integrases, and lysis genes, as well as a putative CRISPR (clustered regularly interspaced short palindromic repeats)-associated endonuclease Cas3 protein and a nearby CRISPR array comprising 4 spacers of 42 or 44 bp. This configuration with several integrase genes and transposases suggests that this large genome fragment was acquired through one or more horizontal gene transfer events resulting from unprecise excision of a former integrated phiCD211 prophage, acquisition of transposons, or both. Further analysis of the CRISPR array revealed that spacer 1 targets known C. difficile phages, more specifically sequences from phages phiCDHM11, phiCDHM13, and phiCDHM14. The spacer is 100% identical to the target sequence, which is localized in the lysogeny module of these phages (17). Spacers 1 and 3 were also detected (100% identity) in several other C. difficile genomes, but spacers 2 and 4 did not match known sequences (Table S2). Of note, phiCDHM11, phiCDHM13, and phiCDHM14 were not detected in the strains carrying phiCD211, further supporting that the CRISPR array in phiCD211 is active and prevents further infection by these phages.

The phage genome also encodes 3 DNA methyltransferases (ORF62, ORF63, and ORF89) and an antirestriction protein (ORF173), which could participate in protection of the phage against host defense mechanisms. We also identified a death-on-curing protein (ORF149) of the Fic/DOC family (36) that provides further evidence for the episomal nature of phiCD211, and loss of the prophage should lead to death of the host. The PROKKA + PHASTER analysis also allowed the identification of a putative protein of the AcrB/AcrD/AcrF family, which comprises efflux pumps generally involved in multidrug resistance (37). It is impossible at this point to comment on the impact of the prophage on antimicrobial susceptibility of C. difficile, since we have been unable to transfer the prophage into another host. Curing of the lysogen from its prophage could be attempted, but this approach would be quite challenging due to the presence of the Fic/DOC addiction system.

A putative YyaC spore protease (ORF68) was also identified in phiCD211. YyaC was shown to cleave small acid-soluble proteins (SASPs) in Clostridium acetobutylicum and acts like a germination protease (GPR) in Bacillus subtilis (38). Therefore, the presence of this protein in phiCD211 could affect sporulation and/or germination of the host. Interestingly, this strain sporulates poorly compared to other strains (39).

In summary, the phiCD211 genome possesses several interesting features, including several genes that could affect host physiology and fitness.

We performed a mass spectrometry analysis in order to detect structural proteins composing the phiCD211 virion. We could not propagate the phage on a suitable bacterial host, so high phage titers were not available. Therefore, we concentrated the phage lysate on filter membranes and used this concentrated lysate for direct mass spectrometry analysis. A total of 14 proteins were detected, of which 5 were predicted to be structural proteins, and 3 others are possibly structural proteins as well (ORF8, ORF13, and ORF34) based on their genomic location (Table 1). A single peptide was also detected for ORFs 10, 12, 15, 31, 36, 37, and 40, which are also putative structural proteins based on their genomic location and/or bioinformatic predictions (not shown). However, only the protein hits for which at least two peptides were detected were reported in Table 1.

New software designed to determine the nature of the physical DNA ends and phage genome packaging mechanism was developed. This software, called PhageTerm (26), was used to analyze the raw sequencing data obtained with phiCD211. The analysis revealed the presence of fixed left and right DNA ends with exact 378-bp direct terminal repeats (DTR) (Fig. S1). The DTR was not initially detected in phiCDIF1296T because the episomal prophage was sequenced (23), whereas phiCD211 was sequenced from the linear encapsidated viral genome, which contains the DTR.

We investigated the prevalence of phiCD211 in 2,584 published genomes from various C. difficile strains (19, 40). Table 2 reports the results obtained, along with the multilocus sequence type (MLST) inferred in silico. The phiCD211 and phiCD211-like prophages were detected in 5% of the strains analyzed (149 out of 2,584) and were distributed in 21 different MLSTs. Some MLSTs were more frequently lysogenized with phiCD211-like prophages, such as MLST 49 (13/29 strains, 45% prevalence), MLST 10 (14/43 strains, 33% prevalence), and MLST 58 (8/28 strains, 29% prevalence). This suggests that phiCD211-like prophages are widely distributed in C. difficile.

We next compared the 149 phiCD211-like phage genomes identified in C. difficile isolates to better appreciate their similarity. To do this, we first performed de novo assembly of the corresponding C. difficile genomes using short sequencing reads publicly available and then used the PHASTER detection tool to retrieve genomic regions containing the phiCD211-like phages (25). Using this approach, only 127 complete phiCD211-like genomes could be retrieved. The other 22 phiCD211-like phages detected by PHASTER were located on separate contigs that could not be joined due to assembly problems. Thus, they were not included for the whole-genome comparison because some parts could be missing or be misassembled. We next used Gegenees to compare the 128 genomes (including phiCD211) at the nucleotide level (41). The three large phages 949 (L. lactis; GenBank accession no. HM029250.1), SPβ (B. subtilis; accession no. NC_001884.1), and c-st (C. botulinum; accession no. NC_007581.1) were included in the analysis. Gegenees uses an all-against-all BLASTn algorithm using short 200-bp fragments to assess whole-genome similarity and nucleotide identity level. The comparison matrix (heatmap) revealed the presence of distinctive clusters of ≥4 phages sharing ≥80% nucleotide identity over their entire genome (Fig. 3). Some of these clusters comprised several phages, like group 8, with 16 members. Other smaller groups contained only 2 or 3 members (e.g., groups 7, 9, 11, and 12). Phages that shared less than 80% identity in pairwise comparisons with other members were assigned a different number (e.g., groups 1, 2, 4, 5, etc.). Overall, 35 distinct phage groups were formed. The reference phage phiCD211 shared similarity with groups 16 and 18 but with slightly less than 80% with some members of these two groups, so we decided to assign it its own group (group 17) (Fig. 3). A table version of the Gegenees results as well as detailed information on the 35 phage groups are provided in Tables S3 and S4.

Representative members of the 35 distinct phage groups were also compared on a whole-genome scale using BRIG (Fig. 2) (42). This type of representation allows rapid identification of conserved and divergent regions in the various phage genomes. For example, we can see in Fig. 2 that the genes coding for the transposases (black arrows in the phiCD211 map) are not conserved in all phages. The genes encoding the DNA methylases (ORF62 and ORF63), the putative cell surface protein (ORF50), the ORFs encoding phage repressors possibly involved in lysogeny (ORF133 to ORF137), and the antirestriction protein (ORF173) all vary significantly between phiCD211 and the other phages. Note that because the reference phage in this analysis was phiCD211, genes that were present in other phages but that were absent from phiCD211 could not be visualized in the alignment. Therefore, additional genomic variations between the groups (insertions and deletions) likely exist. As an example, we identified one phiCD211-like phage with a significantly larger genome (phage 11, ERR339803, ∼170-kb). When we used it as the reference genome in a BRIG analysis, other genomic differences were observed (Fig. S2).

We next selected one representative phage from the 11 groups that contained at least 4 members sharing ≥80% overall identity (Tables S3 and S4) and aligned them using Phamerator to get a more detailed view of the regions of homology (Fig. 4). The phiCD211 phage was included as a reference. This analysis enabled us to visualize regions of high similarity and to identify multiple insertions, deletions, and inversions between these groups of phages. The first 16 ORFs, among which some likely code for capsid proteins based on their genomic location (next to the terminase) (22), were highly conserved in all phage groups. Next to these ORFs, there was a region comprised of ORFs coding for tail structural proteins that seemed to diverge between the groups. The regions flanking ORF16 and ORF25, encoding the tail tape measure in phiCD211, seemed to be prone to recombination, which likely explains the variations observed in the intervening region. Likewise, ORF31, encoding a putative regulator of chromosome condensation, showed a high degree of sequence variability between the phage groups despite strong similarity in the flanking regions. This could also represent a hot spot of recombination that led to modular reorganization of phage genes. In agreement with this, a large genomic region has been inverted in phage ERR340162 (phage 75) compared to the other phages (Fig. 4). Also, several ORFs from phage ERR251779 (phage 6) were not conserved in the region comprising ORF19 to ORF80. This further supports the hypothesis that these regions represent hotspots of recombination which led to the exchange of functional modules, including structural genes composing the phage tail.

The C. difficile strain from which phiCD211 was isolated is ATCC 9689, but for simplicity and to be consistent with previous publications, we will refer to that strain as CD211 (9). Bacteria were routinely grown at 37°C in TY broth (3% tryptose, 2% yeast extract, pH 7.4) inside an anaerobic chamber (Coy Laboratories). All media were prereduced overnight under anaerobic atmosphere (5% CO₂, 10% H₂, and 85% N₂).

To induce resident prophages from CD211, early-log-phase bacteria (OD₆₀₀ of 0.1) were treated with 1, 3, and 5 μg/ml of mitomycin C (BioShop Canada), as previously described (9). The OD₆₀₀ was monitored following mitomycin C addition to observe the typical drop in cell density that is characteristic of prophage induction (57). Crude lysates were filtered through 0.45-μm syringe filters (Sarstedt) and stored at 4°C.

Transmission electron microscopy (TEM) was performed as described before (9). Briefly, phage particles were washed 3 times in ammonium acetate buffer and transferred onto 400-mesh Formvar/carbon-coated copper grids. Negative staining was performed with 2% uranyl acetate, and phage particles were observed using a Hitachi H-7500 transmission electron microscope equipped with a Hamamatsu 10-megapixel digital camera controlled with AMT (Advanced Microscopy Techniques) software. The average size of viral capsids and tails was determined from five different images of isolated phage particles.

Phage DNA was isolated from 50 ml of clarified induction lysate using the Qiagen Lambda minikit by following the manufacturer's recommendations. The quality of the DNA preparation (e.g., the absence of smears suggesting partial degradation) was verified by running an aliquot through a 0.8% agarose gel. The absence of contaminating bacterial genomic DNA was assessed by PCR amplification of the bacterial gene triosephosphate isomerase (tpiA) as described before (58). Single-end multiplex libraries (TruSeq DNA Library Prep kit) were created with a multiplexed protocol according to Illumina's specifications. The sequencing was performed on the Illumina HiSeq 2000 platform of the Pasteur Institute (Paris, France), which generated 110-bp reads. The sequencing reads were first scanned to remove the adaptor sequences and then were assembled de novo into a single contig using Velvet (59) (with the following options: Kmer, 107; -short; -fastq; -unused_reads; -read_trkg). We then used the MicroScope workflow (60) to generate an automatic functional annotation for each ORF, which was manually curated afterward. We first performed a standard annotation using BLASTp searches against the NCBI database and InterProScan conserved domain analyses. To get further insights on potential protein functions and relationships with other phages and to improve function assignment and overall annotation quality, we also performed a complementary genome annotation procedure using PROKKA (24) run locally with the most recent PHASTER-curated phage database (25) and using an E value threshold of 10⁻³.

Phage particles from 60 ml of crude phage induction lysate were concentrated to ∼1 ml by ultracentrifugation for 1 h at 100,000 × g (Beckman SW 41ti rotor). The phage sample was further concentrated to a final volume of ∼0.2 ml with Amicon Ultra 0.5-ml 3-kDa columns (EMD Millipore). Following concentration, the sample was digested with trypsin and analyzed by liquid chromatography-tandem mass spectrometry (LC-MS/MS) at the proteomic and mass spectrometry platform of the Cancer Research Pavilion at the Faculty of Medicine of the Université de Sherbrooke (Sherbrooke, Québec, Canada). Only the proteins for which at least two peptides were detected were kept in the final table.

The multilocus sequence types (MLST) (61) were inferred from raw sequences of 2,584 published strains (19, 40). To detect the presence of phiCD211 in other C. difficile strains, we mapped the raw sequencing reads of each strain on the sequence of the phage using Bowtie2 (27). A strain was considered to carry phiCD211 or a closely related prophage when the overall genome coverage was ≥78%. This value was determined after conducting preliminary analyses with different cutoff values. Generally, ≥80% gave robust results, but several additional phages could be detected in C. difficile strains if we slightly lowered the cutoff to 78% instead, so we decided to include them to expand the number of phages in our analyses.

The packaging mode and the nature of the phiCD211 genome termini were determined using PhageTerm software, version 1.0.8, with default options. PhageTerm uses randomly fragmented sequencing reads and the assembled phage genome to calculate coverage values in order to determine information on the termini (26).

Whole-genome nucleotide similarity analysis and generation of heatplots were performed using the Gegenees program, v2.2.1 (41) (with the following options: fragment size, 500; step size, 500; no identity threshold). The analysis was performed without threshold for the filtration of nonconserved sequences. Excluding this threshold means that the values shown in Fig. 3 are based on both genome size and genome conservation. A cluster was defined by at least 4 phage members with an overall pairwise sequence identity of ≥80%. This value was chosen based on our first screen using Bowtie2 that used a threshold coverage of 80%. Phage genomic maps and alignments were generated using the Phamerator program (62) with default options. Overall identity and E value were used to group similar gene products into phamilies (where pham indicates a family of homologous proteins with an identity threshold of ≥32.5% and an E value of ≤10⁻⁵⁰). Phamerator also performed a pairwise DNA alignment between the genomes and used a spectrum of colors for shading areas to indicate the level of similarity, from purple for the highest similarity to red for the lowest similarity and no color when there was an absence of DNA similarity using BLASTn and a threshold value of 10⁻⁴. Prior to performing the alignments, the phage genomes were reverse complemented when necessary and manually reorganized to set the start position at the terminase gene. Arbitrary numbers and colors were automatically assigned to each family (phams) by the program. Multiple-genome comparisons were also performed using BLAST Ring Image Generator (BRIG). The local similarity of each phage representative was calculated based on BLASTn high-scoring pairs and plotted against the reference phage, phiCD211 (42).

The complete genome sequence of phiCD211 initially was deposited in the European Nucleotide Archive under accession no. LN681537 and in the National Center for Biotechnology Information (NCBI) database under the accession no. NC_029048.1. This entry was later replaced with an updated version of the genome, starting at the physical genome termini determined with PhageTerm (25) and including the direct terminal repeats detected in the encapsidated viral particles (accession no. NC_029048.2).