Molecular Evolutionary Characteristics of the 2019 Novel Coronavirus (SARS-CoV-2) Contracted by Tunisian Citizens : Comparison and Relationship to Other Human and Animal Coronaviruses Based on Spike Glycoprotein-Coding Gene Sequences Analysis

A database of 115 spike glycoprotein-coding gene sequences comprised 15 sequences of each of SARS-CoV-1, SARS-CoV-2, MERS-CoV, HCoV-229E, HCoV-NL63, HCoV-OC43, and HCoV-HKU-1 coronavirus species infecting humans, and 10 accessions of other coronaviruses infecting animals (bats, civet, pangolin, tiger, bovine, dromedary/camel) were used in this study (Table 1). With regard to SARS-CoV-2, the sampling was based on nasopharyngeal swab from Tunisian travellers ¹³ who for many of them were at the origin of the spread of the virus and the outbreak of the epidemic in the country as well as from people permanently resident in Tunisia.

CLUSTAL X 2.1 ¹⁴ software with default settings was used to conduct sequence multiple alignment. Phylogeny inference was carried out by the use of the algorithm Maximum Likelihood (ML) implemented in MEGA-X version 10.0.5 software ¹⁵. Based on the results obtained from best fit model test incorporated and conducted in MEGA-X, inferred phylogenetic relationships were achieved under the assumption of the substitution model GTR (General Time Reversible) coupled to a discrete Gamma distribution (+G) with five rate categories ¹⁶. The substitution model parameters estimated were (i) base frequencies: f(A) = 0.271; f(T) = 0.360; f(C) = 0.181; f(G) = 0.187, (ii) substitution rates: r(AT) = 0.074; r(AC) = 0.045; r(AG) = 0.094; r(TA) = 0.056; r(TC) = 0.131; r(TG) = 0.023; r(CA) = 0.067; r(CT) = 0.261; r(CG) = 0.036; r(GT) = 0.044; r(GA) = 0.135; r(GC) = 0.035, and (iii) transition/transversion ratios: R = 1.55. The Bayesian Information Criterion value (BIC = 102303.728) with GTR+G (+G = 1.25) model was the lowest among the 24 models tested. The branching of the tree was statistically evaluated using bootstrap analysis with 1,000 resamplings of the data. Moreover, in determining phylogenetic network inference, split neighborNet method, based on NJ distance method and incorrected P-distances algorithm contained in SplitsTree4 program ¹⁷, were used. Distance in the evolutionary concept is the smallest amount of evolution that could have occurred. Each changed base pair is considered as one change; whereas, all unchanged base pairs are counted zero.

In order to detect possible recombination signatures in all 115 accessions involved in this study, three different methods were utilized: RDP4.97 package, RECCO, and GARD. RDP4.97 package ¹⁸ put together a number of published recombination detection methods into a single set of tools: RDP ¹⁹, GENECONV ²⁰, BOOTSCAN ²¹, MAXCHI ²², CHIMAERA ²³, SISCAN ²⁴, and 3SEQ ²⁵. Although BOOTSCAN, RDP, and SISCAN are phylogeny-based methods, GENCONV, MAXCHI, CHIMAERA, and 3SEQ are substitution-based methods. In all cases, defaults variables were used. Moreover, these algorithms consider the sequences as linear and set statistical significance at a value of p < 0.05 with Bonferroni correction for multiple comparisons. Only events predicted by more than half of the methods are considered as significant. RECCO, the algorithm performed and developed by Maydt and Lengauer ²⁶ as a fast, simple and sensitive method for detecting recombination in a collection of aligned sequences and locating potential recombination breaking points, is based on cost minimization. This method has only two tunable factors: recombination and mutation cost. In practice the only parameter considered is α, representing the cost of mutation relative to recombination. When α undergoes a change from 0 to 1, the cost of mutation weighted by α increases, and the cost for recombination weighted by 1- α decreases. More specifically, the element α controls the ambiguity between mutation and recombination. GARD (Genetic Algorithm for Recombination Detection) has the purpose to screen a multiple sequence analysis aiming at detecting recombination signals and inferring a unique evolutionary history for each determined block situated on the left and on the right sides of the breaking point ^{27, 28}. Used as a basic component of GARD, the small-sample corrected Akaike Information Criterion (AICc) ²⁹ is essential for breakpoint assignments score. Detected breakpoints are then evaluated for significance through KH test ³⁰ of the HyPhy package ³¹.

Seeking for selection signatures occurring in 115 sequences of spike glycoprotein-coding gene partitioned into clusters and subclusters, three approaches were explored : site, branch, and gene. Precisely, the first approach consisted of finding individual sites under positive/negative selection. Hence, four different methods i.e., the Single-Likelihood Ancestor Counting (SLAC), Fixed Effects Likelihood (FEL), Random Effects Likelihood (REL) ³², and Fast Unbiased Bayesian Approximation (FUBAR) ³³ were utilized. They are codon-based maximum likelihood methods used to estimate dN/dS ratio at every codon in the alignment. In addition, Internal Fixed Effects Likelihood (IFEL), an appropriate model for investigating whether sequences sampled from a population were subject to selective constraints at the population level (i.e., along internal branches), as well as Mixed Effects Model of Episodic Selection (MEME) ³⁴ aiming at looking for evidence of both diversifying, and importantly, episodic selection at individual sites, were also used. The second procedure comprised two methods: aBSREL, and GA-branch. The former model i.e., adaptive Branch-site Random Effects Likelihood (aBSREL) ³⁵ was suited to test if positive selection occurred on a proportion of branches. In contrast, the latter i.e., GA-branch (Genetic Algorithm-branch) ³⁶ is a codon-based genetic algorithm which can divide all branches of the phylogenetic tree setting out non-recombinant input into classes as stated by dN/dS ratio values. Various inferred models provided confidence intervals on dN/dS for each branch. The third process dealt with Branch-site Unrestricted Statistical Test for Episodic Diversification (BUSTED) whose principle was to provide a gene-wide test for positive selection and thus determine whether a gene underwent positive selection at at least one site on at least one branch. BUSTED utilized a codon model with three rate classes, constrained as ω₁ ≤ ω₂ ≤ 1 ≤ ω₃. It involved two models namely unconstrained (alternative) and constrained (null hypothesis). If the null hypothesis is rejected, there is, subsequently, evidence that at least one site has experienced, at least some of the time, positive selection on the foreground branches ³⁷. Last, the Partitioning Approach for Robust Inference of Selection (PARRIS) ³⁸ enlarges classic codon-based likelihood ratio tests to reveal if a part of sites in the alignment underwent positive selection (ω ˃ 1). The model accomplished tests in three phases: (i) fitting an optimal model to estimate relative branch lengths, (ii) fitting the null hypothesis (M1) with three synonymous rate classes, and (iii) fitting the alternative model (M2) with three synonymous rate classes. It is worth mentioning that all selection algorithms used and previously indicated, were made available in the web server of HyPhy package implemented at http://www.datamonkey.org ^{39, 40, 41}

The Maximum Likelihood (ML) algorithm incorporated in MEGA-X software allowed the reconstruction of a phylogenetic tree where 115 spike glycoprotein-coding gene sequences were split into four distinct groups (Figure 1). Although, group I comprised three components (SARS-CoV-1, SARS-CoV-2, and SARS-CoVs infecting exclusively animals), group II encompassed all MERS-CoV lineages, group III gathered all human coronaviruses lineages of the species HCoV/OC43, and HCoV/HKU1, and finally group IV which was formed by all human coronaviruses lineages of the species HCoV/229E, and HCoV/NL63. It is worthwhile to note that three members of the animals SARS-CoVs have integrated the human coronaviruses clusters. It is indeed about SARS-CoV-2/tiger lineage which became member of human SARS-CoV-2 species cluster. Similarly, MERS-CoV/dromedary-camel became a component of MERS-CoV species clade, and finally SARS-CoV/bovine which was part of HCoV/OC43 species group. Howbeit, phylogeny inference and its significance, using ML algorithm implemented in MEGA-X software, may become poorly relevant if reticulate events (recombination, lateral gene transfer, gene duplication, reassortment) take place. For this reason, network phylogeny inference, highlighting the presence of conflicting signals and expressing the real phylogenetic relationships among all coronaviruses lineages used in this study, was adopted as shown in Figure 2. Herein, the produced network, using SplitTree4 software, showed almost the same topology with, however, a few differences which were twofold: (i) a clear-cut delineation of subgroups within the respective clusters. Each subgroup contained specifically all lineages of a distinct coronavirus species (including the new integrated animal lineage) as mentioned above and shown in Figure 1, (ii) a reshuffling of group I. Indeed, in subgroup I, the accession AY508724 was no longer member of SARS-CoV-1 species and moved to subgroup III becoming part of it. In similar fashion, the Tunisian accessions MT499217, MT955170, and MT559038 of SARS-CoV-2 species (subgroup II) shifted to subgroup III (Figure 2, Table 1). More explicitly and in furthering networked phylogenetic relationships among members of each identified subcluster, various analyses using SplitTree4 software, were carried out. Five out of eight coronavirus species groupings (i.e., reshuffled animal SARS-CoVs, MERS-CoV, HCoV-229E, HCoV-NL63, and HCoV-OC43 subgroups), showed network topologies containing conflicting signals represented by boxes and edges. Regarding animals SARS-CoVs subgroup, networked phylogenetic relationships showed that 11 accessions fell into seven different clades. Interestingly, three Tunisian SARS-CoV-2 lineages (MT559038, MT499217, and MT955170) were much closer to SARS-CoV/bat RatG13 (MN996532), and to SARS-CoV/pangolin (MT040333) than to any other SARS-CoV/bat lineages (MG772933, MG772934, DQ022305, and KF367457) nor civet lineage (AY686863) either. Likewise, the shifted lineage of SARS-CoV-1 (AY508724), although it was genetically close to civet lineage, it remained distant from all the other aforementioned (Figure 3). Concerning MERS-CoV lineages, they, phylogenetically, split into eight distinct clades. In addition, the Qatari accession (MN507638) wasn’t so different since it joined other MERS-CoV lineages from Saudi Arabia and formed a homogenous phylogroup (Group II). Seemingly, this group was less impacted by the reticulate factors effects compared to all other groups as shown by the absence of conflicting signals (boxes). Furthermore, the MERS-CoV lineage (KF917527) isolated from camel formed alone a distinct phylogroup (group V) (Figure 4). Referring to HCoV-229E species, two major clades differing geographically, were set out. There was indeed on the one hand, the German accessions (NC_002645, and AF304460) and on the other hand, the American ones distributed in six distinct groups. Besides, the depicted network showed heavy incompatibilites with a noticeable acuity in the USA cluster (Figure 5). With respect to HCoV/NL63 species, inferred networked phylogeny described three major groups and their split-up was overly controlled by the geographic origin of the sampled lineages particularly those belonging to the groups I (from USA) and II (from China) albeit the third group contained lineages from both regions as well as from South Korea. In addition and notoriously, the conflicting signals have spread along the whole network (Figure 6). As regards and analogously, 16 accessions of HCoV-OC43 species segregated into four distinct clades and their split-up was somewhat under the influence of the geographic origin of the lineages. In fact, the Mexican lineage was detached from the American ones as it was even more the lineage SARS-CoV/bovine (Group II) which was completely apart. It is noteworthy that striking reticulations congregated in the network extending from group I to group IV (Figure 7).

Figure 1.The circle form of a phylogenetic tree produced by the maximum likelihood algorithm under assumption of the general-time reversible (GTR) substitution model coupled to a discrete gamma distribution (+G) option of MEGAX software 15. Four major groups were delineated. The numbers above the branches indicate the bootstrap confidence value. The scale bar shows the number of substitution per nucleotide. Thick black lines indicate the demarcation line between phylogenetic groups. Double headed-arrows show the demarcation line between coronavirus species groups. Single headed-arrows represent the animal coronavirus species having shifted and integrated the human coronavirus phylogroups.

Figure 2.Networked phylogenetic relationships among 115 coronavirus species accessions which split into four clear-cut distinct clusters subdivided into subgroups (except cluster II).

Figure 3.NeighborNet for 11 sequences of spike glycoprotein-coding gene of assembled coronaviruses lineages of SARS-CoV-1, SARS-CoV-2, and animal SARS-CoVs belonging to cluster I (subgroup III). The networked relationships indicate the presence of reticulate events. Boxes imply likelihood of recombination. The phylogenetic network constructed, using SplitsTree4 software, delineated seven distinct clusters. The scale bar shows the number of substitution per nucleotide.

Figure 4.NeighborNet for 16 sequences of spike glycoprotein-coding gene of coronaviruses lineages of MERS/CoV belonging to cluster II. The networked relationships indicate the presence of reticulate events. Boxes imply likelihood of recombination. The phylogenetic network, constructed, using SplitsTree4 software, delineated eight distinct clusters. The scale bar shows the number of substitution per nucleotide.

Figure 5.NeighborNet for 15 sequences of spike glycoprotein-coding gene of coronaviruses lineages of HcoV/229E belonging to cluster III (subgroup I). The networked relationships indicate the presence of reticulate events. Boxes imply likelihood of recombination. The phylogenetic network constructed, using SplitsTree4 software, delineated seven distinct clusters. The scale bar shows the number of substitution per nucleotide.

Figure 6.NeighborNet for 15 sequences of spike glycoprotein-coding gene of coronaviruses lineages of HCoV/NL63 belonging to cluster III (subgroup II). The networked relationships indicate the presence of reticulate events. Boxes imply likelihood of recombination. The phylogenetic network constructed, using SplitsTree4 software, delineated seven distinct clusters. The scale bar shows the number of substitution per nucleotide.

Figure 7.NeighborNet for 16 sequences of spike glycoprotein-coding gene of coronaviruses lineages of HCoV/OC43 belonging to cluster IV (subgroup I). The networked relationships indicate the presence of reticulate events. Boxes imply likelihood of recombination. The phylogenetic network constructed, using SplitsTree4 software, delineated seven distinct clusters. The scale bar shows the number of substitution per nucleotide.

Pairwise nucleotide identity comparison among 115 sequences subdivided in clusters and subclusters (eight subdivisions) (according to NeighborNet phylogeny inference performed in SplitTree4 software) of spike glycoprotein-coding gene of coronavirusspecies infecting humans (SARS-CoV-1, SARS-CoV-2, MERS-CoV, HCoV/229E, HCoV/NL63, HCoV/OC43, and HCoV/HKU1) and a mixed group containing various species of animals particularly bat, civet, camel, pangolin, tiger, and bovine, was carried out and the analyses resulted in the following percentage intervals : 89.99-100.00, 99.86-100.00, 26.76-99.92, 99.40-99.97, 40.35-100.00, 35.22-99.97, 31.12-100.00, and 54.47-100.00 ; and the number of most distant and closest couple of accessions were as follows : 9-16, 1-15, 1-1, 2-3, 4-2, 1-3, 4-2, and 1-4, for members of the subgroups I, II, and III (cluster I), cluster II, subgroups I, and II (cluster III), and subgroups I, and II (cluster IV), respectively (Table 1S). From this analysis, it emerged that the sequences of the lineages belonging to Subgroups I, and II (Cluster I), and Cluster II, were highly similar. In contrast, a lower similarity was observed in the sequences of the lineages affiliated to the subgroups I, and II of the clusters III, and IV. Undoubtedly, the mixed group (subgroup III, cluster I) comprising heterogenous members, showed the lowest similarity in the sequences of spike glycoprotein-coding gene.

Although mutations are changes in nucleotide sequences due to errors in replication or repair, substitutions are mutations that have passed through the filter of selection. Substitution often involve amino acids with similar chemical characteristics supporting two evolutionary principles : (i) mutations are rare events, (ii) most dramatic changes are removed by natural selection. The analysis of both number and the type of substitutions, that have occurred during the evolution, are of central importance for the study of molecular evolution. In determining the substitution pattern and rates in the sequence of spke glycoprotein-coding gene of 115 lineages studied here, a substitution matrix was estimated for each type of subdivision as described earlier. Accordingly, it was revealed for subgroup I (cluster I) that rates of different transitional substitutions oscillated from 11.87 to 18.75; while, rates of transversional substitutions ranged from 3.76 to 5.93. Regarding subgroup II (cluster I), rates of transitional substitutions varied from 10.96 to 18.42 ; whereas, rates of transversional substitutions spanned from 3.85 to 6.46. Concerning subgroup III (cluster I), although rates of transitional substitutions extended from 9.51 to 26.03, rates of transversional substitutions fluctuated from 2.77 to 8.10. With regard to cluster II, whilst the span of rates of transitional substitutions ranged from 4.92 to 43.28, that of rates of transversional substitutions varied from 1.78 to 3.14. With respect to subgroup I (cluster III), as the extent of rates of transitional substitutions was from 7.12 to 36.08, that of rates of transversional substitutions was from 2.65 to 5.13. As regards, in subgroup II (cluster III), the expanse of rates of transitional substitutions was from 7.51 to 45.11, that of rates of transversional substitutions was from 1.46 to 3.31. Relating to subgroup I (cluster IV), rates of transitional substitutions oscillated from 6.49 to 39.25 ; while, rates of transversional substitutions fluctuated from 2.19 to 4.86. Referring to subgroup II (cluster IV), rates of transitional substitutions expanded from 9.77 to 28.64; whereas, rates of transversional substitutions ranged from 2.37 to 6.94 (Table 2). Therefore, it came out from this analysis that the highest rate of transitional substitutions regarded the residues T and C had a peack of 45.11 (subgroup II, cluster III). In contrast, the lowest rate of transitional substitutions concerned the residues G and A, had a downward peak of 4.92 (cluster II). On the other hand, the highest rate of transversional substitutions related to the residues T and A, culminated at 8.10 (subgroup III, cluster I). In opposition, the lowest rate of transversional substitutions dealt with the residues C and A, and C and G having the value of 1.46 (subgroup II, cluster III) (Table 2).

Generally, recombination is consisting of an exchange of genetic material in segmented as well as in unsegmented virus genome. It allows the virus to acquire new traits. In order to obtain reliable and relevant results, three bioinformatics methods were used : RDP package, RECCO, and GARD. The objective of the use of RDP package version 4.97 was threefold : (i) to identify the recombinant coronaviruses, (ii) to determine the different site recombination locations, (iii) to recognize the potential parentals. Accordingly, seven lineages having the accession numbers: KF530112, from USA, MK334043, MK334044, MK334046, MK334047 and JX104161, from China, and MG772808, from South Korea belonging to the coronavirus species HCoV/NL63, were recombinants ; the number of recombination sites were as follows : 5, 5, 3, 2, 3, 4, and 4, respectively. Furthermore, there was no geographical impediment with regard to the donor of genetic material. In other words, both major and minor parents could have any geographical origin (Table 3). Besides, only one lineage of HCoV/HKU1 species from USA was recombinant: MK167038 and only two recombination sites were detected. Here too, the geographical origin (USA or China) of the donors did not also constitute a restraint whatsoever (Table 3). Last, RDP package permitted the detection of breakpoints in the sequences of three coronaviruses infecting bat affiliated to cluster I (subgroup III) : SARS-CoV/SL-CoV-ZC45 (MG772933), SARS-CoV/SL-CoV-ZXC21 (MG772934), and SARS-CoV/HKU3-1 (DQ022305). Although the former two underwent recombination in four distinct sites, the latter recombined in only two sites. It was pointed out also that, in the vast majority, both major and minor parents had intraspecies origins. However and interestingly, the Tunisian SARS-CoV-2 lineage MT499217 was a possible donor (major parent) of genetic material for the bat lineages : SARS-CoV/SL-CoV-ZC45 (MG772933), and SARS-CoV/SL-CoV-ZXC21 (MG772934) (Table 3). This possibility was not so surprising since the genetic connection between SARS-CoV-2 and coronaviruses infecting bats was established by the phylogenetic relationships among lineages of cluster I (subgroup III) as shown in Figure 3 and reported elsewhere ⁴². Still, this observation needs to be ascertained by undertaking further analyses involving other genes constituting the SARS-CoV-2 genome in relation to SARS-CoV/bat genome. In contrast, RECCO algorithm demonstrated that only two and different coronaviruses species infecting bats i.e., SARS-CoV/HKU3-1 (DQ022305), and SARS-CoV/WIV1 (KF367457), were recombinants. Likewise, there was an additional recombinant: SARS-CoV/PCOV_GX-P4L isolated from pangolin (Figure 8a). Whereas the former lineage recombined in 23 sites, most of them (9) had a size as long as three residues, the latter two possessed nine and 14 breaking points, most of them (5 and 8), were 17 and 57 residues long, respectively (Table 2S). On the other hand, RECCO was able to detect recombination signals in the sequences of representatives of human coronaviruses : HCoV/229E (2 isolates) (Figure 8b), HCoV/NL63 (5 isolates) (Figure 8c), HCoV/OC43 (1 isolate) (Figure 8d), and HCoV/HKU1 (1 isolate) (Figure 8e). Furthermore, RECCO pointed out that in most cases, the size of the swapped segment exceeded three residues reaching even a size of 1,001 residues for the lineage MK334044 of the species HCoV/NL63. As regards, GARD and despite the fact that it wasn’t able to identify precisely which accession was recombinant by contrast to RDP4.97 and RECCO, it dealt with codon-based alignment and succeeded in placing breakpoints in various positions in the sequences. As a matter of fact and based on statistics, GARD determined 3, 1, 4, 5, 5, and 6 breaking points in aligned sequences of cluster I (subgroup III), cluster II, cluster III (subgroups I, and II), and cluster IV (subgroups I, and II), respectively (Table 4). Moreover and based on the use of Kishino-Hasegawa test (KH test) for significance estimation, it compared phylogenies constructed on both sides of the breakpoints and demonstrated that the topology was incongruent (significance at p = 0.1, p = 0.05, and p = 0.01) at the positions : 206 ^{cluster I (subgroup III)}, 1049 and 2260 ^{cluster III (subgroups I and II respectively)} ,448 ^{cluster IV (subgroup I)}, and 1677, and 2136 ^{cluster IV (subgroup II)} in aligned sequences of different lineages (Table 4). It is worth mentioning that while significant topological incongruence denotes that different tree topologies took place between portions delimited by recombination breakpoints, KH-insignificant breakpoints arise most frequently due variation in branch lengths among fragments. This could be due to some processes of recombination or to other activities such as spatial rate variation, heterotachy, etc. Finally and by way of conclusion, all three statistics-based methods used in this study (RDP package, RECCO, and GARD) failed to detect any recombination event in the sequences of Severe Acute Respiratory Syndrome coronaviruses 1 and 2 (SARS-CoV-1 and SARS-CoV-2).

Figure 8.Graphs displaying potential recombination breakpoints illustrated by downward peaks in aligned sequences of spike glycoprotein-coding gene of coronaviruses lineages belonging to a) cluster I (subgroup III) (SARS-CoV/bats and pangolin), b) cluster III (subgroup I) (HCoV/229E), c) cluster III (subgroup II) (HCoV/NL63), d) cluster IV (subgroup I) (HCoV/OC43) and e) cluster IV (subgroup II) (HCoV/HKU1) determined by RECCO algorithm.

This evolutionary mechanism, defined as the existence of genetic elements that encode for the same function resulting in sequence redundancy, was searched in all sequences of spike glycoprotein-coding gene considering each subdivison separately as was the case for previous analyses. Based on the reconstruction of a newick tree performed by MEGA-X algorithm, search was done in all branching points of the topology. The search for duplication events was accomplished by determining the position of the root on a branch(es) that generated the minimum number of duplication events using an unrooted gene tree (newick) analysis. As a result, although gene duplication events were revealed in all nodes of each tree of seven subdivisions ^{cluster I (subgroups I and II) cluster II cluster III (subgroups I and II) and cluster IV (subgroups I and II)}, they were located in only two branching points (the rooted node and the branching point of two Tunisian SARS-CoV-2 lineages (MT955170, and MT499217) belonging to subgroup III of cluster I (Figure 9). Thus, from this analysis, it was possible to conclude that the more highly similar the sequences, the greater the possibilities of detecting gene duplication events. In point of fact, the pairwise nucleotide identity comparison produced for subgroup III (cluster I) (a highly heterogenous subgroup), an interval with the greatest difference in percentages i.e., 26.76%-99.92% (Table 1S).

Figure 9.Cladogram depicting gene duplication events in the rooted ancestor as well as in the node of two Tunisian SARS-CoV-2 lineages (shown by two arrows). Identification was done by searching for all branching points in the topology with at least one species that is present in both subtrees of the branching point. An unrooted gene tree was used for the analysis such that the search for duplication events was performed by finding the placement of the root on a branch or branches that produced the minimum number of duplication events. Evolutionary analyses were conducted in MEGA-X.

Natural selection, briefly defined as the unequal survival and reproduction of hereditary material due to environmental forces resulting in the preservation of favorable adaptations, was screened in aligned sequences of all eight described subdivisions. To find selection signature over sequences, site-specific models (SLAC, REL, FEL, IFEL, FUBAR, and MEME) and branch-specific models (aBSREL and GA-branch), were used. As indicated in Table 5 for the first category of models (site-specific), whereas purifying selection was constantly detected in the subdivisions : cluster I (subgroup III), cluster II, cluster III (subgroups I and II), and cluster IV (subgroups I and II), adaptive selection was irregularly detected. More specifically, regardless to the subdivision tested, the models IFEL (6/8), FEL (5/8), and MEME (5/8) were more suited to detecting positive selection signatures than SLAC (1/8), FUBAR (1/8), and REL (3/8) models. The highest number of positively and negatively selected sites were revealed in the sequences of HCoV/HKU1 and HCoV/NL63 isolates with 206 and 195 sites, respectively and imparted by REL model. Contrariwise, only a single positively selected site was detected through the sequences of MERS/CoV (IFEL, and FUBAR), HCoV/229E (REL) and the mixed group (cluster I, subgroup III) (REL) isolates. In like manner, only a single negatively selected site was found in the sequences of MERS/CoV isolates (IFEL). Overall, a significant number of sites underwent positive selection as reaveled by MEME model particularly in the sequences of the mixed group (cluster I, subgroup III) (46 sites) and HCoV/229 (43 sites) isolates to a lesser extent in the sequences of HCoV/NL63 (17 sites), HCoV/OC43 (34 sites), and HCoV/HKU1 (34 sites) isolates. At all events, most sites were under purifying selection. Outstandingly, compared to all the other models, the counting method SLAC was a lesser performer in detecting positive as well as negative selection signatures. More importantly, all previously described models failed to detect neither positive nor negative selection in the sequences of SARS-CoV-1 (cluster I, subgroup I) and SARS-CoV-2 (cluster I, subgroup II) isolates (Table 5). Relating to the second category of models (branch-specific), with aBSREL algorithm and after correcting for multiple testing, significance was assessed using the Likelihood Ratio Test (LRT) at a threshold of p ≤ 0.05 resulting in the evidence of episodic diversifying selection on one branch (out of 18), one branch (out of 21), five branches (out of 21), one branch (out of 23), seven branches (out of 23), and five branches (out of 17) of the subdivisions representing the mixed group, MERS-CoV, HCoV/229E, HCoV/NL63, HCoV/OC43, and HCoV/HKU1 species, respectively. On the contrary, no evidence of episodic diversifying selection was found in the sequences of the remaining subdivisions representing SARS-CoV-1, and SARS-CoV-2 species (Table 6). In parallel and in order to shed light onto the lineage-specific nature of the selective constraints exerted on each branch of the reconstructed trees, GA-branch algorithm was employed. On the basis of the use of the statistical approach employing small-sample AIC (c-AIC) best scores, the branches were grouped into two classes for the subdivisions of SARS-CoV-1 and SARS-CoV-2 species and the mixed group, four classes for the subdivisions of MERS-CoV, HCoV/229E, and HCoV/OC43 species, and last, five classes for the subdivisions of HCoV/NL63, and HCoV/HKU1 species with the support of 275, 1876, 7288, 2831, 1526, 1395, 1871, and 2548 tested models at 95% confidence set, respectively (Table 7). Withal, although the lineages of the subdivisions of SARS-CoV-1, and SARS-CoV-2 were not under selection, those of the mixed group, MERS-CoV, and HCoV/229E species experienced predominately purifying selection ; indeed, the total percentages of the branches of the tree where dN/dS <1 were 82%, 53%, and 75%, respectively. In contrast, members of the remaining subdivisions i.e., HCoV/NL63, HCoV/OC43, and HCoV/HKU1 species were predominately under adaptive selection with the following rates : 99%, 58%, and 65%, respectively (Table 7). In other respects, BUSTED with synonymous rate variation found evidence (LRT with p = 0.00 ≤ 0.05) of gene-wide episodic diversifying selection only in the selected branches of the phylogeny of the clade representative of the mixed group (cluster I, subgroup III). As a consequence, there was evidence that at least one site on at least one test branch had experienced diversifying selection. At the end, PARRIS model was able to detect positive selection in only the clades HCoV/NL63 and HCoV/HKU1 species (Table 3S).