Computational protocol: Ancient Recombination Events between Human Herpes Simplex Viruses

Similar protocols

Protocol publication

[…] We removed low quality bases from raw reads using Trimmomatic (LEADING: 30, TRAILING 30, SLIDINGWINDOW: 4:30, MINLEN: 40; ). We then mapped trimmed reads to a HSV-2 reference genome from which we had removed the TRL and TRS regions (NC_001798) using BWA-MEM (). We finally sorted mapping files and removed duplicates using the SortSam and MarkDuplicates tools from Picard ( Coverage plots were generated with Geneious (). Coverage was uneven and usually lower in fast evolving regions, e.g., IRL and IRS (, online, p. 21). We generated consensus sequences using Geneious by calling bases at positions covered by more than 20 reads and for which more than 95% of the reads agreed. We also annotated putative coding sequences with Geneious using a conservative similarity threshold of 60% to coding sequences annotated in the HSV-2 reference genome (NC_001798).We collected all complete and near complete HSV-2 and HSV-1 genome sequences for which geographic location was known as well as the chimpanzee alpha-1 herpesvirus (JQ360576) and cercopithecine herpesvirus 2 (NC_006560) genome sequences. We discarded 7 HSV-2 (KU310662-8) and 5 HSV-1 (KU310657-61) genome sequences which originated from a single unpublished study and exhibited considerable variation in many coding sequence lengths; the final data set comprised 42 HSV-2 and 42 HSV-1 genome sequences (, online, p. 2–5). We also annotated 34 unannotated HSV-2 genomes using Geneious, as abovementioned.We aligned the extracted consensus and reference sequences using MAFFT v7 (), resulting in an alignment comprising 104 sequences and 200,070 positions (Whole genome 1; see Data Availability). We explored Whole genome 1 for evidence of recombination using RDP4 (). We set the automated search to use five recombination detection methods (RDP, GENECONV, MaxChi, BootScan, and SiScan) and validated recombination events where >2 methods agreed. We particularly focused on interspecies recombination events, as we initially suspected recombination may have occurred with the chimpanzee alpha-1 herpesvirus ().To confirm the detection of recombinant fragments, we extracted the respective ORF and performed ML analyses on the nonrecombinant and recombinant parts of the alignments (UL15 nonrecombinant, UL15 recombinant, UL29 nonrecombinant, UL29 recombinant, UL30 nonrecombinant, UL30 recombinant, UL39 nonrecombinant, and UL39 recombinant; see Data Availability). For this, we first selected the best model of nucleotide substitution using jModelTest v2.1.4 and the Bayesian information criterion; the selected models were HKY + I+G (UL15), F81 + G (UL15 recombinant), HKY + G (UL29), HKY + G (UL29 recombinant), GTR + I+G (UL30), HKY + G (UL30 recombinant), GTR + I+G (UL39), and HKY + G (UL39 recombinant; ). We then ran phylogenetic analyses per se using PhyML v3 (). Tree topology, branch lengths and all parameters of the substitution model were optimized. Tree search started from a BioNJ tree and new topologies were generated using an algorithm combining nearest neighbor interchange and subtree pruning and regrafting. Branch robustness was estimated using Shimodaira-Hasegawa approximate likelihood ratio tests (SH-aLRT).To perform phylogenomic analyses, we discarded HSV-1 and the cercopithecine herpesvirus 2 sequences from Whole genome 1, selected conserved blocks therein using Gblocks (), as implemented in SeaView v4 (), and finally excluded the interspecific recombinant regions as well as the TRL and TRS regions. The resulting alignment included 61 sequences and 135,445 positions (Whole genome 2; see Data Availability). We performed phylogenetic analyses on Whole genome 2 and on a reduced version thereof excluding the chimpanzee alpha-1 herpesvirus (to avoid any long branch attraction effect; Whole genome 3) in both ML and Bayesian frameworks. We also performed ML analyses on the ten longest nonrecombinant subalignments that we could identify by running a breakpoint distribution analysis with RDP4 (Nonrecombinant 1–10, 4,654–12,640 bp; see Data availability). These ML analyses were all performed as described above, using GTR + I+G as a model of nucleotide substitution.Bayesian Markov chain Monte Carlo analyses were run with BEAST v1.8.2 (). To determine whether a strict molecular clock model could be applied, we analyzed the clocklikeness of the ML trees using TempEst (). The coefficient of variation was <0.03 for the tree derived from Whole genome 2 and only 0.13 for the tree derived from Whole genome 3. We considered these low values as compatible with a strict clock model. We followed the recent suggestion that HSV-2 is itself the result of a cross-species transmission event from the panine lineage and calibrated the clocks by either describing the root age with a normal distribution with mean 1.6 My (million years) and standard deviation 0.1 My (95% confidence interval: 1.4–1.8 My; ) (Whole genome 2); or using the substitution rate derived from the analysis of Whole genome 2 (7 × 10−2 substitutions site−1 My−1; 95% highest posterior density: 6–8 × 10−2 substitutions site−1 My−1) to inform the clock rate with a normal distribution of mean 7 × 10−2 substitutions site−1 My−1 and standard deviation 7 × 10−3 substitutions site−1 My−1 (95% CI: 6–8 × 10−2 substitutions site−1 My−1). In both cases, the overall model also included a tree shape component, for which we used a coalescent model assuming a constant population size. We used Tracer v1.6 to check that runs had converged and that the posterior was adequately sampled. We then used LogCombiner and TreeAnnotator to combine trees and generate a maximum credibility tree summarizing the posterior tree sample (both softwares are distributed with BEAST). Branch robustness was estimated with their posterior probability. All trees were annotated and prepared for publication using iTOL v3 ().Finally, to control for possible biases induced by our reference-based assembly method, we performed a de novo assembly of all HSV-2v genomes with Tadpole 35.82 (written by Brian Bushnell), as implemented in Geneious. Most contigs were relatively short and we therefore had to use reference mapping to further determine longer contiguous consensus sequences. We could identify the same recombinant and nonrecombinant sequences and the trees determined through ML analyses of amended versions of Whole genome 2 and Whole genome 3 were very similar to those generated using HSV-2v consensus sequences obtained via reference-based assembly (, online, p. 22–23; Whole genome de novo 2 and 3; see Data availability). […]

Pipeline specifications

Software tools Trimmomatic, BWA, Picard, Geneious, MAFFT, RDP4, jModelTest, PhyML, Gblocks, SeaView, BEAST, TempEst, iTOL
Application Phylogenetics
Organisms Human alphaherpesvirus 1, Homo sapiens, Human alphaherpesvirus 2, Viruses, Human poliovirus 1 Mahoney