Members of Coronaviridae family of viruses have started to gain a substantial interest due to their potential role as causative agents of emerging infections in humans (Fehr and Perlman, 2015) . This was manifested by the 2002-2003 SARS outbreak, 2012 MERS outbreak, and the current coronavirus disease 2019 pandemic, the first documented coronavirus pandemic, which can be viewed as the full-blown consequence of coronavirus threat (Cherry and Krogstad, 2004; Lu and Liu, 2012; Peiris et al., 2003) .
The causative agent of this unprecedented pandemic is an enveloped RNA virus of the subfamily Betacoronavirinae . Similar to other RNA viruses, severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2), is presumed to have a relatively high mutation rate that is mostly related to its RNAdependent RNA polymerase, with minimal proofreading activity (Duffy et al., 2008; Sevajol et al., 2014) . In addition, the high frequency of recombination in coronaviruses augments its genetic diversity and its ability of cross-species transmission (Su et al., 2016; Woo et al., 2009) .
The aforementioned features are accompanied by ubiquitous presence of coronaviruses in various animal reservoirs (Guan et al., 2003) . Thus, cross-species transmission, including spread to humans seems an inevitable outcome (Graham and Baric, 2010; Woo et al., 2009) . This is mainly related to human, ecologic and economic factors, which explain the increased frequency of zoonosis (Delabouglise et al., 2017; Karesh et al., 2012; Morse et al., 2012) .
The transcripts of SARS-CoV-2 include nine sub-genomic RNAs, of which one structural protein, spike glycoprotein, being responsible for attachment of the virus to its cellular receptor (angiotensin-converting enzyme 2 [ACE2]) (Fehr and Perlman, 2015) . Host proteases' cleavage of the spike glycoprotein is essential for virion entry into the target cells (Ou et al., 2020) . The receptor-binding domain (RBD) in the S1 subunit binds ACE2 and facilitates fusion with host cell membrane (Tai et al., 2020) . For the S2 domain of the spike glycoprotein, its function facilitates fusion of the viral and host cell membranes (Xia et al., 2020) . is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint
The copyright holder for this this version posted August 25, 2020 . . https://doi.org/10.1101 Studying the SARS-CoV-2 S gene attracts a special attention, particularly from an immunologic and evolutionary points of view (Chen et al., 2020; Korber et al., 2020; Robson, 2020) . The spike gene is the subject of an immune selective pressure, and antibodies against its protein product can inhibit the viral entry into the target cells (Korber et al., 2020; Walls et al., 2020) . The selective forces directed against the S gene can increase genetic variability in the region, which can be used to infer the evolutionary relationships between viral sequences in a shorter time, compared to use of less variable regions (e.g. RNA-dependent RNA polymerase gene (RdRp)), where mutations occur, but appear to be more costly (Duffy et al., 2008; Moya et al., 2004; Pachetti et al., 2020; Robson, 2020) .
Genetic variability in the S gene can be demonstrated by continuous emergence of mutations, that were reported at a global level (Korber et al., 2020) . Some of these mutations appeared to have a significant epidemiologic value, with the replacement of aspartic acid by glycine at position 614 of the spike glycoprotein (D614G), which is associated with a higher viral shedding and increased infectivity (Korber et al., 2020; Maitra et al., 2020; Zhang et al., 2020) . This mutation currently appears to be dominating the pandemic (Grubaugh et al., 2020) . However, the clinical effect of such mutation is yet to be fully determined (Eaaswarkhanth et al., 2020; Kim et al., 2020b; Korber et al., 2020) . Other mutations in the S gene have also been reported, with the most frequent including: D936Y/H, P1263L, and L5F (Korber et al., 2020; Lokman et al., 2020) .
Similar to other RNA viruses, SARS-CoV-2 can be the subject of phylogenetic analysis due to its high evolutionary rate, and the application of molecular clock analysis might be of value to determine the timing of introductions of large clusters that imply networks of transmission (Duffy et al., 2008; Forster et al., 2020; Pybus and Rambaut, 2009 ). State-of-the-art methods for phylogeny construction include maximum likelihood and Bayesian tools (Anisimova et al., 2013) . is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint
The copyright holder for this this version posted August 25, 2020 . . https://doi.org/10.1101 affected early on during the course of COVID-19 pandemic, with an overwhelming number of cases in some countries (e.g. Iran) (Karamouzian and Madani, 2020; Sawaya et al., 2020) . The first confirmed cases of COVID-19 in the MENA dated back to February 2020 and were reported in UAE, Iran, and Egypt (Daw et al., 2020; Karamouzian and Madani, 2020; Mehtar et al., 2020) . The total number of diagnosed cases of COVID-19 in the MENA exceeded 1,175,000 with more than 32,000 deaths reported as a result of the disease, as of July 25, 2020 (Worldometer, 2020) .
Special attention to COVID-19 infections is needed in the countries of the MENA region, where political and economic factors might lead to devastating effects on the countries affected by the current pandemic (Karamouzian and Madani, 2020; Sawaya et al., 2020) . Particular attention should be paid to countries like Yemen, Syria and Libya, where the ongoing instabilities can result in underreporting of COVID-19 cases and heavy burden on their health-care systems (Da'ar et al., 2020; Daw, 2020; Karamouzian and Madani, 2020; Sawaya et al., 2020) .
The aims of this study included an attempt to phylogenetically analyze S gene sequences and to analyze the spike gene mutation patterns in the MENA region. In addition, we aimed to characterize the temporal changes of D614G mutation spread in the region. is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint
The copyright holder for this this version posted August 25, 2020. . https://doi.org/10.1101/2020.08.24.20176792 doi: medRxiv preprint
All SARS-CoV-2 sequences from the MENA countries, were retrieved from the global science initiative and primary source for genomic data of influenza viruses (GISAID) (Elbe and Buckland-Merrett, 2017). We also downloaded the following sequence metadata if available: date of sequence collection, age, gender, city of collection together with country of sequence collection. The sequences were then aligned to the reference SARS-CoV-2 sequence (accession number: NC_045512) and alignment was conducted using multiple alignment program for amino acid or nucleotide sequences (MAFFT v.7) (Rozewicki et al., 2019) . The MENA sequences that did not contain the complete S region were filtered out. In addition, we removed the sequences that contained indels, the nucleotide ambiguity (N); while other ambiguities were retained. The sequences that contained stop codons were removed as well. Each sequence header was also edited to include data in the following order: country of collection, collection date in days starting from January 5, 2020 (the date of reference sequence collection), city, accession number, gender, and age. The final dataset included 553 MENA S nucleotide sequences that were collected during January 2020 until June 2020.
Analysis of the full MENA SARS-CoV-2 S gene sequences was conducted in Molecular Evolutionary Genetics Analysis software (MEGA6) (Tamura et al., 2013) .
Visual inspection of the aligned MENA amino acid sequences was done, and mutations were identified based on comparison to the reference SARS-CoV-2 sequence (accession number: NC_045512), which was considered as the wild-type.
Amino acids that were translated from codons containing ambiguous bases (e.g. R, Y), were excluded from mutation analysis.
The whole MENA S gene dataset was analyzed phylogenetically using the maximum likelihood (ML) approach in PhyML v3, with selection of the best nucleotide substitution model using Smart Model Selection (SMS), and depending . CC-BY-NC-ND 4.0 International license It is made available under a perpetuity.
is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint
The copyright holder for this this version posted August 25, 2020. . https://doi.org/10.1101/2020.08.24.20176792 doi: medRxiv preprint on Akaike Information Criterion (AIK) (Guindon and Gascuel, 2003; Lefort et al., 2017) . The model which yielded the smallest AIC was the general time-reversible plus invariant sites (GTR + I) nucleotide substitution model with an estimated proportion of invariable sites of 0.625. The estimation of nodal support in the ML tree was based on the approximate Likelihood Ratio Test Shimodaira-Hasegawa like (aLRT-SH) with 0.90 as the statistical significance level (Anisimova et al., 2011) . The ML analysis was repeated ten times and the ML tree with the highest likelihood was retained for final analysis, and determination of the MENA phylogenetic clusters was done by examining the ML tree from root to tips looking for branches with aLRT-SH ≥ 0.90, with large clusters having ≥ 15 sequences.
For the large phylogenetic clusters (containing ≥ 15 sequences and identified using ML analysis), tMRCAs were estimated using the Bayesian Markov chain Monte Carlo (MCMC) method implemented in BEAST v1.8.4 (Drummond et al.,
nucleotide substitution model with discrete gamma-distributed rate heterogeneity, uncorrelated relaxed clock model with a normally-distributed rate prior (initial and mean values of 0.0068, standard deviation=0.0008), and a Bayesian skyline tree density model (Tang et al., 2020) . For each large phylogenetic cluster, one run with 200 million chain length was performed. Samples of trees and parameters were collected every 20,000 steps after discarding a burn-in of 20%, and convergence was analyzed in Tracer v1.6.0 (Rambaut et al., 2015) . The runs were accepted based on effective sample sizes (ESS) of ≥200 and convergence in the trace file. The maximum clade credibility (MCC) trees were assembled using TreeAnnotator in BEAST and were visualized using FigTree (Rambaut, 2012) .
Chi-squared test (χ 2 test) was used to detect differences between the D614 and D614G groups in relation to gender and region (Middle East vs. North Africa).
Mann-Whitney U test (M-W) was used to assess the difference between the D614 and D614G groups in relation to age. Linear-by-linear test for association (LBL) was . CC-BY-NC-ND 4.0 International license It is made available under a perpetuity.
is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint
The copyright holder for this this version posted August 25, 2020. . https://doi.org/10.1101/2020.08.24.20176792 doi: medRxiv preprint used to assess the temporal changes in D614G prevalence. The statistical significance for all aforementioned tests was considered for p<0.050.
A complete list of the MENA SARS-CoV-2 sequence epi accession numbers that were analyzed in this study is provided in (Appendix S1). These sequences are available publicly for registered users of GISAID (Shu and McCauley, 2017) .
. CC-BY-NC-ND 4.0 International license It is made available under a perpetuity.
is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint
The copyright holder for this this version posted August 25, 2020. . https://doi.org/10.1101/2020.08.24.20176792 doi: medRxiv preprint
The total number of MENA SARS-CoV-2 S gene sequences that were included in final analysis was 553, distributed as follows: Oman (n=159), KSA (n=140), Egypt (Table 1) . is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint
The copyright holder for this this version posted August 25, 2020. . https://doi.org/10.1101/2020.08.24.20176792 doi: medRxiv preprint
A total 55 unique non-synonymous mutations in the S gene were detected as compared to the reference SARS-CoV-2 genome. Eight mutations were identified in spike receptor binding domain (SRD), compared to 21 mutations in S2 glycoprotein domain and 26 in other S regions. The most frequent mutation detected in the whole S region was D614G (n=435), followed by Q677H (n=8), and V6F (n=5). The majority of mutations were detected sporadically (n=43, 78.2%, Table 2 ). The highest number of unique S gene mutations (including D614G) was noticed in Oman (n=16), followed by Egypt (n=15), Bahrain (n=9), and KSA (n=6, Table 2 ). 
Oman=1 S13I Oman=1 S12F Egypt=1 1 KSA: Kingdom of Saudi Arabia, 2 UAE: United Arab Emirates, 3 Others: Mutations in S gene regions other than Spike receptor binding domain and S2 glycoprotein, 4 N: The D614G mutation which dominated the sequences were analyzed separately in the main manuscript.
. CC-BY-NC-ND 4.0 International license It is made available under a perpetuity.
is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint
The copyright holder for this this version posted August 25, 2020. Figure   1 ). In addition, no statistical difference was found in analysis per country upon comparing the two variants based on age, age group, or gender. is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint
The copyright holder for this this version posted August 25, 2020. . https://doi.org/10.1101/2020.08.24.20176792 doi: medRxiv preprint
in the whole MENA region as a single unit revealed an increasing prevalence of D614G from 63.0% in January 2020 to reach 98.5% in June 2020 (p<0.001; LBL, Figure 2 ). The same pattern was detected upon comparing the first three months of 2020, compared to April, May and June 2020 (59.5% vs. 90.7%; p<0.001; χ 2 test).
SARS-CoV-2: Severe acute respiratory syndrome coronavirus 2, LBL: Linear-by-linear test for association.
. CC-BY-NC-ND 4.0 International license It is made available under a perpetuity.
is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint
The copyright holder for this this version posted August 25, 2020. . https://doi.org/10.1101/2020.08.24.20176792 doi: medRxiv preprint
To assess the possible presence of phylogenetic clusters in the MENA, ML analysis was conducted. The constructed ML tree showed a star-shaped pattern with short internal branches and long terminal branches (Figure 3) . A total of 13 phylogenetic clusters (aLRT-SH ≥ 0.9) were determined; eight of which included sequences from a single MENA country and five clusters contained sequences collected in more than one MENA country (Appendix S2). Five clusters contained two sequences, and two large clusters were identified, each containing 26 MENA sequences. The highest percentage of clustering sequences was found in Iran (n=3/5, 60.0%), followed by KSA (n=39/149, 27.9%), and Tunisia (n=2/8, 25.0%, Figure 4 ). The overall proportion of phylogenetic clustering was 15.4% (n=85/553). is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint
The copyright holder for this this version posted August 25, 2020. . https://doi.org/10.1101/2020.08.24.20176792 doi: medRxiv preprint 
Bayesian phylogenetic analysis was conducted on the two large clusters identified previously using the ML approach. One Egyptian sequence was removed from each cluster due to the lack of exact collection date. This resulted in analysis of two clusters, each containing 25 sequences. The first cluster contained 14 Saudi sequences, ten Omani sequences and a single Egyptian sequence, with a range of sequence collection between February 13 and May 11.
The median estimate for tMRCA for this cluster having the D614G mutation was is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint
The copyright holder for this this version posted August 25, 2020 . . https://doi.org/10.1101 6.46×10 -3 substitutions/site/year (s/s/y) for the first cluster (95% HPD: 4.87×10 -3 -8.03×10 -3 s/s/y), and 6.50×10 -3 s/s/y for the second cluster (95% HPD: 4.91×10 -3 -8.03×10 -3 s/s/y). The upper MCC tree with sequences having the D614G mutation, while the lower MCC tree represents the D614 cluster. The terminal branches are colored based on country of collection. Internal branches with posterior values ≥ 0.70 are shown in red.
Timeline is represented in day/month/2020. is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint
The copyright holder for this this version posted August 25, 2020 . . https://doi.org/10.1101 
In this study, phylogenetic analysis tools were utilized to assess origins, spread and mutations of SRAS-CoV-2 in the MENA. Phylogeny construction can help to formulate hypotheses regarding the spread of certain taxa having a common origin (Ciccozzi et al., 2019; Pybus and Rambaut, 2009 ). In addition, molecular clock analysis can help to establish a timeline for origins of monophyletic clades (Jenkins et al., 2002; Nasir and Caetano-Anolles, 2015) . Phylogenetic analysis of the MENA S gene SARS-CoV-2 sequences showed a relatively low level of phylogenetic clustering (15%), which hints to a large number of virus introductions into the region. In addition, molecular clock analysis suggests an early introduction of the virus into the MENA which might have been circulating in the region from early February 2020 or even earlier, with subsequent spread into large networks of virus transmission. This estimate of an early virus introduction is supported by the close proximity in time of official reporting of confirmed COVID-19 cases in the region (Karamouzian and Madani, 2020) .
In this study, no evidence of distinct SARS-CoV-2 genetic variants was found.
Plausible explanations might be related to the use of sub-genomic part of the genome (the S gene) rather than utilizing the whole genome. The rationale behind selecting the S region was for two reasons: first, the variability of this region is expected to be higher than other parts of the genome (e.g. RdRp, where mutations are more costly) (Agostini et al., 2018; Shannon et al., 2020) . Second, mutations in the S gene can have significant impact particularly for vaccine development and utility of neutralizing antibodies (Lokman et al., 2020) . The absence of distinct SARS-CoV-2 genetic variants in this study does not provide a conclusive evidence of its genuine absence from the region. These two genetic variants (named L and S lineages) were reported previously, however, a recent report by MacLean et al. carefully discussed the potential pitfalls of such premature conclusions (MacLean et al., 2020; Tang et al., 2020) .
For the estimated evolutionary rate of the two large MENA clusters identified in this study, we based the rate prior selection on the previous finding by (Giovanetti et al., 2020) . This estimate appears higher than other estimates for SARS-CoV-2 and . CC-BY-NC-ND 4.0 International license It is made available under a perpetuity.
is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint
The copyright holder for this this version posted August 25, 2020. . https://doi.org/10.1101/2020.08.24.20176792 doi: medRxiv preprint should be interpreted with caution based on our selection of a strong prior.
However, the rate estimate might appear plausible, since it represents the S gene, rather than the whole genome. For ML analysis, the MENA sequences yielded a starlike phylogeny suggesting a recent growing epidemic (Colijn and Plazzotta, 2018) .
The major result of this study was the demonstration of a temporal shift of SARS-CoV-2 from D614 into D614G variant, which dominated the most recent sequences collected in the region. Such trend was revealed at the global level by Korber et al. , and our results indicated a similar pattern in the MENA (Korber et al., 2020) . In the aforementioned comprehensive study, Korber et al. estimated the global prevalence of D614G at 71.0%, whereas our estimate in the MENA was 78.7%, which appears reasonable, bearing in mind the protracted duration of sequence collection in this study. The explanation for such an observation is most likely related to the association of D614G with a higher viral load and subsequent higher quantities of the virus shed by infected individuals, which increases the likelihood of infection by such a mutant, although an early founder effect of this variant cannot be ruled out (Deng et al., 2020; Farkas et al., 2020; Yurkovetskiy et al., 2020; Zhang et al., 2020) .
Whether this variant can have an effect on severity and outcome of COVID-19 is yet to be fully determined (Becerra-Flores and Cardozo, 2020; Eaaswarkhanth et al., 2020; Korber et al., 2020) . This mutation appeared in all MENA countries, except in Qatar and Iran, which might be related to the low number of sequences from these two countries that were found in GISAID, and the early time of sequence collection (less than 10 sequences from each country were found, dating back to March, 2020).
The emergence of D614G and its increasing prevalence have been reported by several published papers and preprints including a report from North Africa by Laamarti et al., albeit with a fewer number of sequences than the one analyzed in the current study (Gong et al., 2020; Kim et al., 2020b; Laamarti et al., 2020; Maitra et al., 2020) .
Other mutations that were found in the study included Q677H (found only in Egypt), and L5F found in three different countries (Oman, Egypt, and Morocco). The L5F mutation is located in the signal peptide domain of the spike glycoprotein and might be related to recurring sequencing errors (Korber et al., 2020; N. De Maio, . CC-BY-NC-ND 4.0 International license It is made available under a perpetuity.
is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint
The copyright holder for this this version posted August 25, 2020 . . https://doi.org/10.1101 .08.24.20176792 doi: medRxiv preprint 2020 . Nevertheless, its appearance in different studies warrants further investigation to determine its significance (Korber et al., 2020) . The functional importance of Q677H has not been determined yet despite a previous report describing its occurrence (Kim et al., 2020a) .
Limitations of this study should be clearly stated and taken into consideration.
The most obvious caveat in the study was sampling bias. In spite of reporting COVID-19 in all MENA countries, the following countries did not have S sequences submitted to GISAID: Syria, Libya, Yemen, Sudan and Palestine (Iraq had partial sequences that did not include the S gene). In addition, bias was observed for timing of sequence collection. Furthermore, only two countries (Oman and KSA) had more than 100 sequences available for analysis. Another point that should be considered is related to the molecular clock analysis, where we used a strong informative prior which may have affected our tMRCA estimates for dating the origins of the two large phylogenetic clusters. Sequencing errors should also be taken into account, which can partly explain some sporadic mutations that were found in this study.
In the current study, we demonstrated that the D614G variant of SARS-CoV-2 appears to be taking over COVID-19 epidemic in the MENA, similar to what have been reported in other regions around the globe. Local transmission of SARS-CoV-2 might have been established earlier than previously thought, and this illustrates the importance of vigilant surveillance in such conditions of outbreaks by novel viruses. The mutational patterns of SARS-CoV-2 should be closely monitored as the virus seems to be heading into an endemicity in the human population, particularly in relation to mutations' potential impact on passive and active immunization. is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint
The copyright holder for this this version posted August 25, 2020. . https://doi.org/10.1101/2020.08.24.20176792 doi: medRxiv preprint
Appendix S1: A complete list of the MENA SARS-CoV-2 sequence epi accession numbers that were analyzed in this study. is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint
The copyright holder for this this version posted August 25, 2020. is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint
The copyright holder for this this version posted August 25, 2020. . https://doi.org/10.1101/2020.08.24.20176792 doi: medRxiv preprint . CC-BY-NC-ND 4.0 International license It is made available under a perpetuity.
is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint
The copyright holder for this this version posted August 25, 2020. . https://doi.org/10.1101/2020.08.24.20176792 doi: medRxiv preprint . CC-BY-NC-ND 4.0 International license It is made available under a perpetuity.
is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint
The copyright holder for this this version posted August 25, 2020 . . https://doi.org/10.1101 sensitivity to antibodies that target the receptor binding domain. bioRxiv, .2007 .187757. Zhang, L., Jackson, C.B., Mou, H., Ojha, A., Rangarajan, E.S., Izard, T., et al., 2020 The D614G mutation in the SARS-CoV-2 spike protein reduces S1 shedding and increases infectivity. bioRxiv.
. CC-BY-NC-ND 4.0 International license It is made available under a perpetuity.
is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint
The copyright holder for this this version posted August 25, 2020 . . https://doi.org/10.1101 