Sequencing of bovine herpesvirus 4 v.test strain reveals important genome features

Background Bovine herpesvirus 4 (BoHV-4) is a useful model for the human pathogenic gammaherpesviruses Epstein-Barr virus and Kaposi's Sarcoma-associated Herpesvirus. Although genome manipulations of this virus have been greatly facilitated by the cloning of the BoHV-4 V.test strain as a Bacterial Artificial Chromosome (BAC), the lack of a complete genome sequence for this strain limits its experimental use. Methods In this study, we have determined the complete sequence of BoHV-4 V.test strain by a pyrosequencing approach. Results The long unique coding region (LUR) consists of 108,241 bp encoding at least 79 open reading frames and is flanked by several polyrepetitive DNA units (prDNA). As previously suggested, we showed that the prDNA unit located at the left prDNA-LUR junction (prDNA-G) differs from the other prDNA units (prDNA-inner). Namely, the prDNA-G unit lacks the conserved pac-2 cleavage and packaging signal in its right terminal region. Based on the mechanisms of cleavage and packaging of herpesvirus genomes, this feature implies that only genomes bearing left and right end prDNA units are encapsulated into virions. Conclusions In this study, we have determined the complete genome sequence of the BAC-cloned BoHV-4 V.test strain and identified genome organization features that could be important in other herpesviruses.


Background
Gammaherpesviruses are archetypal persistent viruses which are ubiquitous in both human and animal populations. The human gammaherpesviruses, Epstein-Barr virus (EBV) and Kaposi's Sarcoma-associated Herpesvirus (KSHV), infect respectively some 90% [1] and 30% [2] of human populations and cause several cancers [2,3]. Although much effort has been invested on these viruses, studies of EBV or KSHV are difficult to perform directly because these viruses show limited lytic growth in vitro and have no well-established in vivo infection model. Related animal gammaherpesviruses are therefore an important source of information.
Bovine herpesvirus 4 (BoHV-4) belongs to the Gammaherpesvirinae subfamily, and to the Rhadinovirus genus [4]. Similarly to its human counterparts, BoHV-4 was found to be widespread in all bovine populations and to persist in the vast majority of individuals as a lifelong, asymptomatic infection [5]. However, in some circumstances, BoHV-4 has been associated with various clinical symptoms such as skin lesions, respiratory diseases, metritis, malignant catarrhal fever or tumors [5].
The virus was first isolated in Europe by Bartha et al. from calves with respiratory diseases [6] and later in North America by Mohanty et al. [7]. Besides cattle, BoHV-4 has also been detected in a variety of ruminants. In particular, BoHV-4 seems to be highly prevalent among wild African buffalo (Syncerus caffer) which could be considered as the natural reservoir of the virus [8][9][10]. Overall, more than 40 BoHV-4 strains have been isolated across the world. These strains can be classified in three groups: the European strains (or Movar 33/63like strains), the American strains (or DN 599-like strains) and the African buffalo strains [9].
It is estimated that the taurine and buffalo strains diverged around 730,000 years ago [9] and that the European and North American clades diverged around 260,000 years ago [9]. The genome of the BoHV-4 66-p-347 North American strain has entirely been sequenced [11]. However, the BAC-cloned reference strain V.test [12,13] belongs to the European clade [9,14]. Previous studies suggested that the BoHV-4 V-test strain contains regions of high dissimilarity compared to the BoHV-4 66p-347 strain. Indeed, the nucleotide identity between the two strains has been previously measured to be as low as 88% on the BORFB2 region [11]. However, the lack of a complete genomic sequence for the V.test strain prevents from drawing a general view concerning this divergence level. Therefore, the low quality of the genomic information hampers the use of the BAC-cloned BoHV-4 V.test strain as a good model for studying gammaherpesvirus biology. In this study, we have determined the genomic sequence of the BoHV-4 V.test strain and analyzed its overall differences with the available sequence of the BoHV-4 66-p-347 strain [11,15]. The results obtained highlighted important differences between BoHV-4 66-p-347 and V.test strains. Moreover complete sequencing of the BoHV-4 V.test strain also revealed genome features potentially important in other herpesviruses.

BAC sequencing
BAC DNA was purified using Qiagen large-construct kit as described by the manufacturer. The complete BAC cloned viral genome of BoHV-4 V.test strain was determined by pyrosequencing using the 454 GS FLX Titanium (Roche) high-throughput sequencer and resulted in 48,967 reads of an average read length of 265 nucleotides and a total of 12,997,275 bases. A targeted ABI-Sanger sequencing of fragments of the prDNA region was also conducted using the primers listed in Table 1. The raw 454 data has been deposited in the NCBI Sequence Read Archive (SRA) database with accession number SRA037246.

BoHV-4 genome LUR assembly
The reads were de novo assembled with gsAssembler (Roche), where the E. coli genome was used as a contaminant to filter out cellular reads [16]. The filtering removed 1,167 contaminant cellular reads. The de novo assembly yielded 11 contigs which were subsequently BLASTed against 66-p-347's long unique region (LUR) and polyrepetitive DNA (prDNA) -accession numbers NC_002665 and AF092919-to define their relative positions [17]. Contigs were assembled into a large scaffold using two previously published V.test sequences (accession numbers Z46380 and Z46385 [18]) overlapping contig borders. A careful comparison of the bordering contigs with the previously sequenced fragments showed a high percent identity (> 99.99%). After verification of the quality of the assembly, the BAC sequence was removed and the genome sequence was annotated as detailed hereunder.

BoHV-4 genome prDNA assembly
The prDNA was determined by a hybrid 454/ABI-Sanger strategy where 17 ABI-Sanger fragments of prDNA were de novo assembled with the 454 reads. Briefly, in order to correctly assemble the prDNA and to disentangle different prDNA units, this second de novo assembly was optimized for highly repetitive segments using MIRA [19]. 454 reads and quality information were extracted from the raw .sff file with 'sff_extract'. The base-calling and quality-calling for Sanger sequences were inferred from the .ab1 raw chromatogram files using 'phred' [20,21] and the sequences were qualitytrimmed using 'lucy' [22]. MIRA assembler (v 3.2.0) was used to build an assembly of the V.test genome with the following flags and options: "-job = denovo, genome, accurate, sanger, 454 -highlyrepetitive -AS:klrs = no 454_SETTINGS -AS:urdcm = 1.1:ardml = 100''. This assembly yielded a very large contig containing a complete prDNA unit, and a second contig containing an incomplete unit bearing the prDNA/prDNA junction. The complete prDNA unit was extracted from the first contig and identified as being the last prDNA unit before the LUR junction and noted prDNA-G following Bublot et al. [14]. By analysing the contig bearing the prDNA/prDNA junction in GAP4 [23], we determined a 518 bp fragment of the prDNA-inner unit (as named by Bublot et al. [14]) bordered on the left by lower read qualities and coverage, and on the right by the beginning of a new prDNA unit. This end was joined to the beginning (2,089 bp) of the prDNA-G unit in order to obtain a complete prDNA-inner unit (2,607 bp). We verified that this complete unit was compatible with previously published information [14].

BoHV-4 genome annotation
All Open Reading Frames (ORFs) from all 6 frames were retrieved from the complete genomic sequence and matched against the Conserved Domain Database [24] using the position-specific scoring matrices (PSSM) based Reverse PSI-BLAST [25]. For all ORFs sharing the same STOP and containing a PSSM match, the smallest ORF containing the largest PSSM match was retained. 59 ORFs were thus considered evolutionarily conserved and were annotated with the corresponding matching conserved domains. Out of the 79 CDS from the previously published 66-p-347 strain, all 59 ORFs matched previously annotated 66-p-347 ORFs. The 20 remaining CDS were added by similarity to this strain and were annotated as such. Repeat segments and special features were annotated according to 66-p-347 if they were present in V.test. The complete genome sequence containing the LUR, prDNA-G and prDNA-inner were annotated and submitted to GenBank with respective accession numbers: JN133502, JN133503 and JN133504.
Comparative genomics analysis of 66-p-347 and V.test The LUR and prDNA sequences of the 66-p-347 strain were joined into a complete genome (accession numbers NC_002665 and AF092919) and aligned against the joined LUR and prDNA-inner V.test sequences with ClustalW 2.0.10 [26]. Percent divergence, percent insertions and deletions, and percent G+C content were computed (i) along the alignment on a 100 bp sliding window of step 3 bp and (ii) on all individually aligned proteins. Analyses and figures were conducted using R [27] and the seqinr [28] package in combination with ad hoc programs written in Python and using the Biopython libraries [29,30].

RT-PCR analysis
These experiments were performed as described elsewhere [31]. Briefly, subconfluent monolayers of MDBK cells were infected with BoHV4 V.test strain at a m.o.i. of 1 PFU/cell. 18 hours after infection, cytoplasmic RNA was extracted, purified and treated for RT-PCR. The cDNA products were amplified by PCR using specific primers listed in Table 1.

BAC sequencing and genome assembly
Pyrosequencing of herpesviral genomes is often limited by the high concentration of contaminating cellular DNA [32]. We therefore prepared the BoHV-4 V.test strain DNA from BAC maintained genomes and sequenced it using a high-throughput pyrosequencing approach [16]. This yielded 48,967 reads among which 47,800 were BoHV-4 specific (> 97% of the reads). After assembly, the mean genome coverage was of the order of 96×. In comparison to the whole genome sequencing of another herpesvirus based on DNA isolated from virus particles, which exhibited a 13× average base pair coverage [32], our strategy showed a more than 7-fold increase. This is probably mainly due to the high proportion of viral to cellular reads present in our dataset. Indeed, only 1,167 Escherichia coli contaminant reads had to be discarded from the data, indicating less than 2.38% of contaminated reads, compared to the previously reported 62.72% contaminating cellular reads in [32]. Our sequencing strategy based on a BAC cloning approach, thus revealed itself very powerful in terms of contamination and subsequent coverage.

V.test genome analysis and comparison to other BoHV-4 strains
The BoHV-4 genome has a B-type structure consisting of a long unique region (LUR) flanked by several polyrepetitive DNA units (prDNA). We assembled the complete LUR of the V.test strain BoHV-4 genome into a 108,241 bp sequence. The average G+C content is of 41.21%. This value as well as the G+C% variation observed on Figure 1 is in agreement with previously reported results on the 66-p-347 strain, namely on the high G+C content of R2a region [11]. The observed-toexpected CpG ratio is of 0.225 on the LUR and is compatible with the value measured on Bos taurus (0.234) [33] suggesting (i) a high degree of methylation of CpG nucleotides and (ii) similar methylation mechanisms acting on the viral and cellular genome. As expected, the nucleotide identity between our assembled genome and previously published V.test strain sequence data was of 99.55% in average, falling into the ranges of comparison between 454 and Sanger sequencing [34].
Compared to the 66-p-347 strain, the V.test strain had previously shown divergence up to 12% on the region surrounding BORFB2 (ORF 16, v-Bcl-2) [11]. However, the lack of a complete genomic sequence for the V.test strain prevented from drawing a general conclusion concerning this divergence level. Compared to 66-p-347 strain, the overall V.test nucleotide identity is high (99.1%), but shows a large variability at the genome level ( Figure 1). As expected, the repetitive regions contained in the LUR (R1, R2a and R2b) exhibit a high nucleotide divergence, up to more than 40%, as well as large gaps ( Figure 1). This indicates that the very high divergence levels seem confined to specific repetitive genomic regions. However, some rather high divergence levels were also identified in other regions ( Figure 1) and namely in ORF-containing regions such as ORF 10, Bo5, ORF 57, and ORF 68 region. We also note a Figure 1 Map of the BoHV-4 V.test strain genome and divergence with the 66-p-347 strain sequence. The LUR of both strains have been aligned. Genome features are represented in the upper part as grey and red oriented arrows. Red arrows represent genes with an in-frame STOP codon, an early Methionine or a high divergence level in the V.test strain compared to 66-p-347. Dark (resp. light) grey arrows represent genes with (resp. without) an evolutionarily conserved domain (see Methods). The exons of spliced genes are indicated under the given gene as thin light-blue lines. Percent divergence is shown as a black-filled curve, percent insertions and deletions are shown as a blue-filled curve. Percent G +C content is shown as a thin green curve, with the mean G+C content drawn as a thin horizontal green line. These percentages are measured in a 100 bp window sliding 3 bp. Repeat regions (R1, R2a, R2b) are depicted as hatched areas. The oriLyt region is mapped as a light-grey area within R2b and the conserved quasi-palindromic motif in the oriLyt region is indicated by a small vertical arrow. large deletion and a high divergence at the beginning of the LUR compared to the 66-p-347 strain. Overall, these differences in protein-coding region as well as in repetitive regions that bear predicted microRNA coding sequences [35] will require specific experiments to identify possible links with observed phenotypic differences between strains.

Conserved protein-coding genes
In order to develop an ab initio approach of gene annotation, we extracted all possible ORFs in all 6 frames from the complete genomic sequence of the BoHV-4 V. test strain. On each of these ORFs, we ran a Reverse PSI-BLAST [25] against all protein domains from the Conserved Domain Database [24]. ORFs containing an evolutionarily conserved domain were defined as the smallest ORF containing the longest CDD match (see Methods). This approach revealed 59 ORFs containing a conserved CDD domain (  [36] such as the BoHV-4 Bo17 gene that encodes a homologue of the cellular core 2 beta -1,6-N-acetylglucosaminyl-transferase M [37]. These results will deserve further studies to identify the evolutionary history responsible for these observations.

Non-conserved protein-coding genes
The remaining 20 annotated ORFs were determined by similarity to the 66-p-347 strain, and correspond for most of them to ORFs unique to BoHV-4 as described previously (Table 3) [11]. Some of these ORFs, however, contain odd characteristics that needed to be investigated (Figure 2, Additional file 1 Figures S1-9). Indeed Bo1, Bo6, Bo7, Bo12 and Bo13 genes of the BoHV-4 V.test strain present inframe STOP codons. Bo5 presents rather high divergency levels and large insertions/deletions (> 5% of its coding sequence as shown in Figure 2) compared to the genomic sequence of the 66-p-347 strain. Moreover, ORFs 36, 67.5 and 75, which bear an evolutionary conserved domain, present late methionines compared to the 66-p-347 annotation. Indeed, in ORF 36 (see Additional file 1 Figures S5), the smallest ORF containing an evolutionary conserved domain is slightly shorter than the one annotated in 66-p-347 and there is no evidence that the previously annotated   Figures S9), leading to the absence of the 66-p-347 annotated methionine in the V.test strain. All these annotated genes requested therefore an investigation of their transcription in mRNA products. As these sequence properties could be specific to the BAC clone of the BoHV-4 V.test strain, we investigated the transcription of these genes on MDBK cells infected with the BoHV-4 V.test WT strain as described in the methods. The primers used are described in Table 1 and highlighted in Additional file 1. For all couple of primers, cDNA from BoHV-4-infected MDBK cells gave rise to the expected PCR products (Figure 3). The absence of contaminant viral DNA in the mRNA preparations was confirmed by a lack of PCR product without reverse transcriptase. The size of the Bo5 RT-PCR product was also consistent with its known mRNA splicing (868 bp rather than 1140 bp). Moreover, the sequences of these RT-PCR products were in agreement with the BoHV-4 V.test sequence derived from our BAC cloned genome (data not shown). Therefore, we can conclude that all these coding sequences are transcribed during BoHV-4 infection of MDBK cells. However, further investigation is needed to determine the presence of proteins and ensure their accurate annotation.

BoHV-4 V.test replication origin
A large region containing the potential lytic replication origin (oriLyt) of the BoHV-4 66-p-347 strain was  determined by Zimmermann et al [11]. Based on this information, we mapped this region on the V.test genome ( Figure 1). This region contains Bo12, the R2b region and partially overlaps with Bo11. Compared to the 66-p-347 strain sequence, the corresponding region in the V.test genome is highly divergent (Figure 4). Although this region shows high divergence rates, we expected the replication origin to be conserved between the two BoHV-4 strains. Previous work on other herpesviruses has identified in oriLyt the presence of palindromic motifs essential for viral replication [38][39][40]. When we compared the potential region containing oriLyt in the two strains, a single conserved palindromic region was observed (AATCCAGGCCCCTGATTGGTA-GATTGC TGAAAGCCAATCAGGGGCCTGGATT, Figure 4). Interestingly, this region forms a perfect hairpin structure ( Figure 4B) that resembles DNA structures formed at other herpesvirus origins [41,42] and may therefore represent a common secondary structure used by all herpesvirus family members during the initiation of DNA replication. In the future, this structure will be tested as a candidate for an essential oriLyt replication motif.

BoHV-4 V.test polyrepetitive DNA
In the BAC clone, previous restriction profiles had determined a hypermolar prDNA band indicating that the BAC contained several prDNA units [12]. Therefore, the major pitfall in the assembly of the BoHV-4 V.test strain was the determination of the prDNA sequence. Indeed, (i) the higher per base coverage on this region due to repetition of prDNA units, (ii) the high GC content, along with (iii) the presence of several long repeats within the prDNA and (iv) the variability observed between prDNA units [14] made it extremely difficult to resolve and assemble with pyrosequencing data alone. Interestingly, it has been shown for several rhadinoviruses that the left junction between the prDNA and the LUR is the site of genome rearrangements and that sequences of the prDNA are found within the first base pairs of the LUR. These properties make this region very difficult to sequence [43][44][45][46][47]. Therefore, we adopted a hybrid strategy consisting in adding some ABI-Sanger reads (with the primers described in Table 1) to guide the 454 assembly on the prDNA region (see Methods).
Bublot, et al. [14] described the different prDNA unit variants present in BoHV-4 V.test, and namely the differences between prDNA units. Firstly, the prDNA units vary according to the number of repetitions of a~200 bp Pst-I bordered fragment. Secondly, the last prDNA before the prDNA/LUR junction (prDNA-G) displays a different ending than the inner prDNA units [14]. Our method allowed us to disentangle the repeats and to assemble a contig containing a whole prDNA unit (2,440 bp) along with the left prDNA-LUR junction. This prDNA unit, corresponding to prDNA-G following Bublot et al. [14], was extracted from the contig and annotated. A second contig from this hybrid assembly yielded the prDNA/prDNA junction. The presence of the prDNA/prDNA junction in our assembly confirmed the presence of at least two prDNA units in our BAC clone and allowed us to build a complete prDNA-inner unit (see Methods). The assembled prDNA-G and inner prDNA units have sizes of 2,440 bp and 2,607 bp respectively. Both these units are in agreement with their previously published restriction maps [14].
Specifically, we showed that, comparatively to the 66p-347 strain, the V.test prDNA-inner unit presents several indels including two large indels in the repetitive PstI region ( Figure 5). This PstI-rich repetitive region seems to be the one presenting the most variation as it also presents comparatively large differences between prDNA units within the same strain. Indeed, Bublot et al. [14] roughly determined the size of the V.test major prDNA-inner unit to be around 2,650 bp due to the presence of 4 repetitions of the two small PstI bordered fragments.  In the prDNA-G unit, we established that these two small PstI-bordered fragments make up a fragment of 186 bp and that these are indeed repeated 4 times ( Figure 6). In the prDNA-inner unit, we determined that the last PstI-bordered fragment is actually a variation of the 186 bp fragment where the inner Pst-I site is slightly modified ( Figure 6). Therefore, the rough 200 bp size discrepancy between the prDNA-G (2,440 bp) and the prDNA-inner units (2,607 bp) is due to the presence of a slightly modified repetition of the previous segment. These results are compatible with the restriction profiles presented in Bublot et al. [14] as detailed by the positions of several restriction sites on Figure 6.
In addition to the variations in the PstI-bordered repetitions, one of the major differences between the prDNA-inner units and the prDNA-G lies in their 5' end. Indeed, the prDNA-inner contains a conserved pac-2 cleavage/packaging signal in its right terminal region, which is not the case of prDNA-G ( Figure 6). Both units however, possess a conserved pac-1 cleavage/packaging signal in their left terminal region. Interestingly, the pac-1 and pac-2 cleavage and packaging signals show a good conservation between 66-p-347 and V.test's inner units, despite the presence of these signals in a repeated region bearing high divergence levels. Broll et al. [15] have determined, by transient cleavage/packaging assay, that a single prDNA unit is sufficient for cleavage and packaging. However, from the absence of a conserved pac-2 motif in the prDNA-G, we suggest that, even if a single inner prDNA unit is indeed sufficient for cleavage and packaging, the prDNA-G alone would not suffice. This would therefore indicate that two prDNA units at least are necessary in the context of naturally occurring BoHV-4 genomes for correct cleavage and packaging. The packaging of herpesvirus genomes is still not fully understood, however, detailed studies in herpes simplex virus type 1 (HSV-1), human and murine cytomegaloviruses (HCMV and MCMV) have highlighted the roles of the major conserved motifs and suggested the following general mechanism by which concatemers are cleaved and packaged [48][49][50]. Firstly, the T-box of the pac-2 signal is essential for the cleavage that initiates DNA packaging. Cleavage occurs at a fixed distance from the pac-2 T-box, and the resulting end that contains the pac-2 GC-box and other cis acting elements is inserted into the procapsid. Packaging is therefore directional and proceeds from pac-2 towards the pac-1 terminus [48]. A second cleavage event, directed by pac-1, then terminates DNA packaging. If we apply this model to BoHV-4, the divergence of the pac-2 signal in prDNA-G, namely the absence of a T-box, indicates that it is not a functional pac-2 initiation signal. As the genome packaging is directional from pac-2 to pac-1 (therefore, from the right to the left end of the genome), the lack of a pac-2 initiation signal in prDNA-G ensures that no packaging would lead to a remaining concatemer lacking a left end prDNA. This would therefore guarantee that genomes bearing at least one left and one right end prDNA unit are encapsulated into virions. This model and its implications will require further investigations in the future.

Conclusions
BAC-cloning of the BoHV-4 V.test strain has greatly facilitated the use of this virus as a model for human pathogenic gammaherpesviruses. However, until now, the complete genome sequence of this strain was unavailable. In this study, we have determined the complete genome sequence of the BoHV-4 V.test strain. In comparison with the previously sequenced 66-p-347 strain, we identified important differences in 9 potential open reading frames. Moreover, sequence analyses allowed us to identify genome features that are potentially important for viral replication. All together, these results should have implications for the study of BoHV-4 and herpesviruses in general. Figure 6 The prDNA-G unit does not present a complete pac-2 cleavage/packaging signal. Alignment of the prDNA units from V.test strain (prDNA-inner above and prDNA-G below). The differences observed in the alignment are highlighted in grey. The cleavage/packaging signals pac-1 and pac-2 are represented in boxes, and their composing C-rich, G-rich, GC-rich and T-rich units are indicated. PstI, EcoRI, SstII, BamHI restriction sites (here PstI) are represented in coloured font. One modified PstI restriction site in the prDNA-inner is also highlighted to indicate the divergence between the fragments composing both units.

Additional material
Additional file 1: Alignments of the nucleotide and predicted amino acid sequences of Bo1 ( Figure S1), Bo5 ( Figure S2), Bo6 ( Figure S3), Bo7 ( Figure S4), ORF36 ( Figure S5), ORF67.5 ( Figure S6), Bo12 ( Figure S7), Bo13 ( Figure S8) and ORF75 ( Figure S9) of BoHV-4 V. test and 66-p-347 strains. Nucleotide sequences aligned at the aminoacid level are represented for BoHV-4 V.test (red) and 66-p-347 strains (blue). Mismatching residues are highlighted in a shaded grey box. The predicted amino-acid sequences are respectively drawn for V.test and for 66-p-347 above and below the nucleotide sequences. The STOP codons are highlighted by small colored boxes. The annotated Methionine are highlighted in bold font. In the Bo5 sequence, introns are represented by boxes. Positions of the specific primers used in Figure 3 are underlined.