The Diversification of Zika Virus: Are There Two Distinct Lineages?

Zika virus (ZIKV) has caused explosive epidemics in the Pacific and the Americas, posing a serious threat to public health. Conventional opinion advocates that ZIKV evolved into two distinct lineages, namely, African and Asian. Descendants of this latter lineage dispersed globally causing major epidemics. However, based on shared amino acid replacements and phylogenetic analyses, it was recently contentiously proposed that the Asian lineage was a direct descendant of the African lineage. To address this contentious issue, we reconstructed a phylogenetic tree of ZIKV using the method based on shared amino acid replacements and found that ZIKV evolved into two distinct lineages. This supports the conventional phylogenetic divergence pattern of ZIKV. Evidence of recombination and sequencing errors was identified among the large collection of ZIKV. As such problematic sequences could confound the phylogenetic analyses, they were removed. Bayesian phylogenetic analyses using the improved sequence data enabled estimates for the divergence time in the past of the African and Asian lineages of ∼180 years ago. Moreover, we found that the Asian lineage viruses did not evolve at an elevated rate. Our findings provide additional support for the conventional opinion that the Asian lineage of ZIKV diverged from the African lineage.

Zika virus (ZIKV), a mosquito-borne flavivirus, was first discovered in the serum of a rhesus monkey caged in the Zika Forest of Uganda in 1947 (Dick et al. 1952). Since then, ZIKV was reported sporadically in both Africa and Asia (Haddow et al. 2012; Shen et al. 2016). Recently, ZIKV has undergone a rapid geographic expansion. In 2007, a ZIKV outbreak was reported in Yap Island, Federated States of Micronesia (Duffy et al. 2009). In 2013–2014, French Polynesia in the Southern Pacific experienced a ZIKV outbreak (Cao-Lormeau et al. 2014, 2016). Since the report of autochthonous transmission of ZIKV in Brazil in May 2015 (Kindhauser et al. 2016), ZIKV has spread rapidly throughout the Americas, among which Brazil has the highest number of reported cases of ZIKV. Although ZIKV infection usually leads to no symptom or mild illness (Goodman et al. 2016), ZIKV can also cause a range of neurological disorders, including Guillain–Barré syndrome and microcephaly (Petersen et al. 2016; Krauer et al. 2017). But the risk of microcephaly due to ZIKV infection in pregnancy is relatively low (Cauchemez et al. 2016).

New Evolutionary Tree of ZIKV?

Numerous conventional phylogenetic analyses of ZIKV genomes reveal the presence of two main viral lineages, that is, African and Asian lineages (Haddow et al. 2012; Faria et al. 2016; Liu et al. 2016; Pettersson et al. 2016; Shen et al. 2016; Wang et al. 2016; Zhu et al. 2016) (fig. 1supplementary fig. S1, Supplementary Material online). It should be noted that phylogenetic analyses using E and NS5 genes reveal three major lineages of ZIKV; an additional lineage has been circulating in Africa (designated African II lineage) (Gong et al. 2016; Shen et al. 2016). Unfortunately, no genome sequence has been available for strains of African II lineages to date. Therefore, we still follow the traditional two-lineage nomenclature, Asian and African (actually African I in the new nomenclature) lineages (Gong et al. 2016). However, Yokoyama and Starmer (2017) inferred a rooted phylogenetic tree using a method based on shared amino acid replacements but without any technical detail (fig. 1). They suggested that MR766, the strain isolated in Uganda in 1947, is basal to all the other ZIKV strains and the Asian and African lineages cannot be readily distinguished. Statistical tests reveal that the topologies between the conventional tree and the Yokoyama and Starmer’s (YS) tree are significantly different (P < 0.05 for the Shimodaira–Hasegawa [SH] test and P < 0.05 for the AU test).

. 1.

—Different classifications of ZIKV. The phylogenetic tree (A) was inferred using maximum likelihood method. The conventional Asian and African lineages are indicated by yellow and blue colors. The numbers on the branches represent the bootstrap values. The tree proposed by Yokoyama and Starmer is shown in the right panel (B). The Asian lineage is indicated by yellow color. The topology of the tree (C) was inferred based on the shared amino acid replacements (D). The conventional Asian and African lineages are indicated by yellow and blue colors. For the shared amino acid replacements, the mutations were labeled in colored boxes. The numbers indicate the residue positions in the polyprotein of ZIKV. Strain information (country, year, accession No.): Ada1 (Senegal, 2001, KF383118), Ada2 (Senegal, 2001, KF383119), Aaf1 (African Republic, 1968, KF383115), Aaf2 (Central African Republic, 1976, KF268950), Aaf3 (Central African Republic, 1979, KF268948), Aaf4 (Senegal, 1984, HQ234501), Aae-Malaysia (Malaysia, 1966, KX377336), Alu1 (Senegal, 1968, KF383116), Alu2 (Senegal, 1997, KF383117), Aop (African Republic, 1980, KF268949), MR766 (Uganda, 1947, LC_002520), and Nigeria (Nigeria, 1968, HQ234500).

The YS tree was reconstructed “by establishing the tree topology based on the shared amino acid replacements in their polyproteins” (Yokoyama and Starmer 2017) and by first “evaluating the number of shared amino acid changes for each pair of sequences” (personal communication with Prof. Yokoyama). Yokoyama and Starmer (2017) proposed that “the ancestral Asian strain diverged from the common ancestor of Ada1 and Ada2 because of one shared amino acid replacement at position 3,040 in the polyprotein” (Yokoyama and Starmer 2017). On close inspection of the data (supplementary figure 2 in Yokoyama and Starmer 2017), no mutation was indicated in MR766 and thus they used MR766 as the ancestral strain, essentially the outgroup. However, either of these assumptions could be invalid, because: 1) One shared replacement might be caused by convergent evolution (see Chen et al. 2016 for cases in influenza virus) or reversion (Han 2017). Actually, the residue 3,040 (residue 3,044 in fig. 1) in the Asian strain Aae-Malaysia, Ada1, and Ada2 is likely to represent the ancestral state, rather than a shared amino acid replacement as identified by Yokoyama and Starmer (2017), because the residue in these strains shares the same amino acid as Spondweni virus; 2) The Asian strain shared more replacements with other African strains than with Ada1 and Ada2 (supplementary figure 2 in Yokoyama and Starmer 2017; fig. 1); 3) Relative date of sampling is not necessarily an indication of chronological order of appearance, for example, sequences of the 1918 Spanish flu cluster within “modern” flu viruses (Worobey et al. 2014). Therefore, MR766 could not be directly used as the outgroup.

We re-examined the phylogenetic relationships among African and Asian ZIKV strains following the principle of shared amino acid replacements (see Materials and Methods for details). Our phylogenetic tree (fig. 1) had a very similar topology to the previously published conventional phylogenetic trees and was therefore significantly different from the one proposed by Yokoyama and Starmer (2017). The phylogenetic positions of Ada1, Ada2, and MR766 differ from previous analyses (Shen et al. 2016; Yokoyama and Starmer 2017), possibly because the shared replacements are limited in number and might not be phylogenetically informative. Moreover, we also performed a Bayes Factor (BF) comparison by estimating the ratio of the marginal likelihood of MR766, Ada1, and Ada2 forming a monophyletic group (the conventional tree, fig. 1 mean marginal likelihood [ln] = −36,628.07) to the marginal likelihood of MR766, Ada1, and Ada2 forming a paraphyletic group (the YS tree, fig. 1 mean marginal likelihood [ln] = −36,639.33), using the Stepping Stone sampling method. Apparently, the BF > 100 suggests that the conventional tree is more strongly supported by the data. Therefore, these findings provide no evidence for the Asian lineage being a continuation of the African lineage but do support the conventional concept of an African ancestral strain of ZIKV diverging to produce distinct African and Asian lineages.

Evolutionary Stability or Experimental Errors?

Yokoyama and Starmer (2017) observed that MR766 does not exhibit “a much shorter branch length than many other African strains” in the conventional phylogenetic analysis (fig. 1) and thus disputed the power of conventional phylogenetic methods. However, numerous previous phylogenetic analyses produced the short branches of sequences from archived samples (e.g., a case of HIV of 1960 in Worobey et al. 2008 and a case of Yersinia pestis from victims of the Black Death in Bos et al. 2011). We believe the absence of short length branches might be explained by reasons other than the choice of phylogenetic methods. Strains MR776 and Ada1, sampled 54 years apart, cluster together with very short length, which seemingly implies a very low evolutionary rate of virus in this cluster. However, similar “evolutionary stasis” was observed in influenza virus and is most likely to be caused by laboratory contamination or other experimental errors (Worobey 2008).

To identify problematic sequences, we performed a comprehensive phylogenetic analysis of 519 ZIKV genomes. We looked for closely related sequences with noticeably longer or shorter branch lengths than might be expected, arguing that they could have arisen through experimental errors (Worobey 2008). We identified ten MR766 strains that clustered together. However, these strains showed significant genetic variation which could have arisen due to their high passage level under laboratory conditions. For example, strains HQ234498, KU955594, and KU963573 had been passaged >150 times. Moreover, MR766 strain DQ859059 did not cluster with the other MR766 strains. Instead, it clustered with a strain isolated in Uganda in 1962. Accordingly, these problematic sequences were not included in our analyses.

Recombination also confounds phylogenetic analyses, because different regions have different evolutionary histories (Han and Worobey 2011). Therefore, we performed recombination analyses within ZIKV genomes. We identified a total of nine sequences of hybrid origin, among which eight were supported by at least three recombination detection statistics (Table 1). It is unclear whether these sequences represent authentic recombinants, because recombinants can also arise through experimental or sequencing errors (Han and Worobey 2011). Indeed, the observations that potential parental strain for some recombinants is MR766 and the isolate dates between recombinants and their potential parents are tens of years apart make experimental errors a plausible explanation. Nevertheless, these sequences of possible hybrid origin should be excluded from phylogenetic analyses.

Table 1

Recombination Analysis of ZIKV Genomes

Recombinant	Event	Major Parent	Minor Parent	R	G	B	M	C	S	T
1968_Senegal_ArD7117	1	1984_Senegal_41525-DAK	1968_Central African_ArB1362	+	+	+	+	+	−	+
1997_Senegal_ArD128000	1	1968_Central African_ArB1362	1968_Senegal_ArD7117	+	+	+	+	−	−	+
	2	1968_Senegal_ArD7117	2001_Senegal_ArD158095	+	+	+	+	+	+	−
	3	1984_Senegal_41525-DAK	2001_Senegal_ArD158095	+	+	+	+	−	−	−
2000_Senegal_ArD142623	1	2001_Senegal_ArD157995	Unknown	−	−	−	−	−	−	+
	2	Unknown	1984_Senegal_ArD41519	+	+	+	+	−	−	−
	3	Unknown	1984_Senegal_ArD41519	+	+	+	+	+	+	−
	4	1984_Senegal_ArD41519	2001_Senegal_ArD158095	+	+	+	−	−	−	+
	5	1947_Uganda_MR766	Unknown	−	+	−	+	−	+	−
	6	Unknown	1997_Senegal_ArD128000	−	−	−	+	+	−	+
	7	1997_Senegal_ArD128000	Unknown	−	−	−	+	−	+	−
2001_Senegal_ArD157995 (Ada2)	1	1947_Uganda_MR766	Unknown	+	+	+	+	+	−	+
	2	1947_Uganda_MR766	2000_Senegal_ArD142623	+	+	+	+	−	+	+
	3	2001_Senegal_ArD158095	1997_Senegal_ArD128000	+	+	+	+	+	+	−
2001_Senegal_ArD158095	1	2001_Senegal_ArD157995	1947_Uganda_MR766	+	−	−	+	+	+	+
2016_Singapore_SG-010	1	2016_Singapore_SG-030	Unknown	−	−	−	−	−	−	+
2016_Singapore_SG-030	1	2016_Singapore_SG-001	Unknown	+	+	+	+	+	+	+
2016_Singapore_SG-042	1	2016_Singapore_SG-030	Unknown	−	−	+	+	−	+	+
2016_Singapore_SG-047	1	2016_China_ZK-YN001	1947_Uganda_MR766	+	+	+	−	−	+	−

Note.—Recombination detection method: R, RDP; G, GENECONV; B, Bootscan; M, MaxChi; C, Chimaera; S, SiScan; T, 3Seq.

Elevated Evolutionary Rate of Asian Lineage?

To further explore the evolutionary time scale and rate of ZIKV, we performed Bayesian phylogenetic analyses using a clean data set, which excludes the potential problematic and recombinant sequences. The Bayesian phylogenetic analyses support the conventional classification of ZIKV. The median time of the most recent common ancestor (TMRCA) of ZIKV was estimated to be 1,837.8 (95% highest posterior density [95% HPD]: 1,776.5–1,887.7) (fig. 2), much earlier than 1930–1945 estimated by Yokoyama and Starmer (2017). Previous studies suggested that the Asian lineage seems to have a higher evolutionary rate (Liu et al. 2016; Yokoyama and Starmer 2017). However, we did not find significant difference in the median evolutionary rates among the whole ZIKV population (6.1 [95% HPD: 4.5–7.9] × 10−4 substitutions per site per year) and the Asian lineage (7.6 [95% HPD: 5.7–9.6] × 10−4 substitutions per site per year) as claimed by Yokoyama and Starmer (2017) and several other studies (Liu et al. 2016; Pettersson et al. 2016) (fig. 2). Unfortunately, due to the limited number of African strains, we could not directly estimate the evolutionary rate. Our findings provide no evidence for an elevated rate for Asian lineages. The difference in evolutionary rates observed in previous studies might be caused by the inclusion of more pandemic ZIKV strains. Indeed, the higher evolutionary rate within outbreaks is not unexpected, given sequences within outbreaks were sampled in a relatively short time frame and slightly deleterious mutations were not removed by purifying selection (Duchene et al. 2014; Aiewsakun and Katzourakis 2016; Holmes et al. 2016).

. 2.

—Evolutionary timescales and rates of ZIKV. (A) MCC tree of ZIKV. The tips correspond to the isolate years. The Asian and African lineages are highlighted in blue and orange, respectively. (B) The 95% HPD for the distribution of evolutionary rates of ZIKV. The distribution of all ZIKV and Asian lineage is highlighted in gray and blue, respectively.

Our time-scaled tree also provides novel insights into the evolutionary history of functionally important amino acids. An amino acid substitution (A188V in nonstructural protein 1 [NS1] protein) is associated with increased infectivity of the ZIKV strain responsible for the epidemic in the Americas in the mosquito vector Aedes aegypti (Liu et al. 2017). The A188V substitution was thought to contribute to the recent emergence of ZIKV in the South Pacific and South America (Liu et al. 2017). However, our Bayesian phylogenetic analysis shows that the 188 V arose in Southeast Asia in between 1,998.9 (95%HPD: 1,994.1–2,003.0) and 2,003.6 (95%HPD: 2,000.1–2,006.6), which suggests the substitution occurred earlier than expected. The residue 188 of NS1 protein in strains isolated in South Asia after 2012 are all Valine, suggesting the A188V might have been fixed in the viral population of South Asia. Given the strain with 188 A was isolated in Philippines in 2012, the virus variants with 188 V and 188 A might have been cocirculating for some time. Although it remains unclear what drove the fixation of A188V, it appears to be not directly related to the emergence of the recent ZIKV epidemic, as proposed recently (Han 2017).

Materials and Methods

Phylogenetic Analyses Based on Conventional Methods and Shared Amino Acid Replacements

Genome sequences of ZIKV and Spondweni virus were retrieved from GenBank. The accession numbers and isolation information of those viruses are in supplementary figure S2, Supplementary Material online. The nucleotide sequences were aligned using MAFFT (Katoh and Standley 2013) and manually edited. The conventional phylogenetic analyses were performed based on both nucleotide and protein sequences using maximum likelihood method with GTR + I + Γ substitution model for DNA and JTT +I + Γ substitution model for protein implemented in MEGA7 (Kumar et al. 2016). The substitution model was selected using jModelTest 2 (Darriba et al. 2012). Support for the tree was assessed using 1,000 bootstrap replicates. The evolutionary tree of ZIKV was also inferred using a method based on shared amino acid replacements. Polyprotein of Spondweni virus was used as the outgroup to identify the shared amino acid replacements in ZIKV. Only the residues where Spondweni virus and some ZIKV strains are the same were considered, because in this case the residue of Spondweni virus is likely to reflect the ancestral state. For each pair of ZIKV, the number of shared amino acid replacements was counted. The tree was inferred by finding the sequences with maximum shared replacements and clustering together sequentially.

Statistical Tests of Tree Topologies

We compared whether the two competing tree topologies (our conventional tree vs. the YS tree) are significantly different using both the SH test and the AU test (Shimodaira 2002) implemented by PAUP* 4.0a. Moreover, to explore which topology is a better fit to the data, we estimated BF, the ratio of the marginal likelihood of one topology to the marginal likelihood of the competing topology, using the Stepping Stone sampling method (Morey et al. 2016) implemented in MrBayes 3.2.6 (Ronquist 2012). The GTR + I + Γ substitution model was used. A total of ten million generations were run, with diagnostics frequency set to once every 5,000 generations.

Large-Scale Phylogenetic and Recombination Analyses

All currently available ZIKV genome sequences (a total of 519) were retrieved from GenBank. The nucleotide sequences were aligned using MAFFT (Katoh and Standley 2013) and manually edited. Due to high computational cost, we reconstructed a comprehensive phylogenetic analysis of all ZIKV genomes using an approximately maximum likelihood method implemented in FastTree 2 (Price et al. 2010) with GTR and CAT model. Recombination events within ZIKV genomes were scanned using 3SEQ method implemented in RDP 4 (Martin et al. 2015). The recombination events were confirmed by other recombination detection methods (RDP, GENECONV, Bootscan, MaxChi, Chimaera, and SiScan) and only the events detected by more than three methods were considered to be positive ones.

Bayesian Phylogenetic Inference

Bayesian phylogenetic analyses of ZIKV and Asian lineage were reconstructed using BEASTv.1.8.2 (Drummond et al. 2012). Following Liu et al. (2016), we accommodated variation in evolutionary rate among lineages using a lognormal distribution in an uncorrelated relaxed molecular clock model and modelled changes in effective population size through time using the Bayesian Skyline model (Drummond et al. 2005). The GTR + I + Γ substitution model was used. MCMC chains were run for 200 million generations, sampling every 1,000 generations. The convergence was indicated by effective Sample Size (ESS) > 200. The parameters (median and 95% HPD intervals) were summarized after discarding a 10% burn-in using Tracer. The TreeAnnotator was used to infer maximum clade credibility (MCC) trees.

Supplementary Material

Supplementary data are available at Genome Biology and Evolution online.

Click here for additional data file.