Risk assessment and genomic characterization of Zika virus in China and its surrounding areas.

BACKGROUND: Zika virus (ZIKV) has emerged as a global pathogen causing significant public health concerns. China has reported several imported cases where ZIKV were carried by travelers who frequently travel between China and ZIKV-endemic regions. To fully characterize the ZIKV strains isolated from the cases reported in China and assess the risk of ZIKV transmission in China, comprehensive phylogenetic and genetic analyses were performed both on all ZIKV sequences of China and on a group of scientifically selected ZIKV sequences reported in some of the top interested destinations for Chinese travelers.
METHODS: ZIKV genomic sequences were retrieved from the National Center for Biotechnology Information database through stratified sampling. Recombination event detection, maximum likelihood (ML) phylogenetic analysis, molecular clock analysis, selection pressure analysis, and amino acid substitution analysis were used to reconstruct the epidemiology and molecular transmission of ZIKV.
RESULTS: The present study investigated 18 ZIKV sequences from China and 70 sequences from 16 selected countries. Recombination events rarely happens in all ZIKV Asian lineage. ZIKV genomes were generally undergone episodic positive selection (17 sites), and only one site was under pervasive positive selection. All ZIKV imported into China were Asian lineage and were assigned into two clusters: Venezuela-origin (cluster A) and Samoa-origin cluster (cluster B) with common ancestor from French Polynesia. The time of most recent common ancestors of Cluster A dated to approximately 2013/11 (95% highest posterior density [HPD] 2013/06, 2014/03) and cluster B dated to 2014/08 (95% HPD 2014/02, 2015/01). Cluster B is more variable than Cluster A in comparison with other clusters, but no varied site of biological significance was revealed. ZIKV strains in Southeast Asia countries are independent from strains in America epidemics.
CONCLUSIONS: The genetic evolution of ZIKV is conservative. There are two independent introductions of ZIKV into China and China is in danger of autochthonous transmission of ZIKV because of high-risk surrounding areas. Southeast Asia areas have high risk of originating the next large-scale epidemic ZIKV strains.

Introduction

Zika virus (ZIKV) is a single-stranded, positive-sense RNA virus in the family Flaviviridae, genus Flavivirus. While it is transmitted among humans mainly through mosquito-bites (Aedes aegypti and Aedes albopictus), it can also be spread through maternal-fetal transmissions, blood transfusions, and sexual activities. In humans, although ZIKV infections typically cause mild and self-limiting febrile illness, it is established that it can also cause severe neurological complications (eg, Guillain-Barré syndrome) and adverse fetal outcomes (eg, microcephaly).[ ZIKV was first isolated in 1947 from a sentinel rhesus monkey in the Zika forest in Uganda,[ but it has become a significant global concern in 2015, when it broke out in Brazil,[ and subsequently spread throughout America. According to the latest update of WHO, as of February 2018, ZIKV has been reported in 86 countries and territories all over the world.[

Previous genetic studies have revealed that ZIKV has evolved into three distinct genotypes: West African (Nigerian cluster and Senegal cluster), East African (Uganda cluster), and Asian.[ However, Asian lineage is the only genotype known to cause ZIKV outbreaks, whereas no outbreak of African lineage isolates has been reported.[ While research efforts have mainly focused on ZIKV in Americas, little is known about the genetic relationships among ZIKV strains in Asia, although China and many Southeast Asian countries are in ideal climate conditions of tropical and sub-tropical weather and abundant in competent mosquito vectors for ZIKV transmission. China has experienced a series of imported cases since 2016. From February 2016 to May 2017, 29 ZIKV infection cases were imported into China.[ Singapore has experienced an outbreak of ZIKV in August, 2016,[ similarly Malaysia, the Philippines, India, Thailand, and Indonesia all reported evidence of ZIKV circulation and transmission before 2016.[

Up to now, much remains unknown about ZIKV epidemiology and evolution, in part owing to a lack of genomic data. However, due to the recent epidemics and enabling technologies such as rapid and full-length sequencing, substantial amount of ZIKV genomes can be obtained in a short period of time. China still under high risk of importation and potential localization of ZIKV. So, here we give a detailed description of all the ZIKV isolates and illustrate their introduction path in China. A comprehensive phylogenetic and genetic analysis on isolates in China and surrounding areas was performed. By exploring phylogenetic relationships, selective pressure, and recombination events of ZIKV genomes, much valuable information can be provided for exploring the genetic evolution of ZIKV in China and assessing the risk of ZIKV transmission in the future.

Methods

Sample collection

Selected ZIKV complete/coding sequences were retrieved from the National Center for Biotechnology Information (NCBI) Virus Variation Resource database[ (hereinafter called NCBI database) as of August 2018. In selecting the ZIKV sequences, all sequences reported from China and only these selected from Zika-endemic areas that are popular Chinese tourist or business destinations were considered. The country list with the number of passengers imported was obtained from China National Bureau of Statistics (http://www.stats.gov.cn/tjsj/), and the top 15 countries with over 5 million passengers arriving China in 2016 were checked for their level of ZIKV-exposing conditions according to WHO Classification Table of ZIKV.[ Countries that belonged to Categories 1 and 2 were defined as ZIKV-popular areas, and the ZIKV strains from these countries were retrieved from the NCBI database. In addition, ZIKV strains in Brazil, USA, French Polynesia, and Yap Island where ZIKV outbreak has been reported, and those in Venezuela, Samoa, Suriname, and Guatemala where ZIKV in China are imported from, were also collected. ZIKV strains in Uganda, Senegal, and Nigeria were collected as African lineage representatives. Only full-length sequences were included for more informative molecular epidemiological analyses.[ For convenience of calculations, the number of ZIKV strains was stratify sampled according to the total number of strains in database. If the total number was less than 5, all strains were included; if the total number was between 5 and 10, 5 strains were included; if the total number of strains was between 10 and 20, half of strains were included; if the total number of strains were more than 20, 10 strains were included. The set of virus sequences included were isolated at different time points. The detail of sampling process can be found in Supplementary Data 1.

Sequence alignment

All sequences were aligned together using MUltiple Sequence Comparison by Log-Expectation (MUSCLE) implemented in Molecular Evolutionary Genetics Analysis (MEGA) 7.0 (Center for Evolutionary Medicine and Informatics, Philadelphia, PA, USA)[ and both ends were trimmed by manual editing. Alignment was also carried out by MUSCLE based on codons for selection pressure analysis and the non-synonymous (dN) and synonymous (dS) substitution rates analysis.

Recombination detection

To prevent potential bias during phylogenetic inference due to recombination, sequences of all 88 ZIKV genomes were first analyzed with Recombination Detection Program version 4 (RDP4 program)[ (University of Cape Town, South Africa) that incorporates seven methods: RDP, Gene Conversion, Chimaera, Maximum Chi-squared, Bootscan, Sister-Scanning, and three Sequences, to uncover evidence for recombination events. Only events with P values ≤0.01 that were confirmed by four or more methods were considered to be recombination events.

ML phylogenetic analysis

The suitability of substitution models for our ZIKV alignment was assessed using jModelTest 2.1.10 (GitHub, USA),[ which performed a statistical model selection procedure based on the Akaike information criterion. This identified the best fitting substitution model General Time Reversible with gamma-distributed rate variation across sites (GTR + Γ) for ML phylogenetic analysis. A ML phylogeny was estimated from this alignment using PhyML 3.0 (ATGC bioinformatics platform, France)[ under a general time-reversible nucleotide substitution model, with a gamma-distributed among site rate variation (GTR + Γ), as determined by jModelTest2. Statistical support for nodes of the ML phylogeny was assessed using a bootstrap approach with 100 replicates. Support values for all nodes were embedded in the phylogenetic tree files.

Molecular clock and phylogeographic analyses

Temporal evolutionary signal in ML tree was evaluated by using TempEst v1.5.1 (Institute of Evolutionary Biology, USA),[ which plots sample collection dates against root-to-tip genetic distances obtained from the ML phylogeny tree. A linear regression between sample collection dates and root-to-tip genetic distances obtained from the ML phylogeny indicated that the feasibility of a molecular clock approach. The plot indicated that the data set contained sufficient temporal signal for molecular clock analysis. Molecular clock phylogenies were estimated using the Bayesian Markov Chain Monte Carlo (MCMC) approach implemented in BEAST v1.10.1 (Institute of Evolutionary Biology).[ We computed five independent runs of 100 million MCMC steps, sampling parameters, and trees every 10,000 steps. An uncorrelated lognormal relaxed molecular clock model[ and a Bayesian skyline coalescent model[ were used; previous studies have demonstrated this combination to be the best fitting model combination for ZIKV in the Americas.[ In a previous study, different substitution models produced similar results and divergence date estimates are robust among different combinations of prior distributions, molecular clock models, and coalescent models.[ In each run, a SRD06 substitution model[ was used, which employs an Hasegawa, Kishino, and Yano (HKY) nucleotide substitution model, a gamma-distribution among site rate variation (HKY + Γ) and a codon position partition (positions [1 + 2] vs. position 3). The program Tracer v1.7.1 (Institute of Evolutionary Biology)[ was used to evaluate MCMC chain convergence and to compute marginal posterior distributions of parameters, after removal of 10% of the chain as burn-in. The program Logcombiner (Institute of Evolutionary Biology) was used to combine and sub-sample posterior tree distributions after a 10% burn-in. TreeAnnotator (Institute of Evolutionary Biology) was used to generate a summary maximum clade credibility (MCC) tree from the posterior distribution of trees (after removal of MCMC burn-in of 10%). The MCC phylogeny tree was drawn using Figtree (Institute of Evolutionary Biology).

Selection pressure analyses

The selection pressures acting on the transmission of ZIKV was inferred by using Datamonkey web server (http://www.datamonkey.org/, USA).[ To test if positive selection has occurred on a proportion of branches, adaptive branch-site random effects likelihood (aBSREL) method[ was used. In the scenario where all branches were tested, P values at each branch were corrected for multiple testing by using the Holm-Bonferroni correction. Significance was assessed using the likelihood ratio test at a threshold of P ≤ 0.05. In addition, we estimated the dN and dS substitution rates on a per-codon site basis for each genome. Sites evolving under positive selection under a proportion of branches were detected by using the mixed effects model of evolution (MEME).[ Sites under positive selection in the entire phylogeny were detected by using single likelihood ancestor counting (SLAC)[ and Fast, Unconstrained Bayesian AppRoximation (FUBAR)[ methods. The significant level was set at P ≤ 0.05 and posterior probabilities >0.9 in FUBAR.

Amino acid substitutions analyses

MEGA7.0 was used to analyze the variation of amino acid between cluster A and B in China, as well as reference clusters in Singapore, French Polynesia, and Malaysia. The structure of protein was modeled by Swiss-Model.[ SuperPose Version 1.0 (University of Alberta, Canada) was used to compare the protein structures and get root mean squared deviations (RMSD) value. The structures are considered very similar when RMSD < 1.1 Å.[

Results

ZIKV isolates selection

A total of 88 ZIKV sequences were included for analyses [see Supplementary Data S1, S2], including 18 from China and 70 from selected countries/regions. It should be mentioned that there are 27 sequences from China in the NCBI database, but nine of them were excluded after checked their sources and sequences from one identical case were removed. The remained 18 sequences are imported from Samoa (9), Venezuela (5), Cambodia (2), Guatemala (1), and Bangladesh (1). The distribution map of the included ZIKV sequences and the introduction path of ZIKV into China are shown in Figure 1.

Figure 1

Map where samples were collected in each of the colored areas. The introduction path of ZIKV in China are shown in dotted line. ZIKV: Zika virus.

Recombination events in ZIKV strains

The RDP analysis suggests four recombination events in ZIKV genomes, showed in Table 1. The recombinant strains are all from Senegal and the major parent strains are from Uganda and Senegal. The results suggest that ZIKV African lineage may have experienced several recombination events in the early stage of transmission. However, after it has spread out from Africa and come to form the Asian lineage, the recombination events rarely happen.

Table 1

Recombination events detected in all ZIKV genome sequences.

Polygenetic analyses

The ML analysis of complete coding regions of the 88 sequences indicates distinct Asian and African lineages. African lineage consists of three independent clusters: Senegal, Uganda, and Nigeria. Asian lineage has caused all the recorded outbreak or endemic circulation of ZIKV in recent decades. All strains imported into China belong to the Asian lineage [Figure S1]. TempEst plot shows a positive correlation between genetic divergence and sampling time in all sequences (R2 = 0.649) [Figure S2A]. The correlation becomes stronger within the Asian lineage (R2 = 0.996) [Figure S2B]. The plot indicated that the data set contained sufficient temporal signal for molecular clock analysis.

The evolutionary rate of ZIKV genomes is estimated to be 8.61 × 10−4 substitutions/site per year (95% highest posterior density [HPD] interval 6.34 × 10−4–1.11 × 10−3 substitutions/site per year) by using Bayesian molecular clock phylogenies. As shown in Figure 2, most of the ZIKV Asian lineage isolates that cause recent outbreaks since 2015 are originated from the French Polynesia outbreak lineage in 2013, except Singapore and Thailand. ZIKV isolates imported into China were assigned into two clusters, cluster A and cluster B (posteriori probability [PP] = 1.0). Cluster A is denoted as Venezuela-origin cluster since it consisted of strains all imported from Venezuela into Guangdong, China. Cluster B is denoted as Samoa-origin cluster as most of strains were imported from Samoa into Zhejiang or Guangdong, China, except for two, KY967711 and MF03611 that were imported from Cambodia to Beijing, China. There are three discrete ZIKV isolates that cannot be aggregated into clade: KU744693 from Venezuela to Jiangxi, China; MF593625 from Guatemala to Henan, China, and KY328290 from Bangladesh to Yunnan, China. ZIKV genomes from an outbreak in Singapore 2016 cluster within a distinct, well-supported clade, and are independent from French Polynesia and American strains. The timing of ZIKV introduction into Singapore seems to be parallel with the introduction of ZIKV into French Polynesia. Interestingly, the ZIKV isolates in Brazil and USA are interspersed among isolates from elsewhere in the Americas and do not form distinct clusters, likely implying multiple introductions of ZIKV into the Americas. ZIKV isolates in Venezuela, Suriname, and Guatemala was likely introduce from Brazil. On the other hand, ZIKV isolates in Southeast Asia countries, such as Malaysia, Thailand, Indonesia, and Philippines are not linked to the recent outbreaks in Brazil and USA after 2015. Malaysia clade was perhaps the oldest ZIKV Asian lineages and may be the ancestor of most Asian strains that have recently caused outbreaks.

Figure 2

The maximum clade credibility tree of all included Zika virus genome sequence. Node labels are posterior probabilities indicating support for the node.

The time of most recent common ancestors (TMRCA) of Venezuela-origin cluster (cluster A) in China dated to approximately 2013/11 (95% HPD 2013/06, 2014/03) and Samoa-origin cluster (cluster B) dated to 2014/08 (95% HPD 2014/02, 2015/01). This indicates that there were two independent introductions of ZIKV into China. More specifically, the Venezuela-origin cluster in China was likely transmitted from Brazil in 2013/03, and this is supported by the evidence in Figure 2 suggesting that the isolates in Venezuela were all from Brazil branch. The Samoa-origin cluster in China was introduced from French Polynesia in 2012/05, which appears to be parallel with the introduction of ZIKV in Brazil from French Polynesia. The TMRCA of Singapore clade was estimated in 2015/01 (95% HPD 2014/07, 2015/05), which is before the first ZIKV infection case reported in Singapore (2016/05). The divergence time of Singapore cluster (2011/01, 95% HPD 2009/05, 2012/07) is earlier than French Polynesia cluster (2012/12, 95% HPD 2012/06, 2013/05). ZIKV genomes sampled in Malaysia form a robust monophyletic cluster (PP = 1.0) and share a common ancestor with other ZIKV Asian strains that circulated in recent decades. The Malaysia cluster diverged at a much earlier stage (1957/10, 95% HPD 1949/01, 1964/02) of the Asian lineage.

Evidence of episodic diversifying selection is detected in three out of all 165 branches in the phylogeny of 88 sequences of ZIKV genome (KU744693_CHINA_2016_02, KF383116_SENEGAL_1968_01, and KF383118_SENEGAL_2001_01) by aBSREL; however, no episodic diversifying selection is identified on internal nodes. As shown in Table 2, there is evidence of episodic positive/diversifying selection at 17 sites as detected by MEME. However, no sites are detected by SLAC. One positively selected site (amino acid site 2449, posterior probability = 0.932, belong to polyprotein NS4B) is detected by FUBAR. The positive selection sites are distributed throughout polyprotein of ZIKV, but most of them located on non-structural proteins, especially on NS3 (5 sites), NS5 (4 sites), and NS4B (3 sites).

Table 2

Amino acid sites under episodic positive/diversifying selection detected by MEME.

Amino acid substitutions in ZIKV clusters

Genetic variation can provide important insights into ZIKV phylogeny and pathogeny. Amino acid substitutions in ZIKV cluster A and cluster B in China, Singapore cluster, French Polynesia cluster and Malaysia cluster were explored in the analyses. As seen in Table 3, there are eight amino acid sites varied between cluster A and cluster B in China, and they are located throughout ZIKV genome. It is also noted that half (4/8) of the different sites are located in NS5 protein. After compared the two structure models of NS5 protein predicted by Swiss-Model according to two sequences different between these four sites, the structures are similar (RMSD of all atoms = 1.13) [Figure 3]. In comparison with Singapore cluster, French Polynesia cluster, and Malaysia cluster, Samoa-origin cluster (B) in China is more variable than Venezuela-origin cluster (A).

Table 3

The sites of amino acid substitution in different ZIKV clusters.

Figure 3

The modeled structure of NS5 protein with amino acid substitutions between two clusters in China.

Discussion

Although ZIKV was first isolated in Uganda in Africa, cases of ZIKV infection were only detected sporadically during the next decades. Only 14 human cases were reported before April 2007.[ The initial ancestor of ZIKV and the potential transmission routes remain elusive mostly due to potential sampling bias[ and inconclusive adaptation signal of ZIKV from mosquito to human likely as a result of limited availability of isolates.[ ZIKV has been imported into China since 2016; however, no comprehensive phylogenetic and genetic analyses have been performed based on all ZIKV sequences in China. Since the accuracy of molecular evolution of the virus significantly depends on genome sequences included in the study, we performed a scientific sampling method and suitable sampling size on ZIKV sequences to help fill the gap in understanding ZIKV transmissions in China.

The evolutionary rate of all ZIKV genomes (8.61 × 10−4 substitutions/site per year) is similar to the rate only in Asian lineage (1.05 × 10−3 substitutions/site per year), which is consistent with other study,[ and the rate is comparable to that of other Flaviviruses.[ Based on the result of polygenetic analysis, ZIKV isolates imported into China are assigned into two clusters, that agrees well with previous study.[ However, there are three sequences (KU7744693, MF593625, and KY328290) that do not appear to belong to the two clusters. KU7744693 is isolated from the first imported case in China from Venezuela.[ It was imported into Jiangxi, that is a likely reason why it does not cluster with other strains in Venezuela-origin cluster which were all imported into Guangdong. MF593625 is a strain imported from Guatemala and it shares a common ancestor with strains in Guatemala. It is remarkable to note that that strain KY328290 was originally imported into Yunnan, China from Bangladesh in 2014.[ It is undetected until retesting the ice-stored serum samples retrospectively in 2016 and reported as the first recorded ZIKV imported case in China. KY328290 share a common ancestor with the ZIKV strain in Thailand in September 2010. Although the case is imported from Bangladesh, no ZIKV strain in Bangladesh was included in the NCBI database. It was confirmed as an Asian-lineage transition-state intermediate between low-prevalence and explosive pandemic strains. This strain indicates that China may have imported ZIKV long before 2016, so, the risk of autochthonous transmission of ZIKV in China has increased.

ZIKV was first reported in Asia in Malaysia in 1966[ and it may be the ancestor of all ZIKV Asian strains. While ZIKV is widespread throughout Thailand, very few strains were reported. This may reflect the limited capacity of ZIKV detection in Thailand.[ The outbreak lineage in Singapore probably originated within Asia.[ It is consistent with our polygenetic analyses that outbreak in Singapore in 2016 is not linked to the South American epidemics. There are two divergent ZIKV strains in Singapore. The MF988734 isolate is a strain imported into Singapore from Cuba, and it is distinctly different from autochthonous strains circulating in Singapore.[ The KY241788 isolate is from the first confirmed case of ZIKV infection in Singapore in May 2016, who had traveled from Brazil. ZIKV in French Polynesia was introduced from Southeast Asia and then transmitted into South America causing pandemic spread of Zika disease. The strain that responsible for the epidemic of ZIKV in Yap Island in 2007 and the Cambodian cases were also most likely originated in Southeast Asia.[ Therefore, it is likely that ZIKV may have established latent circulation in Southeast Asia without causing any noticeable outbreak, while still remains as a potential source of strains that can cause outbreaks, which should be considered seriously.

Selection pressure analyses in all included ZIKV sequences shows that three positive selections have occurred on all branches, but no positive selection has been identified on internal nodes. However, this may be not interpreted as ZIKV strains experienced positive selection in the three branches. Because tip branches are measured over a short sampling period, therefore may include a higher proportion of mutations that have not yet been removed through purifying selection. MEME aims to detect the individual sites under episodic positive/diversifying selection. But SLAC and FUBAR were used to detect sites under pervasive positive selection. The different result between MEME and SLAC/FUBAR indicates that ZIKV genomes were undergone generally episodic positive selection but rarely pervasive positive selection. The result shows that ZIKV genome is generally conserved, which is in line with other study.[ Although there are eight amino acid substitutions in China cluster A and B, none of them were reported of significant biological influence. The amino acid substitution S109N, K709R, and N3144S were tested and showed non-significant enhance of ZIKV infection.[ M2634V substitution was reported of negligible contribution to ZIKV pathogenesis.[ The biological significance of other substitutions remains to be clarified.

Based on the polygenetic results in our study, we can conclude that the selection of sequences via scientific sampling in this study is a well representative of all ZIKV sequences. The computational redundancy can be lessened dramatically by scientifically sampling of sequences and without reducing accuracy. It is reasonable to included only one of these sequences isolated from one same case, although these samples were sequenced by different institutes using different methods, few differences were observed in the genomes.[ But there are spaces for improvements to this study in the future. On one hand, the number of sequences in some countries/regions is too little to draw any conclusions. For example, only one sequence was reported in each of the countries/regions, such as the Philippines, Indonesia, and Micronesia and none was reported in India and Samoa. This may reflect differences in surveillance intensity among different countries. On the other hand, the genomes isolated in China used in this study were all retrieved from NCBI database; however, it may be better that all samples are re-collected and re-sequenced to provide more information and avoid possible confounding factors that may affect the results.

In conclusion, the genetic evolution of ZIKV is conservative. There are two independent introductions of ZIKV into China and China are in danger of autochthonous transmission of ZIKV because of high-risk surrounding areas. Southeast Asia areas have high risk of originating the next large-scale epidemic ZIKV strains.

Funding

This study was supported by grants from the Guangdong Provincial Science and Technology Project (2018B020241002, 2018B020207013), the National Science and Technology Major Project (2018ZX10101002), and the National Natural Science Foundation of China (31872499, 81501744).

Conflicts of interest

None.