Introduction
============

The Zika virus (ZIKV) is a member of family Flaviviridae, genus *Flavivirus*, and is transmitted among humans primarily by *Aedes* mosquitos. The single-stranded positive-sense RNA virus has a length of approximately 10,800 nucleotides with a 5′- and 3′-UTR. It encodes a single polyprotein that undergoes post-translational processing to yield capsid (C), membrane precursor (prM), envelope (E), and seven nonstructural (NS) proteins (NS1, NS2A, NS2B, NS3, NS4A, NS4B, and NS5) ([@msw270-B7]; [@msw270-B11]). The ZIKV (MR766) was first isolated in 1947 from a sentinel monkey (*Macaca mulatta*) in the Zika forest of Uganda ([@msw270-B12]). Although dispersing within the African Continent, some strains also spread to Malaysia and then to Philippine and Micronesia, forming a new Asian lineage ([@msw270-B23]; [@msw270-B19]). After sporadic human infections in Africa and Asia, ZIKV outbreaks were reported in Micronesia ([@msw270-B15]) in 2007 and French Polynesia ([@msw270-B37]) in 2013--2014. It is suspected that ZIKV reached Brazil during the Polynesian epidemic ([@msw270-B18]) and has spread rapidly throughout the Americas ([@msw270-B8]; [@msw270-B18]).

About 80% of ZIKV infections in humans are asymptomatic and when they appear the symptoms are usually mild, including rash, fever, red eyes, joint pain, and headache ([@msw270-B15]; [@msw270-B36]; [@msw270-B30]). This description of ZIKV infection had to be changed drastically in the summer of 2015 when the association between ZIKV and microcephaly, Guillan--Barre syndrome and other serious brain disorders was found in the Americas ([@msw270-B4]; [@msw270-B35]; [@msw270-B44]) and retrospectively in French Polynesia ([@msw270-B6]). Because of these findings and the authentication as a direct causative agent by the US Centers for Disease Control and Prevention (CDC) ([@msw270-B41]), ZIKV has become an emerging infectious public threat ([@msw270-B4]; [@msw270-B6]; [@msw270-B8]; [@msw270-B30]; [@msw270-B35]; [@msw270-B41]; [@msw270-B44]). Indeed, using human stem cells ([@msw270-B20]; [@msw270-B48]) and mice ([@msw270-B10]; [@msw270-B28]; [@msw270-B34]), ZIKV is shown to target neural stem cells (neurospheres), brain organoids, placenta, fetal brains and to cause microcephaly; ZIKV has also been found in human fetal brains ([@msw270-B35]). Curiously, only a portion (1--13%) of mothers infected with ZIKV have babies with microcephaly and other brain dysfunctions ([@msw270-B6]; [@msw270-B25]).

To understand the molecular mechanisms, by which ZIKV causes these serious medical problems, one crucial step is to identify mutations that induce cytopathicity and pathogenicity. One effective way of predicting these mutations is to study how they have accumulated in different ZIKV strains and their possible functions. Expansions of ZIKV strains to various geographic regions and potentially critical mutations have been studied primarily by constructing phylogenetic trees of ZIKV and other flaviviruses ([@msw270-B15]; [@msw270-B23]; [@msw270-B5]; [@msw270-B19]; [@msw270-B8]; [@msw270-B16]; [@msw270-B18]; [@msw270-B45]; [@msw270-B56]). These evolutionary trees show that 1) the strains isolated from Asia, Pacific Islands and the Americas (collectively called Asian strains) and their ancestral African strains form two evolutionarily distinct groups and 2) the numbers of nucleotide substitutions per site (or branch lengths) of various strains measured from the ZIKV ancestor are very similar. However, at least the latter result is in error because MR766 basically stopped evolving in 1947 and American strains accumulated additional mutations for the next 69 years.

As described below, the currently known evolutionary trees of different ZIKV strains have been constructed by basically assuming that various strains evolved with similar speeds ([@msw270-B18]). When this assumption is removed, drastically different tree topologies of ZIKV strains emerge. These new evolutionary trees allow us to properly evaluate the numbers of nucleotide (or amino acid) substitutions per site per year (or evolutionary rates of mutant substitution) for different strains and also identify potentially critical mutations in the Asian strains under the proper historical context.

The goal of this paper is to find a biologically acceptable evolutionary tree of ZIKV strains and try to identify critical mutations. For this purpose, we consider some representative ZIKV with complete genome sequences and exclude partially sequenced African strains and many completely sequenced American strains without losing generality. Throughout the paper, the strains collected from humans are denoted by the countries of their origin and those isolated from mosquitos are indicated by abbreviated species names. The roots of these trees were determined using the Spondweni virus (SM-6V-1) as the outgroup (indicated by stars, [supplementary table S1](#sup1){ref-type="supplementary-material"}, [Supplementary Material](#sup1){ref-type="supplementary-material"} online).

Results and Discussion
======================

Conventional Phylogenetic Trees of ZIKV Strains
-----------------------------------------------

Before finding, a new evolutionary tree of ZIKV strains, those of 9 Asian and 11 African ZIKV strains were constructed first by applying the maximum likelihood (ML), neighbour-joining (NJ), and maximum parsimony (MP) methods directly to their nucleotide sequences using the MEGA6 web server ([@msw270-B47]). The resulting three evolutionary trees show virtually identical topologies and the NJ and ML trees have similar branch lengths ([fig. 1*A*](#msw270-F1){ref-type="fig"} and [Supplementary figs. S1*A* and *B*](http://mbe.oxfordjournals.org/lookup/suppl/doi:10.1093/molbev/msw270/-/DC1), [Supplementary Material](#sup1){ref-type="supplementary-material"} online). These ML and NJ trees were constructed using the Jukes and Cantor (JC) model (see Materials and Methods). When the gamma distances with the JC model were used ([@msw270-B39]), the tree topologies do not differ significantly. For example, the NJ trees with the two distance measures show only two minor differences: MR766 and Ada1 switched their phylogenetic positions in group 1, whereas Alu2 and Aaf4 switched their phylogenetic positions in group 3. Fig. 1.Phylogenetic trees of 20 representative ZIKV strains. Numbers in front of ZIKV strains indicate the year of isolation. The root was determined using Spondweni virus (SM-6V-1) as the outgroup. (*A*) The NJ tree, which was obtained by applying the Neighbour-Joining (NJ)-method directly to the sequence data ([@msw270-B47]). Numbers at various nodes indicate clustering percent supports generated by 1000 bootstrap analyses. Blue and orange rectangles indicate groups 1--3 of African strains and Asian strains which were distinguished using the NJ method, respectively. (*B*) The PE tree, where the tree topology was determined first and then branch lengths were determined using the ML method ([@msw270-B52]).

Another phylogenetic tree was obtained using the Bayesian Evolutionary Analyses Sampling Trees (BEAST v1.8.4) method ([@msw270-B14]). Running 100 × 10^6^ generations with the substitution model (HKY + three codon positions), the clock model (strict) and the tree model (constant), the minimum effective sample size (ESS) of 1795 was obtained ([supplementary fig. S1](#sup1){ref-type="supplementary-material"}*C*, [Supplementary Material](#sup1){ref-type="supplementary-material"} online). Despite assimilating the "time-stamped data," which certainly show a shorter branch length for MR766, the BEAST tree has a very similar topology as those of the ML, MP, and NJ trees ([supplementary fig. S1](#sup1){ref-type="supplementary-material"}*C*, [Supplementary Material](#sup1){ref-type="supplementary-material"} online). The BEAST trees based on the uncorrelated lognormal clock (300 × 10^6^ generations, ESS \> 376) and the uncorrelated exponential clock (300 × 10^6^ generations, ESS \> 326) are very similar to the BEAST tree with the strict clock model, except that Alu2 and Nigeria are switched within group 3.

All of these trees show that the ancestor of African strains MR766, Ada1, and Ada2 (group 1) split from the common ancestor of Aaf1, Aaf2, Aaf3, and Aop (group 2) and Alu1, Alu2, Aaf4, and Nigeria (group 3). These groupings are depicted in the blue rectangle ([fig. 1*A*](#msw270-F1){ref-type="fig"}). Supported by high bootstrap values (the NJ, ML, and MP trees) and high posterior probabilities (the BEAST tree), these evolutionary steps are, statistically speaking, highly reliable and show that MR766 appeared only recently, which is consistent with the currently known evolutionary trees of ZIKV strains ([@msw270-B23]; [@msw270-B19]; [@msw270-B16]; [@msw270-B18]; [@msw270-B45]; [@msw270-B56]). In reality, however, MR766 was isolated in 1947 and has a low passage history ([@msw270-B23]), it is expected to have a much shorter branch length than many other African strains.

To explain these "conventional" tree topologies and similar branch lengths of different strains measured from the ZIKV ancestor, we have to assume that the evolutionary rates of nucleotide substitution (or mutation rates) had to be rigorously adjusted among them so that their branch lengths become very similar at various times of virus isolation---a highly improbable proposition.

A New Evolutionary Tree
-----------------------

To implement the variable sampling history and possibly different evolutionary rates of ZIKV strains, a rooted phylogenetic tree of the same 20 ZIKV strains was constructed first by obtaining its topology based on the shared amino acid replacements (or parallel changes) in their entire polyproteins ([supplementary figs. S2 and S3](#sup1){ref-type="supplementary-material"}, [Supplementary Material](#sup1){ref-type="supplementary-material"} online) and then evaluating the numbers of nucleotide substitutions per site for all branches using the ML-method ([@msw270-B52]). In the following, this new evolutionary tree of ZIKV is referred to simply as the parallel evolution (PE) tree ([fig. 1*B*](#msw270-F1){ref-type="fig"}).

The PE tree differs dramatically from the traditional trees ([fig. 1*A*](#msw270-F1){ref-type="fig"}, [supplementary fig. S1](#sup1){ref-type="supplementary-material"}*A*--*C*, [Supplementary Material](#sup1){ref-type="supplementary-material"} online) and exhibits three major characteristics. First, as expected, MR766 is most distantly related to the other ZIKV strains and consequently the African and Asian strains cannot be clearly distinguished; moreover, MR766 has a very short branch length after it split from the common ancestor of the other ZIKV strains. Second, the ancestors of four African lineages (MR766, Ada1, Ada2, and the rest) and Aae-Malaysia diverged from each other about the same time, strongly suggesting that ZIKV evolution experienced a bottleneck. Third, the branch lengths of different strains, exemplified by Aae-Malaysia and Ada2, which were isolated in 1966 and 2001, respectively, do not reflect the times of their isolation, showing that the evolutionary rates of nucleotide substitution are extremely variable among different strains. However, the PE tree and the conventional trees have one common feature---the American strains can be traced back to the ancestors of Polynesia, Thailand, Cambodia, Micronesia (or Philippine) and Aae-Malaysia, in that order.

In constructing the PE tree, it was assumed that the ancestral Asian strain diverged from the common ancestor of Ada1 and Ada2 because of one shared amino acid replacement at position 3040 in the polyprotein ([supplementary fig.S2](#sup1){ref-type="supplementary-material"}, [Supplementary Material](#sup1){ref-type="supplementary-material"} online). Recently, considering the E and NS5 regions separately, it has been claimed that ZIKV strains should be distinguished into two distinct lineages---lineage 1 includes the African strains we are considering ([fig. 1*A* and *B*](#msw270-F1){ref-type="fig"}), whereas lineage 2 contains ZIKV strains sampled from Senegal and Cote d'Ivoire ([@msw270-B45]). Again, the trees based on the E and NS5 regions show that MR766 appeared only recently and are unrealistic. However, comparing these and the ZIKV sequences considered here ([supplementary table S1](#sup1){ref-type="supplementary-material"}, [Supplementary Material](#sup1){ref-type="supplementary-material"} online), we can identify four shared amino acid changes: mutations from I to V (I → V) (at position 3035 in the polyprotein), N → S (at 3040) and R → Q (at 3041) that are shared by Ada1, Ada2, lineage 2 and Asian strains and also Q → H (at 3150) is shared by lineage 2 and African group 3. Hence, the Asian strains seem to have evolved from the common ancestor of Ada1, Ada2, and the lineage 2 strains. To check the effect of this assumption on the final PE tree, another tree was obtained by assuming that the Asian and African ancestors took independent paths after splitting from MR766. This alternate tree ([supplementary fig. S4](#sup1){ref-type="supplementary-material"}*A*, [Supplementary Material](#sup1){ref-type="supplementary-material"} online) is practically indistinguishable from the PE tree ([fig. 1*B*](#msw270-F1){ref-type="fig"}) and the three major characteristics of ZIKV evolution still hold.

The conventional trees, represented by the NJ tree ([fig. 1*A*](#msw270-F1){ref-type="fig"}), can be easily explained from the PE tree ([fig. 1*B*](#msw270-F1){ref-type="fig"}) by considering the highly variable evolutionary rates of ZIKV strains and the methods by which phylogenetic trees are constructed. After all, these trees are constructed basically by searching for strains with smallest sequence differences and clustered together sequentially. Hence, because of their shortest branch lengths measured from the ZIKV ancestor ([fig. 1*B*](#msw270-F1){ref-type="fig"}), MR766, Ada1, and Ada2 appear to belong to the same evolutionary group in the conventional tree (group 1). By the same token, groups 2 and 3 form separate groups and group 1 seem to be more closely related to group 2 than to group 3. Finally, because of the extremely long branch length of Aae-Malaysia, the three groups of African strains and Asian strains seem to form two distinct evolutionary groups in the conventional trees.

Mutation Rates
--------------

To evaluate the evolutionary rate of nucleotide substitution of a ZIKV strain, we need to know the duration between its appearance and the time of sampling (or divergence time). Since the branch length of MR766 is virtually zero (4 × 10^−6^/site), the PE tree ([fig. 1*B*](#msw270-F1){ref-type="fig"}) suggests that the ZIKV ancestor appeared in 1947, the year of its sampling. However, as described below, MR766 also accumulated some mutations and the PE tree, particularly earlier stages of ZIKV evolution, must be interpreted with caution. In fact, when another tree was constructed considering the numbers of nucleotide substitution per "codon" ([@msw270-B52]), the branch lengths of MR766 and the common ancestor of the rest of ZIKV strains are about the same (1 × 10^−4^/codon, [supplementary fig. S4](#sup1){ref-type="supplementary-material"}*B*, [Supplementary Material](#sup1){ref-type="supplementary-material"} online). Hence, it is more likely that the Asian and African strains, MR766 excluded, started to differentiate at around 1947 and the ZIKV ancestor appeared earlier than that (see Origin of ZIKV).

Divergence times (*T*~1~--*T*~9~) and evolutionary rates of nucleotide substitution (*a*~1~--*a*~19~) of ZIKV strains were evaluated by assuming that Ada1, Ada2, Aaf1, and Aae-Malaysia diverged from each other in April of 1947 and setting several short branch lengths to zero ([fig. 2*A*](#msw270-F2){ref-type="fig"}). The initial evolutionary steps of the Asian strains were studied by considering Aae-Malaysia and Micronesia together with MR766 and Aaf1 ([fig. 2*B*](#msw270-F2){ref-type="fig"}). Micronesia, rather than Philippine, was used here because their ancestors diverged from Aae-Malaysia at about the same time, but Micronesia was sampled five years earlier than Philippine ([supplementary table S1](#sup1){ref-type="supplementary-material"}, [Supplementary Material](#sup1){ref-type="supplementary-material"} online). Fig. 2.Estimation of divergence times and evolutionary rates of ZIKV strains. (*A*) A modified PE tree with divergence times (*T*~1~--*T*~9~) and evolutionary rates of nucleotide substitution per site per year (*a*~1~--*a*~19~). Three groups of strains are distinguished by blue, black and red branches with evolutionary rates of 0.05, 0.2--1.5 and 2.0--5.1 × 10 ^−^ ^3^/site/year, respectively. (*B*) A phylogenetic relationship among MR766 (strain 0), Aae-Malaysia (strain 1), Micronesia (strain 3), and Aaf1 (strain 12). *d*~i,j~s are the branch lengths between pairs of strains *i* and *j* (*i*, *j* = strains 0, 1, 3, and 12). Evolutionary rates *a*~1~, *a*~3~, and *a*~12~ are followed by divergence times (in parentheses).

Using 10,176 comparable nucleotide sites, the branch length for Aae-Malaysia (*L*~1~) was evaluated to be 0.081/site from (*d*~0,1~ − *d*~0,12~ + *d*~1,12~)/2, where *d*~i,j~s are the branch lengths between strains *i* and *j* (*i*, *j* = MR766 \[strain 0\], Aae-Malaysia \[strain 1\], and Aaf1 \[strain 12\]) ([fig. 2*B*](#msw270-F2){ref-type="fig"}, [supplementary table S1](#sup1){ref-type="supplementary-material"}, [Supplementary Material](#sup1){ref-type="supplementary-material"} online). These nucleotide substitutions occurred during the period between April 947 and July 1966, or 19.25 years, and therefore *a*~1~ = *L*~1~/19.25 = 4.2 ± 0.15 × 10^−3^/site/year. Next, the branch lengths between node O and Micronesia (*L*~3~) and between node O and Malaysia (*L*~1,1~) were estimated to be 0.030/site and 0.013/site, respectively, and *T*~1 ~=~ ~19.25 × (*L*~1~ − *L*~1,1~)/*L*~1 ~=~ ~16.2 years and *a*~3~ = *L*~3~/(60.17 − *T*~1~) = 0.7 ± 0.04 × 10^−3^/site/year (see Materials and Methods). Similarly, using appropriate formulae ([supplementary table S2](#sup1){ref-type="supplementary-material"}, [Supplementary Material](#sup1){ref-type="supplementary-material"} online), *T*~2~--*T*~9~ and the remaining evolutionary rates were evaluated.

The results reveal five characteristics of ZIKV evolution. First of all, the divergence times for the 5 major evolutionary steps of Asian strains vary between 4 and 32 years, whereas those for the 4 major steps of groups 2 and 3 of African strains range between 4 and 12 years ([table 1](#msw270-T1){ref-type="table"}). Hence, it is likely that the ancestor of the currently known American strains split from the ancestors of Polynesia, Thailand, Cambodia, Micronesia (or Philippine), and Aae-Malaysia roughly in 2013, 2007, 1999, 1995, and 1963, respectively. In African group 2 strains, Aop appeared roughly in 1961 and Afa2 and Aaf3 in 1973, whereas Alu1 split from the ancestor of Aaf1 in 1957 and the other members of group 3 appeared in 1965. Table 1.Divergence Times and the Origins of Various ZIKV Strains.Divergence Time (in years)Evolutionary Rate (×10^−3^/site/year)Strain (year of origin)*T*~1~ = 16.2*a*~3~ = 0.68Philippine and Microneasia (1963)*T*~2~ = 32.0*a*~4~ = 0.58Cambodia (1995)*T*~3~ = 4.3*a*~5~ = 0.67Thailand (1999)*T*~4~ = 7.9*a*~6~ = 0.53Polynesia (2007)*T*~5~ = 6.2*a*~7~ = 1.46Brazil-1, Haiti-1 and Suriname-1 (2013)*T*~6~ = 9.9*a*~13~ = 3.79Alu1 (1957)*T*~7~ = 8.0*a*~14~ = 4.39Nigeria, Alu2 and Aaf4 (1965)*T*~8~ = 3.7*a*~17~ = 0.92Aop (1951)*T*~9~ = 11.6*a*~18~ = 2.18Aaf2 and Aaf3 (1962)[^1]

Second, excluding MR766, the 19 strains can be classified into three groups based on their evolutionary rates: 1) Ada2 (0.05 × 10^−3^/site/year), 2) Aae-Malaysia, Aaf1, Aaf2, Alu1, and Nigeria (2.1--4.4 × 10^−3^/site/year), and 3) the remaining 13 strains (0.2--1.5 × 10^−3^/site/year) ([table 2](#msw270-T2){ref-type="table"}, [fig. 2*A*](#msw270-F2){ref-type="fig"}). Hence, the evolutionary rates of ZIKV strains differ by two orders of magnitudes and are far more variable than previously thought ([@msw270-B18]). The extremely variable evolutionary rates among different strains seem to be the major reason why the BEAST tree with the "time-stamped data" results in a typical traditional tree. Third, as expected from the functional constraints of protein-coding genes in DNA and RNA genomes ([@msw270-B22]; [@msw270-B38]), the evolutionary rate tends to be highest at the third position of a codon, followed by the first and second positions of a codon, in that order ([table 2](#msw270-T2){ref-type="table"}). Fourth, in the rapidly changed Aae-Malaysia, mutations occurred uniformly throughout the protein coding-regions; for example, the largest difference between the evolutionary rates of capsid (3.4 × 10^−3/^site/year) and propeptide (6.0 × 10^−3/^site/year) segments do not differ significantly. However, the evolutionary rates of amino acid replacement for NS2B and NS4A (∼0.4 × 10^−3^/site/year) and capsid (2.7 × 10^−3/^site/year) differ significantly (*P* \< 0.05) and the evolutionary rate for E (0.8 × 10^−3^/site/year) does not differ significantly from the average rate for the polyprotein (1.0 × 10^−3^/site/year) ([supplementary table S3](#sup1){ref-type="supplementary-material"}, [Supplementary Material](#sup1){ref-type="supplementary-material"} online). These characteristics differ significantly from human immunodeficiency viruses (HIV), in which envelope genes (∼1.0 × 10^−2^/site/year) evolved much faster than the rest of the genome (∼1.0 × 10^−3^/site/year) ([@msw270-B29]; [@msw270-B53]). Finally, Brazil-2 suggests that the evolutionary rates of nucleotide substitution at the 5′- and 3′-UTR are statistically indistinguishable from those of the coding region ([supplementary table S3](#sup1){ref-type="supplementary-material"}, [Supplementary Material](#sup1){ref-type="supplementary-material"} online). Table 2.Evolutionary Rates of Nucleotide and Amino Acid Substitutions (×10 ^−^ ^3^/site/year).ZIKVDivergence Time (years)Position Within a CodonAmino AcidFirstSecondThirdAll1. Aae-Malaysia19.251.4 ± 0.140.3 ± 0.07[\*\*](#msw270-TF3){ref-type="table-fn"}12.6± 0.504.2 ± 0.151.0 ± 0.122. Philippine16.00.3 ± 0.070.1 ± 0.04[\*\*](#msw270-TF3){ref-type="table-fn"}1.5 ± 0.170.7 ± 0.060.2 ± 0.053. Micronesia44.00.3 ± 0.040.1 ± 0.03[\*\*](#msw270-TF3){ref-type="table-fn"}1.9 ± 0.120.7 ± 0.040.2 ± 0.044. Cambodia15.10.3 ± 0.070.1 ± 0.03[\*\*](#msw270-TF3){ref-type="table-fn"}1.5 ± 0.170.6 ± 0.060.1 ± 0.055. Thailand14.80.3 ± 0.080.1 ± 0.05[\*](#msw270-TF2){ref-type="table-fn"}(\*\*)1.6 ± 0.180.7 ± 0.070.3 ± 0.076. Polynesia6.20.2 ± 0.090.05 ± 0.05[\*\*](#msw270-TF3){ref-type="table-fn"}1.4 ± 0.260.5 ± 0.090.1 ± 0.077. Brazil-11.30.9 ± 0.460(\*\*)3.3 ± 0.871.5 ± 0.331.4 ± 0.568. Haiti-11.10.8 ± 0.470(\*\*)2.8 ± 0.871.3 ± 0.341.4 ± 0.609. Suriname-11.90.6 ± 0.320.2 ± 0.16(\*\*)3.5 ± 0.741.4 ± 0.270.5 ± 0.2710. Ada1540.2 ± 0.030.2 ± 0.03(\*)0.3 ± 0.040.2 ± 0.020.3 ± 0.0411. Ada2540.02±0.010.01±0.01(\*\*)0.1±0.020.05±0.010.01±0.0112. Aaf1210.6 ± 0.090.2 ± 0.05[\*\*](#msw270-TF3){ref-type="table-fn"}5.7 ± 0.302.1 ± 0.100.4 ± 0.0713. Alu132.11.3 ± 0.190.3 ± 0.08[\*\*](#msw270-TF3){ref-type="table-fn"}10.6 ±0.563.8 ± 0.190.7 ± 0.1314. Nigeria3.11.6 ± 0.361.3 ± 0.33(\*\*)10.7 ±0.964.4 ± 0.352.0 ± 0.4015. Alu232.10.4 ± 0.060.1 ± 0.03[\*\*](#msw270-TF3){ref-type="table-fn"}0.8 ± 0.090.4 ± 0.050.2 ± 0.0516. Aaf419.10.1 ± 0.040.0 ± 0.02[\*](#msw270-TF2){ref-type="table-fn"}(\*\*)1.4 ± 0.150.5 ± 0.050.1 ± 0.0317. Aop19.40.3 ± 0.070.0 ± 0.01[\*\*](#msw270-TF3){ref-type="table-fn"}2.5 ± 0.200.9 ± 0.070.1 ± 0.0418. Aaf23.80.8 ± 0.250.2 ± 0.11[\*](#msw270-TF2){ref-type="table-fn"}(\*\*)5.7 ± 0.672.2 ± 0.240.5 ± 0.2119. Aaf36.80.5 ± 0.140.1 ± 0.06[\*\*](#msw270-TF3){ref-type="table-fn"}3.2 ± 0.371.2 ± 0.130.3 ± 0.12[^2][^3][^4]

Origin of ZIKV
--------------

To study when the ancestors of various ZIKV strains appeared, the nucleotide sequences of MR766 (strain 0) and Spondweni virus (strain 57) were compared with a designated third strain and the length of the third strain-specific branch (*L*~3rd~) was evaluated. For each third strain, the time span between the year of its divergence from MR766 (the time of its appearance) and the year of its isolation in nature was determined by dividing *L*~3rd~ values by their evolutionary rates (*a*~3rd~) and then the year of its actual appearance was estimated by subtracting this *L*~3rd~/*a*~3rd~ value from the year it was sampled. With the exception of Ada1, the branch lengths (*L*~3rd~) and the resulting estimates for the origin of ZIKV are not affected much by using either the gamma distance with the JC model or the standard JC model ([supplementary table S4](#sup1){ref-type="supplementary-material"}, [Supplementary Material](#sup1){ref-type="supplementary-material"} online).

The results show that Ada2 with the slowest evolutionary rate (*a*~3rd~ = 0.05 × 10^*−*3^/site/year) has a *L*~3rd~/*a*~3rd~ of 56 years, implying that the Ada2 ancestor appeared in 1945 ([supplementary table S4](#sup1){ref-type="supplementary-material"}, [Supplementary Material](#sup1){ref-type="supplementary-material"} online). It also suggests that the ancestor of Ada1 appeared in 1928--1941. On the other hand, Nigeria with the fastest evolutionary rate (*a*~3rd~* *= 4.4 × 10 ^*−*^ ^3^/site/year) has *L*~3rd~/*a*~3rd~ of 8.8 years and suggests that its ancestor appeared in 1959. Likewise, the evolutionary rates for Aae-Malaysia, Aaf1, and Aop lineages suggest that the ZIKV ancestor appeared sometimes during 1952--1958. Since the ZIKV ancestor appeared prior to the isolation of MR766, the latter results are unrealistic. These results show two things: 1) the mutation rate in ZIKV was initially relatively low, probably less than 0.2 × 10 ^*−*^ ^3^/site/year, but it increased significantly in many descendant lineages and 2) ZIKV appeared roughly in 1930--1945.

Accumulation of New Mutations in Asian Strains
----------------------------------------------

Using the Phylogenetic Analyses by Maximum Likelihood (PAML) with Jones--Taylor--Thornton (JTT), Whelan and Goldman (WAG), and Dayhoff models of amino acid substitution ([@msw270-B52]), amino acid changes at all branches of the PE tree ([fig. 1*B*](#msw270-F1){ref-type="fig"}) were identified. The results consistently show that the ancestors of Aae-Malaysia, Philippine, Cambodia, Thailand, Polynesia, and 24 American strains have accumulated 33, 14, 3, 5, 2, and 1 new amino acid changes, respectively ([supplementary table S5](#sup1){ref-type="supplementary-material"}, [Supplementary Material](#sup1){ref-type="supplementary-material"} online). In this inference, mutants that are polymorphic among Asian strains or those found in the Spondweni virus are excluded because it is possible that the mutant forms are actually ancestral forms. None of these mutations identified here was tested in the previous site-directed mutagenesis experiments using different flaviviruses ([@msw270-B27]; [@msw270-B2]; [@msw270-B31]; [@msw270-B9]; [@msw270-B55]; [@msw270-B13]; [@msw270-B51]; [@msw270-B57]; [@msw270-B24]) and some of these mutations could be responsible for causing serious medical problems in human populations. Recent physicochemical analyses suggest that these mutations have a low chance of causing medical problems ([@msw270-B18]). However, we have little knowledge about the individual and nonadditive effects of these mutations of ZIKV function, which need to be determined by mutagenesis experiments.

RNA Structures
--------------

Within the 5′-UTR, six mutations occurred in Aae-Malaysia and have been maintained virtually in all Asian strains ([supplementary table S6](#sup1){ref-type="supplementary-material"}, [Supplementary Material](#sup1){ref-type="supplementary-material"} online). The RNA-folding structures of the stem-loop (SL-A) of Aae-Malaysia, predicted by *mfold* modeling ([@msw270-B58]) are similar to that of MR766 ([supplementary fig. S5](#sup1){ref-type="supplementary-material"}*A* and *B*, [Supplementary Material](#sup1){ref-type="supplementary-material"} online) and no critical mutation can be identified. The 3′-UTR of ZIKV folds into five major SL structures (SL-I, SL-II, SL-III, SL-DB, and 3′-SL) ([@msw270-B56]). In this region, 16 mutations can be found, 13 of which occurred in Aae-Malaysia and the other three in Philippine have been maintained in all of their descendants ([supplementary table S6](#sup1){ref-type="supplementary-material"}, [Supplementary Material](#sup1){ref-type="supplementary-material"} online). The predicted RNA structures of Aae-Malaysia are either identical (SL-II) or very similar (SL-III) to the structure of MR766 ([supplementary fig. S5](#sup1){ref-type="supplementary-material"}*C* and *D*, [Supplementary Material](#sup1){ref-type="supplementary-material"} online, respectively) and the mutations in these regions do not seem to modify the function of the Asian strains.

Mutations in the rest of the 3′-UTR have distinct possibilities of changing ZIKV functions. The RNA structure of the SL-I of MR766 consists of 24 nucleotides (positions 35 through 58, ΔG = −12.9 Kcal/mol, [fig. 3*A*](#msw270-F3){ref-type="fig"}). Compared with this, Aae-Malaysia with a single mutation (u42c) and Philippine with an additional mutation (u9c) both generate more stable RNA structures with 24-nucleotides (Δ*G* = −15.3 Kcal/mol) or those with 58-nucleotides (Δ*G* = −14.8 to −15.6 Kcal/mol) (common segments in blue strips, [fig. 3*A*](#msw270-F3){ref-type="fig"}). Functional importance of this region for flaviviral replication is indicated by chimeric mutants of Dengue type 1 virus (DENV1) ([@msw270-B46]). Mutations (c257u, u258c, and a266g) make the upper stems of the SL-DB of Aae-Malaysia (and Philippine) longer and more stable (Δ*G* = −34 to − 33 Kcal/mol) than that of MR766 (Δ*G* = −27.3 Kcal/mol) (blue strips, [fig. 3*B*](#msw270-F3){ref-type="fig"}). Similarly, compared with the 3′-SL structure of MR766 (Δ*G* = −53.4 Kcal/mol), mutations (u424g, u426c, and c427u) of Aae-Malaysia make the bottom stem more stable (Δ*G* = −58.4 Kcal/mol) ([fig. 3*C*](#msw270-F3){ref-type="fig"}). This is the exact segment of DENV2 chimeric mutants that grew well in monkey kidney cells but have restricted growth in cultured mosquito cells ([@msw270-B54]). Fig. 3.*mfold* RNA structure models of MR 766 and Asian strains. (*A*) The SL-I structures. (*B*) The SL-DB structures. (*C*) The 3′-SL structures. The comparable segments of the African and Asian strains are shown by blue strips, whereas the mutated sites are highlighted by yellow rectangles. The SL-I, SL-DB, and 3′-SL structures were modeled considering nucleotide sites 1--60, 237--322, and 332--428, respectively, and minor loops for the latter two are not shown.

The UTR of flaviviruses are known to be required not only for such basic functions as viability, replication and translation ([@msw270-B40]; [@msw270-B42]; [@msw270-B50]) but also for inhibiting innate immune pathways ([@msw270-B32]) and small interference RNA- and micro RNA-induced RNAi pathways ([@msw270-B43]; [@msw270-B42]), all of which lead to virus-induced cytopathicity and pathogenicity. The effects of the subgenomic RNA products, which include 3′-SL, at the 3′-UTR of flaviviruses on the viral function and host adaptation have been studied using chimeric mutants that recomposed of the corresponding segments of different strains ([@msw270-B33]; [@msw270-B54]; [@msw270-B1]; [@msw270-B46]; [@msw270-B3]; [@msw270-B49]).

These observations strongly suggest that the seven critical mutations (u42c, c257u, u258c and a266g, u424g, u426c, and c427u), which occurred in Aae-Malaysia prior to 1963, may be causing the currently known medical problems. It seems paradoxical that the association between ZIKV infection and microcephaly and other serious brain disorders has not been found prior to the 2013--2014 French Polynesian epidemic ([@msw270-B30]). However, it can be easily imagined that the ZIKV-disease associations were simply missed in earlier days, as we witnessed in the original survey of the French Polynesian epidemic. The geographic regions of ZIKV outbreaks coincide with the areas at risk of DENV and Chikungunya virus (CHIKV) outbreaks ([@msw270-B5]; [@msw270-B45]). These viruses may cause various medical symptoms by interacting with each other ([@msw270-B17]; [@msw270-B21]). Before, we can fully study such interactions in disease developments, it seems critical to clarify the role of each virus in causing a specific disease, like microcephaly. Fortunately, some animal models for ZIKV infection have been developed ([@msw270-B10]; [@msw270-B28]; [@msw270-B34]). At this stage of ZIKV research, the seven point mutations in the 3′-UTR, possibly together with some of the mutations that occurred in the polyproteins of Asian strains, offer a set of specific testable hypotheses (rather than using chimeric mutants) to examine how ZIKV started causing a variety of problems in brain development.

Materials and Methods
=====================

Construction of Phylogenetic Trees of ZIKV Strains
--------------------------------------------------

The "conventional" rooted evolutionary trees of 20 representative strains were constructed by applying ML, neighbour-joining (NJ), and maximum parsimony (MP) methods directly to their nucleotide sequence data together with the outgroup, Spondweni virus (SM-6V-1), using the MEGA6 web server ([@msw270-B47]). In constructing the ML and NJ trees, branch lengths (*d*) were estimated from the proportion of different nucleotides per site (*p*~1~) using the relationship *d* = −(3/4) ln\[1 -- (4/3) *p*~1~\] under the JC model ([@msw270-B26]) and the gamma distribution with the JC model; the standard error of *d* was calculated from \[{9*p*~1~(1 -- *p*~1~)}/{(3 -- 4*p*~1~)^2^*n*}\]^1/2^, where *n* is the number of nucleotide sites compared. The number of amino acid replacements per site (*K*) was estimated from the proportion of different amino acids per site (*p*~2~), *K* = ln(1 -- *p*~2~), and the standard error was calculated from \[*p*~2~/(1 − *p*~2~)*n*\]^1/2^, where *n* is the number of amino acid compared ([@msw270-B39]). The Bayesian Evolutionary Analyses Sampling Trees (BEAST) analyses were also used to incorporate the variable times of virus sampling (or "time-stamped data") ([@msw270-B14]). For each run, we tried to ensure and achieve the ESS of more than 200.

The new phylogenetic trees were obtained first by establishing the tree topology based on the shared amino acid replacements in their polyproteins ([supplementary figs. S2 and S3](#sup1){ref-type="supplementary-material"}, [Supplementary Material](#sup1){ref-type="supplementary-material"} online) and then by evaluating the numbers of nucleotide substitutions per site for all branches applying the ML method ([@msw270-B52]). Using the ML-based Bayesian method ([@msw270-B52]) with JTT, WAG, and Dayhoff models of amino acid substitution, ancestral amino acids at all nodes of the PE tree ([fig. 2*A*](#msw270-F2){ref-type="fig"}) were inferred and amino acid replacements have been identified.

Estimation of Divergence Times and the Evolutionary Rates
---------------------------------------------------------

The evolutionary rate of nucleotide substitution of each ZIKV strain was estimated by dividing its branch length by the time span that the virus was allowed to evolve. Divergence time estimates were based on the ML tree determined from the numbers of nucleotide substitution per codon ([supplementary fig. S4](#sup1){ref-type="supplementary-material"}*B*, [Supplementary Material](#sup1){ref-type="supplementary-material"} online). It was assumed that the ancestors of Ada1, Ada2, Aaf1, and Malaysia diverged from each other in April 1947. The divergence times (*T*~1~--*T*~9~) and evolutionary rates (*a*~1~--*a*~19~) ([fig. 2*A*](#msw270-F2){ref-type="fig"}) were estimated in the following fashion.

For Asian strains, the branch length for Aae-Malaysia (*L*~1~) was first evaluated from (*d*~0,1~ − *d*~0,12~ + *d*~1,12~)/2, where *d*~0,1~, *d*~0,12~, and *d*~1,12~ are the branch lengths between MR766 (strain 0) and Aae-Malaysia (strain 1), between MR766 and Aaf1 (strain 12) and between Aae-Malaysia and Adf1, respectively ([fig. 2*B*](#msw270-F2){ref-type="fig"}). Aae-Malaysia appeared in April 1947 and sampled in July 1966, or the divergence time of 19.25 years, and therefore *a*~1~ is given by *L*~1~/19.25 ([supplementary table S2](#sup1){ref-type="supplementary-material"}, [Supplementary Material](#sup1){ref-type="supplementary-material"} online). *L*~1~ can also be evaluated using other formulae such as (*d*~1,10~ + *d*~1,12~ -- *d*~10,12~)/2 and (*d*~1,11~ + *d*~1,12~ -- *d*~11,12~)/2. However, some values of *T*~6~--*T*~9~ become negative when the *L*~1~ values are estimated using these formulae and thus they were not used.

To evaluate the evolutionary rates of Phillipine, Micronesia, and Cambodia, we need to estimate *T*~1~ ([fig. 2*A*](#msw270-F2){ref-type="fig"}). The branch lengths between node O and Micronesia (*L*~3~) and between node O and Aae-Malaysia (*L*~1,1~) were estimated from (*d*~1,3~ − *d*~1,12~ + *d*~3,12~)/2 and (*d*~1,3~ + *d*~1,12~ -- *d*~3,12~)/2, respectively ([fig. 2*B*](#msw270-F2){ref-type="fig"}). In this calculation, Micronesia, rather than Philippine or Cambodia, was used because Micronesia was sampled earlier than the other two strains, which minimize the possible effect of random genetic drift on the analyses. Then, using relationships *a*~1~ (19.25 − *T*~1~) = *L*~1,1~ and *a*~3~ (60.17 − *T*~1~) = *L*~3~, we obtain *T*~1 ~=~ ~19.25 (*L*~1~ − *L*~1,1~)/*L*~1~ and *a*~3~ = *L*~3~/(60.17 − *T*~1~) ([supplementary table S2](#sup1){ref-type="supplementary-material"}, [Supplementary Material](#sup1){ref-type="supplementary-material"} online). Similarly, pairs of *T*~2~ and *a*~4~, *T*~3~ and *a*~5~, *T*~4~ and *a*~6~ and finally *T*~5~ and *a*~7~ were evaluated together. *T*~2~ can be evaluated using the branch length of Philippine or Cambodia; however, the branch leading to Cambodia was used to assure positive values for *T*~3~, *T*~4~, and *T*~5~. Knowing these divergence times, *a*~2~, *a*~8~, and *a*~9~ can also be evaluated ([supplementary table S2](#sup1){ref-type="supplementary-material"}, [Supplementary Material](#sup1){ref-type="supplementary-material"} online).

For African strains, *a*~10~ (Ada1) and *a*~11~ (Ada2) were estimated from \[(d~10,11~ + d~10,12~ − d~11,12~)/2\]/54 and \[(d~10,11~ − d~10,12~ + d~11,12~)/2\]/54, respectively ([supplementary table S2](#sup1){ref-type="supplementary-material"}, [Supplementary Material](#sup1){ref-type="supplementary-material"} online). For African group 2 and 3 strains ([fig. 2*B*](#msw270-F2){ref-type="fig"}), *a*~12~ = \[(−d~10,11~ + d~10,12~ + d~11,12~)/2\]/21 was first evaluated and then, using relationships *a*~12~ (21 − *T*~6~) = *L*~12,1~ and *a*~13~ (21 -- *T*~6~) = *L*~13~, *T*~6~ and *a*~13~ were evaluated from 21(*L*~13~ − *L*~12,1~)/*L*~13~ and *L*~13~/(21 -- *T*~6~), respectively, where *L*~13~ = (−*d*~11,12~ + *d*~11,13~ + *d*~12,13~)/2 and *L*~12,1~ = (*d*~11,12~ − *d*~11,13~ + *d*~12,13~)/2. Continuing the process, *T*~7~ and *a*~14~, *T*~8~, and *a*~17~ and finally *T*~9~ and *a*~18~ were evaluated together and, eventually, *a*~15~, *a*~16~, and *a*~19~ were determined. In evaluating all of these evolutionary rates, the branch lengths at all positions of each codon were used so that the effects of sampling biases are minimized.

Inferring Amino Acid and Nucleotide Substitutions and RNA-Folding Modeling
--------------------------------------------------------------------------

Using the PAML with JTT, WAG, and Dayhoff substitution models ([@msw270-B52]), amino acid replacements at all nodes of the PE tree ([fig. 2*A*](#msw270-F2){ref-type="fig"}) were inferred initially by considering the 21 strains (indicated by \*, [supplementary table S1](#sup1){ref-type="supplementary-material"}, [Supplementary Material](#sup1){ref-type="supplementary-material"} online) and the maintenance of these mutations in their descendants were studied considering all of the remaining 36 strains. Because of the restriction on the currently available sequence data, the patterns of mutation accumulation in the 5′- and 3′-UTR in Asian strains were constrained to 20 and 13 strains, respectively (indicated by ^†^ and ^‡^, respectively, [supplementary table S1](#sup1){ref-type="supplementary-material"}, [Supplementary Material](#sup1){ref-type="supplementary-material"} online).

RNA-folding structures of a stem-loop (SL-A) in the 5′-UTR and those of (SL-I, SL-II, SL-III, SL-DB, and 3′-SL) in the 3′-UTR were evaluated using the *mfold* web server ( <http://mfold.rna.albany.edu>) ([@msw270-B58]). In modeling SL-I, SL-II, SL-III, SL-DB, and 3′-SL, nucleotide sites 1--106 for the 5′-UTR and 1--60, 101--170, 180--235, 237--322, and 332--428 are considered, respectively.

Supplementary Material
======================

[Supplementary data](#sup1){ref-type="supplementary-material"} are available at *Molecular Biology and Evolution* online.

Supplementary Material
======================

###### 

Click here for additional data file.

We thank Ruth Yokoyama and two anonymus reviewers for their helpful comments as well as Michael Derry and Darin Melton at Emory University for resolving IT issues. This work was supported partially by National fInstitutes of Health (EY016400) and Emory University (S.Y.).

**Associate editor:** Takashi Gojobori

[^1]: 10,176 comparable nucleotide sites were considered.

[^2]: The evolutionary rate at the second position is smaller than that at the first position (*P* \< 0.05).

[^3]: The evolutionary rate at the second position is smaller than those at the first and third positions (*P* \< 0.01).

[^4]: (\*) and (\*\*) The evolutionary rate at the second position is smaller than that at the third position (*P* \< 0.05 and *P* \< 0.01, respectively).
