Differential Infection Patterns and Recent Evolutionary Origins of Equine Hepaciviruses in Donkeys

Donkey sera (n = 829) were collected in five European countries (Germany, Spain, Italy, Bulgaria, and France), as well as in Asia (Israel), Africa (Kenya), and Latin America (Costa Rica and Mexico), between 1974 and 2016 (Table 1). For three countries (France, Germany, and Italy), sampling was conducted in multiple years, and details of the annual sample characteristics in these countries are displayed in Table 2. In addition, 53 mule sera were sampled in Bulgaria in 2015. All 882 donkey and mule sera were analyzed for the presence of antibodies against the viral NS3 domain by a luciferase immunoprecipitation system (LIPS) (24, 25). Three sampling sites (Israel, Kenya, and Costa Rica) showed no serologic evidence for EqHV infection, whereas all other countries yielded positive test results (Fig. 1A). As shown in Table 1 and Fig. 1B, seroprevalence rates ranged between 8.1 and 10.7% in Germany, Spain, and Mexico. Seroprevalence rates in Italy and Bulgaria were significantly higher at 40.0 to 56.7% (corrected χ² = 62.8 and χ² = 109.1 [P < 0.0001] for Italy and Bulgaria compared to all other countries, respectively). Furthermore, within a specific country the seroprevalence rates varied between sampling years and hinted at the occurrence of focal EqHV epidemics, e.g., leading to 100% of EqHV-seropositive animals in Italy in 2015 (Table 2). However, the underlying factors responsible for the variations in seroprevalence are unknown. LIPS signal intensities from seropositive donkeys were comparable to those from seropositive horses, suggesting validity of the assay used for testing (Fig. 1B). Female donkeys were significantly more likely to be seropositive than male donkeys (35.0 versus 28.0%; corrected χ² = 4.1 [P = 0.044]; risk ratio, 1.25 [lower and upper bounds, 1.01 to 1.54]; Table 1). Seroprevalence increased significantly with animal age from 20.7% in young animals (0 to 5 years of age) to 55.5% in older animals (25 to 30 years) (Fig. 1C).

To allow sensitive molecular detection of EqHV genetic variants in donkeys, all samples were tested using two different nested reverse transcription-PCR (RT-PCR) assays. The first assay targeted specifically the EqHV 5′-untranslated region (5′ UTR) commonly used for HV detection (26) and a second assay targeted the NS3 domain that is more conserved among diverse HVs than the 5′ UTR (11). One donkey from France (sampled in 1979, age and gender unknown), one donkey from Bulgaria (sampled in 2015, a 10-year-old male), and one mule from Bulgaria (sampled in 2015, a 16-year-old female) tested positive for EqHV RNA using the 5′-UTR-based assay (0.3% of all 882 donkey and mule sera). No additional specimens tested positive for HVs using the NS3-based assay, arguing against infection of donkeys with diverse HVs beyond EqHV.

Our data enabled comparisons of EqHV infection patterns between donkeys and horses. First, viral loads, which are a quantitative marker of virus replication, were similar between equine host species infected with EqHV. Viral loads in the RNA-positive specimens from this study ranged from 8.4 × 10⁵ to 3.7 × 10⁷ genome copies/ml of serum, as determined by strain-specific real-time RT-PCR assays. These viral loads were similar to viral loads observed in horses (20, 23, 25), suggesting similar infection intensities in both equine species. Furthermore, the detection of viral RNA at comparable loads in the French sera sampled in 1979 and Bulgarian sera sampled in 2015 implicated suitability of the non-recently sampled specimens for viral RNA detection.

Viral clearance is typically delayed in HV infection, including infection with EqHV in horses (7). In our study, three serial individual specimens taken at different time points over 2 weeks (May-June 1979) were available from the RNA-positive donkey sampled in France. All three specimens, as well as both individual specimens from Bulgaria, tested positive both for EqHV antibodies (indicated as red dots in Fig. 1B) and RNA, providing evidence against immediate antibody-mediated EqHV clearance in donkeys. However, the co-occurrence of viremia and antibodies as a sign of delayed clearance was apparently much lower in donkeys at 1.1% (3 of 278 antibody-positive animals) than in horses at 2 to 30% (17, 20, 24, 25, 27, 28).

Predominantly acute resolving infections were compatible with a generally lower RNA detection rate in donkeys than in horses. Combining all available data from previous studies on horses (17, 20–25, 27, 28), 148 of 2,172 horses tested positive for EqHV RNA (6.8% [range, 0.9 to 35.5%]) compared to only 3 of 1,047 donkeys or mules when combining the data from this study with previous studies (17, 20–22) (0.3%; corrected χ² = 65.9, P < 0.00001). The low number of RNA-positive donkeys could not be explained by a putatively low exposure of donkeys to EqHV, since seroprevalence in donkeys was high at 28.3% (278 of 982 donkeys combining this and the only previous serological study [17]), although still significantly lower than in horses at 34.9% (469 of 1,343 horses from all previous studies performing serological analyses; corrected χ² = 11.1, P < 0.0009). The EqHV seroprevalence increased with the age of donkeys, which was comparable to a study on EqHV in German horses (25) but in contrast to another study on EqHV in Japanese horses (28). Finally, female donkeys were more likely to be seropositive for EqHV than male donkeys. A similar distribution was not observed for horses in two previous studies, one showing no gender-associated differences and another one showing a higher EqHV burden in male horses (25, 28).

Next, we investigated the clinical relevance of EqHV infection in donkeys by determination of aspartate aminotransferase (AST; reference value <536 U/liter), gamma-glutamyl transferase (γGGT; <69 U/liter), and glutamate dehydrogenase (GLDH; <8.2 U/liter) levels in serum of all Bulgarian donkeys (n = 201) as markers of liver damage. As depicted in Fig. 1D, liver enzymes concentrations were mainly within the reference range (29) and were comparable between the seropositive and seronegative groups, including the RNA-positive animals (indicated in color in Fig. 1D), which is in line with the reported subclinical course of infection in horses.

The full viral polyprotein genes were determined for all donkey EqHV strains, including those from the three serial bleedings from the French donkeys and those from the Bulgarian donkey and mule. The polyprotein genes encompassed 8,832 nucleotides from the French donkey EqHV strain, as well as 8,835 and 8,841 nucleotides from the Bulgarian donkey and mule, respectively. Polyprotein length and organization was identical in all cases to that observed before in EqHV from horses with the presence of all typical domains in the order C-E1-E2-p7-NS2-NS3-NS4A/NS4B-NS5A/NS5B. Maximum-likelihood (ML) phylogenetic reconstructions based on the complete polyprotein gene were highly robust, as suggested by high bootstrap support for clusters based on 1,000 replicates. In these ML phylogenetic reconstructions, the novel donkey HVs from France and Bulgaria formed two distinct viral lineages that were not monophyletic. In addition, these donkey HV lineages were interspersed between EqHV from horses and did not cluster in sister relationships to EqHV strains from horses (Fig. 2A). The close phylogenetic relationship between EqHV strains from horses and from donkeys or mules was compatible, with a narrow genetic distance of only 1.8 to 2.5% of the translated polyprotein genes of these strains. Of note, even upon inclusion of the novel donkey viruses, the EqHV patristic distance was only 6.2% on amino acid level in the translated polyprotein gene compared to 33.1% within HCV (calculated using 189 genotype 1 to 7 reference sequences from the Los Alamos National Laboratory [http://hcv.lanl.gov ]). However, most of the previous studies on EqHV in horses characterized only short regions of the viral genome. Therefore, we repeated ML reconstructions using different data sets aiming at inclusion of the complete available EqHV genetic diversity without losing too much genetic information. As expected, statistical support for grouping of basal and intermediate nodes was low for the partial NS3 (helicase/protease) and NS5B (RNA-dependent RNA polymerase) domains commonly analyzed in EqHV studies. However, these reconstructions resulted in similar phylogenies, as shown for the complete polyprotein sequences with regard to the phylogenetic relationships between EqHV strains from donkeys and horses (Fig. 2B to D).

To investigate whether potential cross-species transmission was associated with molecular adaptation, we tested for differential selection among horse and donkey EqHV lineages using codon substitution models that allow for various nonsynonymous/synonymous substitution rate ratios (dN/dS) among branches (30). Branches leading to the two common ancestors of donkey EqHV strains showed a lower dN/dS ratio (0.02) compared to the dN/dS ratio among all other branches in the complete genome data set (0.04), indicating no detectable episodic adaptive signal underlying the transmission of EqHV strains from horses to donkeys. Identical results were obtained for the data set encompassing the full NS3, for which a larger number of horse EqHV sequences were available (Fig. 2B), with a dN/dS ratio of 0.0036 in branches leading to donkey EqHV strains compared to 0.0128 among other branches. An analysis using BUSTED confirmed the absence of any signal of gene-wide episodic diversifying selection along the branches leading to the two donkey clades. A FUBAR analysis to identify site-specific selection only indicated two positively selected sites in the complete polyprotein evolutionary history, which do not appear to be related to equine-to-donkey adaptation because the donkey viruses do not share a particular amino acid residue on those positions. In conclusion, the dN/dS ratios suggested that no host adaptation is needed for the mutual infection of horses and donkeys with EqHV.

In order to determine EqHV intrahost evolutionary patterns, the complete polyprotein gene sequences of the three serial bleedings available from the French donkey were analyzed. Intrahost variability within this viral gene spanning 8,832 nucleotides was 0.17% (15 substitutions) over 2 weeks (between 23 May and 6 June 1979; Fig. 3A). Similar to HCV, most mutations and in particular the majority of nonsynonymous mutations occurred in the antigenic E2 envelope protein (31), consistent with immune pressure influencing EqHV evolution in the infected animal. However, the majority of the observed mutations did not map to the N-terminal hypervariable E2 region described for HCV (31) but accumulated in the C-terminal region of the E2 gene. To investigate whether indeed EqHV generally differs from HCV in the distribution of nonsynonymous mutations in the E2 gene, the homologous domains of 12 EqHV strains infecting horses were analyzed. As shown in Fig. 3B, EqHV strains infecting horses were similar to HCV in that 72 nonsynonymous mutations accumulated in the N-terminal region of E2, compared to only 28 nonsynonymous mutations in the C-terminal region. The different pattern observed in the EqHV-infected donkey is thus likely due to the small data set available, but potential differences of genomic variability among EqHV hosts cannot be excluded at this point. Finally, reversion of two mutations was detected across the serial bleedings (in the viral E2 and NS2 domains, Fig. 3A), which again is similar to intrahost evolution patterns observed in HCV (32). The predicted similarities in EqHV and HCV evolutionary patterns in combination with the low EqHV patristic distance suggested a limited time of EqHV evolution in equines compared to HCV in humans.

However, viral evolution may be limited by noncoding constraints such as genome-scale ordered RNA structures (GORS). Albeit the level of predicted mean folding energy differences (MFEDs; a measure of GORS) across the polyprotein-coding region was slightly higher in donkey EqHV strains than the mean MFEDs within horse EqHV strains, the overall EqHV MFED patterns showed similarities between both equine species in terms of the presence and the extent of predicted stem-loops (Fig. 3C). The overall levels of MFEDs ranging up to 11.5% in our analyses were comparable to previous analyses of EqHV (15, 24) and higher than the 8.5% described before for human HCV (33), which may imply a stronger impact of GORS on EqHV than on HCV evolution (6). However, it seems unlikely that GORS alone can account for the drastic differences between EqHV and HCV genetic diversity.

The donkey HVs sequences from 1979 represent the oldest EqHV strains described so far. In order to investigate whether these sequences could serve to calibrate the molecular clock of EqHV evolution, root-to-tip distances were analyzed as a function of sampling time. To further investigate whether the temporal signal in EqHV was potentially influenced by evolutionary pressure, root-to-tip distances were compared for complete polyprotein gene trees comprising only nonsynonymous (NS) or synonymous (S) substitutions (Fig. 3D). The complete polyprotein-based tree, as well as the trees with branch lengths reestimated in either NS or S substitutions lacked a molecular clock signal, as visualized by plotting root-to-tip divergence against year of sampling (Fig. 3E). Of note, a lack of temporal signal upon inclusion of the 1979 donkey EqHV strains was consistent with the apical phylogenetic position of two EqHV strains sampled from horses in 1997 and 1998 (20) (shown in cyan in Fig. 2C and D). Unfortunately, only a partial NS3 sequence is available for the 1997 EqHV and only a partial NS5B sequence for the 1998 EqHV strain, preventing their inclusion in our temporal analyses.

Donkey sera were collected based on availability in France, Germany, Spain, Italy, Bulgaria, Israel, Kenya, Mexico, and Costa Rica from 1974 to 2016. Animal sera were stored at −20 or −80°C prior to analysis. In addition, 53 mule samples were collected in Bulgaria in 2015. Samples were either collected as part of routine examinations (Germany, Italy, Costa Rica, and France) or under permits issued by the responsible authorities. The permit numbers were as follows: Mexico, SICUAE FMVZ-UNAM F. García-Lacy 12042013; Kenya, IACUC 2015.8; Spain, BOJA55-20/2012; Israel, KSVM-VTH/5_2013; Italy, protocol 45/2013/CEISA/COM; and Bulgaria, FVM 15/15. Host designations were assessed for all EqHV RNA-positive specimens from France from 1979 by characterization of the mitochondrial COI gene as described before (49).

All samples were analyzed for the presence of anti-NS3 antibodies by the previously described LIPS (24). Briefly, sera were diluted 1:10 in buffer A and incubated for 1 h on a rotary shaker. Renilla-NS3 fusion proteins were expressed in Cos1 cells, and 10⁷ relative light units (RLU) were added per well to the diluted sera in a 96-well plate. After incubation for 1 h on a rotary shaker, antibody-antigen complexes were immunoprecipitated by A/G beads and the RLU were determined. Each sample was measured in duplicate wells. The cutoff was calculated by the mean values of wells containing only buffer A, the Renilla-NS3 fusion protein and A/G beads plus three standard deviations as described previously (24). A positive control containing anti-EqHV antibody-positive horse serum was included in each run.

For the detection of hepaciviral RNA a hemi-nested RT-PCR assay targeting the 5′ UTR was developed based on all available EqHV 5′-UTR sequences. The primer sequences were as follows: HCV-F150, GSWSCYYCYAGGICCMCCCC; HCV-R371, CTCRTGIISYAIGGTCTACRAGRCC; and HCV-R342, GGIGCICTCGCAAGCRYGCCYATCA (I = inosine, S = C/G, W = A/T, Y = C/T, M = A/C, and R = A/G). The limits of detection were determined as the number in probit analyses conducted with SPSS V23 (IBM, Ehningen, Germany) using eight replicates per RNA concentration as described previously (26). The 95% lower limit of detection of the EqHV 5′-UTR assay was 5.7 × 10² RNA copies per reaction (range, 3.8 × 10² to 1.2 × 10³), which was well below the commonly observed viral loads in EqHV-infected horses (25). The HV NS3-based assay was described previously (11). Cross-tables were calculated using EpiInfo V7 (http://www.cdc.gov/epiinfo/index.html ) and an online tool (http://quantpsy.org/chisq/chisq.htm ). Sequencing of the complete EqHV polyprotein genes was performed by amplifying genome-spanning islets with degenerate broadly reactive oligonucleotides as described previously (11). Viral loads were determined by strain-specific quantitative real-time RT-PCR (oligonucleotide sequences available upon request) with photometrically quantified in vitro cRNA transcripts used for calculation of the standard curve as described previously (11).

Statistical analyses were done using SPSS V23 (IBM). Sequences were aligned with MAFFT (Geneious 6.1.8). Maximum-likelihood phylogenetic analyses were calculated in MEGA6 (50) and RAxML (51) using a general time reversible model with a discrete gamma distribution and a proportion of invariable sites and 1,000 bootstrap replicates. To estimate branch lengths in synonymous and nonsynonymous substitutions per site, a codon substitution model was applied in HypHy (52) that allows for branch-specific synonymous and nonsynonymous substitution rates (53). PAML (54) was used to fit a codon substitution model that allowed for a different nonsynonymous/synonymous substitution rate ratio (ω) on the branches leading to the two donkey HV common ancestors compared to the ω on the remaining branches (30). In addition, we used BUSTED (55) to search for gene-wide evidence of episodic positive selection along the branches leading to the donkey virus clades, and FUBAR (56) to identify site-specific selection patterns, both implemented in HypHy. Root-to-tip divergence was plotted against sampling time using TempEst (57). Mean folding energy differences were calculated using SSE V1.2 as described previously (12).

All polyprotein gene sequences generated in this study were submitted to GenBank under accession numbers KT880191 to KT880193 and KX421286 to KX421287 .