Genomic evolution, recombination, and inter-strain diversity of chelonid alphaherpesvirus 5 from Florida and Hawaii green sea turtles with fibropapillomatosis

DNA isolation, amplification and sequencing

Biopsies of skin tumors were obtained from green sea turtles that had stranded with fibropapillomatosis (Table 1). For HI samples, turtles found dead or humanely euthanized moribund individuals are classified as diagnostic specimens and are considered exempt according to the Animal Welfare Act. Therefore, the Chair of the US Geological Survey, National Wildlife Health Center Animal Care and Use Committee deemed it unnecessary to review or approve sampling. Florida tumor samples were collected according to Marine Turtle Permit MTP-16-233: Breitbart lab under protocol T-IS0001253 approved by the Institutional Animal Care and Use Committee, Research Integrity and Compliance, University of South Florida. A consent permit was issued from the Florida Fish and Wildlife Conservation Commission to transport marine turtle specimens out of the State.

Sections of tumors were excised, placed in cryovials, and stored at −70 °C. DNA was extracted from approximately 10 mg of tumor tissue using DNeasy Blood and Tissue Kit (Qiagen, Valencia, CA, USA). DNA, quantified with the Qubit double-stranded (ds) DNA High Sensitivity (HS) Assay Kit (ThermoFisher Scientific, Grand Isle, NY, USA), were then normalized to 20–45 ng/µl, with one or both concentrations used for LR PCR.

The reference ChHV5 genome (sequenced by Ackermann et al., 2012; GenBank HQ870327) used for primer design in this study exhibited the typical gene order for a Class D genome of the Scutavirus genus within the subfamily Alphaherpesvirinae. Overall, this genome could be divided into unique long sequence (UL; 101,152 bp), unique short sequence (US; 13,319 bp) and inverted repeats that flanked the US (IRS; 8,831 bp each) (Thiry et al., 2005; Ackermann et al., 2012). In this study, LR PCR primers were designed from the ChHV5 reference genome using Geneious v 6.1.8 (Biomatters Ltd.) to amplify overlapping fragments of approximately 5 kb (kilobase pairs). Given difficulties in culturing and purifying the virus, using targeted ChHV5 LR PCR primers was important in order to avoid host DNA, which would overwhelm a shotgun sequencing approach, especially in the case of latent viruses present at low copy numbers. Primer sequences and their genomic coordinates are provided in Table S1.

Long-range PCRs were initially performed using the Qiagen LongRange PCR Kit (Qiagen, Valencia, CA) following the manufacturer protocol. Reactions containing 0.4 µM of each primer, 500 µM of dNTPs, 1U of LongRange PCR Enzyme Mix, 1X LongRange PCR Buffer, and 1 µl of 45 ng/ul DNA template in a 25 µl reaction were run on a Bio-Rad T100™ thermal cycler as follows: 93 °C for 3 min, followed by 35 cycles of (93 °C for 15 s, 62 °C for 30 s, and 68 °C for 5 min). For LR PCR primer pairs that failed to amplify in a few samples, the TaKaRa LA Taq polymerase kit (ClonTech, Palo Alto, CA), was used. Reactions containing 0.2 µM of each primer, 200 µM of dNTPs, 0.75 U of PrimeSTAR GXL DNA Polymerase, 1X PrimeSTAR GXL Buffer, and 2.5 µl of DNA template in a 25 µl reaction, were run on a Bio-Rad T100™ thermal cycler as follows: 93 °C for 3 min, followed by 30 cycles of (98 °C for 10 s, 60 °C for 10 s, and 68°for 5 min). PCRs that did not successfully amplify using the standard TaKaRa LA protocol were optimized according to the manufacturer’s troubleshooting suggestions, which included decreasing the template concentration from 2 ng/µl to 0.5 ng/µl, plus lengthening the last two phases of thermocycling to 60 °C for 15 s, and 68°for 9 min. LR PCR products were verified by gel electrophoresis. For products with multiple bands resulting from nonspecific priming, the band of the expected size was excised from the gel and purified using the UltraClean^® GelSpin^® DNA Extraction Kit (MoBio, Carlsbad, CA, USA). PCR products were quantified using the Qubit^® double-stranded (ds) DNA High Sensitivity (HS) Assay Kit (ThermoFisher Scientific, Grand Isle, NY, USA), then normalized to 1 ng/µl. Amplicons from each sample were then pooled in equal volumes, barcoded, and sequenced simultaneously. Multiple sequencing runs were performed as we explored different platforms and included additional samples (Table S2).

Two samples (FL_F5 and HI_21553) were run separately on an Ion Torrent Personal Genome Machine (PGM™; Life Technologies, Foster City, CA, USA). For each sample, 20 LR PCR products were analyzed separately with the Bioanalyzer 2100 (Agilent Technologies, Santa Clara, CA, USA) using the DNA1000 reagent kit. For each sample, PCR products were pooled in equimolar concentrations and frozen at −20 °C. Libraries were made from pooled PCR products for each sample following the manufacturer’s protocol for the Ion Xpress Plus Fragment Library Kit (Publication # 4471989, Revision N). Briefly, for both libraries, 50 ng of DNA was fragmented using the Ion Shear™ Plus reagents and purified with Agencourt^® AMPure Beads (Beckman Coulter, Indianapolis, IN, USA). Fragment sizes were determined using the Bioanalyzer 2100 (Agilent Technologies, Clara, CA, USA) with the HS DNA kit. The Ion Torrent adapters were ligated using DNA ligase, nick-translated and purified with Agencourt^® AMPure Beads (Beckman Coulter, Indianapolis, IN, USA). For the FL_F5 and HI_21553 libraries, 200 and 100 base read libraries, respectively, were size selected with the E-gel^® SizeSelect™ 2% gel system (http://www.invitrogen.com). Both libraries were quantified using quantitative PCR (qPCR) using the Ion Library TaqMan^® Quantitation Kit (Life Technologies). Within 24 h, the libraries were diluted and placed on the Ion OneTouch™ system (Life Technologies) for template preparation for the emulsion and enrichment of OT2 200 Ion Sphere™ particles (ISP) onto each library. The recovered template-positive Ion PGM template OT2 200 ISPs and library were cleaned with the OneTouch™ emulsion system (Life Technologies) following manufacturer protocols. Enrichment of the ISP was confirmed using the Qubit^® 2.0 fluorometric analysis and the IonSphere™ Quality Control Kit. The final libraries were prepared for sequencing using the Ion PGM™ Sequencing 200 Kit v2 (Publication no. MAN0007273) using the Ion 318 v2 and Ion 316 semiconductor chips for the FL and HI samples, respectively.

For Illumina MiSeq^® runs, cleaned PCR-products were quantified using the Qubit^® double-stranded (ds)DNA HS Assay Kit (ThermoFisher Scientific, Grand Isle, NY, USA) and normalized to 0.2 ng/µl using 10 mM Tris, pH 8.5. Libraries were prepared using 1 ng total input DNA and the Nextera^® XT DNA Library Prep Kit following the Illumina Library Preparation Guide (Illumina, San Diego, CA, USA). The only modification to the protocol was to the PCR clean-up step, in which we used 0.5x Agencourt^® AMPure XP beads (Beckman Coulter, Indianapolis, IN, USA) to select for a larger insert size. Final library size was validated using the BioAnalyzer HS DNA Kit (Agilent Technologies, Santa Clara, CA, USA). Prepared libraries were quantified using the Qubit^® dsDNA HS Assay Kit (ThermoFisher Scientific, Grand Isle, NY) and normalized to 4 nM using 10 mM Tris, pH 8.5. The pooled libraries were run at 15 pM with a 5% PhiX174 spike on the Illumina MiSeq^® with a 600 cycle v3 reagent kit (Illumina, San Diego, CA, USA).

Data analysis

Genomic re-sequencing using long-range PCR and high-throughput sequencing

The amplification strategy was applied to three tumor samples from FL and five from HI (Table 1). Thirty-three of 40 LR PCR primer sets (82.5%) amplified a product of the expected length from at least one sample (Table S2), resulting in the resequencing of ∼73 and 57% of the FL and HI genomes respectively. Variable recovery of amplification products can be seen in the uneven depth of coverage across the genome (Fig. 1), such that not all ORFs were completely recovered for all samples. Regardless, over 75% of predicted protein coding genes were recovered from at least two samples: 7 of 11 ORFs in the US region, and 59 of 76 ORFs from the UL region of the alphaherpesvirus genome (Fig. 1, Table 2, see Morrison et al., 2018 for metadata and sequence alignments). At least two unique consensus sequences were found for every gene, with a maximum of five (F-US1, F-US3A, HP37 and HP38) and an average of 3.2 distinct sequences per gene (Table 2, Unique). For some primer pairs, (e.g., ChHV5-LSC-02 and 03, HI/FL, Table S1), we saw consistent geographic differences in primer success, suggesting sequence or structural polymorphisms between the two clades that were not ascertainable using this strategy.

Genome-wide single nucleotide polymorphisms confirm broad divergence of Atlantic and Pacific strains

Most sequence variation was between regional ChHV5 strains rather than within a location, as expected based on previous phylogenetic analyses (Herbst et al., 2004; Greenblatt et al., 2005a; Greenblatt et al., 2005b; Ene et al., 2005; Patrício et al., 2012; Rodenbusch et al., 2012; Alfaro-Núňez et al., 2014). A total of 1,001 fixed differences was observed between Hawaiian and Floridian samples across the 66 sequenced proteins. All variable sites in over half of examined genes (35/66) were fixed differences between the two geographic strains (Table 3: Fixed), while few variable sites were detected within strains (Table 2: π FL, π HI). ORFs that exhibited the highest number of fixed differences included those of a DNA packing protein (F-UL17), a major capsid protein (F-UL19), a single-stranded binding protein (F-UL29), DNA polymerase (F-UL30), and a hypothetical protein (HP) with DNA- and RNA- binding domains (HP37; Table 3; Fixed).

The average number of nucleotide differences per site between geographic groups (Dxy; Nei, 1987) varied across ORFs, from 0.002 (F-UL04 nuclear protein and HP33) to 0.067 (HP17), averaging 0.015 (Table 3, Dxy; Fig. 2). Some proteins had notably greater levels of geographic differentiation than average. For example, proteins with inter-region Dxy were more than twice the average including the Ser/Arg-rich protein F-US2, the viral capsid proteins F-UL17 and F-UL18, plus two HPs (HP17 and HP38, Table 3, Dxy, range 3.3–6.7%; Fig. 2).

Divergence by protein function

Protein divergence rates can indicate, for example, the severity of functional constraint, and the occurrence of diversifying selection in different environments. We evaluated protein divergence in terms of the estimated nonsynonymous substitution rate (Ka) as well as the relative rate of nonsynonymous to synonymous substitutions (ω, Yang et al., 2000). Relative rates of nonsynonymous substitution less than, equal to, or greater than 1 indicate that proteins are evolving under purifying selection, neutrality, or positive selection respectively (averaged across sites and branches; Yang et al., 2000). As expected, these two metrics were highly correlated (Pearson’s r = 0.948; see Table 3 Ka∕Ks and ω).

Estimates of protein divergence ranged from 0 (F-UL01, F-UL04, F-UL20, F-UL33 HP33, HP36) to 5.3% (HP38) in the 66 ChHV5 protein-coding genes examined. Several of the most conserved genes are involved with core functions, such as DNA replication, DNA processing and packaging, surface and membrane proteins, and capsid assembly and structure (see McGeoch, Rixon & Davison, 2006). For example, several core proteins were nearly identical at the amino acid level between Floridian and Hawaiian strains, including glycoprotein L (F-UL01), a gene involved in DNA cleavage/packaging (F-UL33), a nuclear protein (F-UL04), an egress protein (F-UL20; Table 3, Ka), along with two hypothetical proteins with unknown function (HP33 and HP36). Given the conserved nature of HP33, HP36 and F-UL04, further investigation of their biological role is warranted. Many of the genes that appear conserved among ChHV5 variants, such as the DNA cleavage/packaging protein (F-UL33), glycoprotein L (F-UL01), and the F-UL20 egress protein (Table 3, Ka), are also highly conserved among strains of the alphaherpesvirus HHV-1 (Szpara et al., 2014). However, there were some core genes that were not conserved at the amino acid level between the two geographic ChHV5 strains, including F-UL17 and F-UL15B, which are both involved with DNA processing and packaging (Table 3, Ka). The most divergent proteins at the amino acid level included the Ser/Arg-rich protein F-US2 and the hypothetical proteins HP17 and HP38 (Table 3, Ka).

Glycoproteins coat virions and play a major role in viral spread because of their critical role in virus-cell adhesion and recognition (Szpara et al., 2011). Functional constraints on glycoproteins vary, as glycoproteins were among both the most conserved and divergent genes among strains of HHV-1 (Szpara et al., 2014). Similarly, in six of the seven known glycoproteins we sequenced for both geographic ChHV5 strains, a few were among the most conserved genes sequenced (glycoproteins L and M), while others were less conserved at the amino acid level, with a maximum divergence of 1.12% (glycoprotein D, F-US4; Table 3, Ka), which is approximately half the maximum divergence estimate for HHV-1 glycoproteins (2.3%, Szpara et al., 2014).

Most of the ORFs, including the glycoproteins and many of the genes with higher levels of inter-region divergence mentioned above, had ω values less than one, indicating the predominance of synonymous substitutions and purifying selection (Table 3, ω). However, several proteins had ω estimates greater than one (F-UL03, F-lec2, plus the hypothetical proteins HP14, HP22, and HP38; Table 3, ω), which may be indicative of positive selection. For most of these genes, the overall level of divergence was low (Dxy ≤1%, Table 3), making ω challenging to estimate. However, divergence at HP38 (ICP4b) was among the highest observed (4.8%, Table 3, Dxy), and both synonymous and nonsynonymous substitutions were more numerous than in other genes with ω values greater than one.

Interestingly, the ChHV5 genome (Ackermann et al., 2012) contained four genes that are atypical for alphaherpesvirus genomes: two C-type lectin-like domain superfamily (F-lec1, F-lec2), a viral sialtransferase (F-sial), and an orthologue to the mouse cytomegalovirus gene M04 (F-M04). These genes were hypothesized to play a role in pathogenesis or immunodeviation and steps were taken in that study to verify their proper annotation. It is therefore worth noting that sequences for three of the four ORFs (F-lec1, F-lec2, F-sial) were recovered in this study and they revealed levels of evolutionary conservation comparable to other ChHV5 ORFs. This provides strong confirmation that these ‘atypical’ genes are not artifacts or insertions in a single lineage, but conserved and functional components of the ChHV5 genetic repertoire. However, the F-lec2 gene had an ω of slightly greater than one (1.064, Table 3), suggesting more relaxed evolutionary constraints relative to most of the examined ORFs.

Our sequencing produced full length sequence reads for 22 of the 38 potential coding sequences designated as hypothetical proteins by Ackermann et al. (2012). The lack of disruptive nucleotide changes and evolutionary rates comparable to ORFs of known function support the coding potential for the majority of these HPs. However, stop codons were detected for all FL variants but none of the HI variants within the coding regions for three HPs: HP25 (amino acid 151); HP31 (amino acid 150), and HP38 (amino acid 174).

We used the Shimodaira–Hasegawa (SH) test (Shimodaira & Hasagawa, 1999) implemented in PAML4.2 to test the concordance of phylogenies estimated for each ORF with a global topology determined below. This provides a statistical basis for identifying loci that exhibit novel phylogenetic patterns, without a presumption of the underlying cause (such as recombination, gene conversion, or major changes in evolutionary rate). Only ORFs that contained overlapping sequence data for all eight variants, i.e., ORFs for which a full locus-specific tree could be compared to the whole-genome phylogeny, were used. Out of 27 SH tests, three adjacent ORFs (F-UL28, F-UL29 and F-UL30) individually approached significance. Given the low number of variable sites used to perform this test, its discriminatory power is likely also low. The genomic clustering of these discordant ORFs and the fact that the nature of the discordance was the same in each case, i.e., that sample FL_USF clustered with HI samples (Fig. S1), adds weight to the biological significance.

Genome-scale phylogeny compared to existing phylogeographic structure

We performed a genome-scale phylogenetic analysis to investigate possible changes between older and more recent HI samples and to generate a single best estimate of topology and branch lengths for evolutionary analysis of individual ORFs. Regions suspected to be recombinant based on the evidence described below, plus those with little sequence data, were masked (see Table 3; Masked).

Sequences for eight ChHV5 variants plus the reference genomic sequences covering 55 complete protein-coding genes were concatenated into an alignment that was 72,828 bp in length. Bayesian phylogenetic analysis produced a well-supported topology in which each variant was unique (Fig. 3). The FL variants formed a well-supported clade containing two similar variants (FL_F5 and FL_G5), while the third variant (FL_USF) was more differentiated. One Hawaiian variant sampled in Maui in 2011 (HI_21610) was slightly differentiated from the other Hawaiian variants, which were minimally differentiated from each other.

We analyzed these newly sequenced FL and HI ChHV5 variants in the phylogeographic context of those identified previously by Herbst et al. (2004). A concatenated 6,280 bp alignment was constructed from partial F-UL30 (DNA polymerase; 2,019 bp), partial F-UL27 (glycoprotein B; 1,691 bp), plus the Amplicon IV (1,214 bp) and Amplicon V (1,356) ChHV5 genomic regions examined in Herbst et al. (2004). A Bayesian phylogenetic analysis of the alignment grouped HI and FL variants together (Fig. 4). Two of the newly sequenced Hawaiian variants, HI_21611 and HI_21553, were identical to the HA variant of Herbst et al. (2004) in this genome region. Similar to the whole genome phylogeny, HI_12379 and HI_12354 were slightly differentiated from the Hawaiian variants. The newly sequenced Floridian variants were closely related to the FL_B variant of Herbst et al. (2004), yet were not identical, differing at 1-5 bases. The FL_G5 variant was intermediate between the Herbst et al., FL_B variant and the other two newly sequenced FL variants, FL_USF and FL_F5. The FL_D variant (Herbst et al., 2004), that originated from a loggerhead sea turtle, was intermediate between the Hawaiian and Floridian variants and was basal to the clade containing all FL variants, yet was quite divergent (note long branch length leading to FL_D). There were 98 fixed differences between the regional clades, and the FL_D sequence matched the Hawaiian variants instead of the Floridian variants at 48 sites (see partial alignment, Fig. 4B).

Recombination is a common feature of ChHV5 evolution

Several lines of evidence reveal recombination to be an important mode of ChHV5 evolution, as in other alphaherpesviruses (Thiry et al., 2005; Szpara et al., 2014). First, unambiguous examples of recombination between FL and HI strains were identified. Analysis of ORF-specific phylogenies showed discordant placement of HI_21610 for ORFs in a window between approximately 50–55 kb of the reference genome. Visual inspection of the genome alignment in this region revealed many shared polymorphisms between HI_21610 and FL variants. For example, in the DNA sequence alignment for packing tegument protein F-UL17 variants, the HI_21610 variant matched the Hawaiian regional pattern of substitutions (fixed differences) for a portion of the alignment (bases 1–830), then matched the Floridian regional pattern for the second portion of the gene (Fig. 5A). The shared substitutions between HI_21610 and FL strains were also apparent in the amino acid alignment, in which 10 of the 27 variable amino acids matched FL variants for HI-21610 (Fig. 5B). The tests for recombination of the HI_21610 variant in this genomic region (50,146–68,418 bp) were highly significant (Table 4).

Additional lines of evidence suggest that recombination with unsampled variants has occurred. The first is the distribution of Tajima’s D calculated for ORFs along the genome, revealing genomic regions with either positive or negative values of this statistic. Strongly negative D values indicate an excess of low-frequency polymorphisms. Indeed, examination of alignments in regions with negative Tajima’s D frequently coincided with stretches of private alleles not found in any other sample (Fig. S2). We therefore plotted the number of private alleles for each sample in 100-bp windows along the genome, transformed as a Z-score in which the mean and standard deviation are from all such windows across all sequenced samples. Values are color-coded to produce a heatmap in which the highest rates of private alleles are red and the lowest values are dark gray (Fig. S3). Several long tracts of private alleles were evident, for example, from positions 14,000–21,000 for sample FL_USF. This region corresponds to several genes in the US region of the genome (F-US1, F-US2, and F-US3A and B, F-US4) that also had high proportions of segregating sites and negative Tajima’s D values (Table 3). Another concentration of singletons for the FL_USF variant occurred in the F-UL07 through F-UL09 ORFs, and two significant break points were detected (FL_USF: 42347–47743; Table 2). The HI_12354 variant had many singletons in the F-UL28 through F-UL30, plus the F-UL38 through F-UL52 ORFs.