Multiscale analysis for patterns of Zika virus genotype emergence, spread, and consequence

Abstract

The question of how Zika virus (ZIKV) changed from a seemingly mild virus to a human pathogen capable of microcephaly and sexual transmission remains unanswered. The unexpected emergence of ZIKV's pathogenicity and capacity for sexual transmission may be due to genetic changes, and future changes in phenotype may continue to occur as the virus expands its geographic range. Alternatively, the sheer size of the 2015-16 epidemic may have brought attention to a pre-existing virulent ZIKV phenotype in a highly susceptible population. Thus, it is important to identify patterns of genetic change that may yield a better understanding of ZIKV emergence and evolution. However, because ZIKV has an RNA genome and a polymerase incapable of proofreading, it undergoes rapid mutation which makes it difficult to identify combinations of mutations associated with viral emergence. As next generation sequencing technology has allowed whole genome consensus and variant sequence data to be generated for numerous virus samples, the task of analyzing these genomes for patterns of mutation has become more complex. However, understanding which combinations of mutations spread widely and become established in new geographic regions versus those that disappear relatively quickly is essential for defining the trajectory of an ongoing epidemic. In this study, multiscale analysis of the wealth of genomic data generated over the course of the epidemic combined with in vivo laboratory data allowed trends in mutations and outbreak trajectory to be assessed. Mutations were detected throughout the genome via deep sequencing, and many variants appeared in multiple samples and in some cases become consensus. Similarly, amino acids that were previously consensus in pre-outbreak samples were detected as low frequency variants in epidemic strains. Protein structural models indicate that most of the mutations associated with the epidemic transmission occur on the exposed surface of viral proteins. At the macroscale level, consensus data was organized into large and interactive databases to allow the spread of individual mutations and combinations of mutations to be visualized and assessed for temporal and geographical patterns. Thus, the use of multiscale modeling for identifying mutations or combinations of mutations that impact epidemic transmission and phenotypic impact can aid the formation of hypotheses which can then be tested using reverse genetics.

Fig 1. Overview of dispersal patterns of nonsynonymous mutations prior to and during the 2015–2016 ZIKV epidemic.

Mutations that are detected in 5 or more genome sequences are shown on the x axis and are listed in order of number of sequences with the mutation (high to low; 398–5). The y axis represents data from 408 ZIKV genomes and due to the number of genomes the genome names are not shown in this condensed version of the figure; see Figs 2 and 3. Cells are colored according to year in which a given mutation was first detected in a ZIKV genome sequence, as shown in legend located in the lower right corner of the figure. Mutation data from Asian strain genomes are condensed to allow mutations to be grouped according to co-occurrence, prevalence, and persistence. S1 and S2 Figs include interactive versions of the data spreadsheet enabling individual cells to be queried by cursor position (as shown in brown inset adjacent to mutation queried). Readout includes: x mutation information, y sample identification, z # corresponding to color code (i.e. 0 no color, data missing), year mutation was first detected, or in cases where more than one mutation occurs at the residue “different than described” for example mutations at nts 3050, 3458, 6026, 8006, and 10265). X-axis labels include nucleotide location: codon: protein: codon mutation: number of genomes containing that mutation.

Fig 2. Overview of the dispersal patterns of nonsynonymous mutations within all sequences from Cluster 0.

Subsection of the graph from Fig 1 showing higher resolution of mutation patterns within the Cluster 0. Cluster 0 is defined by NS5 mutation Thr114Met (marked by a star) first detected in 2006. Mutation accession numbers, year, host type, country and number of mutations detected in genomes are listed on the y axis. Genomes listed in the y axis are listed in order that best illustrates mutation clusters. Y-axis labels include: Accession number. year. species. geographic location the mutation was identified: total number of mutations identified in that genome compared to the reference genome (Malaysia 1966).

Fig 3. Overview of dispersal patterns of nonsynonymous mutations within sequences from Cluster 1.

Subsections of the graph from Fig 1 showing higher resolution of mutation patterns. Cluster 1 is defined by NS5 mutation Thr/Met114Val (Fig 4E) first detected in 2014 (marked by a star). Fragments of clusters 1.1–1.4 with genomes listed in the y axis are selected to show the locations of sequences from microcephaly cases.

Fig 4. Epidemic-associated mutations tend to occur in exposed regions of Zika proteins.

Protein “hydrophobicity surface” representation of Zika proteins Pr, E, NS1, NS3, and NS5 using Kyte-Doolittle scale [47] and the coloring scheme from blue (hydrophilic) to red (hydrophobic), with white color representing neutral residues. For a given mutation, a difference in hydrophobicity score between initial residue and the mutation is indicated by a change in the color intensity used as text background. For example, light red for lower and dark red for higher hydrophobicity, and light blue for lower and dark blue for higher hydrophilicity. Exposed regions of listed mutations are circled and colored in green. Plots were constructed using Chimera software. (A) Model of the Pr protein with epidemic mutations first detected in 2013: Val1Ala, Ser17Asn. (B) E protein with high abundance epidemic mutation Val473Met (observed in 387 sequences; see Table 1) first detected in 2013, and two mutations, Asp393Glu and Thr487Met, first observed in pre-emergent samples from Thailand in 2006, Cambodia in 2010, Philippines in 2012 and French Polynesia in 2013. (C) Homodimer conformation of the NS1 protein (chain B colored in yellow). Mutations Met349Val, Gly100Ala, and Arg324Trp are from different subclusters: 1.1, 1.2, and 1.3, respectively (Fig 1, Table 2). Mutation Tyr122His first detected in Haiti 2014 and Brazil 2015, and a group of three exposed mutations and in close vicinity to each other with significant amino acid property change observed: Trp98Gly, Lys146Glu, and Val264Met. (D) NS3 protein showing 4 microcephaly associated mutations Thr566Ala, Thr567Ile, Met572Leu and Tyr584His, with predicted location in the interface with NS5, and His355Tyr mutation first detected in Haiti 2014, Brazil 2015, and not detected in French Polynesia. (E) NS5 protein showing critical outbreak mutation 114Val (forms Cluster 1; Fig 1 and Table 4), Ile322Val mutation with predicted location in the interface with NS3, and epidemic mutations Arg525Cys, Thr833Ala, and Asp878Glu from different subclusters: 1.2, 1.3, and 1.4, respectively (Fig 1 and Table 4).