Sequence analysis of malacoherpesvirus proteins: Pan-herpesvirus capsid module and replication enzymes with an ancient connection to "Megavirales".

The order Herpesvirales includes animal viruses with large double-strand DNA genomes replicating in the nucleus. The main capsid protein in the best-studied family Herpesviridae contains a domain with HK97-like fold related to bacteriophage head proteins, and several virion maturation factors are also homologous between phages and herpesviruses. The origin of herpesvirus DNA replication proteins is less well understood. While analyzing the genomes of herpesviruses in the family Malacohepresviridae, we identified nearly 30 families of proteins conserved in other herpesviruses, including several phage-related domains in morphogenetic proteins. Herpesvirus DNA replication factors have complex evolutionary history: some are related to cellular proteins, but others are closer to homologs from large nucleocytoplasmic DNA viruses. Phylogenetic analyses suggest that the core replication machinery of herpesviruses may have been recruited from the same pool as in the case of other large DNA viruses of eukaryotes. Published by Elsevier Inc.

Introduction

The order Herpesvirales consists of animal viruses with large linear double-strand (ds) DNA genomes of 125–145 kbp, encoding between 70 and 200 proteins. Herpesviruses are characterized by the following phenotypic attributes (Davison et al., 2009, Pellett et al., 2011): the mature particle contains the genome within a T = 16 icosahedron capsid, composed of 162 capsomers, of which one (the portal) has a distinct structure, specialized for mediating virus DNA entry and exit from the capsid; the capsid is surrounded by a proteinaceous tegument and further wrapped in a lipid envelope that contains several virus-encoded transmembrane proteins; the viral genome contains directed or inverted DNA repeats; the virus DNA replication mechanism generates head-to-tail concatemers that are later cleaved and individually packed into pre-made capsids.

These characteristic molecular and morphological features have been used to assign members to the order Herpesvirales even in the absence of the genome sequence information (reviewed in Davison, 2010). In the genomic era, however, the state-of-the-art taxonomy of herpesviruses is based on phylogenetic analyses of virus nucleotide and amino acid sequences (Davison et al., 2009, Davison, 2010, Pellett et al., 2011). For certain subsets of herpesviruses, additional phylogenetic signal can be extracted from the information about gene synteny and by estimating the minimal amount of rearrangements required to convert the gene order between two species (Hannenhalli et al., 1995, Bourque and Pevzner, 2002, Larget et al., 2005).

Analysis of these and other molecular features resulted in establishing three families within Herpesvirales that are currently recognized by the International Committee for Taxonomy of Viruses (ICTV; Davison et al., 2009; Pellett et al., 2011). The Herpesviridae family (GenBank TaxID 10292, with 73 completely sequenced virus genomes as of May 2017) includes viruses that infect amniotes; typically, a pair of virus species from this family shares dozens of protein-coding genes recognized by high amino acid sequence similarity (Davison, 2010). The family Alloherpesviridae (TaxID 548682, 8 completely sequenced genomes) includes viruses of fishes and amphibians. These viruses have been reported to share only 13 core orthologous genes with each other, suggesting that this group is more divergent than the family Herpesviridae (van Beurden and Engelsma, 2012). The third family, Malacoherpesviridae (TaxID 548685), includes two ICTV-approved members, Ostreid herpesvirus 1 (OsHV-1) and Haliotid herpesvirus 1 (HaHV-1), which infect bivalve and gastropod molluscs, respectively. In addition to the herpesvirus-like morphology of their virions and capsids, as well as the genomic repeat features, these two viruses were reported to have 39 orthologous genes in common (Savin et al., 2010), some of which are related to enzymes that play essential roles in other herpesviruses, for example, the large subunit of virus terminase ATPase and family B replicative DNA polymerase. However, no gene products with similarity to herpesvirus structural proteins have been reported in Malacoherpesviridae. Interestingly, a putative virus related to mollusc herpesviruses has been detected by sequence database searches, among the contigs co-assembled with the portions of the draft genome of a cephalochordate Branchiostoma floridae (Savin et al., 2010). Taken together, these genomic data suggest that the order Herpesvirales is characterized by higher genomic diversity, broader host range, and, likely, deeper phylogeny than was assumed just a decade ago.

In contrast to our understanding of the phylogenetic relationships within Herpesivirales, the evolutionary origin of this order as a whole remains uncertain. The picture is clearer in the case of structural proteins. One domain, “the floor domain”, in the large multi-domain capsid proteins of vertebrate herpesviruses, whose structure has been determined by a cryoEM approach, is structurally similar to bacteriophage capsid proteins with the HK97 fold (Baker et al., 2005, Yu et al., 2017). Moreover, mechanistic similarities in the maturation pathways of these diverse virus particles have been noted (Casjens and King, 1975, Booy et al., 1991, Davison, 2002, Mettenleiter et al., 2009, Veesler and Johnson, 2012), and sequence comparisons have established the ancestral relationship of the capsid maturation protease between phages and Herpesviridae (Cheng et al., 2004, Liu and Mushegian, 2004). That relationship is further supported by the demonstration that the high-resolution structure of phage protease is folded very similarly to the herpesvirus assemblin proteases, whose structure has been determined earlier (Fokine and Rossmann, 2016). Along with the long-known monophyly of the large (ATPase) subunit of phage and virus terminases (Mitchell et al., 2002), these observations point at the common origin of the entire head/capsid formation module in Herpesviridae and in a subset of tailed viruses of bacteria and archaea.

The state of affairs is different for genes involved in virus DNA replication and expression. For several established families of eukaryotic dsDNA viruses, which, unlike herpesviruses, replicate partly or fully in the cytoplasm (Nucleocytoplasmic Large DNA Viruses, or NCLDV), considerable evidence of the monophyletic origin has been derived from comparative genome analyses (Iyer et al., 2006, Yutin et al., 2009, Koonin and Yutin, 2010). The assembly of large DNA viruses consisting of NCLDV and several related lineages has been proposed as a candidate order “Megavirales” (Colson et al., 2013), not yet accepted by ICTV. Large dsDNA viruses replicating in the nucleus, such as herpesviruses and baculoviruses, have several replication enzymes homologous to those of NCLDV (see below), but the evolutionary scenario linking them together has not been established.

In this work, we analyzed the gene repertoire of malacoherpesviruses, using predicted open reading frames (ORFs) of HaHV-1 (synonym Abalone herpesvirus Victoria/AUS/2009) as the starting point. Our study characterized the enlarged repertoire of conserved malacoherpesvirus genes, and for many of them we identified similar sequences in other viruses or in the genomes of cellular organisms. Analysis of sequence similarities and phylogenetic inference on those gene families expand the common gene core of Herpesvirales, support the hypothesis of the common ancestry of their morphogenetic module, and suggest that the herpesvirus replication module may have been recruited from the same ancient gene pool as the replication genes of other large DNA viruses.

Methods

Sequence database searches were mostly performed in July-December of 2016, except for the results presented in Fig. 1 and the accompanying text, for which the searches were repeated in June-September of 2017. Protein queries representing every predicted ORF in HaHV-1 (GenBank Taxonomy ID 860344) were used to query the nucleotide sequence databases (NT and dbEST) and the non-redundant protein sequence database (NR) at NCBI. Searches were done with the BLAST family of programs (Altschul et al., 1997) with the “Composition-based statistics” option. The HHPred server (Söding et al., 2005; Alva et al., 2016) was used to compare profile-Hidden Markov Models (HMM) of virus proteins to profile-HMMs built from the entries in the NCBI CDD database (Marchler-Bauer et al., 2017), with multiple sequence alignment generation method set to PSI-BLAST, re-aligning the results with MAC set to 0.3, and the remaining options set to default values. Multiple sequence alignments were computed using the MUSCLE program (Edgar, 2004), joining closely related (family level and below) sequences, and combining the alignments obtained at that first step using the -profile option of MUSCLE. When distantly related proteins were to be aligned with structural constraints, the PROMALS3D server was used (Pei and Grishin, 2014). Spatial structures were visualized with the PyMOL legacy build 0_99rc6 (https://sourceforge.net/projects/pymol/files/Legacy/).

Fig. 1

Multiple sequence alignments of malacoherpesvirus morphogenesis proteins. GenBank accession numbers as well as ICTV-approved, sometimes truncated, but not italicized, virus species names are shown before each sequence. PDB ID numbers, when available, are also shown. The color codes of conserved amino acids are as follows: bold type and yellow shade, residues with bulky hydrophobic side chains (I, L, V, M, F, Y, W); red type, residues with negatively charged or amido side chains (D, E, N, Q); blue type, residues with positively charged side chains (K, R); white type and red shade, residues with turn- or kink-prone side chains (A, G, P, S). Other conserved residue colors are ad hoc. In the secondary structure lines, h and s indicate helix and strand, respectively. A. Floor and buttress domains of the major capsid proteins. The number of the leftmost amino acid residue in each block, counting from the N-terminus of the precursor protein, is shown before each sequence block. In the secondary structure line, l stands for a loop, and different colors match the loop elements highlighted in Additional Files 1 and 2. B. Procapsid proteases. The number of the leftmost amino acid residue in each block, counting from the N-terminus of the precursor protein, is shown before each sequence. Numbers in parentheses indicate distances between conserved sequence regions, in amino acids. Asterisks indicate the amino acids that form the catalytic triad of the protease. The wild-type sequence of protease from phage T4 (Escherichia virus T4) is shown, though two of the three catalytic residues and a non-catalytic methionine residue have been mutated in the sequence of 5jbl.

Phylogenetic inference was done using a local installation of the PhyML program as well as the PhyML 3.0 web server (Guindon et al., 2010) with the LG substitution model, most other parameters estimated from the data, and bootstrap replicates performed to assess the support of the internal partitions in the tree. The iTOL v.3 server (Letunic and Bork, 2016) was employed for tree examination and visualization.

Statistical tests for relaxed selective constraints in the lineages leading to B.floridae and Capitella teleta were conducted using RELAX (Wertheim et al., 2015). For each set of orthologous proteins, corresponding nucleotide sequences were extracted from GenBank, codons were aligned using the PRANK program (Löytynoja, 2014), and for each multiple sequence alignment, two p-values were computed (one for a relaxation signal in the lineage leading to B. floridae and one for a relaxation signal in the lineage leading to C. teleta). These p-values were computed using a chi-square distribution in a likelihood ratio test, as proposed in RELAX. In order to account for multiple hypotheses testing and control the false discovery rate, p-values were adjusted using the Benjamini-Hochberg procedure (Benjamini and Hochberg, 1995), as implemented in the R function p.adjust.

Results and discussion

Nearly half of malacoherpesvirus gene products are conserved

We compared 118 putative proteins encoded by the genome of HaHV-1 to the protein and nucleotide sequence databases, as well as the databases of known conserved protein domains (the CDD database at NCBI). The identities of the best statistically supported matches in viruses and cellular organisms were recorded; some of the matches had low statistical significance, but could be validated by observing the conservation of characteristic sequence motifs. The combined results of these analyses are presented in Table 1.

Table 1

Putative protein functions in malacoherpesviruses and orthologous relationships in four species. Gray shade, two types of repeats within the HaHV-1 genome. Blue shade, gene products with predicted structural and morphogenetic roles. Yellow shade, gene products with predicted roles in genome replication. Bold type, novel predictions.

Image 1

Image 2

Image 3

Image 4

Image 5

Image 6

Nearly 60 proteins in HaHV-1 have homologs in the better-studied malacoherpesvirus, OsHV-1 (and often also in the unclassified, but highly similar, malacoherpesvirus infecting scallops, Chlamys acute necrobiotic virus; data not shown). The BLASTP pairwise sequence identity within malacoherpesvirus sequences was generally between 20% and 40%, suggesting considerable divergence of these viruses. Some amount of synteny between the two virus genomes was observed, usually in blocks of 2–4 genes.

As reported before (Savin et al., 2010), substantial portions of the malacoherpesvirus genomes also show homology to a part of the draft genome assembly of the cephalochordate B. floridae. Most of the herpesvirus-like matches in the B. floridae genome are mapped to a single contig; some of these genes have been predicted by the original genome annotation, while others have not been reported by Savin et al. (2010), but could be identified by matching malacoherpesvirus protein sequences to the translations of the putative intergenic regions. We also detected a collection of homologous genes in a draft genome assembly of another marine invertebrate, annelid C. teleta. In this case too, at least some of the genes appear to be located on the same contig (Table 1).

Taken together, these matches provide evidence for evolutionary conservation at the family level for at least half of the genes in each malacoherpesvirus genome. Each of these herpesvirad genomes, in addition, has genes that were not found in other herpesviruses but are distantly related to gene products of a different evolutionary origin. In more than 35 of the conserved malacoherpesvirus gene products, we identified specific sequence motifs suggestive of their molecular function (Table 1 and see below).

Structural proteins: herpesvirus orthologs and phage connections

The major capsid protein in Herpesviridae (MCP; also known as VP5/UL19 gene product in Human herpesvirus 1, HHV-1, synonym Herpes simplex virus) is a large multidomain protein. Structural data for the entire protein within the assembled capsid, obtained by cryoEM, indicate that VP5 consists of the “floor” domain, which forms the icosahedral shell and is the only part of the protein that interacts with the capsid interior, as well as the “middle” and “top” domains protruding above the surface of the shell (Zhou et al., 2000). More recently, an assembly of additional domains has been recognized in the MCPs of HHV-1 (Hui et al., 2013) and of human cytomegalovirus (Human herpesvirus 5; Yu et al., 2017). A high-resolution X-ray structure is still available only for the top, or tower, domain of HHV-1, which is contiguous in the sequence of VP5, starting approximately at amino acid 500 (Bowman et al., 2003).

Our analysis detected a match between a profile-HMM seeded with the HaHV-1 ORF68 and the probabilistic model specified by the top domain of HHV1 VP5 (PDB ID 1no7; Bowman et al., 2003). The HHPred search initiated with HaHV-1 ORF68 matched the sequence family model of PDB identifier 1no7 with p-value of 10−4. The match is about 100 amino acids long, and is located close to the N-terminus of the proteolytically isolated top domain of VP5, whereas in ORF68 it is found in the middle, approximately between amino acids 570 and 670. In addition, a PSI-BLAST search against the complete NR database at NCBI (mined on June 14th, 2017) was initiated with the full-length ORF68. This database search retrieved the Malacoherpesviridae homologs at the first round, followed at the next round by the putative capsid protein sequence from Macacine betaherpesvirus 3 (synonym Rhesus cytomegalovirus; protein ID APT40174.1; p-value < 10−4) as well as their homologs from several other betaherpesviruses and gammaherpesviruses. The matches covered two-thirds of the protein lengths, starting close to the N-termini and extending to the same region of similarity as above.

To refine the alignment, we used the PROMALS3D approach (Pei and Grishin, 2014), which collects homologs of all input sequences by PSI-BLAST and realigns them using the compatibility of highly-scoring alignment segments, with the additional constraint imposed by the known tertiary structure of the homologous sequences. There were three regions of high sequence similarity between herpesvirus and malacoherpesvirus major capsid proteins. Two of these regions map to the terminal sequence segments of the tower domain, and the third region is located adjacent to the C-terminal end of the tower domain sequence (Fig. 1A).

The structure of the tower domain of herpesvirus MCP is described as a pyramid with a square base (Bowman et al., 2003). Much attention has been given to the structural elements constituting the sides and the apex of the pyramid, as these sites appear to be involved in multiple protein-protein interactions (Hui et al., 2013). Less is known about the base itself, other than that its structure comprises mostly loops and a few helices (Bowman et al., 2003). Our results, however, indicate that the base is made of the best-conserved portion of the MCP sequence; indeed, conserved regions 1 and 2, distal in the sequence as they are, come together in the three-dimensional structure to account for almost the entirety of the base (visualized in Additional files 1 and 2). Conservation of these regions in the two families of Herpesvirales suggests that the foundation of the tower domain plays an important role in maintaining the integrity of the multidomain MCP and, possibly, its interactions with other virion components.

Our efforts to identify a homologous region in Alloherpesviridae, where mass-spectrometry approaches have assigned the major capsid protein role to ORF39 in Ictalurid herpesvirus 1 (Davison and Davison, 1995) and to its ortholog ORF104 in Anguillid herpesvirus 1 (Van Beurden et al., 2011) have produced ambiguous results. When ORF39/ORF104 and their orthologs from other alloherpesviruses were used as HHPred queries to search the available HMMs, no matches were found. In a targeted PSI-BLAST search initiated by these sequences against all virus proteins in the NR database, capsid proteins from Herpesviridae were detected, but the similarities were not statistically significant.

The N-terminal region of MCPs in Herpesviridae was thought to have the HK97-like fold based on low-resolution cryo-EM (Baker et al., 2005, Hui et al., 2013). Very recently, the atomic structure of the entire MCP of HHV-5 has been determined by cryo-EM to 3.9 Å resolution, revealing that a modified HK97 fold is present in the structure, but it is assembled in space from several distal segments of the peptide chain. This can be explained as a series of insertions of several very large domains into the loops of the HK97 “floor” domain. The relative arrangement of the predicted secondary structure elements along much of the lengths of ORF68 and VP5 homologs is similar (Additional file 3), suggesting structural conservation of multiple domains in these proteins; the structural discontinuity of the HK97 domain helps to understand the difficulties of establishing sequence homology in these regions of herpesvirus MCPs. The conserved region 3 (Fig. 1A) corresponds to the “buttress domain”, which appears to play a role in interactions between MCP subunits within the capsomer, as well as between MCP and “triplexes” that stabilize the capsomers (Yu et al., 2017). A better understanding of the role of conserved sequence elements within the buttress domain shown in Fig. 1A awaits the atomic model of this region.

Another connection to the morphogenesis module of bacteriophages is provided by HaHV-1 ORF73, which appears to belong to the herpesvirus capsid maturation protease family (assemblins), previously found to be homologous to capsid proteases from various viruses of bacteria and archaea (Liu and Mushegian, 2004 and Fig. 1B). Interestingly, the configuration of putative catalytic residues in malacoherpesvirus proteases is reminiscent of phage capsid proteases, in that the most C-terminal member of the predicted catalytic triad is an acidic residue as in phage enzymes, rather than by histidine as in assemblins of Herpesviridae. The C-terminal part of the ORF73 product has a reduced compositional complexity; this is similar to the scaffold protein domains, which in some other herpesviruses are likewise compositionally biased and form a fusion with the protease that self-cleaves following the capsid assembly (Newcomb et al., 2000; Veesler et al., 2014). Here again, sequences from Alloherpesviridae are outliers: we did not recognize any sequence homolog of assemblin among alloherpesvirid genes, including ORF78 product of Cyprinid herpesvirus 3 and its orthologs, for which such role has been assigned based on several lines of circumstantial evidence (Michel et al., 2010, Davison et al., 2013).

Other putative structural proteins were also identified in HaHV-1. ORF14 is homologous to the herpesvirus VP23 family, and ORF81 appears to belong to the VP19C family (Table 1). A heterotrimer (“triplex”) of two molecules of VP19C and VP23 stabilizes the adjacent capsomers of VP5 (Newcomb et al., 2003, Okoye et al., 2006). We detected neither the homologs of the accessory assembly factor VP3 or of the capsid outside-surface protein VP26, nor the components of the capsid portal complex that is required for letting the scaffold fragments out of the capsid and for importing virus DNA. Nonetheless, identification of sequence similarities in three structural components of the capsid and in the capsid maturation protease, together with the catalytic subunit of terminase ATPase detected earlier, suggests that the structural module of malacoherpesviruses is mostly or entirely homologous to other herpesviruses, and the ancient connections to ancestors from viruses of prokaryotes are retained in morphogenetic modules of diverse Herpesvirales.

The identities of tegument proteins as well as virion outer surface proteins in malacoherpesviruses remain to be determined. Interestingly, a pair of paralogous proteins (products of ORF25 and ORF102 in HaHV-1) is found in all malacoherpesviruses. Database searches demonstrated that homologs of ORF25 and ORF102 are found in many animal genomes, including diverse aquatic invertebrates and terrestrial arthropods. In addition, an initial BLASTP search seeded with ORF25 recovered a match (p-value < 10−6) to a putative glycoprotein of Hubei myriapoda virus 8, an RNA virus with a negative-strand RNA genome that has been characterized in a large-scale arthropod metagenomic effort (Shi et al., 2016). Further searches suggest that these proteins are related to structural proteins from several other animal viruses with different types of genome organization. Their distant homologs include spike proteins of positive-strand RNA toroviruses and coronaviruses, as well as baculovirus F proteins (see full list of matches in Additional File 4). All these genes encode virus membrane envelope proteins that facilitate virus contact with the host receptor, virus-cell membrane fusion and virus entry into the cell (Belouzard et al., 2012, Baquero et al., 2015, Palgen et al., 2015). More distant relationships (Additional File 4) link a sub-region within these proteins to the envelope proteins of the Ty3-copia-class retrotransposons in insects (virus family Metaviridae), as well as their domesticated homologs of those proteins, such as the Iris gene product in Drosophila (Malik and Henikoff, 2005). Conservation of this region in integrated copies of metaviruses would explain the matches between ORF25/ORF102 family to other animal genomes.

Non-structural proteins: variations within the herpesvirus enzymatic repertoire

Previous analysis has identified a number of malacoherpesvirus genes that contained sequence signatures strongly suggestive of their role in virus genome replication (Davison et al., 2005). These include two subunits of ribonucleoside reductase (ORF22 and ORF31), putative deoxyuridine 5′-triphosphate nucleotidohydrolase ORF34, two related helicases (most likely, DNA helicases) ORF38 and ORF46, as well as family B DNA polymerase ORF92 and ATP-dependent DNA ligase ORF107.

In addition to these identifications, we predict, on the basis of sequence similarity, molecular functions of two putative nucleases, ORF50 and ORF98, and three additional helicase-related proteins, ORF26, ORF36 and ORF94. Of these, the ORF26 product appears to include only the putative non-catalytic C-terminal domain, whereas ORF36 and ORF94 contain the apparently complete sets of helicase catalytic residues in their N-terminal domains, as well as the C-terminal domain homology. In ORF94, the helicase is fused to the archaeo-eukaryotic-type primase (AEP) domain. One more protein that contains an AEP domain is ORF45; it is adjacent to a helicase, ORF46, on the virus chromosome.

Other genes with possible roles in the replication and expression of the virus genome are two small-molecule kinases ORF105 and ORF111, which could play a role in nucleotide salvage and their delivery to the replication complexes, as has been shown for several phages and viruses (Kim et al., 2005, Jeudy et al., 2006, Guevara-Hernandez et al., 2012). One gene product, ORF105, appears to belong to the herpesvirus UL92 family, which in betaherpesviruses and gammaherpesviruses is essential for the expression of the true-late transcripts encoding the virion proteins (Omoto and Mocarski, 2014). Another gene product, ORF86, is erroneously annotated in the database as a methyltransferase, but appears to be a homolog of the TnpB protein, known to occur in the IS200/IS605 class of transposons in prokaryotes, as well as in various eukaryotic transposable elements and NCLDV genomes (Bao and Jurka, 2013); it contains an N-terminal helix-turn-helix domain and RuvC-like putative nuclease domain with inserted Zn ribbon (Bao and Jurka, 2013, Makarova and Koonin, 2015). Finally, ORF27 and ORF82 each encodes a fusion of an HUH-motif endonuclease and an ATPase-helicase domains. Fusions of those domains are common in proteins involved in replication of plasmids, transposons and small DNA viruses (Koonin et al., 2015). None of the latter three ORFs have homologs in other malacoherpesviruses, and it is unclear whether they are involved in virus genome propagation or represent selfish genetic elements parasitizing virus DNA.

Evolutionary origin of Herpesvirales replication module: sampling from the same gene pool as NCLDV?

Gene sharing between viruses and many groups of virus-related mobile elements have been studied recently using the network analysis techniques, with virus genomes represented as nodes and genes shared between pairs of genomes modeled as weighted edges (Iranzo et al., 2016). A giant module in that network comprises Megavirales – a recently proposed order that joins several families of NCLDV, including Poxviridae, Asfarviridae, Iridoviridae, Ascoviridae, Phycodnaviridae, Mimiviridae and Marseilleviridae (Colson et al., 2013). This order is further linked to several groups of viruses with moderate-sized genomes and to cellular mobile elements. Most viruses in that super-module share subsets of structural and morphogenetic proteins, including a double jelly-roll major capsid protein, whereas the sets of genes involved in DNA replication and transcription differ between modules within the super-module.

Herpesviruses and tailed bacteriophages comprise a separate large module in the same network, held together by their own conserved assembly of morphogenetic proteins. Thus, the structure/morphogenesis modules in the established order Herpesvirales and proposed order Megavirales are, for all practical purposes, unrelated. It is not clear, however, whether the genome propagation modules of viruses in the two orders are related or not.

Among the five hallmark proteins shared by all Megavirales, three play essential roles in virus genome replication and expression: family B DNA polymerase, helicase-primase, and viral late transcription factor 3 (Yutin et al., 2009). The orthologs of the first two gene products are also found in all Herpesvirales. To understand the phylogenetic relationships between these two virus groups, we assessed the phylogenetic position of herpesvirus DNA polymerases. A multiple alignment of viral and cellular Family B polymerases was constructed, and the alignment of 364 amino acid positions was used to infer an unrooted maximum-likelihood tree, shown in Fig. 2 (the source file of the tree in Newick format is available as Additional File 5). The known monophyletic groups of eukaryotic and archaeal DNA polymerases are observed in the tree, and most of the virus taxa approved by ICTV are seen as statistically supported clades, though neither of the two orders, proposed Megavirales and approved Herpesvirales, resolved as monophyletic. Instead, both orders were split by a well-supported (bootstrap value 0.88) internal partition in the tree. On the one side of the split, DNA polymerases of Herpesviridae join an assembly of NCLDV families and, at a longer evolutionary distance, eukaryotic DNA polymerases delta. The other side includes DNA polymerases of Alloherpesviridae and Malacoherpesviridae, together with their poxvirus, baculovirus and bacteriophage homologs and an assembly of various eukaryotic and archaeal polymerases.

Fig. 2

Unrooted maximum-likelihood tree of viral and cellular Family B DNA polymerases. The percentage of bootstrap support is shown for all internal partitions where it exceeds 80%. Monophyletic clades are indicated by different colors. Blue, lineages within proposed order Megavirales; red, families of Herpesvirales; black, other viruses and phages; green, Archaea; gray, Eukarya.

Taking such tree topology at face value, one would have to suggest a very ancient split between Herpesviridae and other Herpesvirales, as well as a similarly ancient division between poxviruses and the rest of NCLDV, with an ad hoc hypothesis required to explain the origin of the eukaryotic polymerases delta. However, the topology of the tree shown in Fig. 2 should be interpreted with caution; viral sequences evolve at a high rate, resulting in long branches in the trees and ensuing topological artifacts, the best-studied of which is long branch attraction (Bergsten, 2005). Indeed, most virus branches in the tree are long, and the branches leading to Alloherpesviridae and Malacoherpesviridae DNA polymerases are the longest ones; together with generally high divergence in the family, this may obscure the identity of their true nearest evolutionary neighbors. An earlier study (Yutin and Koonin, 2012) also showed that in the DNA polymerase tree Megavirales do not resolve as a clade unless the tree topologies are pre-constrained. More recently, another phylogeny of DNA polymerases has been presented, in which Megavirales and Herpesvirales are seen as two monophyletic lineages (Kazlauskas et al., 2016), albeit the sampling of herpesviruses was much narrower in that study than in the current case.

Recently, a comprehensive analysis of viral and cellular primases of the archaeo-eukaryal type (Prim-Pol superfamily), including virus homologs, has been published (Burroughs and Aravind, 2016). The results of the analysis of that second hallmark gene of dsDNA viruses are strikingly similar to our analysis of DNA polymerases detailed here. In particular, Prim-Pol proteins of both Herpesvirales and Megavirales resolved as multiple clades, though in that case Alloherpesviridae and Herpesviridae form a clade to the exclusion of Malacoherpesviridae (refer to Fig. 6 in Burroughs and Aravind, 2016, and online material accompanying their paper). All these virus clades are deep in the tree, and, whether coincidentally or not, the closest tree neighbor of Malacoherpesviridae is again Poxviridae, though their joint clade in the Prim-Pol tree is not well-supported.

Notwithstanding the concerns about tree artifacts, there appears to be a common theme in the phylogenies of the two hallmark families from the proposed order Megavirales and approved order Herpesvirales. In both cases, proteins from viruses that belong to the same order tend not to resolve as a single clade, and in both cases subsets of each order are intermingled. These results are compatible with the notion that the DNA polymerase genes in Megavirales as well as in Herpesvirales are roughly of the same evolutionary age and have been drawn from the same source – probably, from the ancestral pool of phage replication enzymes.

Unexpected tree topologies may be caused not only by the noisy data or methodological artifacts, but also by true past events, such as horizontal gene transfer (HGT) and gene loss, both of which are well-documented in viruses. Whereas no credible scenario involving HGT can be proposed for the two hallmark genes, DNA polymerase and Prim-Pol, a transfer between distant clades may cause unusual evolutionary affinities of several other malacoherpesvirus replication proteins. For example, of the two subunits of HaHV-1 ribonucleoside reductase, the large subunit ORF022 is most similar to eukaryotic homologs, and the small subunit ORF037 is more closely related to phage and bacterial homologs; and HaHV-1 ATP-dependent ligase ORF107 is most similar to the homologs from echinoderms and other invertebrates while lacking homologs in the genome of OsHV-1 (Table 1 and data not shown). Interestingly, exactly these trends, i.e., the pattern of frequent HGTs and repeated loss events, have been observed in the orthologs of those same genes in NCLDV (Yutin and Koonin, 2009, Yutin and Koonin, 2012).

On a separate note, a relatively recent HGT event between a malacoherpesvirus and other viruses is suggested by the presence in HaHV-1 of two copies of the uncharacterized conserved domain DUF1335, which thus far was thought to be restricted to White spot syndrome virus, the sole member of the dsDNA virus family Nimaviridae (Table 1).

Expanded host range and conserved gene core of malacoherpesvirids

In this study, we report a novel genome related of malacoherpesviruses, inadvertently sequenced by the genome project of the annelid C. teleta. This, as well as the earlier discovery of a nearly-complete virus-like DNA scaffold co-deposited with the genome of amphioxus B. floridae (Savin et al., 2010), suggest that herpesviruses may be consorting with diverse marine invertebrates. Putative herpesvirus-like isolated from other aquatic Metazoa, such as two species of stony corals and wild-caught Hydra sp., have been reported (Vega Thurber et al., 2008, Grasis et al., 2014, Correa et al., 2016). Caution has been advised in interpreting these results (Houldcroft and Breuer, 2015; Wood-Charlson et al., 2015; Sweet and Bythell, 2016), and our reanalysis of the three publicly available short sequences from these studies annotated as fragments of anthozoan herpesvirus DNA suggests that they are in fact of bacterial and fungal provenance (Additional File 6). The existence of herpesviruses of coelenterates is nevertheless likely, given that particles with herpesvirus-like morphology have been observed in coral Acropora (Correa et al., 2016). Another recent report of a virus with herpes-like virion morphology in blue king crab Paralithodes platypus (Ryazanova et al., 2015) appears to be corroborated by sequence similarity analysis, which suggests that the DNA fragment specifically amplified from the diseased tissue is closely related to alloherpesvirus DNA polymerase (Additional file 6).

All told, close to 50% of all HaHV-1 virus genes have orthologs in the Ostreid herpesvirus 1 genome, and the majority of these genes also have orthologs in either annelid or amphioxus herpesviruses. Thus, at least 30 genes are conserved within the expanded Malacoherpesviridae family, and their sequence conservation models may serve as this family genetic signature. Two-thirds of these genes represent a variety of essential molecular functions, both structural and informational, and the remaining ones are uncharacterized.

Several gene products with the roles in capsid formation are conserved throughout the family and in Herpesvirales, suggesting that the sequences deposited with the two invertebrate genomes may represent non-defective, particle-forming viruses. To study this possibility further, we tested whether the herpes-like contigs found within assembled genomes of Branchiostoma and Capitella show evidence of relaxed selective constraints. We used the RELAX method (Wertheim et al., 2015), which contrasts two probabilistic evolutionary models, one in which relaxation of selection is allowed in a specific branch and the other in which relaxation is not allowed. The test then compares the log-likelihoods of the two competing models. The results of this analysis, summarized in Table 2, show no evidence of relaxed selection in the Capitella and Branchiostoma virus-like genes. This may suggest that these sequences represent bona fide viruses that were sequenced together with their host. If, on the other hand, this DNA is integrated into the host genome (the possibility supported in the case of B. floridae by the recent deposition in the database of shotgun sequence of the gonad-derived genomic DNA of another cephalochordate, B. belcheri, which also includes a complement of herpesvirus genes, including for example the major capsid protein homolog XP_019628941), the integration event must have taken place recently enough so that a change in selective constraints had no time to leave an imprint on the protein-coding sequences.

Table 2

Testing for relaxation of selection in the herpesvirus genes deposited with the assemblies ofand. For the identity of HaHV-1 orthologs in the annelid and cephalochordate, refer to Table 1. Fourteen statistical tests for relaxation were conducted separately, and the p-values are shown before correcting for multiple testing. After correcting for multiple testing, none of the tests was statistically significant.

Gene in HaHV	Gene function	Capitella	Branchiostoma
ORF098	YkaJ-type endonuclease	0.078	0.67
ORF092	Family B DNA polymerase	0.49	0.12
ORF068	Major capsid protein	0.005	0.49
ORF045	Primase small subunit	0.59	0.3
ORF036	Helicase	0.023	0.55
ORF014	Capsid protein VP23 family	0.41	1
ORF030	Terminase large subunit	0.49	0.21

Improved detection of sequence similarity – as reflected by 33 malacoherpesvirus gene homologs in B. floridae detected here, as compared to only 19 reported in an earlier study (Savin et al., 2010) – provides a new lower bound of the size of the common gene complement in herpesviruses. We expect that re-sequencing of virus-like replicons associated with the amphioxus and annelid genomes will reveal more gene products shared with other herpesviruses. Based on these data and the evidence in the literature, the fraction of widely conserved genes (in the first approximation, having homologs beyond the family level) in herpesviruses may be at least 30% and perhaps closer to 50%, similarly to what has been suggested for bacteriophages with similarly-sized DNA genomes (Kristensen et al., 2011, Grazziotin et al., 2017). Given the observations of the putative herpesviruses associated with diverse marine invertebrates, the family name Malacoherpesviridae may be worth replacing with a more inclusive one, perhaps “Ectoherpesviridae”.

Concluding remarks

Recent studies have provided important advances in our understanding of early evolutionary history of eukaryotic viruses. A general principle of virus evolution has emerged, i.e., that the hallmark genes in major groups of viruses of eukaryotes are related, sometimes distantly, to phage genes with similar functions, with additional contribution of genes acquired from mobile elements and host genomes. The genome replication gene modules in many groups of viruses are derived from different phage sources than the structural genes. Putative order Megavirales, all RNA viruses, and certain classes of DNA viruses with smaller genomes represent distinct variations on this general theme (Koonin et al., 2015; Iranzo et al., 2017).

Here we see the same trend in the evolution of the order Herpesvirales, where the morphogenesis and replication gene sets have distinct evolutionary histories. A core set of morphogenetic genes appears to be conserved in at least two families of Herpesvirales and is more distantly related to homologs in phages with HK97-type capsids. In contrast, several essential herpesvirus genome replication factors, such as Family B DNA polymerase and DNA Prim-Pol primase, appear to have an ancient origin and possible ancestral relationship to the hallmark genes of Megavirales, whose structural module is distinct from that of herpesviruses.

Taken together, this evidence indicates that the ancestors of two orders, Herpesvirales and putative Megavirales, may have recruited their genome propagation module from the same source – most likely, from the replication genes in the community of phages and prokaryotic mobile elements. These genes may have been confined to each lineage of eukaryotic viruses ever since; unlike several other enzymes with roles in DNA synthesis in the same virus orders, these essential replication genes do not appear to have been horizontally transferred between these viruses and their hosts, except for occasional integration of the whole virus DNA into the host genome.

Also at an early stage, the two emerging virus lineages have acquired distinct morphogenetic modules, i.e., the jelly-roll type in Megavirales and the HK97 type in Herpesvirales. Curiously, the capsid proteins of the former type are relatively narrowly distributed in the known phages, but are found, either in single jelly-roll or in an internally duplicated form, in eukaryotic viruses that infect nearly all major clades of eukaryotes (Krupovic and Koonin, 2017). The latter type of capsid displays the opposite trend, i.e., it is found in many groups of bacterial and archaeal viruses, but among eukaryotic viruses it has been identified in just one order infecting Metazoa. The reasons for these asymmetries and for the association of a specific capsid type with nuclear vs. nucleocytoplasmic replication cycle are unknown. It is quite likely, however, that better sampling of viruses from diverse hosts and from the environment will improve virus phylogeny and provide better understanding of the emergence of the two major lineages of eukaryotic large dsDNA viruses with their distinct replication strategies.