The neoplastic disease, fibropapillomatosis (FP) is present in sea turtles worldwide (Quackenbush et al., 1998; Lackovich et al., 1999), primarily affecting green sea turtles (Chelonia mydas; Smith & Coates, 1938) and to a lesser extent loggerhead (Caretta caretta; Harshbarger, 1991), olive ridley (Lepidochelys olivacea; Aguirre et al., 1999), Kemp’s ridley (Lepidochelys kempii; Barragan & Sarti, 1994), hawksbill (Eretmochelys imbricate; D’Amato & Moraes-Neto, 2000) and leatherback (Dermochelys coriacea; Huerta et al., 2002) sea turtles. FP manifests as external fibropapillomas or fibromas and internal fibromas, myxofibromas, or fibrosarcomas of low grade malignancy (Herbst, 1994) and when severe, can lead to immunosuppression (Work et al., 2001), secondary bacterial infections (Work et al., 2003), and death (Work et al., 2004). Fibropapillomatosis is now considered a panzootic in green turtles (Williams et al., 1994) where it is a leading cause of strandings in Hawaii (Aguirre et al., 1998; Chaloupka et al., 2008; Work et al., 2004) and Florida (Foley et al., 2005).
Transmission studies implicate an enveloped virus as the likely etiological agent of FP (Herbst et al., 1995) and herpesvirus-like particles are routinely observed via transmission electron microscopy in association with FP tumors (Jacobson et al., 1991). All available evidence suggests that chelonid alphaherpesvirus 5 (ChHV5) is the etiological agent of FP, with molecular studies showing a strong association between tumor presence and ChHV5 DNA (Lachovich et al., 1999; Lu et al., 2000; Quackenbush et al., 2001; Herbst et al., 2004; Ene et al., 2005; Duarte et al., 2012; Patrício et al., 2012). However, attempts to isolate ChHV5 in cell culture have been unsuccessful (Work et al., 2009), until recently when the virus was grown in organotypic skin cultures (Work et al., 2017). In cases such as this where the virus cannot easily be grown in culture, genomic sequencing presents a valuable strategy for investigating various aspects of viral ecology and evolution to supplement viral morphogenesis and other studies looking at pathogenesis of disease. Unfortunately, efforts to sequence the ChHV5 genome (Lu et al., 1999; Greenblatt et al., 2005b; Work et al., 2009) have been hindered by its latent lifestyle (Alfaro-Núňez et al., 2014), which maintains relatively low ChHV5 copy numbers compared to host DNA (Quackenbush et al., 2001) and precludes virion purification strategies to enrich for the viral genome prior to sequencing (Ng et al., 2009).
Obtaining sequence data from ChHV5 is an important goal because geographic differences between ChHV5 strains may be one potential factor in the differential gross manifestation (high prevalence of oral tumors) in Hawaiian green sea turtles (Aguirre et al., 2002) versus Florida green sea turtles (no oral tumors) (Foley et al., 2005), as well as the long-term decline of FP being observed in Hawaii but not Florida (Chaloupka, Balazs & Work, 2009). To date, molecular epizootiological studies of ChHV5 have focused on restricted and often non-overlapping suites of individual genes or small subsets of genes (Ene et al., 2005; Greenblatt et al., 2005a), resulting in a limited picture of the genetic variation among global populations of the virus. Nonetheless, such studies have shown that geographic structure instead of structuring exists among oceanographic regions, with ChHV5 phylogeographic patterns generally reflecting movements of the sea turtle host (Herbst et al., 2004; Ene et al., 2005; Patrício et al., 2012; Ariel et al., 2017), suggesting that ChHV5 has undergone region-specific co-evolution with sea turtle hosts (Jones et al., 2016).
In 2012, the first complete ChHV5 genome was sequenced by screening over 10,000 bacterial artificial chromosome (BAC) clones from the glottis tumor of a Hawaiian green sea turtle (Ackermann et al., 2012). Although the use of BAC clones is not efficient for sequencing large portions of the genome from numerous samples, the availability of this genome provides new avenues to explore the molecular epizootiology of FP that can overcome some of the methodological challenges described above. Here we utilized the available ChHV5 genome as a reference sequence to develop a simple long-range PCR (LR PCR) amplicon resequencing strategy based on the cost-effective Illumina MiSeq platform. Using this strategy, we sequenced partial viral genomes from the tumor tissue of multiple green sea turtles with FP from Florida (FL) and Hawaii (HI). This expanded collection of ChHV5 genomic sequence data, which is the largest presented to date, enabled the assessment of rates and patterns of change in a majority of ChHV5 protein coding genes as well as the opportunity to evaluate the concordance of multiple loci with existing hypotheses of biogeographical strain divergence derived from previous, more limited, datasets. Furthermore, characterization of the variability in particular viral protein coding genes sheds light on the evolutionary mechanisms that shape diversity and provides targets for future studies examining the linkages between geographic viral variants and differences in disease manifestation and dynamics.
Materials & Methods
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.
|Case ID||Sample ID||Sample date||Location||Tumor location||Sampled by||Mappeda||Genome coverageb|
|21610||HI_21610||6∕15∕2011||Kaaawa, Honolulu, Oahu||Flipper||T. Work||92||42|
|21611||HI_21611||6∕20∕2011||Hookipa Beach, Maui||Flipper||T. Work||83||38|
|21533||HI_21533||2∕9∕2011||Hilo, Hawaii||Flipper||T. Work||78||82|
|12354||HI_12354||3∕25∕1996||Launiupoko Beach Park, Maui||Flipper||T. Work||90||72|
|12379||HI_12379||5∕24∕1996||Halananea, Maui||Flipper||T. Work||90||56|
|ALF20110705-01||FL_F5||9∕8∕2012||Indian River Co., FL||N/A||S. Hirama||92||81|
|PRB20120615-01, Uno||FL_G5||11∕8∕2012||Martin Co., FL||Flipper||R. Butts||92||69|
|7118, Gabriel||FL_USF||8∕16∕2007||Sarasota, FL||N/A||L. Byrd||19||68|
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).
Machine-demultiplexed fastq files were imported into CLC Genomics Workbench v. 8 (https://www.qiagenbioinformatics.com) for trimming of exogenous sequences and low-quality base calls. Sequencing adapters were removed using the default scoring algorithm, bases were trimmed on an error probability of 0.05, and at most three ambiguous bases per read were allowed. Reads containing less than 50 bp after trimming were discarded. For Illumina paired-end reads, reads orphaned after the trimming step were discarded.
The combined filtered short-read data for each sample were mapped to the ChHV5 reference genome sequence (HQ878327). Mapping was performed with smalt (http://www.sanger.ac.uk/science/tools/smalt-0) version 0.7.6 with a word size of 13 and a step size of 5. Mappings for separate sequencing runs of the same biological sample were merged in samtools v. 1.3 (Li et al., 2009). Mappings were filtered to exclude those with Phred-scaled quality less than 20. Raw sequences were deposited in GenBank as a BioProject (ID: PRJNA360405), including biosample numbers SAMN06210144, SAMN06210145, SAMN06210146, SAMN06210147, SAMN06210148, SAMN06210149, SAMN06210150. For sample metadata, as well as processed alignment files associated with this project, see (Morrison et al., 2018).
Consensus sequences were generated from sam-formatted alignments using the samtools mpileup function and the view, call, and consensus functions of bcftools version 1.3.1 (https://samtools.github.io/bcftools). A consensus base was chosen at variable sites following the original samtools model (-c switch) rather than the multiallelic model, and IUPAC ambiguity codes were not used. Coverage was computed with the samtools depth function, and all sites in each sample with coverage less than 5× were converted to “N”. Open reading frames (ORFs) were parsed from consensus sequences based on the coordinates in Ackermann et al. (2012) with some minor corrections to ensure full representation of the final codons.
Polymorphism & divergence
For each ORF, sequences were confirmed by translating the nucleotide sequences into amino acid sequences in Geneious using the standard genetic code, then checking for stop codons within reading frames. Additionally, for each ORF, at least one nucleotide sequence was compared to data available in GenBank using the nucleotide BLAST search. Metrics of DNA polymorphism were computed using complete ORF sequences and the DNA Sequence Polymorphism software (DnaSP v5, Librado & Rozas, 2009) and MEGA7 (Kumar, Stecher & Tamura, 2016). The number of unique sequences, the proportion of segregating sites (Ps), and the nucleotide diversity (average number of nucleotide differences per site between two sequences (π; Nei, 1987)), were calculated for all sequences as well as for each regional set of sequences in MEGA7. DNA divergence was estimated in DnaSP by defining regional sequence sets and calculating (1) the average number of nucleotide substitutions per site (DXY, Nei, 1987, equation 10.20), and (2) numbers of synonymous (Ks) and nonsynonymous (Ka) substitutions per site (Nei & Gojobori, 1986, equations 1–3), using the universal genetic code. In order to find outlier clusters along the genome, the number of “private alleles” or “singletons”, i.e., consensus base calls that were unique to only one sample at a given nucleotide position, were extracted directly from the fasta-formatted consensus sequences in 500 bp genomic windows. In order to smooth the inherently stochastic processes involved in mutation, the count of private alleles per 500 bp window was Z-transformed on an approximately continuous distribution, allowing for outlier clusters of private alleles to be visualized.
Each ORF was initially analyzed using the Geneious Tree Builder with Tamura–Nei (Tamura & Nei, 1993) genetic distances and the neighbor-joining tree building method (Saitou & Nei, 1987), and topologies were visually compared among individual ORF phylogenies. Since 11 of the 66 ORFs we sequenced were overlapping (HP14, HP15, HP19, F-UL26.5, HP22, HP23, HP24, HP25, HP26, HP31 and HP38), they were excluded from the alignment used to construct a global phylogenetic hypothesis. We then created a concatenated alignment of the remaining 55 unique ORFs to estimate a global phylogeny. Sequences with less than 50% of the ORF data present, plus those that were potentially recombinant (multiple sources of evidence were used as described below) were filtered from the alignment. Optimal partitioning schemes and nucleotide substitution models were selected following heuristic searches based upon gene and codon position in PartitionFinder (Lanfear et al., 2012). The optimal scheme was selected using the corrected akaike information criterion (AICc; Sullivan & Joyce, 2005).
Phylogenetic analyses of the manually curated, concatenated 55 ORF alignment (72,828 bp) were completed using MrBayes v3.2.6 x64 (Ronquist et al., 2012). First, the AICc favored partitioned dataset was run. Second, Bayesian model jumping was also employed in MrBayes, in which no model is selected, but instead sampling across all 203 possible general time-reversible rate matrices occurs according to their posterior probabilities (Ronquist et al., 2012). Each MrBayes analysis included five separate runs with four heated chains and 100,000 generations, sampling after every 500 generations, and diagnostics sampled every 1,000 generations and a burn-in fraction of 25%. Convergence among samples was monitored using the average standard deviation of split frequencies reported during the run, and by inspecting that the Potential Scale Reduction Factor (PSRF) diagnostic factor approached 1.0 (Ronquist et al., 2012). Phylogenetic trees were edited in FigTree v1.4.2 (Rambaut, 2014).
In order to compare the FL and HI ChHV5 variants sequenced in this study to strains previously documented in these geographic regions, sequence alignments were created for four gene regions previously utilized to determine strains: Amplicon IV and V, partial F-UL28 (glycoprotein B), and partial F-UL30 (DNApol; Herbst et al., 2004; Ene et al., 2005). The four alignments were concatenated into one alignment (6,160 bp in length) that contained 15 variants.
Selection & recombination
Tajima’s D statistic (Tajima, 1989) was calculated in DnaSP v5 (Librado & Rozas, 2009). Under neutral evolution (Kimura, 1983), the number of segregating sites and the average number of nucleotide differences have the same expected value; deviations from equality (a nonzero D) may reflect demography or natural selection (Tajima, 1989).
The relative proportion of nonsynonymous (Ka) versus synonymous substitutions (Ks) is a common measure of the rate at which protein-coding sequences evolve, providing an index of the average selective pressure on a coding region. Since most genes experience predominantly purifying selection over their evolutionary history, the Ka∕Ks ratio is usually less than one (equivalently, Ka − Ks < 0), whereas for coding sequences evolving neutrally this ratio should converge to 1 as samples become increasingly large (Makalowski & Boguski, 1998). A Ka∕Ks ratio significantly greater than one (equivalently, Ka − Ks > 0) reveals the accumulation of coding changes faster than expected under complete neutrality, implying positive selection. However, the latter test lacks power when the number of positively selected sites is small or occur only in certain lineages (Yang & Nielsen, 2002). The Ka∕Ks ratio is often estimated as the single parameter ω (omega) in a likelihood model fit to an explicit phylogenetic topology and mutational model (Yang et al., 2000). We used PAML 4.2 (Yang, 1997; Yang, 2007) to estimate ω for each ORF. We expect ω to be more sensitive than the Ka∕Ks ratio as it is fit to the phylogenetic topology estimated for the virus as a whole and includes an explicit mutational model. A single ω was estimated for each ORF, as we did not explore lineage- or site-specific models for this small number of samples.
For ORFs that had sequence data for all samples, we tested whether individual ORF phylogenies differed from the global phylogeny using the weighted non-parametric Shimodaira–Hasegawa test (WSH; Shimodaira, 1993; Shimodaira, 1998; Shimodaira & Hasagawa, 1999; Buckley et al., 2001), using the codeml program of the PAML 4.2 package (Yang, 1997; Yang, 2007).
Six tests of recombination (RDP, GENECONV, Bootscan, Maxchi, Chimaera, and 3Seq) were implemented in RDP4 (Martin et al., 2015) for genomic regions with high numbers of singletons. For each test, the ChHV5 reference genome (HQ878327) and the FLF5 variant were used to represent each of the two clades, with the exception of one test of a region of FLF5 as potentially recombinant, in which case the FLG5 sequence was chosen to represent the FL type sequence.
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.
|Designationa||GenBanka||Predicted featuresa||Complete sequences||nb||Uniquec||Psd||πe||πe (FL)||πe (HI)||Tajima’s Df|
|F-US2||AHA93306||Ser/Arg-rich protein (pr)||6HI/3FL||891||3||0.116||0.029||0.071||0.000||−1.708|
|F-US3B||AHA93304||Cdk2 cyclin-dep kinase 2||6HI/3FL||942||3||0.058||0.015||0.034||0.000||−1.491|
|F-US3A||AHA93305||US3 protein kinase protein||6HI/3FL||1,152||5||0.056||0.015||0.031||0.000||−1.435|
|F-UL02||AHA93320||Uracil DNA glycosylase||4HI/3FL||842||2||0.015||0.009||0.000||0.000||2.191|
|HP14||AHA93324||Hypothetical Protein (HP)||4HI/3FL||624||2||0.010||0.005||0.000||0.000||2.032|
|F-UL06||AHA93325||Capsid portal pr||4HI/3FL||1,986||2||0.005||0.003||0.000||0.000||2.146|
|F-UL07||AHA93327||Herpes UL7 superfamily||4HI/3FL||915||3||0.061||0.021||0.033||0.000||−0.942|
|F-UL09||AHA93329||Origin binding pr||4HI/1FL||2,477||2||0.015||0.006||0.000||0.000||−1.377|
|F-UL17||AHA93337||DNA packaging tegument pr||5HI/3FL||1,956||3||0.037||0.019||0.000||0.006||1.085|
|F-UL15B||AHA93338||DNA packing pr||6HI/3FL||1,071||2||0.029||0.014||0.000||0.000||1.823|
|HP18||AHA93340||HP, RLSA-like (human)||6HI/3FL||696||3||0.011||0.006||0.001||0.000||1.458|
|HP19||AHA93341||HP, RL6-like (human)||6HI/3FL||483||2||0.025||0.012||0.000||0.000||1.706|
|HP20||AHA93342||HP, virion transactivator||6HI/3FL||501||2||0.010||0.005||0.000||0.000||1.520|
|F-UL18||AHA93343||Viral capsid pr VP23||6HI/3FL||992||3||0.034||0.016||0.003||0.000||1.518|
|F-UL19||AHA93344||Major capsid pr||6HI/3FL||4,152||3||0.015||0.007||0.000||0.000||1.689|
|F-UL20||AHA93345||Egress pr UL20||6HI/3FL||630||2||0.008||0.004||0.000||0.000||1.520|
|F-UL26||AHA93352||Capsid maturation protease||2HI/2FL||1,650||3||0.012||0.007||0.003||0.000||1.452|
|F-UL26.5||AHA93353||Virion scaffolding pr||2HI/2FL||882||3||0.014||0.009||0.002||0.000||1.727|
|HP23||AHA93355||C3 precursor pr||2HI/3FL||870||3||0.013||0.008||0.001||0.000||−1.231|
|F-UL28||AHA93357||DNA cleavage/packaging pr||4HI/3FL||2,253||2||0.015||0.008||0.000||0.000||−1.704|
|F-UL29||AHA93358||Single-stranded binding pr||4HI/3FL||3,588||2||0.012||0.007||0.000||0.000||−1.139|
|F-UL30||AHA93360||DNA polymerase subunit||3HI/3FL||3,453||3||0.020||0.012||0.000||0.000||1.78|
|F-UL31||AHA93361||UL31 Nuclear egress lamina pr||4HI/3FL||924||3||0.018||0.010||0.001||0.000||2.151|
|F-UL33||AHA93362||DNA cleavage/packaging pr||4HI/3FL||336||2||0.003||0.002||0.000||0.000||1.342|
|F-UL35||AHA93364||VP26 capsid pr||4HI/2FL||366||2||0.008||0.004||0.000||0.000||−0.302|
|F-UL38||AHA93367||Capsid shell pr VP19C||4HI/3FL||1,326||3||0.013||0.007||0.001||0.000||−1.001|
|HP25||AHA93368||HP bipartite NLS||4HI/3FL||510||2||0.010||0.006||0.000||0.000||−1.042|
|F-UL41||AHA93369||Tegument host shutoff pr||4HI/3FL||1,179||4||0.014||0.008||0.001||0.000||−1.121|
|F-UL42||AHA93370||DNA polymerase processivity||2HI/3FL||1,076||3||0.008||0.005||0.001||0.000||−1.535|
|F-UL43||AHA93371||Gallid UL43-like pr||2HI/3FL||1,226||3||0.015||0.009||0.001||0.000||−0.617|
|HP27||AHA93373||HP protein IG||2HI/3FL||653||2||0.012||0.007||0.000||0.000||−1.155|
|F-UL52||AEZ68791||UL52 helicase-primase subunit||6HI/3FL||1,372*||3||0.022||0.011||0.000||0.000||−1.346|
|F-lec2||AEZ68793||C-type lectin domain family||6HI/3FL||530||2||0.013||0.007||0.000||0.000||1.601|
|F-sial||ABX60166||Glycosyltransferase family 29||6HI/3FL||963||2||0.012||0.006||0.000||0.000||1.701|
|F-Nec2||AHA93379||Like CD155 and CD112||6HI/3FL||1,740||2||0.010||0.005||0.000||0.000||1.747|
|HP31||AHA93378||HP Tau-tubulin kinase 1||6HI/3FL||525||2||0.011||0.006||0.000||0.000||1.566|
|HP32||AHA93380||HP Immunoglobulin V-set||6HI/3FL||1,686||3||0.011||0.005||0.000||0.000||1.549|
|HP35||AHA93383||HP heptad repeat regions||4HI/3FL||1,440||3||0.019||0.011||0.000||0.042||0.321|
|HP37||AHA93385||HP ICP4 DNA- RNA- binding||4HI/2FL||2,094||5||0.025||0.011||0.001||0.001||0.289|
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.
High throughput sequencing strategies and large nucleotide dataset additions to public databases now provide a means to augment our understanding of herpesvirus evolution (McGeoch, Rixon & Davison, 2006; Abdelgawad et al., 2016). Genome resequencing of spatially segregated ChHV5 variants provides an in-depth examination of strain variation at the DNA and amino acid levels to facilitate comprehensive phylogeographic analyses that were previously not feasible. Through a combination of long-range PCR amplification and amplicon sequencing adaptable to any high-throughput platform, we have developed an approach to successfully sequence the majority of ChHV5 protein coding genes from unenriched host extracts. With current Illumina HiSeq output, for example, hundreds of samples could be multiplexed together, although an optimized pooling strategy that accounts for variable efficiency of long-range PCRs by region merits further exploration. Similar custom strategies are frequently utilized for virus resequencing (Kvisgaard et al., 2013; Tweedy et al., 2016). This method is cost-effective for the examination of genomic evolution of ChHV5 across space and time, and facilitates the identification of both quickly evolving genes that may be useful for molecular epizootiological studies, and more slowly evolving genes that may be of fundamental importance to viral infection. Since there is no clear, practical understanding of how genetic differences between alphaherpesviruses influence pathogenesis, the detailed examination of viral genomes is likely to lead to an increased understanding of the origins, pervasiveness, and flexibility of viral functions, and the factors that have shaped viral evolution through time.
This study documents divergence between Western Atlantic and Eastern Pacific ChHV5 strains from green sea turtles (C. mydas) within a majority of the protein coding genes of ChHV5. Many of these genes are examined for geographic variability for the first time, including seven genes in the US region and 46 genes in the UL region. Similar to previous studies, the majority of genomic variation was observed between geographic locations, with a total of 1,001 fixed differences detected between FL and HI variants. As suggested by previous molecular studies (Quackenbush et al., 1998; Quackenbush et al., 2001, Greenblatt et al., 2005a; Greenblatt et al., 2005b; Duarte et al., 2012; Patrício et al., 2012; Rodenbusch et al., 2012; Monezi et al., 2016), the divergent strains match predominantly regional movement patterns of sea turtle hosts, indicating that FP is geographically specific (Herbst et al., 2004; Ene et al., 2005; Patrício et al., 2012; Jones et al., 2016; Ariel et al., 2017). Variability in rates of sequence divergence among ChHV5 genes was consistent with results noted for other alphaherpesviruses such as HHV-1, where the most divergent proteins were approximately 2.8% divergent (Szpara et al., 2014). However, it is now possible to examine these genes in ChHV5 at a fine scale to better understand how encoded proteins are evolving and what evolutionary forces may shape geographic strain divergence.
The Amplicon IV and V regions of ChHV5 (Herbst et al., 2004) are commonly utilized for phylogenetic analyses. Here we documented that these target regions contain some of the most variable genes (HP17 and F-UL18, respectively), and quantitatively validated their utility for the determination of geographic variants. Interestingly, some of the variability observed in these gene regions was likely introduced through recombination (see below). We also detected several ORFs with ω values greater than one (F-UL03, Flec2, and HP14, 22 and 38) that were accumulating coding substitutions faster than expected under strict neutrality, which may suggest loss of functional constraint. The UL03 gene was found to be conserved in a study of HHV-1 (Karamitros et al., 2016). Divergence estimates for the other four genes (Flec2, and HP14, 22 and 38) are not comparable with studies of other herpesviruses, as these are specific to ChHV5. We did not perform a secondary likelihood ratio test to see if these estimates were significantly greater than one, as we believe the small number of samples lacks statistical power to robustly implement that test. Rather, these ORFs merit further analysis with targeted sequencing in a larger set of samples.
Overall, variability within geographic locations was low, yet in each location, we detected more variants than have been documented previously. Notably, the three samples from FL sequenced for this study were each unique and although all three were most closely related to the FL-B variant, which is common in the Indian River Lagoon (Ene et al., 2005), none matched the four variants (FL A-D) previously documented from FL coastal waters (Herbst et al., 2004; Ene et al., 2005). The approximately 6-kb fragment initially used by Herbst et al. (2004) to determine virus variants does not capture the divergence that one of our newly sequenced variants (FL_USF) exhibits when additional protein coding genes are examined, suggesting that previous studies may have underestimated strain diversity. Given that several genes in the US genomic region distinguished the FL_USF variant, inclusion of sequence data from these genes may better distinguish FL variants. Similarly, each of the samples sequenced from HI was unique, yet genetic differences among these variants were fewer than those among FL variants.
Homologous recombination, the replacement of a genomic region with sequence from a co-infecting genome via homology-mediated strand exchange, is a feature known to influence alphaherpesvirus evolution (Davison, 2000; Thiry et al., 2005; Thiry et al., 2006). Specifically, other alphaherpesviruses, such as HHV-1 (Bowden et al., 2004), bovine herpes virus 1 (BoHV-1; Thiry et al., 2006), equine herpesviruses (EHV-1 EHV-4 and EHV-9; Pagamjav et al., 2005; Abdelgawad et al., 2016), and varicella-zoster virus (Norberg et al., 2015) undergo recombination. The establishment of a latent viral infection may increase the likelihood for co-infection and in turn, recombination (Ma, Azab & Osterrieder, 2013). As such, recombination may be a driving evolutionary force that introduces variation against a background of more conservative evolution known in alphaherpesviruses (Thiry et al., 2005).
The genomic sequencing of geographic strains of ChHV5 from FL and HI green sea turtles presented here has provided the first compelling evidence of recombination among strains of ChHV5; a previously undescribed yet likely important aspect of ChHV5 epizootiology. Based upon visual inspection of sequence alignments, there was clear evidence of recombination between FL and HI geographic strains in a few ORFs (F-UL12 and F-UL17; see Fig. 5B). These examples are particularly clear because multiple sequences from both regions are available and the two regional strains are strongly diverged. Additionally, for several of the variants sequenced for this study, clusters of singleton substitutions were detected that could be more parsimoniously explained by recombination with an unsampled lineage rather than many individual mutations. Genomic regions that contained potential recombinant sequences were suggested by a phylogenetic pattern that did not match the overall topology, strongly bimodal Tajima’s D statistics, and highly clustered distributions of unique alleles. Several specific instances of hypothesized recombination between FL and HI types were subsequently confirmed by highly significant p-values for a suite of recombination tests. Although lineage-specific diversifying selection or lineage-specific mutational hotspots remain formal possibilities, the recognition of recombination as a driving force in the evolution of other alphaherpesviruses, plus the data presented here, point to recombination as the likely cause of high rates of polymorphism observed in certain regions of the genome. It is important to note that clusters of private alleles are not limited to coding regions or enriched for nonsynonymous changes, discounting selection on protein sequence as a possible cause. There is also no obvious indication of mutational hotspots, as there is no tendency for clusters of private alleles to occur in the same region in different samples. The only reasonable explanation we can postulate for the patterns in Fig. S3 is recombination between the viral lineage and another divergent lineage not represented in our sample.
The complex life histories of sea turtles are likely to have important implications to the regional distribution of ChHV5 strains (Patrício et al., 2012). After spending the first few years in oceanic environments, juveniles form aggregations of individuals from multiple nesting colonies and even different species in coastal foraging habitats (Bowen & Karl, 2007). Adults often undergo long-distance migrations (hundreds to thousands of kilometers) to mating areas and nesting beaches (Bowen & Karl, 2007). Phylogeographic analyses involving green sea turtles have confirmed that there is sufficient exchange between ocean basins to prevent long-term isolation and allopatric speciation (Roberts, Schwartz & Karl, 2004; Bourjea et al., 2007). Leatherbacks are tolerant of cold waters and have a distribution spanning from the tropics to the Arctic Circle, showing no evolutionary partitions between the Indo-West Pacific and the Atlantic and no geographic segregation of lineages by ocean basin (Bowen & Karl, 2007). Other more temperate species, such as loggerheads, ridleys and hawksbills are known to congregate at feeding grounds, and it has been noted that the swarm of these species in Brazil may create an important phenomenon for the evolution of sea turtles (Bowen & Karl, 2007). Such overlapping distributions of sea turtle hosts may create opportunity for inter-species transmission and for co-infection by divergent virus strains that are dominant in different sea turtle species. Although interspecies transmission is unusual for most alphaherpesviruses (Davison, 2000), it has been documented for ChHV5 in regions where multiple sea turtle species co-occur (Herbst et al., 2004; Ene et al., 2005; Patrício et al., 2012). In such scenarios, recombination between divergent strains may give rise to novel variants that may differ in virulence. In this way, even slowly evolving viruses, such as ChHV5 (Herbst et al., 2004; this study), may diversify rapidly.
In light of this information, it is interesting to note that in several phylogeographic analyses of ChHV5, highly divergent strains have been documented from FL waters, such as the FL_D haplotype (Herbst et al., 2004; Ene et al., 2005; called FL_3 variant by Patrício et al., 2012). Inconsistent phylogenetic patterns led to the detection of recombination in the FL_D variant, yet the results were discounted as potential artifacts resulting from mis-identifications in GenBank sequences (Patrício et al., 2012). It has also been hypothesized that highly divergent variants such as FL_D were restricted to different species of sea turtles including the Atlantic loggerhead (Herbst et al., 2004). Yet divergent lineages have been reported in Pacific olive ridley turtles, and in a green sea turtle from San Diego (Patrício et al., 2012). Visual inspection of the sequence alignment of the 6,160 bp ChHV5 segment utilized by Herbst et al. (2004) (Fig. 5A) shows that the FL_D variant shares many fixed substitutions with HI strains as well as FL strains (Fig. 5B), suggesting this lineage may be the result of a previous recombination event that involved divergent oceanic strains. As suggested by Jones et al. (2016), our results suggest that further genomic profiling of ChHV5 from different regions is essential for improving our understanding of zootiology and pathogenesis of this virus.
Now that recombination has been clearly demonstrated among ChHV5 strains, future work should examine whether certain genes are more likely to have histories including recombination. For example, a high frequency of substitutions was observed in HP38 (ICP4b) in our dataset, which is consistent with known recombination in the ICP4 gene of equid herpesviruses (Pagamjav et al., 2005). Also, the genomic feature of inverted repeats shared by the subfamily Alphaherpesvirinae allows for segment inversion due to recombination during DNA replication, given adequate sequence similarity (Davison, 2000; Thiry et al., 2005). Efforts to amplify the two approximately 8-kb inverted repeat regions between the US and UL were not successful, in line with previous reports that these regions are recalcitrant to standard PCR protocols given low conservation rates, presence of microsatellite loci, and palindrome stem loops (Zhao, Haqqi & Yadav, 2000; Kieleczawa, 2005; Kieleczawa, 2006; Karamitros et al., 2016). In HHV-1, the inverted repeat regions plus a 400-bp terminal redundancy play a key role in recombination between the long and short genomic segments (McGeoch et al., 1986; McGeoch et al., 1988). These recombination events are thought to be essential for viral replication and in vivo infection (Jenkins & Roizman, 1986; Jenkins & Martin, 1990). Given this importance of the inverted repeat region in other herpesviruses, future studies could strive to recover full genomic sequences from different isolates. Typically, repeat regions in herpesviruses are determined by cloning of purified viral DNA from monolayer cell cultures into plasmids (McGeoch et al., 1991), an approach not currently feasible with ChHV5 that can only be grown with difficulty in organotypic cultures (Work et al., 2017). The addition of long-read high throughput sequencing may prove useful (e.g., Karamitros et al., 2016).
Methodologies that estimate the prevalence of ChHV5 co-infection, a precursor to recombination, as well as those that determine the relationship between virus genotype and disease phenotype (comparative pathology), should be priority topics of future research. The difficulty in isolating and culturing ChHV5 in the laboratory hampers more traditional approaches to investigate recombination rates by capturing concatemers (e.g., Slobedman, Zhang & Simmons, 1999; Meurens et al., 2004) and by assessment of the biological properties of recombinant strains (e.g., Muir, Nichols & Breuer, 2002). However, the next-generation sequencing approach described here should prove useful in this endeavor, as the sequencing of numerous global variants from all sea turtle species will allow for more robust statistical tests for recombination and better characterization of the frequency and implications of recombination events. Through additional ChHV5 genomic sequencing, researchers have access to a greater array of diagnostics. A coordinated effort will allow integration of datasets, including population genetics of the virus and host species, in addition to molecular and epizootiological views of regional and global infection, for a comprehensive understanding of ChHV5 biology, virulence and transmission, potentially leading to more accurate risk assessments and disease mitigation strategies.
The genome re-sequencing of spatially segregated ChHV5 variants presented here allows for an in-depth examination of strain diversity at the nucleotide and amino acid levels. Genome-wide SNPs confirm broad divergence of geographic ChHV5 strains and general conservation at the protein level. However, this resequencing effort has provided the first compelling evidence that recombination is a common feature of ChHV5 evolution, and is likely an important aspect of the epizootiology of fibropapillomatosis. In the future, we suggest surveys of turtles in coastal foraging habitats, where an overlapping distribution of sea turtle species may provide the opportunity for inter-species ChHV5 transmission, co-infection by divergent virus strains, and possibly novel variants that may differ in virulence.
Long-range PCR primers used to amplify ChHV5 genomic sections in green sea turtles
High throughput sequencing runs, reads, and mapping success to the ChHV5 genome
Alternative placement of samples HI_12354 and HI_12379 for three ChHV5 genes with SH tests that approached significance
(A) Topology of the global phylogeny (B) Phylogenies for F-UL28, 29 and 30, showing that the different placement of the HI-12354 and HI_12379 samples.
Estimates of Tajima’s D statistic (Tajima, 1989) for ChHV5 ORFs
(A) Plotted as a histogram (B) Ordered by genomic position.