GENOMIC ANALYSIS OF VIRAL OUTBREAKS

SEQUENCING METHODS

It is essential that patient samples be processed and sequenced in a way that will provide the highest quality data for downstream analysis. For RNA viruses, timing is especially important: degradation can occur quickly in clinical samples and is common, so the time between sample collection and sequencing should be minimized ((18); see (19) for procedures used in an EBOV diagnostic laboratory). In general, all sample processing should consider both sample preservation and researcher safety (20).

Sequencing itself has progressed from technologies tailored to a specific viral sequence to sequence-independent, high-throughput approaches. Amplicon-based sequencing is the most common sequence-dependent method and was used early in the Ebola outbreak (21). In that study, viral genomes were amplified in long (often ≥ 2 kb) overlapping fragments by reverse-transcriptase polymerase chain reaction (RT-PCR) with EBOV-specific primers; these fragments were then sequenced by Sanger sequencing. This method is popular for detecting and studying viruses because it is fast and can be used to amplify very small amounts of material.

The speed and accuracy of amplicon-based sequencing has made it an effective method for on-site sequencing even in remote field settings (22), and the resulting genomes are sufficient for pathogen identification and basic analysis. However, amplicon-based sequencing does have drawbacks. First, Sanger sequencing is not conducive to the deep coverage needed to detect low-frequency variants. Second, designing PCR primers requires prior knowledge of the viral sequence, which introduces bias and precludes metagenomic analysis. Third, it can be difficult to design primers that produce full-length genomes for all samples, given the high sequence diversity of many viruses. Lastly, degraded samples prevent full-length amplicon production necessary to obtaining whole-genome sequences.

High-throughput sequencing platforms resolve many of the issues of amplicon-based approaches. These platforms, which produce short reads, are better able to capture fragmented or partially degraded samples. They also allow for the ultra-deep sequencing needed to detect low-frequency within-host variants (23). Sequence independence is essential for studying outbreaks caused by new or unknown pathogens, and for metagenomic analysis. Instead of virus-specific primers, these methods rely on random priming followed by high-throughput sequencing (14, 16, 24). Combining sequence-independent primer amplification with selective RNase H-based digestion of contaminating RNA (mainly host ribosomal RNA) enables rapid, unbiased deep sequencing of viral samples, as was done during the Ebola virus epidemic (25).

Other high-throughput sequencing approaches can contribute to viral genomic analysis. Hybrid selection has been used to enrich the viral content of sequencing libraries with high host contamination even after RNase H digestion (17), and is an active area of development. Refining this technology will improve viral genomic analysis during outbreaks, when sample quality may be variable. Other potentially useful technologies still in development include long-read sequencing, which could allow for phasing of variants, and technologies optimized for rapid on-site sequencing. These cheap and portable approaches (26) are useful for rapid diagnostics, but have high error rates that may preclude some detailed genomic analysis.

SEQUENCE ASSEMBLY AND ALIGNMENT

After sequencing, care should be given to the assembly and alignment of genomes. Some of the necessary steps and best practices for processing high-throughput sequencing reads are shown in Figure 1 (see also Supplemental Table 1). After completing these steps, reads that do not map to the database of possible viruses (step ❷) can be investigated using one of several taxonomic analysis tools (27–29). Such reads can be de novo assembled and further investigated using a nucleotide or protein homology search. For a detailed example of how these methods were used to discover a novel flavivirus and two novel rhabodviruses, see (30, 31). Alternatively, comprehensive metagenomics pipelines (32, 33) can be used if rapid pathogen identification is the primary goal.

Recombination can affect downstream phylogenetic analysis, so the final sequence alignment should be screened for recombination. Many methods that check for recombination have been compiled into a single software package, RDP4 (34). As described below, there are alternative phylogenetic tools that should be used when recombination is present, but analysis of recombinant viral sequences is still an area of active development.

VARIANT CALLING

Sequence differences between viral genomes mark the evolutionary history and relationships between samples. Single-base substitutions (single nucleotide polymorphisms, or SNPs) are the simplest variants. Given high-quality consensus sequences aligned to a reference, it is relatively easy to manually identify SNPs. However, more complex approaches—such as those implemented in packages like GATK (35) or Samtools (36)—are helpful when samples contain insertions or deletions, or if regions of the genome have poor quality, low coverage, or high diversity.

Individual SNPs should be annotated—classified as nonsense, missense, or intergenic—and located relative to genes and other genomic elements. Many annotation tools are available online, each requiring only a list of SNPs and an annotated reference genome (37–39).

At this stage, it is also useful to identify variants within individual samples (intrahost variants, or iSNVs), indicating the presence of multiple viral quasispecies. Powerful tools exist for calling low-frequency variants in heterogeneous viral populations (40, 41). To avoid calling sequencing errors as iSNVs, we suggest discarding variant calls with fewer than five forward or reverse reads and those where the number of reads differs greatly between the forward and reverse strands ((25), see supplemental methods). Because PCR errors during library construction can introduce false variants, replicate libraries should be prepared and sequenced whenever possible to confirm the presence of within-host variants at comparable frequencies.

CONSTRUCTING A PHYLOGENETIC TREE

Phylogenetic trees can be constructed using maximum likelihood (44, 45) or Bayesian approaches (46). All methods require only a sequence alignment and a nucleotide substitution model. The nucleotide substitution model describes the rate at which one nucleotide is replaced by another, and is used to estimate the evolutionary distance between sequences. The model is used in calculating the likelihoods of various possible phylogenetic trees, and therefore may greatly affect results (47). A general time reversible (GTR) model is often used for phylogenetic analyses because it is the most general and makes no assumptions about substitution rates or base frequencies (48). Alternatively, several groups have written statistical software to compare substitution models for a given dataset (49).

When constructing or reading trees it is important to keep confidence values in mind. Confidence in maximum likelihood trees is commonly represented by bootstrap values (50). Bootstrapping estimates uncertainty by sampling from a dataset with replacement. In this case, the bootstrap value for a node is the proportion of bootstrap trees in which that particular branch topology occurs. Although there is some debate about the accuracy of bootstrap values (51), reporting these values, at least for important nodes, is common practice. Confidence values are built into Bayesian phylogenies and are the posterior probabilities. A Bayesian approach can be thought of as a faster version of a bootstrapped maximum likelihood approach, though the concordance between the two types of confidence values is variable (52).

Both maximum likelihood and Bayesian methods were used to determine the phylogeny of Ebola viruses sequenced during the outbreak (21, 25, 53). These two methods are often used together to check for agreement: major differences in the resulting trees may suggest a complex evolutionary relationship not fully captured by one or more methods.

ROOTING A PHYLOGENETIC TREE

Without further information, a phylogenetic tree will be unrooted: it will show the relationship of branches relative to one another and the overall topology, but it will not identify the base of the tree or the direction of evolution. It thus cannot tell you which samples are ancestors and which are descendants. Since ancestry is very important to determining the origin of an outbreak, it must be determined by rooting the phylogenetic tree.

There are two primary methods of rooting phylogenetic trees: mid-point rooting and outgroup rooting. Mid-point rooting is done by finding the longest tip-to-tip distance in the tree and setting the root half way between these tips. This method assumes that evolutionary rates are constant throughout the tree, meaning the root should be equidistant from all tips (it also assumes contemporaneous sampling). As discussed in the next section, this assumption (known as the molecular clock assumption) is often incorrect. Therefore, viral outbreaks are typically rooted by selecting an outgroup; that is, a set of sequences known to be more distantly related than anything else in the tree. This can be composed of published sequences from previous outbreaks of the same virus, virus sampled from another host species, or a closely related viral species. Outgroup genomes must be distinct from outbreak genomes, but rooting trees using highly divergent sequences can also be problematic (53). If only outbreak sequences are available, it is also possible to use a particularly divergent cluster of outbreak sequences as the outgroup, if one exists.

Once an outgroup is selected, the phylogenetic tree is reconstructed using these additional sequences; the root is the point of divergence between the outgroup sequences and the rest of the tree (Figure 2). In some cases the root of the tree is ambiguous, as in the recent Ebola virus outbreak. Dudas and Rambaut (53) explain how viral substitution rates and linear regression can be used to select themost likely root for a viral outbreak.

ESTIMATING THE START DATE OF AN OUTBREAK

Rooted trees suggest the infection history of an outbreak; outbreak sequences close to the root of the tree (separated by fewer nodes) came earlier in the outbreak. If sampling dates are available, branch lengths can be converted from units of nucleotide substitutions to units of time, which can be used to date the true origin of the outbreak. This is done using a strict molecular clock model, which assumes that nucleotide substitutions accumulate at a constant rate (54). The number of substitutions on each branch of a tree with dated tips can be used to estimate the substitution rate, which then can be used to extrapolate backwards to the date of origin of a particular outbreak strain (55). Maximum likelihood methods (56) can calculate substitution rates given a phylogenetic tree and sampling dates.

While a helpful simplification, the strict molecular clock does not always accurately model real viral evolution; evolutionary rates can vary over time, space, or between different branches. To address this, more flexible models have been developed that allow for variation in the substitution rate over time (57). The BEAST (Bayesian Evolutionary Analysis by Sampling Trees) package (58) implements a Bayesian Markov Chain Monte Carlo (MCMC) method to determine changing substitution rates over time; this is referred to as a relaxed molecular clock. This framework can be used to co-estimate the phylogeny and divergence times given sequence data and sampling dates.

During the Ebola outbreak, BEAST was used to estimate when outbreak viruses split from lineages documented in other outbreaks (53), and to estimate the date of entry of the virus into Sierra Leone from Guinea (25, 59).

DETERMINING THE CAUSE OF AN OUTBREAK

The same phylogenetic methods can be used to determine the type(s) of transmission causing an outbreak (human-human or animal-human), but this analysis requires sequences from appropriate hosts and/or time periods (Figure 3). For example, analysis of Ebola patient samples showed that there was substantial genetic variation between EBOV outbreaks, but limited variation within each outbreak. This suggested that the virus evolves separately in an animal reservoir, and that a single zoonotic transmission was responsible for the start of each outbreak. This hypothesis was supported by the calculation of divergence times calculated by BEAST: the lineages from the two most recent outbreaks diverged from a common ancestor significantly before the start of either outbreak (25).

CHALLENGES IN ESTIMATING THE ORIGIN OF AN OUTBREAK

Although phylogenetic tools have been used successfully to understand many viral outbreaks, significant challenges remain in correctly establishing the origin of an outbreak. A detailed review of current challenges in phylogenetic methods can be found in (60). Here we highlight those challenges particularly relevant for determining the origin of a viral outbreak.

First, although relaxed molecular clocks allow for some rate variation, current models may still fail to capture the full variation in evolutionary rates. For example, an analysis for pandemic HIV-1 group M found that the time to the most recent common ancestor varies significantly when subtypes are analyzed separately compared to jointly, perhaps because closely related viral lineages have different substitution rates (61).

Additionally, the methods described above cannot account for recombination that occurs in many viruses, since the ancestry of these viruses cannot be represented by a simple branching process. Instead, different parts of a single sample’s genome can be the product of different genealogical trees, and are better modeled by a phylogenetic network or ancestral recombination graph (ARG) that allows for complicated evolutionary relationships (62–64). Because many phylogenetic tools cannot account for recombination, it is important to restrict analysis to parts of the viral genome or tree where recombination is limited. Important advances in the field will come from continued development of phylogenetic tools that can incorporate viral recombination.

Lack of data also poses significant barriers to analyzing many viral outbreaks. For example, lack of sampling dates essentially rules out divergence time estimates, and lack of informative outgroup sequences—perhaps due to limited past sequence data, or because a zoonotic reservoir has yet to be determined—prevents accurate rooting of a phylogenetic tree. Non-random sampling over time or space may significantly bias results (60). Finally, understanding the ecological factors leading to an outbreak at a particular place and time requires detailed surveys of the outbreak location, both during and before its start (65, 66). Without detailed epidemiological surveys, it may be impossible to determine the index case of a viral outbreak.

TRANSMISSION RECONSTRUCTION

Reconstructed viral transmission routes can reveal modes of transmission that are important for containment and prevention (6). Depending on the depth of sampling, rooted trees either coarsely or in detail correspond to the transmission history of the virus (79–82) (Figure 4).

Formal methods have been developed to reconstruct transmission chains from genetic data (83); several of these methods combine genetic and epidemiological data (e.g. sampling dates and locations) into a single likelihood function that is used to sample possible transmission trees (79, 84–86). Jombart et al (85, 87) have developed an R package that constructs transmission trees from genetic and any available epidemiological data (Figure 4b).

Within-host genomic data has recently been recognized as an essential component of phylogenetic analysis and transmission tree reconstruction (55, 86, 88). One challenge of incorporating this type of data is that it requires an understanding of the characteristic within-host dynamics for a given virus before it can be used to effectively inform statistical and epidemiological models. Specifically, it is important to know the underlying viral mutation rate and the typical within-host substitution rate, as well as how much diversity is transmitted during an infection event (the bottleneck size).

Because studying within-host dynamics requires both high sequencing depth and longitudinal sampling, limited information exists for most viruses. Within-host studies are most common in well-studied chronic viral infections such as HIV (89, 90), although similar studies in other viruses are beginning to appear (91). In the same vein, the average size of the transmission bottleneck in HIV is known (92), but is still under investigation in most other viruses. Deep sequences of Ebola viruses published during the 2014–2015 outbreak suggest that the bottleneck size is greater than one (18, 93), but more precise estimates are still needed. Within-host viral dynamics studies, along with the development of robust phylogenetic methods that incorporate within-host variation, are a crucial next step in outbreak research.

CALCULATING THE REPRODUCTION NUMBER

The basic reproduction number (R₀)—the number of secondary cases from a single infection—is a useful measure of the infectivity of a pathogen and is usually estimated from epidemiological models. However, this number can also be estimated from a detailed transmission chain or from genomic data (94, 95). This value often frames the discussion about containment for a disease outbreak and can be used to predict outbreak dynamics and eventual size in the presence or absence of various control measures (in the Ebola outbreak: (69, 96)).

Calculation of epidemiological parameters like R₀ is part of the new and growing field of phylodynamics, the study of infectious disease behavior that arises from a combination of evolutionary and epidemiological processes (89). Incorporating epidemiological metadata can enhance genetic analysis and vice versa in outbreak situations.

CHALLENGES IN UNDERSTANDING TRANSMISSION

Although joint evolutionary and epidemiological analysis has greatly advanced the field of outbreak investigation, there are still challenges associated with determining the route and rate of viral spread. Major hurdles include sampling bias and the difficulty of allowing for spatial and temporal complexity in phylodynamic models (60). For example, Lloyd-Smith et al (97) highlight problems with using the same reproduction number for all individuals and the resulting implications for outbreak control.

THE MUTATION AND SUBSTITUTION RATES

The mutation rate is a major determinant of the overall rate of evolutionary change during and between outbreaks; it is the number of genetic mutations that occur per viral genome replication. It is largely determined by a virus’s biological properties, such as the fidelity of its polymerase, the speed at which it replicates its own genome, and whether the genome is RNA or DNA (98). In general, RNA viruses mutate fastest and DNA viruses slowest. The mutation rate must be measured experimentally because natural selection affects the number of mutations identified in genetic data. For mutation rate estimates for specific viruses, see Drake (99) and Drake and Hwang (100).

The mutation rate should not be confused with the viral substitution rate, which is the rate at which nucleotide substitutions accumulate in a viral lineage. This rate is determined by the mutation rate and by other factors, including natural selection and the effective viral population size. This is the rate most commonly discussed in an outbreak situation, both because it can be calculated from sequence data and because it can be used to understand selective pressures on a viral population during an outbreak. For example, it may be useful to compare the substitution rate within an outbreak to that in a zoonotic reservoir. The virus should have the same intrinsic mutation rate in both hosts, so differences in substitution rate could be due to selection.

While calculating the viral mutation rate requires careful experimentation, the substitution rate can be calculated given a phylogenetic tree and sampling dates. This can be done with maximum likelihood methods (56) or Bayesian methods (58). A Bayesian method like BEAST is more statistically rigorous than most maximum likelihood implementations because it can allow the substitution rate to vary between branches, but is computationally intensive (98). For approximate substitution rates for various viruses, see Jenkins et al (101).

The caveat to all of these methods is that they assume all substitutions are fixed in the population. Because many mutations on recent branches are mildly deleterious and will disappear from the population over time, the substitution rate for any tree containing recent samples may be artificially high. This is common during outbreaks, where a majority of viral genomes may be terminal branches.

RATE-BASED TESTS FOR SELECTION

During an outbreak, identification of variants that confer a fitness benefit to the virus is a key concern. Variants that affect fitness are likely under selection in a population. A well-known statistical test for selection is the d_N/d_S (or K_a/K_s or ω) test (Figure 5). This test compares the rates of synonymous and non-synonymous substitutions in some region of the genome. In the absence of selection, these two rates should be the same (corrected for the frequency of each type of site). Non-synonymous substitutions are more likely to affect the resulting protein, and are therefore more likely to influence fitness. Therefore, the ratio of synonymous to non-synonymous changes can indicate selection: d_N/d_S > 1 suggests positive selection, and d_N/d_S < 1 suggests negative or purifying selection.

This test requires only a codon alignment and is often applied to viral sequences to identify domains or genes under selection (102–104). However, it was originally developed to analyze sequences from divergent species, not to detect selection within a single population (105). The d_N/d_S ratio is not applicable within single populations, and the results obtained from this test are often misleading when applied to microbes (106). Therefore, while this ratio can still be used to analyze sufficiently divergent, separately-evolving outbreaks, users should be wary of using this test to detect selection within a single outbreak population.

If data from other outbreaks is available, the d_N/d_S statistic can be used to identify sites or regions under selection. During the Ebola epidemic, several groups found d_N/d_S > 1 in the glycoprotein gene (GP)—or more specifically, in the disordered, mucin-like region of GP—when including sequences from all known Ebola virus outbreaks (18, 107) (Figure 5b). This is unsurprising because GP encodes the envelope protein: as the only surface-exposed protein on the viral particle, GP is the target of host cell antibodies. Because of this biology, it is suspected that GP undergoes diversifying selection—or relaxed purifying selection—as a response to host immune pressure (108).

In general, comparing non-synonymous to synonymous mutations between sites or between timescales (i.e. comparing inter- and intrahost substitutions, as in (109)) can be used to suggest regions under selection. It is possible to identify selection using only synonymous substitutions by looking for regions of excess synonymous constraint in a virus (110) (Figure 5C). In viruses, many protein-coding regions contain overlapping or embedded functional elements. Because any type of substitution may disrupt an overlapping element, these regions are often characterized by an unusually low synonymous substitution rate. This method was used during the Ebola epidemic to find a constrained region in a known editing site.

OTHER TESTS FOR SELECTION

Tests based on mutation frequency spectra and tree topology can also be used to identify selection in single populations. The basic principle behind these methods is that selective pressure on a population leaves a distinct mark on overall genetic diversity and tree symmetry (89, 111, 112). Statistics that test for selection, or non-neutrality, using these principles include Tajima’s D, Fu and Li’s D, and tree-imbalance metrics (113).

Tajima’s D is a statistic that compares the average number of pairwise differences between sequences to the total number of variable sites within the set of sequences (114). A negative value of Tajima’s D indicates an excess of low frequency polymorphisms compared to a neutral model, and may be due to purifying selection or population size expansion. Conversely, a population size reduction (or bottleneck) and balancing selection keep variants at intermediate frequencies, which translates to a positive value for Tajima’s D. Fu and Li’s D (115) applies this idea to the phylogeny of these sequences and compares the number of mutations on newer, external branches to the total number of mutations on all branches. As with Tajima’s D, a basic understanding of population genetics can be used to interpret the result of this test. For example, under purifying selection, there is likely to be an excess of mutations on external branches because deleterious mutations are generally purged from the population before they can be passed on.

Three measures of tree imbalance are commonly used to test for selection: B_I (116), the cherry count (117), and Colless’s tree imbalance index (118). More asymmetry than expected in a phylogenetic tree suggests non-neutral evolution. These methods have been successfully used to understand selection in HIV, influenza, and other viruses (113, 119).

The major drawback of both frequency spectra and tree imbalance methods is that it can be hard to differentiate between selection and epidemiological effects such as changing population size. For example, a very negative Tajima’s D or Fu and Li’s D can be due to exponential growth rather than non-neutral evolution. Drummond and Suchard (113) address this problem by incorporating a demographic model when analyzing three RNA virus datasets, and show that it is possible to use these tests to identify selective pressure on viral populations.

CHALLENGES IN VIRAL EVOLUTION

Detecting selection in viruses is challenging because most statistical methods have been created for the comparison of divergent populations or species, rather than analysis of a single population that may be rapidly evolving and expanding. Additionally, viruses are very biologically diverse and have highly variable mutation and substitution rates. This makes it difficult to use the same selection tests for all viruses. For example, slow-mutating viruses, which usually have low substitution rates, may require extended sampling to achieve the population diversity needed to identify evolutionary trends. Unfortunately, not all viruses and datasets make good subjects for evolutionary analysis, and even when they do, the results may be relatively uninformative or uninteresting.