Genome Analysis of Epsilon Crass-like Phages Revealed That Some of Them Can Be Temperate Phages Encoding Atypical Diversity-Generating Retroelements

1. Introduction

CrAssphages were the first phages that were officially recognized by the International Committee for Taxonomy of Viruses (ICTV, https://ictv.global/taxonomy) despite the fact that they had not been isolated as individual entities [1,2]. CrAssphages have been discovered when analyzing NGS sequencing data using cross-assembly software [3]. Previously, phages of uncultivated prokaryotes were not systematized and incorporated into the modern taxonomic system. The discovery and recognition of the CrAssphages ensured the search and study of new phages among the metagenome assembled genomes (MAGs). Despite the fact that some phages that were found in the metagenomes were subsequently cultivated and isolated, most of them have been studied only in terms of the organization of their genomes [4,5,6].

СrAss-like phages are double-stranded DNA viruses that infect bacteria from the phylum Bacteroidetes [1]. Previously, all CrAssphages were assigned to the order Crassvirales and divided into four groups (families): Alpha, Beta, Gamma and Delta [2]. ICTV recently approved new family names for crAss-like phages: the Intestiviridae, Steigviridae, Crevaviridae, and Suoliviridae families (corresponding to the former Alpha, Beta, Gamma, and Delta groups, respectively). These four families currently include ten subfamilies, 42 genera, and 73 species [6].

The first crAss-like phage that has been isolated from the culture of Bacteroides intestinalis was p-crAssphage or ΦCrAss001 (now Kehishuvirus primaries from the Steigviridae family, former beta crAss-like phage group) [4]. Later, another 30 crAssphages were identified, cultivated, and studied that allowed them to be attributed to 11 different species of the Steigviridae family. Most of them infect Bacteroides intestinalis; three species infect Bacteroides cellulosilyticus and one is specific to Bacteroides tethaiotaomicron [7,8]. In 2021, the ΦCrAss002 phage belonging to the Jahgtovirus secundus species, the Intestiviridae family (former alpha group of сrAss-like phages) was isolated from Bacteroides xylanisolvens [5]. No сrAss-like phages of other families were isolated from the bacterial culture. Many cultured crAss-like phages require long-term cultivation, while even at high titers they do not lyse planktonic cultures and do not form plaques on a lawn of sensitive cells, or plaques are poorly visible [5,6]. It is assumed that members of the Crevaviridae family were observed by electron microscopy of fecal filtrate containing significant amounts of crAss-like phages [2].

The study of crAssphages has shown their important place in the intestinal microbiota and revealed their wide prevalence and great diversity [6,9,10,11]. It has been shown that certain groups of crAss-like phages dominate in industrialized countries, while other crAss-like phages are more common among the rural population of developing countries. The prevalence of crAss-like phages in the urban population is substantially higher than in the rural population [8,12,13]. This fact can be partly explained by the association of Bacteroides with the diet of the urban population [14]. Due to the high prevalence of crAssphages in the human intestine, their monitoring in the environment has been currently proposed to control the level of contamination with human faeces [15,16,17].

Recently, two new groups of Crass-like phages have been identified among the intestinal circularized MAGs [18]. These groups that were named Zeta and Epsilon showed sabstantial genetic distance from other crAss-like phages assigned to four approved families. CrAss-like phages of the Epsilon group do not encode DNA polymerases, unlike other members of the Crassvirales order, which mainly encode the B family DNA polymerase and sometimes A family DNA polymerase (some species of the Steigviridae family and zeta crAssphages) [18]. It has been proposed to separate CrAss-like phages of Zeta and Epsilon groups into new families [6,18]; however, ICTV currently does not assign them to any families.

In this study, a new epsilon crAss-like phage was found in the human intestinal metagenome. It was shown that this phage encode several atypical proteins (DNA primase, virion RNA polymerase, reverse transcriptase, and receptor binding protein) that substantially differ from those found in other crAss-like phages. The presence of the conservative prophage repressor and anti–repressors that belong to the lysogenic-to-lytic cycle switches indicated that the studied phage and several other epsilon crAss-like phages could be temperate.

2. Results

2.1. Analysis of the Phage Genome

Virome shotgun library was constructed from a stool sample from a healthy voluneer (Project: PRJNA1025976; BioSample: SAMN37731570; SRA: SRR26322353). After sequencing and de novo assembly, a 151,256 bp gapless contig was identified using the SPAdes genome assembler. An average coverage was 126 and it was a pseudo-circular MAG that suggested that this is the complete genome of the unknown Caudoviricetes phage. “Cut open” of this MAG was carried out using PhageTerm analysis and the starting point of the genome was predicted with a probability value 1.5 × 10^-22. The name crAssE-Sib was proposed for this phage. The complete genome sequence of crAssE-Sib phage was submitted into the GenBank database with the accession number OR575929.

Genome of crAssE-Sib phage encodes 202 putative ORFs, six of them encode tRNAs (Figure 1). Of the 196 ORFs, 76 encode proteins with predicted functions, which were determined based on the similarity of their amino acid sequences and domain structure with the described proteins. The remaining 122 ORFs were defined as hypothetical.

2.2. Comparative Analysis of the CrAssE-Sib Phage Genome

Comparison of the crAssE-Sib genome using BLASTn revealed the highest sequence identity with members of a new group of epsilon crAss-like phages. From this group, sequences with a length of more than 140,000 bp were taken for further analysis. Previously, sequences of epsilon crAss-like phages of this length have been characterized as pseudo-circular MAGs, indicating that these were complete phage genomes [18]. The analysis of the proteomic tree calculated using ViPTree (Figure 2) confirmed that the crAssE-Sib phage belongs to the epsilon crAss-like group of phages. In addition, comparative ViPTree analysis indicated that epsilon crAss-like phages formed two clades: clade I and clade II (Figure 2). Clade I contains 35 sequences, including crAssE-Sib; clade II combines 30 sequences (Figure 2). All epsilon crAss-like phages have similar genome lengths ranging from ~144 to 159 kb. Close to epsilon crAss-like phages there were phages infecting Bacteroidota: sixteen Cellulophaga phages and one Polaribacter phage. The lengths of their genomes are smaller than those of epsilon crAss-like phages and vary from 42 kb to 149 kb; however, most of genomes are 70–80 kb. Significant differences of epsilon crAss-like phages proteomes from their close relatives and other crAss-like phages allowed us to distinguish them into a separate family.

To clarify the taxonomy of epsilon сrAss-like phages, a matrix of intergenomic similarities was constructed using the Virus Intergenomic Distance Calculator (VIRIDIC) (Figure 3). A comparative genome analysis was performed on the base of intergenomic similarity, aligned genome fraction, and genome length ratio. The obtained results confirmed the division of epsilon crAss-like phages into two subgroups, which represent two proposed genera Epsilonunovirus and Epsilonduovirus in the group of epsilon crAss-like phages. It can be concluded that there are at least 31 different species in genus Epsilonunovirus, including crAssE-Sib, and 25 species in genus Epsilonduovirus. In some cases, individual species are represented by 2–3 sequences (Figure 3).

Comparative genome alignment of the epsilon crAss-like phages obtained using the ViPTree tool indicated high level of gene synteny both within the proposed genera and between them (Figure 4). Sequences with minimal nucleotide similarity (from 0.582 to 0.700) were taken for the analysis. There were 20 and 10 sequences of the Epsilonunovirus (includes crAssE-Sib) and Epsilonduovirus genera, respectively (Figure 4). At the beginning of the genomes, there is a region that is highly homologous between the related genomes. This region demonstrates the most similarity of sequences between the genera Epsilonunovirus and Epsilonduovirus. This region includes the genes encoding the terminase subunits and genes of structural proteins, namely portal, major capsid, and ring proteins. In addition, a significant homology was found in some other proteins of phages from both the Epsilonduovirus and Epsilonduovirus genera (carbohydrate-binding protein, muscle protein, and some proteins of the replicative apparatus).

2.3. Analysis of the crAssE-Sib Proteins

In order to verify the taxonomy of the crAssE-Sib phage and relative phages, phylogeny of the large subunit of terminase and DNA primase along with their most similar orthologs was analyzed. On the maximum likelihood phylogenetic tree of the terminase large subunit (Figure 5), sequences of epsilon crAss-like phages are divided into two clades according to the Epsilonunovirus and Epsilonduovirus genera, as shown in Figure 2 and Figure 3. In addition, they form an outer group to the rest of the crAss-like phages. This fact confirms that epsilon crAss-like phages should be isolated into a separate family. Notably, there are two additional sequences that are closely related but are not members of the proposed genera Epsilonunovirus and Epsilonduovirus. However, these sequences of the terminase large subunit are derived from short (less than 30 kb) nucleotide sequences, and taxonomic analysis of these phages is impossible.

Phylogenetic analysis of the phage primase showed that, as in the case of the large terminase subunit, two clades corresponding to the genera Epsilonunovirus and Epsilonduovirus can be distinguished (Figure 6). The primase sequences of epsilon crAss-like phages are significantly differ from those of other crAss-like phages indicating that epsilon crAss-like phages belong to a separate family.

In addition to the genes encoding the terminaze large subunits and DNA primase, several other genes that substantially differed from those in other crAss-like genomes were found in crAssE-Sib. Thus, the gene encoding the RNA polymerase that distinguished from those in crAss-like phages was revealed in crAssE-Sib (Figure 1, Data S1). HHPRED analysis showed that the crAssE-Sib virion RNA polymerase is similar to that of the phage phi14:2 with 95% probability [19]. This phage with a heterogeneous group of other phages formed a neighboring cluster on the ViPTree dendrogram (Figure 2). The crAssE-Sib virion RNA polymerase was conservative for epsilon crAss-like phages, and according to phylogenetic analysis, this polymerase significantly differed from the RNA polymerases of other crAss-like phages (Figure 7).

In addition, two transposases genes were found in the genome of the crAssE-Sib phage (Figure 1, Data S1). These two-directional genes encode proteins with an aa identity level of 18.2% between themselves and both transposases are members of the RNA-guided endonuclease TnpB protein family. The probable function of TnpB is to maintain a transposon in the genome. Notably, these genes are absent from the genomes of most epsilon сrAss-like phages.

Moreover, the crAssE-Sib phage genome encodes the NinG protein, transcriptional repressor, and two prophage antirepressors that both are the podoviridae phage lytic switch antirepressors (Figure 1, Data S1). These genes were found only in members of the Epsilonunovirus genus. The presence of such proteins indicates a temperate lifestyle of the phage.

2.4. Analysis of the Diversity-Generating Retroelement (DGR)

DGR cassettes are widely presented in the genomes of various Crass-like phages; some of them contain complete DGR cassettes, whereas only DGR fragments or nothing are found in other phages. In the crAssE-Sib phage genome, the reverse transcriptase (RT) gene was found using HPRED analysis. In addition, sequences located next to the RT gene were identified as the target gene (receptor binding protein, RBP), variable repeat (VR) located at the 3′-end of the target gene, template repeat (TR), initiation of the mutagenic homing (IMH) sequence that was found at the 3′-end of the VR, and IMH* sequence that was detectedat the 3′-end of the TR (Figure 8). All these sequences are essential elements of the complete DGR-cassette [20,21,22].

Screening of the available genomes of epsilon crAss-like phages showed that DGR cassettes were found only in a large group of phages from the Epsilonunovirus genus UNO, whereas DGRs were absent in members of Epsilonduovirus genus (Figure 4 and S1). Notably, all genomes of phages from the genus Epsilonunovirus encoded conservative RBPs. All VRs in these RBPs showed a certain amino acid (aa) pattern despite the variability in their aa sequences (Figure 9). Unlike the gene encoding RBP, the RT gene in some phages was fragmented or deleted indicating that retrohoming is not used by some members of the Epsilonunovirus genus. In addition, TRs were conservative in those DGR cassettes that were complete (Figure 9).

Analysis of the RT phylogeny of the crAss-like phages indicated that sequences from the Epsilonunovirus genus formed a monophyletic group that substantially differed from those from the Intestiviridae (alpha) and Suoliviridae (delta) genera as well as zeta crAss-like phages (Figure 10).

Three-dimensional (3D) modeling of RBP that was encoded by the crAssE-Sib phage was performed using Alphafold2 (Figure 11). Comparison of the obtained model using DALI server with structures of similar proteins showed that RBP has a structural similarity with proteins containing the C-type lectin domain (Figure 11 and Figure S2). The highest similarity of the crAssE-Sib RBP was revealed with the TaqVP protein (pdf 5 VF4, DALI Z-score 19.7) that is encoded by the gene located in the DGR cassette found in Thermus aquaticus (Figure 11 and S2) [23]. However, RBP of crAssE-Sib contained an additional N-terminal part that was absent in TaqVP and probably played a role in the protein attachment to the virion. RBP of crAssE-Sib also showed high similarity to the variable protein 1 from Treponema denticola (pdb 2Y3C, DALI Z-score 11.7), which is also hypermutated by DGR [24]. In addition, structural similarities were found with some enzymes, including the human formylglycine-generating enzyme (pdb 6MUJ, DALI Z-score 10.2). On the contrary, the structure of crAssE-Sib RBP differed from that of the target Mtd protein found in the DGR cassette of the BPP-1 phage (Figure 11). In the Mtd protein, the VR region is located on its surface, whereas it is hidden in a cavity in RBP of the crAssE-Sib phage.

3. Discussion

In this study, the complete genome encoding a new phage crAssE-Sib belonging to a group of epsilon crAss-like phages was found in the human virome. A detailed comparative analysis of the crAssE-Sib genome and genomes of available epsilon crAss-like phages indicated that all epsilon crAss-like sequences can be divided into two proposed genera Epsilonunovirus and Epsilonduovirus. This result is confirmed both by analysis of phages proteomes and phylogeny of phage proteins, including signature ones. In addition, the crAssE-Sib genome contained several genes that significantly differed from their orthologs found in members of the Intestiviridae, Steigviridae, Crevaviridae, and Suoliviridae families and zeta-group of crAss-like phages.

Thus, it has been previously suggested that DNA primase is the only protein that is conservative for all crAss-like phages, with the exception of epsilon ones [6]. It has been shown that the DNA primases of crAss-like phages form a monophyletic clade and probably descended from one crAss-like phage ancestor [1,6,18]. In this study, the primase genes were found in genomes of epsilon crAss-like phages; however, the primase of epsilon phages varied and were not part of the monophyletic clade with the primases of other crAss-like phages (Figure 6). We also found the gene encoding RNA polymerase in the crAssE-Sib genome (Figure 7). This finding became possible when the structure of the virion RNA polymerase of the phage phi14:2 had been available [19], since this enzyme has certain similarity with that of crAssE-Sib.

It was shown that the crAssE-Sib genome contains DGR cassette with all essential elements, namely the RT gene, target gene encoding RBP, VR, TR, IMH, and IMH*, and RBP undergoes mutagenic retrohoming in this phage (Figure 9). Despite the fact that the RT gene was detected with high confidence in the crAssE-Sib phage genome, the myDGR program [25] did not find the presence of DGR in crAssE-Sib and other epsilon crAss-like phages genomes. It can be explain by the fact that crAssE-Sib RT contains the GxxxSP motif in its fourth domain, instead of the characteristic GxxxSQ motif found in DGRs. Usually, the GxxxSP motif is detected in most retroviral and retrotransposon RTs and group II intron maturases [26] and finding this motif in RT from the DGR cassette was unexpected. When we changed this motif in RT of crAssE-Sib, myDGR perfectly detected the RT gene and hence DGR. It should be noted that myDGR is widely used to search and study DGRs. The presence of RT with unusual aa motifs in DGRs would prevent the detection and study of these DGRs. So, prevalence and diversity of DGRs in viral and prokaryotic genomes are insufficiently studied.

As for RBD of crAssE-Sib that is a target protein in DGR, no similar phage proteins were found in the Protein Data Bank (PDB). However, despite low sequence similarity, this RBD showed structural similarity with C-lec fold-containing proteins that are the predominant type of target proteins in DGR cassettes [23]. The highest structure similarity was found between crAssE-Sib RBD and the target protein from DGR of Thermus aquaticus (Figure 11 and Figure S2).

Surprisingly, the genes encoding prophage repressors and anti–repressors absent in crAss-like phages were found in the epsilon crAss-like genomes. It has been previously shown that antirepression system is associated with the switching of the life cycle in the temperate podoviridae phages [27,28]. There are two prophage anti-repressor genes and one transcription repressor gene in the crAssE-Sib phage genome (Figure 1). The alignment of these two prophage antirepressor sequences indicated their significant identity, which suggests that these antirepressors are probably paralogs (Figure S3). Analysis of other crAss-like phage genomes showed that these genes are unique to the Epsilonunovirus genus. Within this genus, one prophage antirepressor gene and one transcriptional repressor gene are conservative.

Attempts were made to detect prophage sequences related to epsilon crAss-like phages in the bacterial genomes. We found two incomplete prophage sequences in chromosomes of Parabacteroides distasonis BFG-238 and Parabacteroides merdae CL06T03C08. Comparative analysis of their proteomes with crAssE-Sib (genus Epsilonunovirus) and OLQV010 (genus Epsilonduovirus) suggested that these prophage sequences belong to unknown genera of temperate Parabacteroides phages. This supposed genus could be the third one in the epsilon crAss-like phages family. Notably, the prophage sequence from the Parabacteroides distasonis BFG-238 strain demonstrates close levels of similarity with the Epsilonunovirus and Epsilonduovirus genera, like these genera among themselves (Figure S4).

Based on sequence similarity of CRISPR spacers in microbial genomes and sequences of crAss-like phage genomes, it has been previously supposed that epsilon crAss-like phages infect bacteria from the Bacteroides, Parabacteroides, and Porphyromonas genera (class Bacteroidia) [9,18]. So, we can assume that the crAssE-Sib phage also infects bacteria from one of these genera.

Our search for the integrase gene in the epsilon crAss-like phages genomes using FASTER and IntegronFinder 2.0 failed. A number of genes involved in recombination processes have been discovered, and the gene encoding the NinG protein is one of them [29]. However, it is impossible to say whether the epsilon crAss-like phage genomes can be inserted into the host genome. It has been known that some temperate phages may not integrate their genomes into the host genome, but save them in the cell as episomes or linear genomes, and that the genes of such phage genomes are expressed along with bacterial genes [30,31,32]. Possibly, the lysogeny of epsilon crAss-like phages occurs by this mechanism. When crAss-like phages from the Intestiviridae, Steigviridae, Crevaviridae, and Suoliviridae families were studied, it was suggested that these crAss-like phages are lytic and incapable of lysogeny. However, these phages could be pseudolysogenic since they are capable of long-term persistence in the host culture, and the probable mechanisms of pseudolysogeny can be the carrier state or chronic infection hibernation inside infected cells [4,6,33,34]. As for members of the genus Epsilonunovirus, we probably found the true temperate crAss-like phages. Finding of prophages that are relative to this genus can confirm this assumption.

4. Materials and Methods

4.1. Phage Genome Sequencing

Isolation of viral DNA from a fecal sample was described in detail previously [35]. Briefly, the sample was clarified by repeated centrifugation and treated with DNase I, 5 U (Thermo Fisher Scientific, MA, USA) and Proteinase K 100 μg/ml (Thermo Fisher Scientific, MA, USA). Then, phenol-chloroform extraction and subsequent ethanol precipitation were used to purofy DNA. The resulting DNA was diluted in 50 µl of TE buffer and DNA concentration was measureed using Qubit 4.0 (Thermo Fisher Scientific, Massachusetts, USA). Virome shotgun library was constructed using the NEB Next Ultra DNA library prep kit (New England Biolabs, Ipswich, Massachusetts, USA). For sequencing, a MiSeq Benchtop Sequencer (Illumina Inc., San Diego, CA, USA) and a MiSeq Reagent Kit 2 × 250 v.2 (Illumina Inc., CA, USA) were used. After de novo assembly of the genome using the SPAdes software V.3.15.2 [36], the resultant single contig was further analyzed.

This work was approved by the Local Ethics Committee of the Center for Personalized Medicine, Novosibirsk (protocol #2, 12 Febriary 2019). The written consent of healthy volunteer was obtained in accordance to guidelines of the Helsinki Ethics Committee.

4.2. Phage Genome Analysis

With the help of PhageTerm [37] (https://galaxy.pasteur.fr, access date 25 August 2023), the termini of the genome that was obtained in a “pseudo-circular” form were determined. Online RAST server v. 2.0 [38] was used for the identification of ORFs and possible tRNA genes within the genome. The prediction of the tRNA genes was additionally performed by tRNAscan-SE V 2.0 [39]. Genome function annotation was done by InterProScan (www.ebi.ac.uk/interpro/search/sequence/), HHPpred (toolkit.tuebingen.mpg.de/tools/hhpred), HHblits (toolkit.tuebingen.mpg.de/tools/hhblits) and BLASTP (NCBI, Bethesda, MD, USA). Predictions on the architecture of the domain and functional motives were made using CD-BLAST (NCBI, Bethesda, MD, USA).

4.3. Comparative Phage Analysis

Identification of nucleotide sequences similar to the crAssE-Sib phage genome was performed using BLASTN. Comparative analysis of crAssE-Sib phage proteome was performed by ViPTree version 3.7 web server (https://www.genome.jp/viptree, accessed on 3 August 23) with default parameters [40]. ViPTree analysis based on genome-wide sequence similarities computed by tBLASTx. The intergenomic similarity was determined using the VIRIDIC tool (http://rhea.icbm.uni-oldenburg.de/VIRIDIC) [41]. The genera were predicted using thresholds of the intergenomic identity 70% [42,43].

4.4. Analysis of Proteins

Amino acid sequences with significant similarity to the studied ones were downloaded from GenBank using crAssE-Sib phage proteins as query in BLAST search against NCBI’s protein database. Protein sequences were aligned using the M-Coffee method in the T-Coffee program [44]. The levels of amino acid identity of proteins were calculated by the BioEdit 7.2.5 program [45]. Maximum Likelihood phylogenetic trees were generated by the IQ-tree program v. 2.0.5 [46]; the best fit substitution model according to ModelFinder [47] was used. Branch supports were assessed using 1000 ultrafast bootstrap replicates [48]. The resulting trees were midpoint rooted and visualized in FigTree v1.4.1. Analysis of the genes included in DGR was carried out by the myDGR program [25]. The search for phage integrase genes was performed using programs IntegronFinder 2.0 [49] and FASTER [50]. Sequence logos representations of amino acid and nucleic acid multiple sequence alignment were calculated by WebLogo version 2.8.2 [51].

4.5. 3D Modeling of Protein Structures

Ribbon and surface representations of proteins were predicted in the AlphaFold2 program available at https://colab.research.google.com/github/sokrypton/ColabFold/blob/main/AlphaFold2.ipynb [52]. The resulting structural models were downloaded and edited using UCSF Chimera molecular visualizer, version 1.17 [53]. 3D-model of the RBP part (aa residues 85−678), which was encoded by the crAssE-Sib phage genome, was modeled using Alphafold2 with a high degree of confidence. The analysis of the models was performed using DALI server [54]