to El-Kafrawy et al., Middle East Respiratory Syndrome Coronavirus: Genetic diversity, phylogeography, infection rates and age associated patterns in imported African and local Arabian dromedary camels Supplementary text:
For initial MERS-CoV genome analyses we amplified three target regions by RT-PCR from 100 imported camel samples. These fragments were located in ORF1a (genome position 10,265 to 11,361 in EMC/2012 reference genome; Genbank Acc No. JX869059); the Spike gene (genome position 23,864 to 24,909); and ORF4b (genome position 26,025 to 27,020). Selection of samples for sequencing was based on availability of sample volume and viral RNA concentration as per initial real-time PCR testing (minimum, CT = 35). All RT-PCR primers and protocols were described before 1 . The three genome regions were selected to generate genomic data before and after known recombination breakpoint 2, 3 , and to include the region within ORF4b that was reported to show specific deletion patterns in African isolates 4 . The concatenated sequences of these fragments covered approximately 3140 nucleotides (~10% of a typical complete MERS-CoV genome with a size of ~30,120 nucleotides). After Sanger sequencing, we selected 22 samples to represent all clades observed in initial phylogenetic analyses using the concatenated fragments. Full genome sequences from 24 local camel samples were also generated.
Full genomes were generated using shotgun Illumina High-Throughput Sequencing (HTS) on all selected positive samples showing CT-values of less than 26. All other samples were first amplified by RT-PCRs as in Corman et al. 1 to amplify the complete genome with overlapping PCR amplicons, followed by PCR amplicon-based HTS. For HTS from clinical sample we used 5µL of extracted RNA for library preparation using the KAPA RNA Hyper Prep kit (Roche Molecular Diagnostics, Basel, Switzerland) according to manufacturer's instructions. Indexed DNA libraries were measured by Qubit dsDNA HS Assay kit (Thermo Fisher Scientific, Karlsruhe, Germany), pooled together at equimolar ratios (up to 20 libraries), and normalized. Sequencing was performed using the 600cycle MiSeq reagent v3 cartridge (Illumina, San Diego, California, U.S.) according to manufacturer's instructions. Genome regions with low coverage (<10) and sequencing gaps were amplified by bridging RT-PCR using MERS-CoV specific primers and sequenced by dye terminator chemistry (Microsynth Seqlab, Göttingen, Germany).
For the amplicon-based approach, we used a hemi-nested RT-PCR assay for amplifying and sequencing of complete genomes from samples with low RNA concentration. The RT-PCR assays were designed to amplify about 800 overlapping base pairs of all known MERS-CoV sequences and were published by our group before 1 . Following amplification, PCR amplicon bands were visualized by agarose gel electrophoresis (2% agarose) with ethidium bromide and purified using KAPA Pure Beads (Roche Molecular Diagnostics). All amplicons from each single sample were pooled and each pool was quantified. Approximately 5 ng of each DNA pool was fragmented, using KAPA Frag Kit (Roche Molecular Diagnostics). Library preparation was done using the KAPA Hyper Prep kit (Roche Molecular Diagnostics) according to the manufacturer's protocol. Indexed DNA libraries were, measured by Qubit dsDNA HS Assay kit (Thermo Fisher Scientific), pooled together at equimolar ratios (up to 50 libraries), and normalized. Sequencing was performed using the 600-cycle MiSeq reagent v3 cartridge (Illumina). The set-up of the instrument and the preparation of the flow cell and were performed according to manufacturer instructions.
The strategy for complete genome assembly comprised three major steps: First MERS-CoV reads in the HTS dataset were identified by using a local BLAST comparison against a database containing all MERS-CoV strains available in Genbank (as of 1st June 2019). Second, all reads matching this database were used for de-novo genome assembly by SPAdes within Geneious 9.1.8. Third, all HTS reads were used for specific reference-mapping against the de-novo assembled consensus sequences and MERS-CoV strain EMC/2012 in Geneious 9.1.8. After performing multiple iterations nearly complete consensus genomes could be generated for all samples. The denovo and the mapping approach did not result in discordant consensus sequences. Nevertheless, genome regions with low coverage (<10 reads), sequencing gaps, and deletions as compared to reference sequences from Genbank were amplified by bridging RT-PCR using MERS-CoV specific primers and sequenced by dye terminator chemistry (Microsynth Seqlab, Göttingen, Germany).
For sequence analysis the complete or partial sequences of the MERS-CoV from this study were aligned with other MERS-CoV strains available in GenBank. The alignments included representatives of MERS-CoV lineages representing MERS-CoV clades A and B 3 and all published clade C (non-A, non-B) MERS-CoV complete genomes, as of 1st Feb 2019 4 . Alignments were made using MAFFT 5 in Geneious R11 (www.geneious.com, Biomatters Ltd, New Zealand). Bayesian phylogeographic reconstructions were made using Beast 1.10.4 6 under assumption of a Hasegawa, Kishino and Yano (HKY) gamma 4 nucleotide substitution model as in Dudas et al. 7 a strict molecular clock, and a calibration of substitution rate based on recorded or database-derived sampling dates of taxa. Geographic locations were reconstructed at all tree nodes based on phylogeographic diffusion in discrete space using trait matrices that incorporate geographic centers of countries of origin (Morocco, Burkina Faso, Nigeria, Egypt, Sudan, Djibouti, Saudi Arabia, Jordan, United Arab Emirates). A Bayesian stochastic search variable selection (BSSVS) procedure was chosen, in order to limit the number of migration rates to only those that are necessary to explain the overall changes of location traits in light of the underlying phylogenetic tree structure 8 . Significance of migration rates is determined by toggled rate indicators for each migration rate in the trait change matrix, followed by Bayes factor analyses of posterior versus prior probability of the rate of interest activated. Phylogeographic rendering as well as Bayes factor analysis was done in SpreaD3 9 . Recombination analyses were made using SimPlot V3.5 10 and RDP V4.95 11 .
For analysis of effective reproductive number of MERS-CoV circulating in Saudi-Arabia we applied Bayesian Birth Death Skyline analyses in Beast2 12 using an alignment covering a non-recombinant part of the MERS-CoV genome (corresponding to position 140 to 15985 in the MERS-CoV/EMC/2012 reference strain) and containing all MERS-CoV sequences of clades A and B as available in GenBank as of 1st April 2019. In accordance with Dudas et al. 7 , we followed the rationale that human primary cases or unique viruses observed in outbreaks represent one-time samples from the ongoing enzootic circulation in camels. We removed all human cases except primary cases and outbreak index cases. For data published up to 2015, we used the same index sequences from human outbreaks as in Dudas et al. 7 . For sequences published in GenBank after 2015, index cases could often not be identified due to a lack of associated scientific publications. We therefore identified clusters of similar humanderived sequences that fell into timeframes compatible with outbreak settings (typically, phylogenetically clustered human samples taken within up to 3 months from each other) and removed all sequences except one representative from these potential human-to-human clusters. More than one sequence per time cluster was only left in the dataset when these sequences were considerably different from each other, suggesting no transmission link between humans (=these sequences were regarded as likely to represent independent index sequences). This curation of the dataset caused reduction of the original alignment from 261 to 178 sequences (refer to supplementary material, Table 2 ).
The following parameter and prior settings were chosen in Beast2: Based on Dudas et al. 7 , the MRCA of the sequences under analysis falls in the end of 2011. As the origin date in BDSKY analysis has to be older than the root age, its lower bound was set to 2372 days (latest date in our sample set, May 2018). The "become uninfectious" rate (the reciprocal of the duration of infection) was set to 0.05, based on Dudas et al. 7 (we note that the true infectious period is probably shorter because Adney et al. 13 , as used for reference in Dudas et al. 7 , infected naïve camels and detected infectious nasal shedding only up to day 7 while the prior setting implicates 20 days being infectious). Priors for the sampling proportion followed the consideration that our alignment contains 178 sequences sampled over 2106 days (=5.76 years). According to a projection in Dudas et al. 7 
There are approximately 14 million camels worldwide (1.6 million in KSA) of which 90% are dromedary camels. In KSA, dromedary camels are indispensible and fulfil several important domestic, commercial and nutritional demands. Camel trade is a large industry and local demands are met by importation of camels from outside KSA. Camels serve community, trading, farming, transportation, and leisure-related purposes and are also sources of milk, meat, wool and leather. The main source of dromedary camels for the Arabian Peninsula is the importation of life camels from the Horn of Africa (Somalia, Ethiopia, Sudan, Kenya and Djibouti), and North Africa (mainly Egypt, Libya) mainly through the ports of Jeddah and, to a much lesser extent, Jizan in Saudi Arabia 14 . In 2013, official records in Saudi Arabia reported the import of 131,932 camels representing more than 70% of the number of animals slaughtered in the country 15 . Live imported camels are distributed across the Kingdom.
The increased urbanization of the human population in KSA has led to changes in patterns of camel husbandry. A report by Abdallah et al. 16 described 6 types of camel farms and husbandry in Saudi Arabia; commercial farms, racing farms, use of camels for leisure for urban dwellers (weekend farms), traditional camel farms, camels for renting and camels for prestige. Camel herds are frequently moved across the Arabian Peninsula and may cross the borders of other Gulf countries for reasons such as grazing, participation in camel races and camel trade fairs. However, there is no official national identification procedure, nor is there an obligatory vaccination campaign focused on camels; it is difficult to quantify the extent of animal movements. Some general patterns are known. For instance, camels from both UAE and Qatar are moved every year during the winter season (November -February) for grazing in the Eastern region of Saudi Arabia. Detailed spatiotemporal studies of camel movement patterns are lacking, but camels are generally moved and traded over wide distances. This is also reflected in scientific data. For instance, a recent study observed little phylogeographic signal in genetic analyses of the modern dromedary population, indicative of extensive gene flow 17 .
The regular movement of camels across many countries may contribute to the transfer, evolution and amplification of MERS-CoV. For instance, a high degree of relatedness between MERS-CoV sequences obtained from UAE, Burydah/KSA and Al-Ahsa in camels and humans was seen 18 . MERS CoV has been detected in camels crossing into the UAE from Saudi Arabia 19, 20 . 