Equine herpesvirus 4 infected domestic horses associated with Sintashta spoke-wheeled chariots around 4,000 years ago

Abstract

Equine viral outbreaks have disrupted the socio-economic life of past human societies up until the late 19th century and continue to be of major concern to the horse industry today. With a seroprevalence of 60-80 per cent, equine herpesvirus 4 (EHV-4) is the most common horse pathogen on the planet. Yet, its evolutionary history remains understudied. Here, we screen the sequenced data of 264 archaeological horse remains to detect the presence of EHV-4. We recover the first ancient EHV-4 genome with 4.2× average depth-of-coverage from a specimen excavated in the Southeastern Urals and dated to the Early Bronze Age period, approximately 3,900 years ago. The recovery of an EHV-4 virus outside the upper respiratory tract not only points to an animal particularly infected but also highlights the importance of post-cranial bones in pathogen characterisation. Bayesian phylogenetic reconstruction provides a minimal time estimate for EHV-4 diversification to around 4,000 years ago, a time when modern domestic horses spread across the Central Asian steppes together with spoke-wheeled Sintashta chariots, or earlier. The analyses also considerably revise the diversification time of the two EHV-4 subclades from the 16th century based solely on modern data to nearly a thousand years ago. Our study paves the way for a robust reconstruction of the history of non-human pathogens and their impact on animal health.

Keywords: equine herpesvirus 4; horse; palaeogenomics; sintashta; virus evolution.

Figure 1.

Sample map and phylogenetic relationships. (A) Location of the samples used in this study, including subclades, sample dates, and numbers of genomes. For samples with geographic information at country-level, the coordinates of the regional or country capital were used. (B) Maximum-likelihood tree and node bootstrap support. Sample names comprise Genbank accession numbers, sampling locations and, in parentheses, sampling dates in years BP with ‘2017’ set as ‘Present’ based on the most recently sampled isolates. The symbols refer to Fig. 1A. EHV-1 was used as an outgroup.

Figure 2.

EHV-4 genome organisation and protein structure predictions. (A) %GC and coverage was calculated within non-overlapping 50 bp and 100 bp sliding windows, respectively. Coverage is given as a percentage based on a maximum of 7.3×. The open reading frame encoding for the tegument protein (tp) UL14 (highlighted in red) was found to have been positively selected along the branch ancestral to subclades I and II (Branch M, Fig. 1B). The grey area on the coverage track shows the drop in coverage represented by the duplicated repeat regions including ORF64–ORF66. cmp: capsid maturation protein; cpp: capsid portal protein; csp: capsid scaffold protein; ct: capsid triplex; dutp: deoxiuridine triphosphatase; gp: glycoprotein; hps: helicase primase subunit; hsp: host shutoff protein; mcp: major capsid protein; nelp: nuclear egress lamina protein; nemp: nuclear egress membrane protein; np: nuclear protein; pcs: DNA polymerase catalytic subunit; obh: origin binding helicase; pps: DNA polymerase processivity subunit; rnr: ribonucleotide reductase; scp: small capsid protein; ssbp: single-stranded DNA binding protein; stk: serine/threonine protein kinase; tk: thymidine kinase; tp: tegument protein; ung: uracil DNA glycosylase. (B) Secondary structures predicted by I-TASSER for UL14 proteins. Amino acid substitutions are shown in light purple (substitutions differing between the ancient and modern versions) and in dark purple (substitutions differing from the most commonly observed version). The confidence score C denotes the quality of predicted models ranging from −5 to 2. The predictions for the most common version and UR17x29 are shown in green and red c-scores, respectively.

Figure 3.

Time-scaled EHV-4 evolutionary history. (A) Bayesian phylogeny under the strict clock and birth–death skyline serial models based on the ancient and modern dataset. Node posterior supports are provided above nodes to the left, while median estimates of divergence times are provided below nodes together with the 95%  CI in square brackets [Height_95 per cent_HPD]. Sequence labels refer to Genbank accession numbers, sampling locations and, in parentheses, sampling dates formatted as years BP with ‘Present’ set to ‘2017’. These dates were used for tip-dating calibration. The symbols refer to Fig. 1A. (B) Posterior distribution of the tMRCA of all EHV-4 strains considered in this study (N1). Time is indicated in calendar years (BCE: Before Common Era, CE: Common Era).