Genomic data offer new insights into epidemiologic and evolutionary dynamics, and sequencing pathogen samples is becoming increasingly routine, providing new insights for a range of issues of public health importance ([@kwx182C1]).  Pathogen genomic data allow us to determine the phylogeny of isolates, which in turn sheds light on the potential transmission networks between the hosts from whom they were collected.  Identifying routes of infection enables the estimation of risk factors for disease transmission, which in turn can inform the implementation of infection control strategies.  As such, inference of transmission trees using genomic data is an increasingly well-studied field ([@kwx182C2]--[@kwx182C8]).  Although low-resolution pathogen typing has been used for some time to discriminate between independent outbreaks ([@kwx182C9]--[@kwx182C11]), whole-genome sequencing provides additional resolution with which genetic distance between identical phenotypes may be ascertained ([@kwx182C12]--[@kwx182C14]).  However, this too has limits.  Studies have shown that while transmission clusters may be identified with genomic data, individual-level transmission routes can rarely be identified with a great degree of certainty ([@kwx182C3], [@kwx182C4]).  Characterizing an infected host by a single pathogen genome (isolation and purification of a single colony for bacteria, or using the consensus sequence for viral pathogens) is common practice yet neglects within-host diversity.  The variation in sampled genetic distances can be large relative to the expected number of mutations between hosts, rendering the number of single nucleotide polymorphisms a rather crude measure of relatedness on an individual level ([@kwx182C15]).  As such, particularly for rapidly evolving pathogens, or those whose mode of transmission is associated with a large and potentially diverse inoculum ("transmission bottleneck"---the number of pathogens transmitted in an infection event), single-genome sampling can cause hosts to appear misleadingly similar or dissimilar. 

Deep sequencing can potentially provide new insights into within-host diversity.  Currently, sequencing a mixed population sample to sufficient depth to identify minor nucleotide variants has mostly been limited to viral samples.  While consensus sequences may appear identical for 2 samples, comparing minor variants can offer additional resolution.  For instance, if the same nucleotide variation is observed at the same locus in pathogen samples from 2 individuals (henceforth referred to as a "shared variant" (SV)), this could be considered strong evidence for direct transmission, particularly if the variant is not observed in any other host.  This naturally relies on the possibility that a pathogen population of size greater than 1 survives the transmission bottleneck; otherwise, each infection must initially be monoclonal, implying that any variation found within distinct hosts must have arisen independently.  There is evidence for larger bottlenecks occurring in viral pathogens such as influenza ([@kwx182C16], [@kwx182C17]) and Ebola ([@kwx182C18]), and it is plausible that this is also the case for bacterial pathogens ([@kwx182C19]). 

The connection between SV presence and direct transmission has previously been suggested.  Gire et al. ([ @kwx182C20]) noted the presence of SVs in Ebola virus samples from individuals who were potentially linked by transmission.  Data collected from 2 influenza A animal transmission studies were used to explore the presence of SVs between hosts, and it was shown that such data were consistent with known contact patterns ([@kwx182C21]).  This study used known contact patterns to identify characteristics of SVs that were more likely to be associated with transmission, allowing variants to be split into those consistent and inconsistent with transmission, minimizing false connections.  Poon et al. ([ @kwx182C22]) identified routes of influenza transmission occurring during a household contact study using both consensus whole-genome sequence data and the presence of SVs.  In the case of bacterial pathogens, the diversity in *Staphylococcus aureus* infections, which can be considerable, has been linked to transmission in a veterinary hospital ([@kwx182C23]). 

Pathogens vary considerably in their bottleneck size, mutation rates, and transmission dynamics.  It remains unclear how methods based on SVs are expected to perform in different regions of this parameter space.  Establishing this is a crucial component of the interpretation of SVs and the value of the approach. 

In the present study, we investigated the predictive power of SVs for identifying transmission routes.  In addition to pathogen genomes, many other sources of data may contribute information towards inference of transmission routes, including temporal and spatial data, contact patterns, and expression of symptoms.  However, here we aimed to examine the information contributed by genomic data alone and, in particular, the additional benefit offered by considering SVs. 

METHODS {#kwx182s2}
=======

We generated infectious disease outbreaks with within-host pathogen evolution under various mutation rates and bottleneck sizes by simulation, using the R (R Foundation for Statistical Computing, Vienna, Austria) package Seedy ([@kwx182C24]).  In contrast to other simulation approaches, this method explicitly simulates within-host evolution and allows sampling of mixed pathogen populations across time.  We expanded upon methods previously used to infer transmission routes using deep sequence data ([@kwx182C20]--[@kwx182C22]), comparing their performance with analogous genetic distance--based approaches.  We additionally proposed hybrid approaches that combine SV information with phylogenetic distance data.  We considered the following approaches: Weighted variant tree: For each host, we weighted potential sources by the number of observed SVs, such that the host sharing the largest number of variants is attributed the largest weight.  Hosts sharing no variants with any other were not assigned a source.  Weighting edges provided an extension to previous approaches ([@kwx182C20]--[@kwx182C22]).Maximum variant tree: For each host, we defined the source to be the individual with whom the largest number of SVs was observed.  Hosts sharing no variants with any other host were not assigned a source. Weighted distance tree: Using consensus sequences, the genetic distance (number of single nucleotide polymorphisms) between isolates was calculated, and potential sources were weighted inversely by this metric.  This approach has been described previously ([@kwx182C3]).Minimum distance tree: Using consensus sequences, the source of a given host was defined to be the carrier of the genetically closest isolate to that of the host.  This approach, with the incorporation of sampling times to provide directionality, has been described previously ([@kwx182C7]).Hybrid weighted tree: First, the weighted variant tree was constructed.  Hosts with no source were then assigned potential sources based upon weighted genetic distance. Hybrid maximum tree: First, the maximum variant tree was constructed.  Hosts with no source were then assigned potential sources based upon minimum genetic distance. 

These 6 simple heuristics by no means comprise an exhaustive list of approaches to identify routes of transmission but are instead a range of readily implemented, distance-based approaches that require neither knowledge of evolutionary dynamics nor infection or sampling times.  As has previously been demonstrated, simple methods based on genomic data alone can provide powerful insights into transmission dynamics ([@kwx182C25]).  Further details of the approaches used here, as well as the metrics used to assess the accuracy of tree reconstruction and reliability of estimated transmission routes, are provided in [Web Appendix 1 (](#sup1){ref-type="supplementary-material"}available at <https://academic.oup.com/aje>).  We additionally applied SV methods to previously published data collected during the Ebola virus outbreak in West Africa in 2014 ([@kwx182C20], [@kwx182C26]). 

RESULTS {#kwx182s3}
=======

Simulation studies {#kwx182s3a}
------------------

As expected, the proportion of cases in which a SV was observed in at least 1 other host increased rapidly with mutation rate and bottleneck size (Figure [1](#kwx182F1){ref-type="fig"}A).  The majority of SVs were observed in exactly 2 individuals, with the proportion shared among larger groups declining rapidly as the size of the group increased (Figure [1](#kwx182F1){ref-type="fig"}B).  For each simulation, we constructed a weighted transmission tree according to the 6 methods outlined previously.  An example simulated outbreak of a pathogen with characteristics similar to *S. aureus* (see [Web Appendix 1](#sup1){ref-type="supplementary-material"}) is shown in Figure [2](#kwx182F2){ref-type="fig"}, along with reconstructions based upon two of these methods.  While many edges are bidirectional and symmetric, asymmetry can occur under most methods due to the lack of commutativity (i.e., even if B is the closest host to A, A may not be the closest host to B). 

![ Summary of genetic variant frequency across simulated outbreaks.  We simulated 10 outbreaks for each combination of 6 bottleneck sizes and 3 mutation rates (180 in total).  A) Mean proportion of cases in the outbreak who shared at least 1 variant with another host.  B) Distribution of shared-variant group size for different bottleneck sizes, with a mutation rate of 5 × 10^−4^ per genome per generation.](kwx182f01){#kwx182F1}

![ Simulated and reconstructed transmission trees.  A simulated tree (A) was generated with bottleneck size 10 and mutation rate 0.001 per genome per generation.  Cases are represented as numbered circles linked by arrows representing transmission routes.  Based on the simulated genetic data from this outbreak, trees were reconstructed according to the maximum variant approach (B) and the minimum genetic distance approach (C) described in "Methods."  Arrows denoting unambiguous routes (weight = 1) are solid, while dashed arrows denote edges with weight \<1.  Networks were plotted with the igraph package in R (Foundation for Statistical Computing) ([@kwx182C38]).](kwx182f02){#kwx182F2}

We used two metrics to assess the reliability of individual estimated transmission routes and compared the different methods described above.  First, we considered the true path distance between inferred transmission pairs.  We found that under the maximum variant tree, the mean path distance was typically less than 2, outperforming the minimum distance approach (Figure [3](#kwx182F3){ref-type="fig"}A).  Second, we examined the mean weight assigned to the true source of each host.  In the case of small (\<5) bottleneck sizes, methods based on SVs perform poorly because the likelihood of a monoclonal infection is high, resulting in most true links being assigned a weight of zero.  Furthermore, those links that are inferred by SVs in small-bottleneck settings identify direct transmission with high confidence.  The hybrid approaches perform best for small bottlenecks, incorporating SV information when available but not relying upon it.  For larger bottlenecks, the distance-based approaches were markedly outperformed by the variant-based approaches (Figure [3](#kwx182F3){ref-type="fig"}B).  As the rate of mutation increases, SV approaches outperform those based on distance alone to an increasing degree as mutation generates increasing diversity in the infecting population ([Web Figure 1A--1C](#sup1){ref-type="supplementary-material"}). 

![ Reliability of estimated transmission routes.  A) A comparison of the minimum distance (black) and the maximum variant (gray) approaches.  The true path distance between estimated transmission pairs gives insight into the extent to which transmission links are misspecified.  A perfect reconstruction would have mean path length of 1; greater values indicate increasing misspecification.  Maximum variant path lengths are averaged over identified transmission pairs---that is, excluding hosts with no shared variants.  B) The mean weight assigned to transmission links for each tree reconstruction method, under a range of scenarios and methodologies.  Results are shown for a mutation rate of 5 × 10^−4^ per genome per generation and were averaged across 10 simulated outbreaks for each scenario.](kwx182f03){#kwx182F3}

In addition to the reliability of individual links, we also considered the overall accuracy of a transmission tree reconstruction.  This was measured by the area under the receiver operating characteristic curve (AUC) statistic.  For small bottlenecks, variant-based methods provide a poor tree reconstruction by this metric (Figure [4](#kwx182F4){ref-type="fig"}).  We note that here we regard assigning no source (e.g., where no SVs exist) as an incorrect classification, leading to values below 0.5, the value expected under random source assignment.  A tight bottleneck leads to little diversity persisting across transmission events, and as such, SVs are rarely observed, leading to a sparsity of informed links across the network.  However, larger bottlenecks led to a considerably better performance under this measure for the variant-based approaches, which exceeded even the weighted distance approach with a sufficiently large bottleneck size and mutation rate (Figure [4](#kwx182F4){ref-type="fig"}, [Web Figure 1D--1F](#sup1){ref-type="supplementary-material"}).  In contrast, distance-based approaches typically declined in accuracy as the bottleneck size increases, for reasons that are well understood ([@kwx182C3], [@kwx182C27])

![ Transmission tree reconstruction accuracy.  The area under the receiver operating characteristic curve (AUC) metric provides an overall measure of network accuracy: Values closer to 1 indicate better overall reconstruction.  Results for a mutation rate of 5 × 10^−4^ per genome per generation are shown here.  The mean AUC across 10 simulations was measured for each scenario.](kwx182f04){#kwx182F4}

We additionally investigated the effect of "mutational hotspots," which can generate potentially confounding homoplasy.  We found that while variant approaches performed less well, they generally continued to outperform distance-based approaches for larger (\>10) bottleneck sizes (see [Web Appendix 1, Web Figures 2 and 3](#sup1){ref-type="supplementary-material"} for further details). 

Ebola virus data {#kwx182s3b}
----------------

We next examined previously published data sets for Ebola in Sierra Leone, for which sequence data are available, and we can determine the presence and properties of intrahost variants.  In order to reduce the risk of counting variant calling errors as true intrahost variants, we identified only variants in which the minor within-host allele frequency was at least 5% (routes estimated under a 1% threshold are shown in [Web Figure 4](#sup1){ref-type="supplementary-material"}).  Figure [5](#kwx182F5){ref-type="fig"} shows the transmission trees reconstructed for each data set under the weighted variant approach, using no epidemiologic information.  In the first data set (Figure [5](#kwx182F5){ref-type="fig"}A, ([@kwx182C26])) 19 of 78 hosts were found to share a variant with at least 1 other individual.  Four pairs of patients shared more than 1 variant (3 pairs with 2 SVs and 1 pair with 4), while 1 additional pair shared 1 unique variant.  Consistent with transmission, each of these pairs originated from the same geographic location, and permutation testing revealed that this geographic similarity was significantly higher than would be expected via random selection (*P* = 0.0075).  Pairs were also temporally clustered; 3 of these links were sampled 2 or fewer days apart, while the remaining 2 were sampled 12 and 22 days apart, which is plausible given serial interval estimates for Ebola virus infection of 15.3 ± 9.3 days ([@kwx182C28]).  Under the minimum distance tree, 2 of these pairs were reproduced, 2 pairs belonged to a much larger group of samples with identical consensus sequences, and 1 pair, differing by a single nucleotide according to consensus sequences, remained unconnected due to the presence of other identical sequences ([Web Figure 5](#sup1){ref-type="supplementary-material"}). 

![ Estimated Ebola transmission routes, Sierra Leone, 2014.  Transmission links between sampled hosts during an Ebola outbreak under the maximum variant approach.  Node colors denote distinct geographic regions to which hosts belong, while arrow color and thickness correspond to the relative weight attributed to each potential transmission event.  Variant detection threshold was 5%.  Cases are shown from the first ([@kwx182C20]) (A) and the second ([@kwx182C26]) (B) data sets.  Further details are provided in [Web Appendix 1](#sup1){ref-type="supplementary-material"}.](kwx182f05){#kwx182F5}

While a consistent result of our simulations was the sharing of variants among small numbers of hosts, rarely more than 2, in the Ebola data collected by Gire et al. ([ @kwx182C20]), 1 variant was shared by 11 hosts.  Samples with this variant are highly clustered geographically (10 of 11 in the same chiefdom, *P* = 0.022) and temporally (observed within an 18-day period), as well as phylogenetically, lending support to this group representing a transmission cluster. 

In the second data set (Park et al. ([ @kwx182C26])) 26 of 150 (for which replicate sequencing and variant calling was performed) shared a variant with at least 1 other host (Figure [5](#kwx182F5){ref-type="fig"}B).  There were 5 pairs of individuals sharing a unique variant.  As before, 1 variant was shared by multiple hosts, but unlike the previous data set these were not geographically or temporally clustered, coming from different regions and spanning several weeks.  Furthermore, while some of these samples clustered on the phylogenetic tree, many fell in different clades ([@kwx182C26]), suggesting the group is unlikely to represent a single transmission cluster but rather multiple transmission events in combination with homoplasy. 

DISCUSSION {#kwx182s4}
==========

We have described some simple methods for reconstructing transmission trees using SVs, testing how well this approach performs for a range of parameters governing the rates of diversification within and between hosts.  We have then applied the methods to data from the 2014 Ebola outbreak to identify links, using the genomic data alone, that are likely to be consistent with transmission given time and location. 

For the great majority of parameter space, excluding only very low mutation rates and tight bottlenecks, these methods outperform genetic-distance comparison methods, which have increasingly been used to identify potential transmission events ([@kwx182C7], [@kwx182C29], [@kwx182C30]).  The limitations of distance-based methods that characterize a single genome are well appreciated.  We note that although, for the purpose of comparison, additional data sources were not included in our inference of transmission routes, incorporating these independently would be a relatively straightforward step with these methods.  Most simply, sampling dates could be used to provide directionality to inferred connections. 

The additional information we derive from SVs can inform the transmission tree in two distinct ways, depending on the region of parameter space.  First, small but nonsingular bottlenecks (e.g., for airborne influenza transmission ([@kwx182C31]) or sexually transmitted HIV ([@kwx182C32])) produce few inferred transmission pairs, but these are highly accurate.  The small bottleneck means that the probability of observing a SV between individuals who are in the same transmission cluster, but not directly linked, is negligible.  Second, SV data for pathogens with larger transmission bottlenecks (e.g., Ebola ([@kwx182C18]), influenza transmitted via contact ([@kwx182C31]), or intravenous drug--associated HIV transmission ([@kwx182C33])) provide good information on the overall tree structure and transmission clusters, but individual links may be more uncertain.  In all cases, higher mutation rates make the observation of SVs more probable, which typically results in better inference of transmission routes. 

A hybrid approach that combines SVs and the sequence of either an individual, sequenced genome, or the consensus offers substantial benefit in the case of small bottleneck sizes (\<5), where we predict a method based on SV alone would struggle.  Because transmission routes are assessed independently of one another, estimated transmission trees frequently comprise several unconnected nodes or clusters.  Such unconnected clusters could be linked to one another if further structure is required, using the weighted distance approach on pooled within-cluster samples.  Here we have simply used genetic distance, which is predicted to be efficient and reliable under the relatively short time scale of an outbreak, but more sophisticated models of sequence evolution could be applied. 

We applied these methods to Ebola data collected from Sierra Leone in 2014.  While the first data set is thought to represent relatively dense coverage of the initial stages of the epidemic in the country, with around 70% of cases sampled ([@kwx182C20], [@kwx182C34]), the later data set comprised a sparser sample.  While this reduces the expected number of linked cases in the data set, the reliability of transmission routes identified via SVs remains largely unaffected ([Web Figures 6 and 7](#sup1){ref-type="supplementary-material"}).  As such, while only relatively few transmission routes were identified in the data sets, this is likely a function of both the proportion of missing data and the relatively low mutation rate of Ebola virus ([@kwx182C20], [@kwx182C35]).  Confidence in the transmission pairs identified was reinforced by investigating temporal and geographic clustering, which proved to be significant, and while the aim of our study was to assess the accuracy of transmission-route identification via genomic data alone, methodology combining spatial and temporal data sources will naturally provide further insight.  Identifying even a small proportion of direct-transmission pairs can be of great interest in studying pathogen-level transmission dynamics as well as conducting outbreak investigations. 

Studying the Ebola data showed that both data sets contained a large group of hosts sharing the same variant, which was rare in all our simulations.  The observation can be explained in at least two ways: recurrent mutation (as might arise through selection) or an anomalously large number of contacts with large bottleneck size (such as might be associated with a funeral-based exposure).  Park et al. ([ @kwx182C26]) suggest that the large group in the second data set likely arose through a combination of patient-to-patient transmission and recurrent mutation.  Subsets of this group do cluster on the phylogenetic tree, and having identified clusters of hosts with SVs, we might partition these groups by genetic background (e.g., ruling out transmission between hosts with \>1-nucleotide difference present across nonvariant sites). 

Sample contamination may be an additional source of error.  Cross-contamination may potentially lead to SVs observed between unlinked hosts.  However, in many settings we do not believe this would present a major concern.  If no minor variants are contained in the contaminating sample, SVs would not link this to the contaminated sample.  Still, it remains important to verify that observed SVs are consistent with transmission and to minimize the risk of contamination as much as possible.  The second Ebola data set highlights the potential for false-positive connections to be identified between spatiotemporally inconsistent pairs.  SVs observed across several samples with differing genetic backgrounds may be indicative of contamination.  While such rules can eliminate obvious sources of error, further work is required to formally evaluate the risk of contamination based on deep sequence data. 

Another deliberate simplification in the present work is the assumption of neutral evolution.  While this is plainly faulty over longer time scales, over the relatively short timescale of an outbreak it is a first approximation, and this is supported by real data from outbreaks ([@kwx182C20], [@kwx182C26]) and even longer periods ([@kwx182C36]) showing evidence of incomplete purifying selection.  Selection may not, however, have as severe an effect on these methods as we might assume.  If a specific variant is maintained through balancing selection, it is likely to be found in multiple hosts and, as a result, will be less informative as to specific transmission links; if several hosts are connected by the same SV, this will be misleading only if no additional variants are observed.  In contrast, diversifying selection by the host immune system is expected to produce the mutational hotspots we have studied here, which again have little impact.  A similar argument can be made that sequencing errors will be less important than expected, because they are likely to be found in just one sample and hence be uninformative as to links to other samples.  A more formal approach to this problem would be to test for selection and down-weight the identified loci from the analysis. 

As yet, there are still few studies in which adequate data have been collected in order to use SVs as a feature to identify transmission routes.  Deep and high-quality sequencing is required to reliably call minor variants, as well as dense sampling of the outbreak population such that the majority of infection sources are included in the study population.  It is likely that such data will become more commonly collected in the near future, for both viral and bacterial pathogens, as the associated sequencing costs fall and the benefits become more evident.  This work demonstrates that deep sequence data can be informative in outbreak investigations and epidemiologic studies, and it should motivate both the wider collection of such data and the further development of methodologies that might accommodate scenarios that stray further from neutral evolutionary dynamics.  It is noticeable that bottleneck size in nature, as opposed to minimal infectious dose, has not received the attention it deserves.  The importance of this parameter for these methods, as well as other factors such as the evolution of virulence ([@kwx182C37]), should motivate further study. 

We have demonstrated the power of deep sequencing data to identify transmission routes with greater resolution than analogous methods using the genome of a single isolate.  While homoplasy and contamination might generate false-positive results and, in some settings, might be relatively common, consideration of additional data sources (temporal, spatial, and similarity of genetic background) can sometimes identify and rule out such cases.  Rigorous collection of epidemiologic data remains a crucial component of outbreak investigation, and combining this with deep sequencing and SV analysis can provide unprecedented insight into individual-level transmission dynamics.  The development of models that suitably incorporate all data sources remains an important goal. 

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

###### 

Click here for additional data file. 

Author affiliations: Center for Communicable Disease Dynamics, Department of Epidemiology, Harvard T.H. Chan School of Public Health, Boston, Massachusetts (Colin J. Worby, Marc Lipsitch, William P. Hanage).  C.J.W. is currently at the Department of Ecology and Evolutionary Biology, Princeton University, Princeton, New Jersey. 

Research reported in this paper was supported by the National Institute of General Medical Sciences of the National Institutes of Health (grant U54GM088558). 

The content is solely the responsibility of the authors and does not necessarily represent the official views of the National Institute of General Medical Sciences or the National Institutes of Health.  The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript. 

Conflict of interest: none declared.