MERS-CoV spillover at the camel-human interface.

Middle East respiratory syndrome coronavirus (MERS-CoV) is a zoonotic virus from camels causing significant mortality and morbidity in humans in the Arabian Peninsula. The epidemiology of the virus remains poorly understood, and while case-based and seroepidemiological studies have been employed extensively throughout the epidemic, viral sequence data have not been utilised to their full potential. Here, we use existing MERS-CoV sequence data to explore its phylodynamics in two of its known major hosts, humans and camels. We employ structured coalescent models to show that long-term MERS-CoV evolution occurs exclusively in camels, whereas humans act as a transient, and ultimately terminal host. By analysing the distribution of human outbreak cluster sizes and zoonotic introduction times, we show that human outbreaks in the Arabian peninsula are driven by seasonally varying zoonotic transfer of viruses from camels. Without heretofore unseen evolution of host tropism, MERS-CoV is unlikely to become endemic in humans.


Introduction
Middle East respiratory syndrome coronavirus (MERS-CoV), endemic in camels in the Arabian Peninsula, is the causative agent of zoonotic infections and limited outbreaks in humans. The virus, first discovered in 2012 (Zaki et al., 2012;Boheemen et al., 2012), has caused more than 2000 infections and over 700 deaths, according to the World Health Organization (WHO) (World Health Organization, 2017). Its epidemiology remains obscure, largely because infections are observed among the most severely affected individuals, such as older males with comorbidities (Assiri et al., 2013a;The WHO MERS-CoV Research Group, 2013). While contact with camels is often reported, other patients do not recall contact with any livestock, suggesting an unobserved community contribution to the outbreak (The WHO MERS-CoV Research Group, 2013). Previous studies on MERS-CoV epidemiology have used serology to identify factors associated with MERS-CoV exposure in potential risk groups (Reusken et al., 2015(Reusken et al., , 2013. Such data have shown high seroprevalence in camels (Müller et al., 2014;Corman et al., 2014;Chu et al., 2014;Reusken et al., 2013Reusken et al., , 2014 and evidence of contact with MERS-CoV in workers with occupational exposure to camels (Reusken et al., 2015;Müller et al., 2015). Separately, epidemiological modelling approaches have been used to look at incidence reports through time, space and across hosts (Cauchemez et al., 2016).
Although such traditional epidemiological approaches yield important clues about exposure patterns and potential for larger outbreaks, much inevitably remains opaque to such approaches due to difficulties in linking cases into transmission clusters in the absence of detailed information. Genomic epidemiology, however, can fill this critical gap and has repeatedly shown the utility of viral sequence data in outbreak scenarios (Liu et al., 2013;Gire et al., 2014;Grubaugh et al., 2017), where data are relatively cheap to produce. These data can stand in for diagnostics and often yield a highly detailed picture of an epidemic when complete genome sequencing is performed consistently and appropriate metadata collected (Dudas et al., 2017). Sequence data can help discriminate between multiple and single source scenarios (Gire et al., 2014;Quick et al., 2015), which are fundamental to quantifying risk (Grubaugh et al., 2017). Sequencing MERS-CoV has been performed as part of initial attempts to link human infections with the camel reservoir (Memish et al., 2014), nosocomial outbreak investigations (Assiri et al., 2013b) and routine surveillance (Park et al., 2015). A large portion of MERS-CoV sequences come from outbreaks within hospitals, where sequence data have been used to determine whether infections were isolated introductions or were part of a larger hospital-associated outbreak (Fagbo et al., 2015). Similar studies on MERS-CoV have taken place at broader geographic scales, such as cities .
It is widely accepted that recorded human MERS-CoV infections are a result of at least several introductions of the virus into humans  and that contact with camels is a major risk factor for developing MERS, per WHO guidelines (World Health Organization, 2016). Previous studies attempting to quantify the actual number of spillover infections have either relied on traditional epidemiological approaches (Cauchemez et al., 2016) or employed methods agnostic to signals of population structure within sequence data (Zhang et al., 2016). Here we use a dataset of 274 MERS-CoV genomes to investigate transmission patterns of the virus between humans and camels.
Here, we use an explicit model of metapopulation structure and migration between discrete subpopulations, referred to here as demes , derived from the structured coalescent (Notohara, 1990). Unlike approaches that model host species as a discrete phylogenetic trait of the virus using continuous-time Markov processes (or simpler, parsimony based, approaches) Lycett et al., 2016), population structure models explicitly incorporate distinct sampling patterns, population dynamics within demes and migration between demes. By estimating independent coalescence rates for MERS-CoV in humans and camels, as well as migration patterns between the two demes, we show that long-term viral evolution of MERS-CoV occurs exclusively in camels. Our results suggest that spillover events into humans are seasonal and might be associated with the calving season in camels. Once human MERS-CoV infections are established, however, we find that MERS-CoV is poor at transmitting between humans. Using Monte Carlo simulations we show that R 0 for MERS-CoV circulating in humans is much lower than the epidemic threshold of 1.0 and that correspondingly the virus has been introduced into humans hundreds of times.

MERS-CoV is predominantly a camel virus
The structured coalescent approach we employ (see Methods) identifies camels as a reservoir host where most of MERS-CoV evolution takes place (Figure 1), while human MERS outbreaks are transient and terminal with respect to long-term evolution of the virus ( Figure S1). Across 174 MERS-CoV genomes collected from humans, we estimate a median of 56 separate camel-to-human cross-species transmissions (95% highest posterior density (HPD): 48-63). While we estimate a median of 3 (95% HPD: 0-12) humanto-camel migrations, the 95% HPD interval includes zero and we find that no such migrations are found in the maximum clade credibility tree ( Figure 1). Correspondingly, we observe substantially higher camel-to-human migration rate estimates than human-to-camel migration rate estimates ( Figure S2). This inference derives from the tree structure wherein human viruses appear as clusters of highly related sequences nested within the diversity seen in camel viruses, which themselves show signicantly higher diversity and less clustering. This manifests as different rates of coalescence with camel viruses showing a scaled effective population size N e τ of 3.49 years (95% HPD: 2.71-4.38) and human viruses showing a scaled effective population of 0.24 years (95% HPD: 0.14-0.34).
We believe that the small number of inferred human-to-camel migrations are induced by the migration rate prior, while some are derived from phylogenetic proximity of human sequences to the apparent "backbone" of the phylogenetic tree. This is most apparent in lineages in early-mid 2013 that lead up to sequences comprising the MERS-CoV clade dominant in 2015, where owing to poor sampling of MERS-CoV genetic diversity from camels the model cannot completely dismiss humans as a potential alternative host. The first sequences of MERS-CoV from camels do not appear in our data until November The vast majority of MERS-CoV evolution is inferred to occur in camels (orange) with human outbreaks (blue) representing evolutionary dead-ends for the virus. Confidence in host assignment is depicted as a colour gradient, with increased uncertainty in host assignment (posterior probabilities close to 0.5) shown as grey. While large clusters of human cases are apparent in the tree, significant contributions to human outbreaks are made by singleton sequences, likely representing recent cross-species transmissions that were caught early.
2013. Our finding of negligible human-to-camel transmission is robust to choice of prior ( Figure S2).
The repeated and asymmetric introductions of short-lived clusters of MERS-CoV sequences from camels into humans leads us to conclude that MERS-CoV epidemiology in humans is dominated by zoonotic transmission (Figure 1 and S1). We observe dense terminal clusters of MERS-CoV circulating in humans that are of no subsequent relevance to the evolution of the virus. These clusters of presumed human-to-human transmission are then embedded within extensive diversity of MERS-CoV lineages inferred to be circulating in camels, a classic pattern of source-sink dynamics. Our analyses recover these results despite sequence data heavily skewed towards non-uniformly sampled human cases and are robust to choice of prior. This suggests that instances of human infection with MERS-CoV are more common than currently thought, with exceedingly short transmission chains mostly limited to primary cases that might be mild and ultimately not detected by surveillance or sequencing.

MERS-CoV shows seasonal introductions
We use the posterior distribution of MERS-CoV introduction events from camels to humans (Figure 1) to model seasonal variation in zoonotic transfer of viruses. We identify four months (April, May, June, July) when the odds of MERS-CoV introductions are increased ( Figure 2) and four when the odds are decreased (August, September, November, December). Camel calving is reported to occur from October to February (Almutairi et al., 2010), with rapidly declining maternal antibody levels in calves within the first weeks after birth (Wernery, 2001). It is possible that MERS-CoV sweeps through each new camel generation once critical mass of susceptibles is reached (Martinez-Bakker et al., 2014), leading to a sharp rise in prevalence of the virus in camels and resulting in increased force of infection into the human population. Strong influx of susceptibles and subsequent sweeping outbreaks in camels may explain evidence of widespread exposure to MERS-CoV in camels from seroepidemiology (Müller et al., 2014;Corman et al., 2014;Chu et al., 2014;Reusken et al., 2013Reusken et al., , 2014. Although we find evidence of seasonality in zoonotic spillover timing, no such relationship exists for sizes of human sequence clusters ( figure 2B). This is entirely expected, since little seasonality in human behaviour that could facilitate MERS-CoV transmission is expected following an introduction. Similarly, we do not observe any trend in human sequence cluster sizes over time ( Figure 2B, Spearman ρ = 0.06, p = 0.68), suggesting that MERS-CoV outbreaks in humans are neither growing nor shrinking in size. This is not surprising either, since MERS-CoV is a camel virus that has to date, experienced little-to-no selective pressure to improve transmissibility between humans.

MERS-CoV is poorly suited for human transmission
Structured coalescent approaches clearly show humans to be a terminal host for MERS-CoV, implying poor transmissibility. However, we wanted to translate this observation into an estimate of the basic reproductive number, R 0 , which is more familiar to epidemiologists and provides insight into epidemic behavior. The parameter R 0 determines expected number of secondary cases in a single infections as well as the distribution of total cases that result from a single introduction event into the human population (Equation 1, Methods). We estimate R 0 along with other relevant parameters via Monte Carlo simulation in two steps. First, we simulate case counts across multiple outbreaks totaling 2000 cases using Equation 1 and then we subsample from each case cluster to simulate sequencing of a fraction of cases. Sequencing simulations take place at different levels of bias, wherein bias enriches sequencing of larger case clusters. This is a particularly pressing issue, since a priori we expect large hospital outbreaks of MERS to be overrepresented in sequence data, whereas sequences from primary cases will be sampled exceedingly rarely. We record the mean, median and standard deviation of sequence cluster sizes in each simulation (left-hand  Figure 3) and extract the subset of Monte Carlo simulations in which these summary statistics fall within the 95% highest posterior density observed in the empirical MERS-CoV data from structured coalescent analyses. We record R 0 values, as well as the number of case clusters (equivalent to number of zoonotic introductions), for these empirically matched simulations. A schematic of this Monte Carlo produre is shown in Figure S3. Generally, higher R 0 results in fewer larger transmission clusters, while lower R 0 results in many smaller transmission clusters ( Figure 3).

panels in
We find that observed phylogenetic patterns of sequence clustering strongly support R 0 below 1.0 (middle panels in Figure 3). For increasing levels of bias mean R 0 values observed in matching simulations are 0.844, 0.730, and 0.683, respectively. While the 95% percentiles for R 0 values are close to 1.0 (0.720-0.995) for the unbiased sequencing simulation (i.e. uniform sequencing efforts, in which every case is equally likely to be sequenced), we also note that increasing levels of bias are considerably more to likely to generate MERS-CoVlike sequence clusters ( Figure 3). Under unbiased sequencing only 0.6% of simulations fit our phylogenetic observations, while 2.7% and 2.7% of simulations fit for bias levels of 2.0 Each row corresponds to a different bias value used to concentrate the subsampling of sequences from cases, and goes from 1 (no bias) to 2, and 3 (increasing levels of bias which make large case clusters to be more likely to be sequenced). Leftmost scatter plots show results of individual Monte Carlo simulations, coloured by the R 0 value used for the simulation. The dotted rectangle identifies the 95% highest posterior density bounds for sequence cluster size mean and standard deviation observed for empirical MERS-CoV data. The distribution of R 0 values found within the dotted rectangle is shown in the middle, on the same y-axis across all levels of bias. Bins falling inside the 95% percentiles are coloured by R 0 , as in the leftmost scatter plot. The distribution of total number of introductions associated with simulations matching MERS-CoV sequence clusters is shown in the plots on the right, on the same y-axis across all levels of bias. Darker shade of grey indicates bins falling within the 95% percentiles. These Monte Carlo simulations indicate R 0 for MERS-CoV is likely to be below 1.0, with biased sequencing and numbers of zoonotic transmissions numbering in the hundreds. and 3.0, respectively. Correspondingly, we estimate 10% support for a model with bias level 1.0, 45% support for a model with bias level 2.0, and 45% support for a model with bias level 3.0. Model averaging would suggest plausible R 0 values between 0.57 and 0.91.
Lower values for R 0 in turn suggest relatively large numbers of zoonotic transfers of viruses into humans (right-hand panels in Figure 3). The median number of cross-species introductions observed in simulations matching empirical data without bias are 339 (95% percentiles 255-431). These numbers jump up to 558 (95% percentiles 424-720) for bias = 2 and 650 (95% percentiles 481-850) for bias = 3 simulations, which as mentioned previously match considerably better to MERS phylogenetic data. Model averaging would suggest plausible numbers of introductions between 299 and 818. Our results suggest a large number of unobserved MERS primary cases. Given this, we also performed simulations where the total number of cases is doubled to 4000 to explore the impact of dramatic underestimation of MERS cases. In these analyses R 0 values tend to decrease even further as numbers of introductions go up, although very few simulations match currently observed MERS-CoV sequence data ( Figure S4).
Overall, our analyses indicate that MERS-CoV is poorly suited for human-to-human transmission, with an estimated R 0 < 0.91 and sequence sampling likely to be biased towards large hospital outbreaks. Given these findings, and the fact that large outbreaks of MERS occurred in hospitals, the combination of frequent spillover of MERS-CoV into humans and occasional outbreak amplification via poor hygiene practices (Assiri et al., 2013b;Chen et al., 2017) appear sufficient to explain observed epidemiological patterns of MERS-CoV.

Recombination shapes MERS-CoV diversity
Recombination has been shown to occur in all genera of coronaviruses, including MERS-CoV (Lai et al., 1985;Makino et al., 1986;Keck et al., 1988;Kottier et al., 1995;Herrewegh et al., 1998). In order to explore the role of recombination in shaping MERS-CoV genetic diversity we used two recombination detection approaches across partitions of taxa corresponding to inferred MERS-CoV clades. Both methods rely on sampling parental and recombinant alleles within the alignment, although each quantifies different signals of recombination. One hallmark of recombination is the ability to carry alleles derived via mutation from one lineage to another, which appear as repeated mutations taking place in the recipient lineage, somewhere else in the tree. The PHI (pairwise homoplasy index) test quantifies the appearance of these excessive repeat mutations (homoplasies) within an alignment (Bruen et al., 2006). Another hallmark of recombination is spatial clustering of alleles along the genome, due to how template switching, the primary mechanism of recombination in RNA viruses, occurs. 3Seq relies on the spatial structure of nucleotide similarities between sequence triplets -two potential parent-donors and one candidate offspring-recipient sequences (Boni et al., 2007).
Both tests can give spurious results in cases of extreme rate heterogeneity and sampling over time (Dudas and Rambaut, 2016), but both tests have not been reported to fail simultaneously. PHI and 3Seq methods consistently identify most of the apparent 'backbone' of the MERS-CoV phylogeny as encompassing sequences with evidence of recombination ( Figure S5). Neither method can identify where in the tree recombination occurred, but each full asterisk in Figure S5 should be interpreted as the minimum partition of data that still captures both donor and recipient alleles involved in a recombination event. This suggests a non-negligible contribution of recombination in shaping existing MERS-CoV diversity. As done previously (Dudas and Rambaut, 2016), we show large numbers of homoplasies in MERS-CoV data ( Figure S6) with some evidence of spatial clustering of such alleles. Although the evolutionary centrality of camel viruses (Figure 1) may be sufficient to argue that camels are the host where MERS-CoV recombines, incidence of MERS-CoV is known to be much higher in camels (Müller et al., 2014;Corman et al., 2014;Figure 4. Recombinant features of MERS-CoV phylogenies. A) Heatmap showing the posterior probability that viruses from trees of different genomic fragments belong to the same introduction event. Genomic fragment 1 represents nucleotide positions 1 to 21000 and genomic fragment 2 represents nucleotide positions 21001 to 29364. The same set of viruses are arrayed on the x-axis and on the y-axis; the x-axis shows identity of these viruses along genomic fragment 1 and the y-axis shows identity of these viruses along genomic fragment 2. Viruses are ordered by their appearance in tree of genomic fragment 2 reduced to just the human tips and coloured by inferred host (blue for human, orange for camel) on the left. Human clusters are largely well-supported as monophyletic and consistent across trees of both genomic fragments. B) Phylogenies derived from genomic fragment 1 (left) and genomic fragment 2 (right), reduced to just the human tips. Taxa descended from the same introduction in fragment 1 are connected to their counterparts in fragment 2 with combinations of line dotting and colour unique to each introduction clade. Branches are coloured by inferred ancestral host state (human in blue, camel in orange). Human clusters in blue change phylogenetic positions between the trees together, with little incongruence within clusters. This is evidence for recombinant viruses generated in the camel reservoir entering human populations. Chu et al., 2014;Reusken et al., 2014;Ali et al., 2017). This provides ideal conditions for co-infection with distinct genotypes, which is a pre-requisite for detectable RNA virus recombination to occur.
Conversely, our results strongly suggest that co-infection of humans with distinct lineages of MERS-CoV should be exceedingly rare. We find little evidence that recombination will interfere with the inference of human outbreak clusters ( Figure 4A). Between the two trees in figure 4B four (out of 54) 'human' clades are expanded where either singleton introductions or two-taxon clades in fragment 2 join other clades in fragment 1. For the reverse comparison there are five 'human' clades (out of 53) in fragment 2 that are expanded. All such clades have low divergence (figure 4B) and thus incongruences in human clades are more likely to be caused by differences in deme assignment rather than actual recombination. And while we observe evidence of distinct phylogenetic trees from different parts of the MERS-CoV genome ( Figure 4B), human clades are minimally affected and large portions of the posterior probability density in both parts of the genome are concentrated in shared clades ( Figure S7). Critically, we observe the same source-sink dynamics between camel and human populations in trees constructed from separate genomic fragments as were observed in the original full genome tree (Figures 1, 4B). Observed departures from strictly clonal evolution suggest that while recombination is an issue for inferring MERS-CoV phylogenies, its effect on the human side of MERS outbreaks is minimal, as expected. MERS-CoV evolution on the reservoir side, though complicated by recombination, is nonetheless still amenable to phylogenetic methods, in part through limited diversity of the virus in camels (see next section). In humans MERS-CoV evolution should be far easier to track as the only detectable and problematic recombinants are more likely to arise within the transmission chain, than through human co-infection with distinct MERS-CoV lineages.

MERS-CoV shows population turnover in camels
Here we attempt to investigate MERS-CoV demographic patterns in the camel reservoir. We supplement camel sequence data with a single earliest sequence from each human cluster, treating viral diversity present in humans as a sentinel sample of MERS-CoV diversity circulating in camels. This removes conflicting demographic signals sampled during human outbreaks, where densely sampled closely related sequences from humans could be misconstrued as evidence of demographic crash in the viral population.
Despite lack of convergence, neither of the two demographic reconstructions show evidence of fluctuations in the relative genetic diversity (N e τ ) of MERS-CoV over time ( Figure 5). We recover a similar demographic trajectory when estimating N e τ over time with a skygrid tree prior across the genome split into ten fragments with independent phylogenetic trees to account for confounding effects of recombination (Figures 5, S8). However, we do note that estimates of relative genetic diversity are relatively low overall with a mean estimate of N e τ of 3.49 years (95% HPD: 2.71-4.38), and consequently MERS-CoV phylogeny resembles a ladder often seen in human influenza A virus phylogenies (Bedford et al., 2011). This estimate of N e τ can be translated into an estimate of effective prevalence N e given an estimate of generation time or serial intervals. Given a 20 day duration of infection in Here, we compare this estimate of effective prevalence to a rough estimate of census prevalence. Given extremely high seroprevalence estimates within camels in Saudi Arabia (Müller et al., 2014;Corman et al., 2014;Chu et al., 2014;Reusken et al., 2013Reusken et al., , 2014, we expect camels to usually be infected within their first year of life. Therefore we can estimate the rough number of camel infections per year as the number of calves produced each year. We find there are 830 000 camels in Saudi Arabia (Abdallah and Faye, 2013) and that female camels in Saudi Arabia have an average fecundity of 45% (Abdallah and Faye, 2013). Thus, we expect 830 000 × 0.50 × 0.45 = 186 750 new calves produced yearly and correspondingly 186 750 new infections every year, which spread over 10 day serial intervals gives an average of N = 5042 infections.
We thus arrive a ratio of N/N e = 39.1 (31.1-50.4). This is less than the equivalent ratio for human measles virus (Bedford et al., 2011) and is in line with the expectation from neutral evolutionary dynamics plus some degree of transmission heterogeneity (Volz et al., 2013) and seasonal troughs in prevalence. Thus, we believe that the ladder-like appearance of the MERS-CoV tree can likely be explained by entirely demographic factors.

MERS-CoV epidemiology
In this study we aimed to understand the drivers of MERS coronavirus transmission in humans and what role the camel reservoir plays in perpetuating the epidemic in the Arabian peninsula by using sequence data collected from both hosts (174 from humans and 100 from camels). We showed that currently existing models of population structure  can identify distinct demographic modes in MERS-CoV genomic data, where viruses continuously circulating in camels repeatedly jump into humans and cause small outbreaks doomed to extinction (Figures 1, S1). This inference succeeds under different choices of priors for unknown demographic parameters ( Figure S2) and in the presence of strong biases in sequence sampling schemes (Figure 3). From sequence data we identify at least 50 zoonotic introductions of MERS-CoV into humans from the reservoir (Figure 1), from which we extrapolate that hundreds more such introductions must have taken place (Figure 3). We also looked at potential seasonality in MERS-CoV spillover into humans. Our analyses indicated a period of three months where the odds of a sequenced spillover event are increased, with timing consistent with an enzootic amongst camel calves ( Figure  2). As a result of our identification of large and asymmetric flow of viral lineages into humans we also find that the basic reproduction number for MERS-CoV in humans is well below the epidemic threshold ( Figure 3).
Strong population structure in viruses often arises through limited gene flow, either due to geography (Dudas et al., 2017), ecology (Smith et al., 2009) or evolutionary forces (Turner et al., 2005;Dudas et al., 2015). On a smaller scale population structure can unveil important details about transmission patterns, such as identifying reservoirs and understanding spillover trends and risk, much as we have done here. When properly understood naturally arising barriers to gene flow can be exploited for more efficient disease control and prevention, as well as risk management.

Transmissibility differences between zoonoses and pandemics
Severe acute respiratory syndrome (SARS) coronavirus, a Betacoronavirus like MERS-CoV, caused a serious epidemic in humans in 2003, with over 8000 cases and nearly 800 deaths. Since MERS-CoV was also able to cause significant pathogenicity in the human host it was inevitable that parallels would be drawn between MERS-CoV and SARS-CoV at the time of MERS discovery in 2012. Although we describe the epidemiology of MERS-CoV from sequence data, indications that MERS-CoV has poor capacity to spread human-to-human existed prior to any sequence data. SARS-CoV swept through the world in a short period of time, but MERS cases trickled slowly and were restricted to the Arabian Peninsula or resulted in self-limiting outbreaks outside of the region, a pattern strongly indicative of repeat zoonotic spillover. Infectious disease surveillance and control measures remain limited, so much like the SARS epidemic in 2003 or the H1N1 pandemic in 2009, zoonotic pathogens with R 0 > 1.0 are probably going to be discovered after spreading beyond the original location of spillover. Even though our results show that MERS-CoV does not appear to present a serious global threat, we would like to highlight that its epidemiology is nonetheless concerning.
MERS-CoV may join the list of pathogens able to jump species barriers but not spread efficiently in the new host, but every system is distinct. Pathogens such as Bacillus anthracis, Andes hantavirus (Martinez et al., 2005), monkeypox (Reed et al., 2004), triple reassortant and H3N2v influenza A viruses (Shinde et al., 2009;Epperson et al., 2013) belong to this list and yet only triple reassortant viruses eventually contributed to a pandemic, due to epidemiology and phylodynamics confined to influenza A virus. When information about a system is lacking either because it is recent or entirely novel or if the system is inherently unpredictable, sequence data can and should be used jointly with whatever data sources are available to provide the unique pathogen perspective on an infectious disease outbreak. Although MERS-CoV seems to cause self-limiting outbreaks in humans with no evidence of worsening outbreaks over time, sustained evolution of the virus in the reservoir and flow of viral lineages into humans provides the substrate for a more transmissible variant of MERS-CoV to possibly emerge. We encourage continued genomic surveillance of MERS-CoV in the camel reservoir and from sporadic human cases.

Sequence data
All MERS-CoV sequences were downloaded from GenBank and accession numbers are given in Table S1. Sequences without accessions were kindly shared by Ali M. Somily, Mazin Barry, Sarah S. Al Subaie, Abdulaziz A. BinSaeed, Fahad A. Alzamil, Waleed Zaher, Theeb Al Qahtani, Khaldoon Al Jerian, Scott J.N. McNabb, Imad A. Al-Jahdali, Ahmed M. Alotaibi, Nahid A. Batarfi, Matthew Cotten, Simon J. Watson, Spela Binter, and Paul Kellam prior to publication. Fragments of some strains submitted to GenBank as separate accessions were assembled into a single sequence. Only sequences covering at least 50% of MERS-CoV genome were kept. Sequences were annotated with available collection dates and hosts, designated as camel or human and aligned with MAFFT (Katoh and Standley, 2013) and edited manually. Protein coding sequences were extracted and concatenated, reducing alignment length from 30130 down to 29364. The final dataset consisted of 174 genomes from human infections and 100 genomes from camel infections (Table S1).

Structured coalescent analyses
For our primary analysis, the MultiTypeTree module ) of BEAST v2.4.3 (Bouckaert et al., 2014 was used to specify a structured coalescent model with two demes -humans and camels. This model estimates four parameters: rate of coalescence in human viruses, rate of coalescence in camel viruses, rate of migration from the human deme to the camel deme and rate of migration from the camel deme to the human deme. Analyses were run on codon position partitioned data with two separate HKY+Γ 4 (Hasegawa et al., 1985;Yang, 1994) nucleotide substitution models specified for codon positions 1+2 and 3. A relaxed molecular clock with branch rates drawn from a lognormal distribution (Drummond et al., 2006) was used to infer the evolutionary rate from date calibrated tips. Default priors were used for all parameters except for migration rates between demes for which an exponential prior with mean 1.0 was used. All analyses were run for 200 million steps across ten independent Markov chains (MCMC runs) and states were sampled every 20 000 steps. Due to the complexity of multitype tree parameter space 50% of states from every analysis were discarded as burn-in, convergence assessed in Tracer v1.6 and states combined using LogCombiner distributed with BEAST v2.4.3 (Bouckaert et al., 2014). Three chains out of ten did not converge and were discarded altogether. This left 70 000 states on which to base posterior inference.
As a secondary analysis to test robustness to choice of prior, we set up an analysis where we increased the mean of the exponential distribution prior for migration rate to 10.0. All other parameters were identical to the primary analysis and as before 10 independent MCMC chains were run. In this case, two chains did not converge. This left 80 000 states on which to base posterior inference.
As an additional secondary analysis, we also considered a scenario where alignments were split into two fragments (fragment 1 comprised of positions 1-21000, fragment 2 of positions 21000-29364), with independent clocks, trees and migration rates, but shared substitution models and deme population sizes. Fragment positions were chosen based on consistent identification of the region around nucleotide 21000 as a probable breakpoint by GARD (Pond et al., 2006) by previous studies into SARS and MERS coronaviruses (Hon et al., 2008;Dudas and Rambaut, 2016). All analyses were set to run for 200 million states, subsampling every 20 000 states. Chains not converging after 200 million states were discarded. 20% of the states were discarded as burn-in, convergence assessed with Tracer 1.6 and combined with LogCombiner. Three chains out of ten did not converge. This left 70 000 states on which to base posterior inference.

Introduction seasonality
We extracted the times of camel-to-human introductions from the posterior distribution of multitype trees. This distribution of introduction times was then discretised as follows: for sample k = 1, 2, . . . , L from the posterior, Z ijk was 1 if there as an introduction in month i and year j and 0 otherwise. We model the variable Y ij = L k=1 Z ijk with the hierarchical model: σ m ∼ Cauchy(0, 2.5); α j ∼ Normal(µ y , σ y ); µ y ∼ Normal(0, 1); σ y ∼ Cauchy(0, 2.5).
Odds ratios of introductions can then be computed for each month i as OR i = exp(β i ).

Epidemiological analyses
As of 10 May 2017, the World Health Organization has been notified of 1952 cases of MERS-CoV. We thus simulate final transmission chain sizes using equation 1 (Lloyd- Smith et al., 2005;Blumberg and Lloyd-Smith, 2013) until we reach an epidemic comprised of 2000 cases. 10 000 simulations were run for 121 uniformly spaced values of R 0 across the range [0.5-1.1] with dispersion parameter ω fixed to 0.1. Each simulation results in a vector of outbreak sizes C, where C i is the size of the i th transmission cluster and K i=1 C i = 2000, where K is the number of clusters generated. We sample from the case cluster size vector according to a multivariate hypergeometric distribution to simulate sequencing (algorithm 1). The resulting sequence cluster size vector S contains K entries, some of which are zero (i.e. case clusters not sequenced), but K i=1 S i = 174 which reflects the number of MERS-CoV sequences used in this study. Since the sampling scheme operates under equi-probable sequencing we also simulated biased sequencing by using concentrated hypergeometric distributions where the probability mass function is squared (bias = 2) or cubed (bias = 3) and then normalized. This makes clusters likely to be 'sequenced' even more likely to be sequenced and vice versa. We performed a smaller set of simulations with 2500 replicates and twice the number of cases, i.e. K i=1 C i = 4000, to explore a dramatically underreported epidemic.
Let R 0 be the basic reproductive number (0 < R 0 < 1), and ω > 0 be the a dispersion parameter, then the probability of observing a stuttering chain (cluster) r of size j is (Blumberg and Lloyd-Smith, 2013) P r(r = j|R 0 , ω) = Γ(ωj + j − 1) Γ(ωj)Γ(j + 1) Although case clusters and their sizes are difficult to infer directly and require detailed epidemiological follow up, sequence data have fewer limitations.Our structured coalescent analyses indicate that every MERS outbreak is a contained cross-species transmission of the virus from camels into humans. The distribution of the number of these cross-species transmissions and their sizes thus contain information about the underlying distribution of case clusters. We employ a Monte Carlo simulation approach to identify simulations where the recovered distribution of sequence cluster sizes fall within the 95% highest posterior density intervals for three summary statistics of MERS-CoV sequence cluster sizes recovered via structured coalescent analyses: mean, median and standard deviation. These values are 2.90-3.70 for mean sequence cluster size, 4.87-6.07 for standard deviation of sequence cluster sizes and a median sequence cluster size of 1.
Data: Array of case cluster sizes in outbreak C = [C 1 , C 2 , . . . , C K ], sequences available M , total outbreak size K i=1 C i . Result: Array of sequence cluster sizes sampled: S = [S 1 , S 2 , . . . , S K ]. Draw S i from a hypergeometric distribution with C i successes, Compute the probability mass function (pmf) for all possible values of S i , p = [p(0) bias , p(1) bias , . . . , p(C i ) bias ] × ( i p bias i ) −1 , where p(·) is the pmf for a hypergeometric distribution with C i successes, K i=1 C i − C i failures after M trials; Draw a sequence cluster size S i from array of potential sequence cluster sizes [0, 1, . . . , C i ] according to p; end Add remaining sequences to last sequence cluster C K = M − S K−1 ; Algorithm 1: Multivariate hypergeometric sampling scheme. Pseudocode describes the multivariate hypergeometric sampling scheme that simulates sequencing. Probability of sequencing a given number of cases from a case cluster depends on cluster size and sequences left (i.e. "sequencing capacity"). The bias parameter determines how probability mass function of the hypergeometric distribution is concentrated.

Demographic inference of MERS-CoV in the reservoir
In order to infer the demographic history of MERS-CoV in camels we used the results of structured coalescent analyses to identify introductions of the virus into humans. The oldest sequence from each cluster introduced into humans was used to represent a random draw from the diversity of MERS-CoV circulating in camels. These sequences were combined with existing sequence data from camels to give us a dataset with minimal demographic signal coming from epidemiological processes in humans. Sequences belonging to the outgroup clade where most of MERS-CoV sequences from Egypt fall were removed out of concern that MERS epidemics in the Arabian Peninsula and Egypt are distinct epidemics with relatively poor sampling in the latter. A flexible skygrid tree prior (Gill et al., 2013) was used to recover estimates of relative genetic diversity (N e τ ) at 50 evenly spaced grid points across six years, ending at the most recent tip in the tree (2015 August) in BEAST v1.8.4 (Drummond et al., 2012), under a relaxed molecular clock with rates drawn from a lognormal distribution (Drummond et al., 2006) and codon position partitioned (positions 1 + 2 and 3) HKY+Γ 4 (Hasegawa et al., 1985;Yang, 1994) nucleotide substitution models. We set up five independent MCMC chains to run for 500 million states, sampling every 50 000 states. This analysis suffered from poor convergence, where two chains converged onto one stationary distribution, two to another and the last chain onto a third stationary distribution, with high effective sample sizes. Demographic trajectories recovered by the two main stationary distributions are very similar and differences between the two appear to be caused by convergence onto subtly different tree topologies. This non-convergence effect may have been masked previously by the use of all available MERS-CoV sequences from humans which may have lead MCMC towards one of the multiple stationary distributions.
To ensure that recombination was not interfering with the skygrid reconstruction we also split our MERS-CoV alignment into ten parts 2937 nucleotides long. These were then used as separate partitions with independent trees and clock rates in BEAST v1.8.4 (Drummond et al., 2012). Nucleotide substitution and relaxed clock models were set up identically to the whole genome skygrid analysis described above (Drummond et al., 2006;Hasegawa et al., 1985;Yang, 1994). Skygrid coalescent tree prior (Gill et al., 2013) was used jointly across all ten partitions for demographic inference. Five MCMC chains were set up, each running for 200 million states, sampling every 20 000 states.

Data availability
Sequence data and all analytical code is publicly available at github.com/blab/structuredmers.
acknowledge the contribution of the following scientists for sharing MERS-CoV sequence   human 2013-06-12 Figure S1. Evolutionary history of MERS-CoV partitioned between camels and humans. This is the same tree as shown in Figure 1, but with contiguous stretches of MERS-CoV evolutionary history split by inferred host: camels (top in orange) and humans (bottom in blue). This visualisation highlights the ephemeral nature of MERS-CoV outbreaks in humans, compared to continuous circulation of the virus in camels. Figure S2. Posterior migration rate estimates for two choices of prior. Negligible flow of MERS-CoV lineages from humans into camels is recovered regardless of prior choice. Plots show the 95% highest posterior density for the estimated migration rate from the human deme into the camel deme (blue) and vice versa (orange). Dotted lines indicate exponential priors specified for migration rates, with mean 1.0 (bottom) or 10.0 (top).  Figure S3. Monte Carlo simulation schematic. Case clusters are simulated according to Equation 1 until an outbreak size of 2000 cases is reached. We sample 174 cases from each simulation to represent sequencing of human MERS cases. 'Sequencing' is carried out by using multivariate hypergeometric sampling, representing sampling cases without replacement to be sequenced. Sequencing simulations take place at three levels of bias: 1.0, where every case is equally likely to be sequenced, and 2.0 and 3.0, where cases from larger clusters are increasingly more likely to be sequenced. The distribution of simulated sequence clusters is summarised by its mean, median and standard deviation. A simulation is considered to match if the mean, median and standard deviation of its sequence cluster sizes falls within the 95% highest posterior density interval of observed MERS-CoV sequence clusters. R 0 values that ultimately generate data matching empirical observations, as well as associated numbers of 'introductions' are retained as estimates. These estimates are summarised in Figure 3. Figure S4. Results of Monte Carlo simulations with vast underestimation of cases. The plot is identical to Figure 3, but instead of 2000 cases, simulations were run with 4000 cases. With more unobserved cases the R 0 values matching observed MERS-CoV sequence clusters can only be smaller, with a corresponding increase in numbers of zoonotic transmissions. However, the numbers of simulations that match MERS-CoV data go down as well. Figure S5. Tests of recombination across MERS-CoV clades. Maximum clade credibility tree of MERS-CoV genomes annotated with results of two recombination detection tests (PHI and 3Seq) applied to descendent sequences of each clade. Both tests identify large portions of existing sequence data as containing signals of recombination. Note that markings do not indicate where recombinations have occurred on the tree, merely the minimum distance in sequence/time space between recombining lineages. Figure S6. MERS-CoV genomes exhibit high numbers of non-clonal loci. Ancestral state reconstruction (right) identifies a large number of sites in which mutations have occurred more than once in the tree (homoplasies, orange) or are reversions (red) from a state arising in an ancestor. Mutations that apparently only occur once in the tree (synapomorphies) are shown in grey. The maximum likelihood phylogeny on the left is coloured by whether sequences were sampled in humans (blue) or camels (orange). Figure S7. Human clade sharing between genomic fragments 1 and 2. Central scatter plot shows the posterior probability of human clades shared between genomic fragments 1 and 2, in their respective trees. Left and bottom scatter plots track the posterior probability of human clades only observed in fragment 2 (left) or fragment 1 (bottom). The cumulative probability of human clades present in either tree are tracked by plots on the right (fragment 2) and top (fragment 1). Most of the probability mass is concentrated within human clades that are present in trees of both genomic fragment 1 and 2 (0.9701 and 0.9474 of all human clades across posteriors, respectively). Figure S8. Skygrid comparison between whole and fragmented genomes. Inferred median N e τ recovered using a skygrid tree prior on whole genome (bottom) and ten genomic fragments with independent trees (left), coloured by time. Dotted line indicates the one-to-one line.