Bacteriophages of Gordonia spp. Display a Spectrum of Diversity and Genetic Relationships

Grouping of Gordonia phages into clusters and subclusters.

The Gordonia phages span considerable genetic diversity, and comparisons with average nucleotide identities (ANIs) and dot plots reveal small groups of phages that are closely related to each other at the nucleotide level, but which are not closely related to other groups (Fig. 1A and ​andB).B). We initially assigned the Gordonia phages to clusters and subclusters using the same parameters as for the mycobacteriophages, forming 17 clusters and 19 singletons. However, this generates several instances in which phages of different clusters share large proportions of their genes and are likely to have common biological features. This is illustrated by the abundance of variable branch lengths in a SplitsTree network-based gene content phylogeny, creating examples where phages share a substantial portion of their gene content, without necessarily having nucleotide span-length matches exceeding 50% of the genomes (e.g., SoilAssassin and Yeezy [Fig. 1C]). We therefore chose to further refine the clustering parameters such that phages would be grouped into the same cluster if they share at least 35% of their genes (sorted into phamilies) with at least one other member of that cluster (Fig. 1 and Table 1). Thus, phages such as Attis and SoilAssassin that share 37 to 40% of their genes with Kita and BaxterFox are grouped together in Cluster CZ, whereas Orchid, Kampe, and PatrickStar (Cluster CX) share 28% of their genes with GMA3, Jumbo, GordTnk2, GordTuk1, and Gmala1 (Cluster DF), and clustered separately. These parameters place each of the phages in this analysis into one and only one cluster, however, it is conceivable that expansion of the collection with future isolates could yield a phage sufficiently similar to phages of two distinct clusters such as to warrant membership in both groups. Such an event would necessitate further refinement of our clustering parameters.

Open in a separate window

FIG 1 

Analyses of complete genome sequences of phages of Gordonia spp. (A) Heat map of average nucleotide identities (ANIs) of Gordonia phages. Pairwise ANIs were calculated using DNA Master, and the heat map was generated using R. (B) Dot plot of complete Gordonia phage sequences. (C) SplitsTree network of shared gene content between Gordonia phages. All genes were grouped into phams using Phamerator. Each phage was scored by the presence/absence of phams, and the distance between each phage was calculated using SplitsTree. The scale bar indicates 0.001 substitution. Phages are colored according to cluster membership.

Using these revised parameters, 65 of the phages can be grouped into a total of 14 clusters, and the other 14 have no close relatives and are designated “singletons” (Table 1). Thirteen of the clusters contain only Gordonia phages, but a group of six Gordonia phages form a subcluster (A15) within the large group of Cluster A phages that previously exclusively contained mycobacteriophages (Table 1). Seven of the Gordonia phage clusters (CQ, CR, CS, CU, CW, CZ, and DF) can be subdivided into subclusters (Table 1). By comparison, when 60 mycobacteriophages had been sequenced, they were grouped into nine clusters (five with subcluster divisions) and five singletons (9), and when 80 were sequenced, they were grouped into 10 clusters (seven with subcluster divisions) and five singletons (36). Thus, at comparable sample sizes, the Gordonia phages are at least as diverse as the mycobacteriophages (28 clusters plus singletons for Gordonia phages versus 15 clusters plus singletons for mycobacteriophages). Genome maps of these phages are included in Fig. S1 to S27 at figshare (https://doi.org/10.6084/m9.figshare.5149663).

A diverse array of systems for immunity and lysogenic maintenance in Gordonia phages.

The Gordonia phages form a variety of plaque types, ranging from clear to turbid, and the majority (53 of the 79) have genomic features associated with temperate phages, including tyrosine-integrase (Int-Y) or serine-integrase (Int-S), parABS partitioning systems, or repressors; only phages in Clusters CR, CS, CT, and DE and singletons Terapin, Ghobes, GMA2, and GMA6 lack these (Table 1). For 27 of the phages encoding a tyrosine-integrase (Int-Y), a putative attP common core sequence (20 to 40 bp) could be identified by similarity to a putative attB common core in the G. terrae 3612 genome (Table 2; see Table S1 at figshare). A total of 13 different attB sites were identified, 10 of which overlap host tRNA genes (Table 2). Interestingly, about half of these phages have the features of integration-dependent immunity systems described previously for mycobacteriophages (37), in which the attP site is located within the repressor gene (located immediately upstream), and the extreme C terminus of the repressor contains an ssrA-like degradation tag (Table S1). In some phages (e.g., Splinter and Vendetta) there are two potential attP sites corresponding to two different attB sites (Table 2; see Table S1 at figshare [https://doi.org/10.6084/m9.figshare.5149663]), but it is unclear if only one or both are used for integration. In GMA1, there is a single putative attP site, but it corresponds to two potential attB sites, suggesting that GMA1 can form double lysogens in which two prophages can be established in the same cell. We could not identify the attP site for about half of the Int-Y phages, but repressors of some phages (e.g., Nyceirae) have putative degradation tags and thus may also use integration-dependent immunity regulation. The attP sites for the serine-integrase (Int-S) phages could not be bioinformatically predicted as their attP and attB sites only have very short segments of sequence identity (38).

Several phages are unusual in having two integrase-like genes located near the centers of their genomes, although some are partial genes (e.g., Twister6 gene 47 and Wizard gene 44) and are likely to be defective (Fig. 2). For example, Wizard and Twister6 have seemingly intact integration-dependent immunity systems (37) but also have an Int-S-like gene that lacks the catalytic N-terminal domain (38) (Fig. 2). However, Utz, Howe, and Bowser all appear to have two intact integrase genes, all of which are Int-Y genes, except for Utz gene 33; Utz gene 31, Howe gene 36, and Bowser gene 37 are closely related and share nucleotide sequence similarity (Fig. 2). Curiously, the attP site identified in Howe is absent from Bowser and Utz, such that even though the Int-Y proteins are related (71% amino acid identity), the integration system may be nonfunctional for lack of an attP site. Bowser’s second Int-Y has features of an integration-dependent immunity component and contains a putative ssrA-like tag at its C terminus, as does the repressor gene (gene 40) (Fig. 2). However, we have not been able to identify a putative attP site either in the canonical locations within the repressor gene or elsewhere, and it is plausible that this integration system evolved to operate in a host other than G. terrae. We also note that the attP site associated with the Howe gene 36 Int-Y can potentially use two different attB sites (Fig. 2 and Table 2; see Table S1 at figshare). Finally, in Clusters CX and DF, the putative integrase genes are located in the middle of the right arm genes and lack the N-terminal domain associated with lambda-like Int-Y’s, and attP sites could not be identified. Although these systems are complex, it is not yet clear which of these recombinases mediate integration and which may be involved in recombination functions other than in phage integration (14). Taken together, the Gordonia phages show an unusual level of complexity of these site-specific recombinases that is not seen in phages of other actinobacterial hosts.

Open in a separate window

FIG 2 

Gordonia phage genomes exhibit multiple integrases. Phamerator map of the integration cassettes of the six phages with double integrases. Integrases are labeled with catalytic residue (Y or S), and attP sites are indicated when known. Wizard gene 44 and Twister6 gene 47 are predicted to be nonfunctional.

Potential for prophage-mediated defense systems.

We have previously described a series of prophage-mediated defense systems encoded by Cluster N mycobacteriophages that defend against heterotypic viral attack (39). The genes conferring defense are located adjacent to the integration/immunity cassette in a region that is highly variable within a set of otherwise closely related genomes (39). Several groups of Gordonia phages share this genomic architecture, notably in Clusters CV, CZ, and DB (see Fig. S7, S10, and S11 at figshare [https://doi.org/10.6084/m9.figshare.5149663]). The phages within Cluster CV have 4 to 10 genes in the region between the minor tail protein genes and the integration/immunity systems, and although the functions of most of these are unknown, there are several with putative functions similar to those in Cluster N phages. These include Blueberry gene 33, which possibly codes for a restriction protein, CaptainKirk2 genes 34 and 35, and Utz genes 34 and 35 coding for toxin-antitoxin (TA) pairs, and genes (e.g., CaptainKirk2 genes 32 and 37) coding for predicted membrane proteins. Likewise, in Cluster CZ phages, genes in this region code for candidate TA systems (e.g., Nymphadora genes 42 and 43 and Zirinka genes 33 and 34) and membrane proteins (e.g., Kita gene 33), and BaxterFox genes 31 and 32 are analogous to the MichelleMyBell viral defense genes 29 and 30. Similarly, in the Cluster DB phages, there are 6 to 8 genes in this location, including genes for putative antitoxins (e.g., Bowser gene 35 and Schwabeltier gene 33) and membrane proteins (e.g., Schwabeltier gene 35); Bowser gene 31 codes for a putative lipase. Several other genomes also code for TA pairs, including Bantam and Eyre (see Table S2 at figshare). Taken together, we predict that the temperate Gordonia phages are replete with novel viral defense systems that warrant further investigation.

Lysis systems in Gordonia phages.

Most mycobacteriophages code for three closely linked gene products involved in lysis—a holin, an endolysin (lysin A), and a mycolylarabinogalactan esterase (lysin B)—although some (e.g., Che12) lack the lysin B gene (40). Because the lysin B enzyme is responsible for cleavage of the mycobacterial outer membrane during lysis and because Gordonia cells contain mycolic acids, it is not surprising that most (all except for Clusters CR, CW, and DE and the singleton BetterKatz) of the Gordonia phages carry lysin B genes. However, there is great variation in where and how the lysis cassettes are organized. In many genomes (e.g., Cluster CR), the endolysin and holin genes are located downstream of the tail genes, a position where they are commonly located in the mycobacteriophages. However, in several phages (e.g., Clusters CV and CZ), they are found atypically within the minor tail protein genes. In a few genomes (e.g., Clusters CV, CZ, DB, DC, and DF), the lysin B genes are near the other lysis functions, but in most, they are displaced from it and are commonly positioned with the right arm genes. Finally, in some genomes (e.g., Cluster DF) the multiple domains found in mycobacteriophage endolysins (41) appear to be split into different genes.

Atypically small Gordonia phage genomes.

The Cluster CW phages have genomes that are much smaller (17.1 to 17.7 kbp) than the smallest of the mycobacteriophage genomes (41 kbp), although phages with similarly sized genomes have been found for other actinobacterial hosts, including Rhodococcus (e.g., RRH1, 14.3 kbp) (42) and Arthrobacter (e.g., Maggie, 15.5 kbp) (61). All of these have siphoviral morphologies, and more than 80% of the genome coding capacity is devoted to the virion structure and assembly genes (Fig. 3). Although the small-genome phages of Arthrobacter, Rhodococcus, and Gordonia are not closely related, some of the gene products have amino acid sequence similarity, including all of the large terminase subunits (31 to 68% amino acid identity [Fig. 3]). Although RRH1 and Maggie are seemingly lytic, the Gordonia Cluster CW phages are temperate and form turbid plaques, and these carry integrase-dependent immunity systems as noted above. Although the repressors of GMA5 and GRU3 have diverged from those of Jeanie and McGonagall, all share an identical 37-bp attP site corresponding to an attB site overlapping a tRNA^Lys gene (BMC27_10365) in G. terrae 3612 (see Fig. S8 at figshare [https://doi.org/10.6084/m9.figshare.5149663]). GMA5 deviates somewhat from the canonical integration-dependent immunity system organizations in that an additional open reading frame (gene 18) of no known function is inserted between the integrase and repressor genes (Fig. 3); a homologue is present in Maggie (gene 20), although it lacks an immunity system (Fig. 3).

Open in a separate window

FIG 3 

Small-genome actinobacteriophages. Shown are genomic maps of the small-genome phages included in the Actinobacteriophage_789 database. The central rule indicates nucleotide position—with large tick marks every 1 kbp—and the colored boxes indicate predicted genes. Genes are colored according to pham membership, with the number on the top of each box representing the pham number, followed by the number of members of the pham in the database in parentheses. Pairwise DNA sequence similarity calculated by BLASTn is shown between adjacent genomes and is spectrum colored, with violet being the most similar and red the least above a threshold E value of 10⁻⁵. Phages infect the following hosts: RRH1, Rhodococcus sp.; Jeanie, G. neofelifaecis; Maggie, Arthrobacter; McGonagall, G. neofelifaecis; GMA5, G. malaquae; and GRU1, G. rubripertincta.

Virion morphology variation in Cluster CZ.

Gordonia phage Cluster CZ is unusual in that it contains phages with distinctly different virion morphologies. All nine have long flexible noncontractile tails, but five (Howe, BaxterFox, Yeezy, Attis, and SoilAssassin) have isometric heads, whereas the other four (BatStarr, Kita, Nymphadora, and Zirinka) have prolate heads (Fig. 4). Genome comparisons show that the two sets of genomes differ in the leftmost eight genes that code for the terminase, portal, capsid maturation protease, and capsid proteins. Interestingly, the portal, protease, and capsid proteins of the prolate CZ phages have sequence similarity (65 to 67% amino acid identity) to homologues in the Cluster I mycobacteriophages, such as Che9c, which also have prolate capsids. This example illustrates a genetic link between the Gordonia and Mycobacterium phages and how the exchange of genes can result in morphological variation among phages with otherwise similar overall genomic architectures.

Open in a separate window

FIG 4 

Gordonia phages in Clusters CZ exhibit multiple virion morphologies. (Top) Genome maps of Gordonia phages Kita and Yeezy and mycobacteriophage Che9c (see Fig. 3 for details) displayed in two tiers. Yeezy and Kita are both Gordonia phages but have different capsid morphologies: Yeezy is isometric, and Kita is prolate. Che9c is a mycobacteriophage that shares capsid structure and assembly genes with Kita and is also prolate. (Bottom) Electron micrographs of Kita, Yeezy, and Che9c. Scale bar = 100 nm.

Relationships of Gordonia phages to phages of other actinobacterial hosts.

It is striking that six of the Gordonia phages (Table 1) fall into Cluster A—which previously contained solely mycobacteriophages—and form a distinct subcluster (A15) (Fig. 5; see Fig. S1 and Table S3 at figshare [https://doi.org/10.6084/m9.figshare.5149663]). The Subcluster A15 Gordonia phages are closely related to each other, and phage KatherineG shares the greatest sequence similarity with Mycobacterium phage Che12, Subcluster A2 (BLASTn score of 77% nucleotide identity spanning 47% of their genome lengths; they share 64% of their genes) (Fig. 5A and ​andB).B). However, the Gordonia Subcluster A15 phages are distant relatives of all other Gordonia phages, and only about 15% of their genes are present in other Gordonia phages. The Gordonia Cluster A phage homologues are distributed broadly among the other Gordonia phages, including representatives of Clusters CQ, CS, CT, CU, CW, CZ, DB, DC, DE, and DF and singletons BritBrat, GMA4, GMA6, and Yvonnetastic.

Open in a separate window

FIG 5 

Relationships between Gordonia and Mycobacterium phages. (A) Genome maps of Gordonia phage KatherineG and mycobacteriophages Phlei and Che12 (see Fig. 3 for details). Stoperators are indicated with a + or − above or below the map to indicate sequence orientation. (B) SplitsTree (60) network representation of shared gene content among the Cluster A phages. Phage groups are colored according to subcluster. The positions of Che12 (A2), Phlei (A13), and KatherineG (A15) are shown. (C) Shared gene phamilies of phages of Mycobacterium, Gordonia, and Arthrobacter. Phamily membership was determined using Phamerator, and the proportions of shared phams were calculated: 24 phams out of 8,294 total are shared between phages of all three hosts, 395 phams are shared between Gordonia and Mycobacterium, 57 phams are shared between Gordonia and Arthrobacter, and 56 phams are shared between Arthrobacter and Mycobacterium. The percentages shown are determined relative to the total number of phams present in each host (the number in boldface in parentheses in each circle).

The relationship of Gordonia Subcluster A15 phages to the mycobacteriophages is distinct and different from those of the other Gordonia phages (Fig. 5C). Of the 2,848 gene phamilies in the Gordonia phages, only 13.9% are shared with the Mycobacterium phages; the proportion of shared phamilies is only 10% if the A15 phages are excluded. Although few genes are shared, this is still greater than the number of phamilies shared with phages of the more distantly related host Arthrobacter (1.9% [Fig. 5C]), although this likely also reflects the differences in the total numbers of phages and phamilies that have been identified. As additional phages of different hosts within this phylogenetic space are characterized, further genetic connections between the phages are anticipated.

Cluster A is the largest group of related actinobacteriophages (192 in the database used here) and contains 16 subclusters (see phagesdb.org for a complete list of all current actinobacteriophage clusters and cluster membership). They are temperate (or recent lytic derivatives of temperate parents) and share overall genome architectures. They also have a common but unusual immunity regulation scheme in which the repressor binds to multiple binding sites (operators and “stoperators” [43]) located intergenically throughout the genomes. The Cluster A genomes all contain either an integration system using a tyrosine-integrase (44) or serine-integrase (45), although some from Subcluster A6 encode parABS partitioning systems (46, 47). The Subcluster A15 Gordonia phages are closely related to each other, all contain a parABS partitioning system, and the repressors are >99% identical to each other, indicating these likely form a homoimmune group. The A15 phages each contain 21 or 22 stoperator sites corresponding to the 13-bp asymmetric consensus sequence 5′-GGGGATTGTCAAG. The sequence conservation among the sites is similar to that reported previously (36, 43, 48), with positions 1, 9, 10, 11, 12, and 13 being invariant and positions 2 to 8 correlating with differences in immune specificity.

Gordonia phages display a spectrum of relatedness.

We noted above that the parameters used for cluster assignments were revised for grouping of the Gordonia phages, reflecting the lack of distinct boundaries between genome types and a spectrum of diversity as described for the mycobacteriophages (4). To explore the relationships between the Gordonia phages and other actinobacteriophages further, we calculated gene content dissimilarity (GCD) for pairwise genome comparisons between all 79 phages, where GCD is the proportion of shared phams relative to each genome’s total number of phams, averaging the two proportions and subtracting from 1 (see Materials and Methods) (49). Values range from 1, where no genes are shared (i.e., the two phages are 100% dissimilar) to 0, where the gene content is identical. The diversity of the Gordonia phages is illustrated by the range and distribution of GCD values (Fig. 6). A high proportion (88.6%) of pairwise comparisons have fewer than 10% shared genes (GCD, >0.9), with all 79 phages involved in at least one of the comparisons (Fig. 6). However, about 3% of the pairwise comparisons have values between 0.3 and 0.7 (i.e., reflecting 70% and 30% shared genes, respectively), and 38 of the 79 phages (47%) participate in at least one of these comparisons, representing 10 of the 28 groups (clusters plus singletons; 35.7%). Thus, although a value of 35% pairwise shared genes (GCD, 0.65) was used to group genomes into clusters, this is an arbitrarily chosen cutoff and not one that reflects a fundamental separation point among genome comparisons. It is helpful to note the cluster assignment parameters used here require 35% shared gene content with only one other phage, but because clusters (and subclusters) themselves can be quite diverse, the GCD values vary enormously within both clusters and subclusters.

Open in a separate window

FIG 6 

Gene content dissimilarity in phage populations. For all plots, the x axis is the ordered-by-magnitude individual pairwise comparisons, the y axis is the gene content dissimilarity (GCD [where 1 and 0 correspond to no shared genes and 100% shared genes, respectively]). (A) GCD in mycobacteriophages. (B) GCD in Gordonia phages. (C) GCD in Arthrobacter phages. (D) GCD in cyanophages. Both Mycobacterium and Gordonia phages exhibit a smooth curve of pairwise GCD values when ordered by magnitude. Cyanobacteriophages also exhibit a continuum of diversity with respect to GCD values; however, the slope of the line exhibits a number of plateaus, likely reflecting the “discrete” lineages described previously. Arthrobacter phages exhibit a large discontinuity, reflective of the few shared genes between clusters. Brackets indicate the interval between 30% and 70% gene content dissimilarity.

Bacteriophages of Arthrobacter provide an informative comparison (Fig. 6). Pairwise comparison of 47 phages of Arthrobacter (61) shows greater delineation between groups that are closely related and those that are unrelated than for the Gordonia phages, (Fig. 6). For example, only 2 of the 1,129 (0.12%) pairwise genome comparisons have GCD values in the 30 to 70% range, contributed by just three (0.5%) of the phages, representing 2 of the 13 groups (16.7%) (11 clusters, two singletons). For the 550 mycobacteriophages in the data set used here, 7% of the pairwise genome comparisons have a GCD in the 30 to 70% range, and 80% of the mycobacteriophages—13 of the 30 groups (43%)—participate in at least one comparison; similar proportions were identified when we repeated the analysis using three random samplings of 80 mycobacteriophages. The parsimonious explanation is that the mycobacteriophages and Gordonia phages are exchanging genes at a higher rate with other phages of the same host than the Arthrobacter phages are. We similarly examined a group of cyanophages as these were reported elsewhere to form discrete lineages (8, 10). These have a more complex GCD pattern (Fig. 6), with 27% of the comparisons within the 30 to 70% bracket (Fig. 6), and 188 phages of 209 cyanophages participate in these comparisons.

We thank Beckie Bortz, Emily Furbee, Sarah Grubb, and the SEA-PHAGES and PHIRE students at the University of Pittsburgh for phage isolation and annotation, Karen Klyczek, Fred Bonilla, and the SEA-PHAGES students at the University of Wisconsin—River Falls, for the isolation, annotation, and preliminary comparative analysis of Jumbo, Bantam, Remus, JSwag, and Strosahl, and Randy DeJong, John Wertz, and the SEA-PHAGES students at Calvin College for the isolation, annotation, and preliminary comparative analysis of Cucurbita, and members of the Hatfull lab for helpful discussions.

This work was supported by funding from National Institutes of Health grant GM116884 and Howard Hughes Medical Institute grant 54308198 to G.F.H. and by National Science Foundation Graduate Research Fellowship grant 1247842 to T.N.M. The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.