Introduction

Peste des petits ruminants (PPR), caused by Small ruminant morbillivrus (SRMV), is a contagious transboundary disease of domestic and wild ruminants1,2. Despite exhaustive vaccination, the disease is endemic across many regions/countries in Africa, Middle East and Asia, where occurrence of frequent disease outbreaks is not uncommon3,4,5,6,7. Currently, the PPR is threatening approximately 80% of the global population of sheep and goats with an estimated loss of USD 2.1 billion per year8.

The SRMV belongs to the genus Morbillivirus within the family Paramyxoviridae. It is a pleomorphic and enveloped virus that carries a negative sense RNA genome9 of variable length, from 15,927 to 16,058 nucleotides (NCBI database). The genome encodes six structural and two non-structural proteins in an order of 3′-N-P/C/V-M-F-HN-L-5′. Non-structural proteins (V and C) are encoded either by alternate open reading frames or mRNA editing in the phosphoprotein (P) gene. Based upon either N gene (255 bp) or F gene (322 bp), four distinct lineages of SRMVs (I-IV) are reported so far. Lineage I-II viruses are mostly reported from West African countries. Lineage III viruses seem restricted to the Middle East and East African countries. Lineage IV viruses have been reported from Asian and African countries1,10. The lineage IV is replacing prevalence of other lineages (i.e. I-III) territories and the occurrence of lineage IV is overwhelming even in Africa. These features demonstrate that lineage IV possess stronger positive selection and host-adaptation potential in a wide spectrum of hosts and geographical areas11,12.

Given the fact that genetic variations within a population of viruses could alter their pathogenicity and host spectrum, viral genetic diversity is considered a key to unleash viral evolution13. Using complete or partial sequencing of single genes (H, N or F), the clustering pattern, genomic and residue characteristics of SRMVs have widely been studied and discussed across the globe5,10,12,14. However, based upon each of these particular genes, the deduced genomic and residue characteristics may not be considered enough to predict ongoing evolutionary patterns across the whole length of the genome. In addition, many aspects of SRMVs evolution, including ancestral strain links, historical and geographic patterns of strain dispersal, divergence and time of origin remain poorly understood. These aspects are important because evolution within a single gene may not necessarily be occurring at the same rate as that of the whole genome15. Also, being RNA viruses, SRMVs are more prone to mutations during acute infection and therefore could present a polymorphic population11. Therefore, genetic diversity driven from consensus sequences of partial genomes could be far from representing the actual polymorphism across the whole length of the genome. Taken together, understanding comparative phylogenomics and evolutionary dynamics by exploiting complete genome sequence data of SRMVs facilitate better elucidate the genetic diversity, trends in its evolution and disease distribution pattern across diverse geographical regions. With this background, complimented by two complete genome sequences from Pakistan we used complete genome sequence data of SRMVs accessible in public database (until October 01, 2019) and analyzed for genetic diversity, phylogenomics and residue characteristics through different bioinformatics tools. We extend the analysis to each of the coding genes and identified potential ancestral relationship among SRMV-lineages reported from different countries during different time-period. In addition, we analyzed coding genes of all reported complete genomes to determine SRMV’s divergence time, and identified another candidate gene to be used as a phylogenetic marker. Together, the outcome will be expected to enhance our understating of phylogenetic and evolutionary dynamics of SRMVs across the globe.

Results

Comparative genome features

The comparative genomic analysis revealed a variable length of genomes as 15927, 15942, 15948, 15954, 15957 and 16058 nucleotides. Most of sequences across the globe had 15948 nucleotides (n = 39) whereas a number of Chinese strains (n = 31) and a Mongolian strain (KY888168) possessed a genome length of 15954 nucleotides. Only a single SRMV strain reported from India (KT270355) carried 15942 nucleotides. One Israeli strain (MF678816) had 15927 nucleotides. Two unusual genome lengths of 15957 (KM089831) and 16058 nucleotides (KM816619) were exclusively reported from China (Table 1). Excluding complete genomes of unusual lengths (MF678816, KM089831 and KM816619) while performing complete genome-specific analysis, the percentage for GC and AT contents was 47% and 53%, respectively. The proportion of GC content was found highest in N gene (50%) followed by P (48%), each of M, F, H (46%), L genes (43%), trailer (41%) and leader region (38%) (Table 2). The study genomes had a 52 nucleotide (nt) long leader in 3ʹ UTR at 107 nt long genome promoter region and a 73 nt long trailer at 5ʹ UTR in 109 nt long anti-genome promoter region. The total length of each of the genes varied across the whole genome: N gene (1578 nt) encoded 526 aa of 58 KDa, P gene (1530 nt) encoded 510 aa of 55 KDa, M gene (1008 nt) encoded 336 aa of 38 KDa, F gene (1641 nt) encoded 546 aa of 59 KDa, H gene (1830 nt) encoded 610 aa of 69 KDa and L gene (6552 nt) encoded 2184 aa of 247 KDa. Although all complete SRMV sequences showed variable genome length, the coding region for each of the genes was the same. The varying genome length was due to insertion of nucleotides in a non-coding region between P and M genes, and between M and F genes (Table 3). However, all genes were separated by similar conserved non-coding intergenic trinucleotide (CTT) except for one intergenic region between L gene and the trailer sequence (CTA).

Table 1 A brief summary of dataset on SRMVs available at public database incluidng under-study Pakistan-originated strains
Full size table
Table 2 A brief descriptions on genome atlas including coding and non-coding regions of so far reported SRMVs worldwide
Full size table
Table 3 A comparative analysis for the coding genes and intergenic regions present in the whole genome of SRMVs reported from different regions of the globe.
Full size table

Percentage identity of nucleotide and comparative residue analysis

We found a varying nucleotide divergence among strains representing different lineages and geographical settings. For instance, a maximum nucleotide divergence (12.7%) was observed among Mongolian, Georgian (lineage IV) and Asian strains (lineage III). This was followed by 11.9% divergence between Pakistani (lineage IV) and other Asian strains (lineage III), and 11.8% divergence between Chinese (lineage IV) and rest of Asian strains (lineage III). As high as 11.5% nucleotide divergence was observed between Asian (lineage II) and African strains (lineage III) of SRMV. Similarly, a total of 11% nucleotide divergence was observed between African strains of lineages II and III. However, a variable divergence (8.5–10.3%) was noticed between SRMVs of lineage I and IV whereas, a divergence of 1.0–4.9% was revealed among strains within lineage IV (Table 4).

Table 4 Percntage nucleotide identities and divergence derived from complete genome consensus sequences of SRMVs strains (lineage I–IV) reported so far in the public database.
Full size table

Comparative residue analysis of different proteins across the entire genome length revealed conserved functional and/or structural motifs; however, few substitutions were noticed in some of the studied strains. A hypervariable region of varying length was observed in each of the SRMV proteins i.e., 423–456 aa in N, 74–111 aa in P, 73–197 aa in M, 6–16 aa in F, 174–179 aa in H and 617–627 aa in L protein. The nuclear export and nuclear localization signal, and RNA binding motifs appeared conserve in N protein of all strains. In P protein, a Soyuz 1 motif was also conserved in all strains except for the consensus sequence of Africa/1994–2012 (lineage III) where a total of six substitutions (L5Q, V10N, E11K, A14E, L16I and F20K) were observed. A serine residue (151S) in the P protein and a cell membrane anchor in the M protein were conserved in all of the SRMV sequences (Table 5). The signal peptide in F protein has previously been reported to be hypervariable (Table 6); however, while comparing SRMVs of different lineages, we proposed a relatively conserved long stretch of residue (1MTRVAILTFLFLFPNVVAC19) (Fig. 1). The cleavage motif (103GRRTRR108) was conserved in the F protein of all sequences. The fusion peptide motif was conserved for 109–133 aa in all SRMV strains except for consensus sequence of China/2013–15 strain where, phenylalanine (F) was replaced by leucine (L) at 1st position of the motif. Substitutions were observed in leucin zipper domain of consensus sequence in lineage II (African/2009–15, V479I), lineage III (African/1994–2012, I463V) and lineage IV (Bangladesh/2008, A464T). All consensus strains from lineage II including China/2013–15, Mongolia/2016, Georgia/2016 and Ethiopia/2010 carried a conserved residue pattern for hydrophobic anchor membrane of F protein; however, two substitutions (A486V and G489S) were observed predominantly in sequences from lineage IV. While comparing residue type and position in the H protein, several substitutions were revealed. For strains within lineage IV, these included a substitution in the N-terminal anchor of an Indian strain (India/2014–16, A41V) and in Georgian strain (Georgia/2016, Y481H). A substitution common to all SRMV strains within lineage III was observed in SLAM binding site where tyrosine (Y) was replaced by phenylalanine (F) at position 553, whereas a substitution in asparagine N-linked glycosylation site (215NVT217) was exclusive to strains reported from Africa during 1994–2012 (Table 7). For the N protein, all functionally and structurally important motifs were conserved in strains representing lineage I-IV.

Table 5 A summarized comparative residue analysis of important domain and motif at NP, P and M proteins of SRMVs for their structural, functional and biologic activities.
Full size table
Table 6 A summarized comparative residue analysis of important domain and motif in the F and H proteins of SRMVs for their structural, functional and biologic activities.
Full size table
Figure 1
figure 1

WebLogo-based diversity and/or conserveness of residues in proposed stretch at fusion protein of so-far reported SRMVs

Full size image
Table 7 A summarized comparative residue analysis of important domain and motif in the L protein of SRMVs for their structural, functional and biologic activities.
Full size table

Estimation of evolutionary and divergence rates

Using a Bayesian coalescent approach, a molecular clock analysis of the whole genome and all coding gene sequences was performed to estimate the mean rate of evolution. Based on this analysis, the mean evolution rate for the complete genome of SRMV was estimated to be 9.953 × 10–4 substitutions per site per year. Best growth model was used for individual SRMV gene dataset to estimate the TMRCA and substitution rate per site per year. A cumulative interpretation of individual gene-based analysis showed rates of evolution in N, P, M, F, H and L genes as 1.1 × 10–3, 1.23 × 10–3, 2.56 × 10–3, 2.01 × 10–3, 1.47 × 10–3 and 9.75 × 10–4 site per year, respectively. The N gene (1.1 × 10–3) showed a lesser evolution rate as compared to other genes whereas it was highest for the L gene (9.75 × 10–4).

Phylogenetic topology based on geographical pattern

Utilizing each of the coding genes, the phylogenetic analysis of SRMV sequences revealed a distinct pattern of clustering according to the geographical locations. However, the complete N gene-based clustering pattern was more authoritative and conclusive followed by L, H, F, M and P genes (Fig. 2a,b). Within lineage II viruses, variations in clustering pattern were related to the reporting period from different regions in the African continent. In contrast, lineage III viruses from Africa showed variations in their clustering pattern on the basis of each of five gene used for analysis. For lineages IV viruses, there were significant variations in clustering pattern for each of the coding gene. The clustering pattern derived from N and L genes was similar to M, P, F and H genes. Since N and L genes-based topology of phylogenetic relationship among geographically distinct strain was found to be more precise and conclusive, the L gene is suggested to be employed in future epidemiological investigations.

Figure 2
figure 2

Individual coding gene-based phylogenetic analysis of so-far reported SRMVs revealed mismatching for monophyletic clustering of strains.

Full size image

Based on the analysis of the complete N gene dataset, we proposed a geography and timeline based-classification of SRMV strains within lineage IV. A substantial analysis revealed a 6% and 2% nucleotide divegernce as a considerable cut-off criterion for the clasification of SRMV lineages and sub-lineages, respectively. Further analysis identified a total of six sub-clades (a-f) where sub-clade “a” represented strains from India, Turkey and Israel during 1994–2017, sub-clade “b” contained Chinese strains reported in 2007–08, sub-clade “c” represented strains from Africa and Georgia during 2008–2016, sub-clade “d” had Chinese strains reported during 2013–2016, sub-clade “e” possessed strains reported from India during 2014–2016, and sub-clade “f” was exclusive to Pakistan-originated strains which were reported in 2015 (Fig. 3).

Figure 3
figure 3

The complete N gene-based intra-lineage classification of strains within lineage IV.

Full size image

Nucleotide diversity and selective pressure analysis

The average nucleotide diversity (Pi-value) was 0.04889 for complete genome of all SRMV strains. With a variance (0.00001) and standard deviation (0.002) for haplotype diversity (Hd = 1.000), the average nucleotide differences among all haplotypes was found to be k = 788.690. A total of 5891 mutations were observed in DnaSP analysis, where 10831 were monomorphic and 5117 were polymorphic. The polymorphic mutations consisted of 1311 singleton variable sites with 3806 parsimony informative sites. While an assessment for neutrality, the Tajima’s D value was found to be negative for all genes (p > 0.10). The reliability of the analysis, as determined by HKA test, was found to be 6.078 (X-square value) in T = 6.732 (divergence time) at a significant level (p = 0.0131) (Table 8). An analysis of the genetic diversity within the coding genes across the whole length of the genome revealed an occurrence of hotspot event (300 nt window size per ten nt overlapping steps) between 5ʹ UTR of M gene and 3ʹ UTR of F gene (Fig. 4). The nucleotide diversity across the coding genes of nucleotide sequence haplotypes was found to be highest in H gene (0.05171) followed by P (0.04527), F (0.04409), N (0.04068), L (0.03982) and M (0.03931) genes. On the other hands, the haplotype diversity (Hd) was observed to be higher in L gene (0.996) followed by F (0.926), P (0.921), H (0.920), N (0.909) and M (0.900) genes (Table 8).

Table 8 A brief description on genome polymorphism for secletion sites in the complete genome and each of the coding regions in SRMVs.
Full size table
Figure 4
figure 4

Nucleotide diversity plotamong whole genome sequences of SRMVs derived from DnaSP.

Full size image

Datamonkey output for selective pressure analysis across CDS regions is summarized in Table 9. Although none of the gene carried a mean dN-dS greater than 1 at p < 0.05, it was highest for P gene (0.44679) followed by H (0.20017), N (0.12168), F (0.10253), L (0.08976) and M (0.06601) genes. At p < 0.05, analyzing through different algorithmic approaches (SLAC, FEL, IFEL, REL and MEME), revealed that the L gene showed a highest positive selection sites (96) followed by N (27), F (21), P (16), H (12) and M (2) genes. The plots against codon positions for individual genes were drawn using SLAC statistical approach based on dN-dS value (Fig. 5).

Table 9 Data Monkey analysis based brief summary of positive and negative substitution sites in each of the coding gene of so far reported SRMVs.
Full size table
Figure 5
figure 5

Differences in codon position, synonymous and non-synonymous substitutions (dN-dS values) for each of individual genes.

Full size image

Recombination analysis

Lying between 5ʹ UTR of the M gene (3406–4888 bp) and 3ʹ UTR of the F gene (4892–7306 bp), apparently a putative recombination event was observed in the complete genome (4607–5425 nts) of Pakistan-origin strain of SRMV but is attributed to bioinformatics errors. Therefore, no recombination was found in the current study.

Discussion

We presented a comparative genetic, phylogenomic and evolutionary analysis of SRMV strains reported so far in public database. Whole genome sequences and open reading frames (ORFs) of individual genes of representative strains were used in subsequent higher-resolution bioinformatic analysis. This is because a specific gene might not evolve at the same rate as does the whole genome15 and, therefore, can provide precise information on viral evolutionary dynamics and necessary epidemiological investigations in future16. While considering the “rule of six” for whole gnome atlas, comparative complete genome analysis revealed a varying length of complete genome suggesting the potential of the virus to evolve over a period of time. A few sequences showed unsusual lengths (e.g., MF678816; 15927 bp, KM089831; 15957 bp and KM816619; 16058 bp) where, for each of these sequences, a nucleotide insertion/deletion was observed in the noncoding region between the M and F genes17,18. Interestetingly, each of these sequeunce was deriven from the next generation sequencing approach and, therefore, such an unusual length may correspond to the sequencing errors. Owing to the fact that all paramyxovirus including SRMV follow a polyhexameric genome length for the effective replication in host cells19, SRMV sequences erroneously not following the “rule of six” in genome atlas were excluded from the specific analysis.

Comparative residue analysis of viral proteins showed several conserved motifs20,21. Among these, the N protein had three conserved motifs. These included export signal, nuclear localization signal and RNA binding motif. The first two are considered responsible for transport of the N protein to nucleus of host cell, while the third one was believed to be involved in interaction of N-N monomers of RNA during genomic RNA binding and N-N self-interaction20. Developing polymerase complex with N and L proteins, the P protein plays a significant role in virus replication and RNA biosynthesis22. The protein contains a variable N-terminus whereas C-terminus is believed to be the most conserved, and is required for the interaction with L protein in synthesis of polymerase complex23. The Soyuz 1 motif and presence of 151S residue, responsible for viral transcription via altering its phosphorylation status24, were found in all study-included strains21. The M protein is a core organizer of viral morphogenesis and has the ability to interact with other proteins for maturation of viral progeny25. For all of the investigated strains, this protein carried a previously known residue pattern21 for late domain or cell membrane anchor, which has a known role for localization of cell membrane and budding activity26. An unusually long and GC rich non-coding region was observed between 3ʹ UTR-M and 5ʹ UTR-F genes in studied SRMV sequences. While no biological or functional significance is warranted, a previous study has suggested an up- and/or down-regulation of these proteins to differences in their lengths and therefore may alter cyto-pathogenicity and survival fitness of the virus in nature27.

Three motifs were also noticed in the F protein as signal peptide, cleavage site (responsible for virulence and adaptation in the environment) and a leucine zipper domain. These are known to be involved in maintenance of protein tertiary structure20,22. Since the signal peptide motif was located in a variable region28, we performed a comparative analysis to investigate the conserveness of specific residue at a specific position among all reported strains from different geographies and proposed a stretch of consensus residues at the global level. The H protein is considered responsible for attachment of the virus to host cell membrane via cleavage of sialic acid residue in cellular glycoprotein29. As observed in the current study, the protein has a hydrophobic domain at the N-terminus that acts as a signal peptide to anchor the protein into the membrane20. The findings of SLAM receptor binding sites during the analysis highlight the epitheliotropic and lymphotropic nature of SRMVs30. Herein, a high number of glycosylation sites were found in the N protein, which plays a major role in protein translocation31. The large protein (L) contributes in viral replication, transcription and polyadenylation using different domains that were found to be conserved in this study. Domain I, II and III are considered responsible for polymerase and kinase activity where GDDD and QGDNQ residues carry a prime significance32 and, as observed in a previous experimental study33, any substitution in these residues can abolish the polymerase activity of the L protein. Two highly conserved hinge regions were also observed in a pattern typically corresponding to established hinge regions of other closely related morbilliviruses34. Taken together, the potential influence of these substitutions in the functionality of corresponding proteins is scarce and, therefore, requires future investigations to determine impact of these variations in conserved domains. 

The phylogenetic analysis, either based upon complete genome or each of the complete coding genes, showed a clustering pattern according to distinct geographical setting and time-period e.g., strains clustered within a distinct clade represented same country of origin within a specific time period. Therefore, while presenting a global perspective, a clustering and subsequent sub-clade grouping is proposed in the current study as an imporved and updated version of previous proposal35. This is simply becuase the previous classification proposal was limited to sub-grouping of Indian strains along with a few of those reported from the Middle East and Africa. Not only that the said proposal excluded strains reported from China and Georgia but also did not represent a well-defined evolutionary cut-off for the lowest taxonomic node (sub-lineage or sub-grouping). In addition to that, Kumar et al.35 have classified the strains into clades and subclades which contradicts previously proposed standard classification criteria for the lowest taxonomic node or sub-grouping of the viruses within a lineage or genotype36,37. Though such a classification may provide some pre-liminary assessment exclusively for Indian-origin strains, a limited geographic-pattern based classification may raise controversies for SRMV classification at global scale. Therefore, these are considered unreliable to present molecular epidemiology of SRMVs worldwide. Indeed, with a substantial increase in the number of SRMV sequences in future, following a uniform classification criterion such as presented in the current study (IVa, IVb, IVc, IVd, IVe and IVf), is necessary for a more precise clustering at the lowest taxonomic node. While comparing different coding genes (P, M, F, H and L) of SRMV strains (Fig. 2a,b), minor differences were observed in the clustering pattern indicating an influence of nucleotides in genetic diversity of SRMVs. Nevertheless, the N gene-based topography was closer to those of the L gene (RNA-dependent RNA polymerase) and complete genome sequences. Thus it (L gene) could be employed alternatively for a precise evolutionary relationship of SRMV strains originating from different geographical regions. This is important because, considering SRMV a member of the family Paramyxoviridae, L protein is now considered as a standard criterion for classification of some of the closely related members of the sub-family Avualvirinae38. The observed topology of the N gene revealed evolutionary dynamics of circulating SRMV strains consistent with observations made previously10,20. Therefore, it is suggested that complete N and L gene-based phylogeny analysis can provide an accurate evolutionary relationship of the circulating strains in particular geographical settings10,38, especially for those regions where full-genomes have not yet been reported or have limited resources.

Nucleotide diversity analysis was used to unleash the genomic variation (polymorphism) within a given dataset39 where a substitution rate is considered a prime parameter to elucidate virus evolution over a period of time. The average number of pairwise nucleotide difference among the whole genome of all SRMV sequences was found to be 788.690 with a diversity in nucleotide sequences (0.04889 ± S.D. 0.00468) and haplotype variance (0.00001). Gained observations correspond to distinct features of RNA viruses where there is a lack of proofreading activity by reverse transcriptase40. In contrast to previous observations14, a lower diversity in nucleotide and haplotype variance, and nucleotide difference in the current study may largely be ascribed to inclusion of a smaller number of complete nucleotide sequences than those employed in the current study (n = 37 vs n = 68). In addition to this, evidenced by significant nucleotide diversity over a period of time (p < 0.05), the HKA test outcome indicated an ongoing evolution or adaptation of virus in the environment.

The DnaSP based nucleotide diversity analysis revealed higher diversity in the H gene than others of SRMVs. Owing to significant roles in attachment and subsequent genome replication, the gene has been proposed to assess the evolutionary relationship of SRMV strains41. Though it ascertains further research, the substitutions in the H gene may have an influence on host adaptability and pathogenicity to susceptible host such as observed previously for SRMV14 and influenza virus42. A diverse nucleotide hotspot was obsereved between 5ʹ UTR of M and 3ʹ UTR of the F genes in the whole genome. This aligns with observations made previously where a hotspot was identified at similar position between M and F genes14, highlighting potential variations in the genome size and corresponding substitutions43 in each of the gene. An influence of these spontaneous mutations in genome was assessed by employing Tajima’s D statistics that showed a non-significant negative value for all coding genes in DnaSP analysis, suggesting a lack of influence of spontaneous mutations on the fitness of individual virus. Such observations suggest positive selection among coding region of sequences with a lower level of sequence diversity and an excess of low-frequency variants reflecting the role of natural selection in SRMV genomes. Contrary to current study findings where analysis showed negative value for each of the coding genes, positive values in F and H genes has previously been suggested14.

The non-synonymous/synonymous rate (ω = dN-dS) is an important indicator of selective pressure at the protein level where ω = 1 means neutral mutations, ω < 1 correspond to purifying selection while ω > 1 indicates diversifying positive pressure44. Herein, as reported in a previous study14, the dN-dS plot for each protein showed value not more than 1 indicating a slow genetic evolution of SRMV. Indeed, such a comparison of rates of synonymous and non-synonymous mutations provides an understanding towards the mechanisms of molecular sequence evolution. The positive selection sites were found in all coding genes (N, P, M, F, H and L) using different statistical approaches. Though these sites were found to be non-significant with a ratio less than 1 by Tajima’s D statistics, it seldom happens in structural domains of genome. However, the impact of such positive selection sites with lower level of sequence diversity may cause the emergence of variants44. According to the neutral theory of molecular evolution, such type of molecular variations, which arise via spontaneous mutations, has no influence on individual’s fitness45. However, the biological significance of these sites still remains unknown and needs to be explored in future.

The occurrence of recombination events is considered a significant source of genetic diversity for RNA viruses46. Occurrence of recombination in negative-sense RNA viruses is extremely rare, still analysis for the detection of recombination event/s is recommended as a standard component of every phylogenetic analysis to serve an important quality-control function to weed out laboratory and analytical errors47. In the current study, no recombination event was found.

Materials and Methods

Complete genome sequencing of SRMVs from Pakistan and dataset information

The complete genome sequencing of two SRMV isolates [KY967609 (SRMV/Faisalabad/UVAS/Pak/2015) and KY967610 (SRMV/Layyah/UVAS/Pak/2015)] was performed as per primers and protocols described previously5. Later, including these two strains, a total of 75 whole genome sequences of SRMVs were accessed (https://www.ncbi.nlm.nih.gov/, October 01, 2019) and processed for subsequent bioinformatic analysis. Among these 75 SRMV sequences, four were attenuated vaccine strains (KJ867542, KF727981; HQ197753, X74443) and were excluded from the dataset used in the current study. Furtermore, given the “rule of six” genome atlas or polyhexameric genome length, 03 sequences including MF678816 (15927 bp), KM089831 (15957 bp) and KM816619 (16058 bp) were also excluded from comparative whole genome-specific analysis. However, owing to length of coding region comparable to each of the protein of SRMV, only the coding regions of these sequences were included and processed further in comparative genomic and residue analsyis. All essential information related to whole genome sequences of study-included strains is presented in Table 1.

Comparative genomic analysis

The complete genome (15954 bp) dataset was aligned to equal length using ClustalW methods in BioEdit version 5.0.648 and, based upon nucleotide number and position across the whole length of the genome, different genomic features were compared among all SRMV sequences. The consensus sequences were made for those SRMV sequences that had a highest nucleotide similarity and were originated from similar geographical regions. Nucleotide identity and divergence among all consensus whole genome sequences of lineages I-IV was assessed by Pairwise Sequence Comparisons (PASC) analysis in MEGA version 6.0649. The conserved domains, functional and structural motif/s, trans-membrane regions and unique substitutions in open reading frames were predicted using ORF Finder (http://www.ncbi.nlm.nih.gov/gorf/gorf.html), Conserved Domain Prediction tool (http://www.ncbi.nlm.nih.gov/Structure/cdd/wrpsb.cgi) and HMMTOP program (http://www.enzim.hu/hmmtop/index.php). The potential N-glycosylation sites (N-X-T/S, where X denoted any residue except a Proline) were predicted by NetNGlyc 1.0 server (http://www.cbs.dtu.dk/services/NetNGlyc) and accepted if the G-score was 0.5. Similarly, the diversity and/or conserveness of residues at important but hypervariable motif/s were analysed through WebLogo version 3.1 (accessible at http://weblogo.threeplusone.com/create.cgi).

Estimation of evolutionary and divergence dates

Using a Bayesian Markov Chain Monte Carlo (MCMC) approach implemented in Bayesian evolutionary analysis sampling trees (BEAST) software package version 1.8.050, the molecular evolutionary and divergence rates were co-estimated for complete genome and individual genes. For each dataset, a total of three independent runs of MCMC were conducted under a strict molecular clock model, using the Hasegawa–Kishino–Yano model of sequence evolution with a proportion of invariant sites and gamma distributed rate heterogeneity (HKY + I + C) with partitions into codon positions, and the remaining default parameters in the prior’s panel. For each gene, the MCMC run was 36107 steps long and the posterior probability distribution of the chains was sampled every 1000 steps. Convergence was assessed on the basis of an effective sampling size after 10% burn-in using Tracer software, version 1.5 (http://tree.bio.ed.ac.uk/software/tracer/). The estimations were the mean values obtained for the three runs. The mean time of the most recent common ancestor (TMRCA) and the 95% CI were calculated, and the best-fitting models were selected by a Bayes factor using marginal likelihoods implemented in Tracer51.

Phylogeography-based reconstruction of evolutionary tree

A reliability of a gene for molecular epidemiology was assessed by comparing all coding genes (N, P, M, F, H and L) extracted from whole genome sequence of SRMV and aligned separately by ClustalW methods incorporated in the BioEdit version 5.0.648. The phylogenetic trees were constructed by neighbour-joining method with best-fit substitution model for each set of sequences using MEGA version 6.0649. A 1000 replication bootstrap value was adjusted to better elucidate the probability and reliability of clustering of isolates or any change in their clustering pattern.

Nucleotide diversity and natural selective pressure analysis

Based upon variable sites for mutations, and average numbers of pairwise nucleotide differences, the nucleotide diversity among coding sequences (CDS) of complete genome sequences was assessed for genomic polymorphism by DnaSP version 5.10.01 (accessible at http://www.ub.es/dnasp). The departure from neutrality in all isolate’s sequences was tested by Tajima’s D statistical method52. Divergence time in nucleotide diversity was estimated by a direct statistical model (HKA test). Data-monkey adaptive evolution server (http://www.datamonkey.org/) was used to evaluate synonymous (dS) and non-synonymous (dN) substitution rate per codon among CDS of all sequences53. Later, the positive and negative selection sites under natural selection were determined through six different genetic algorithms including Single Likelihood Ancestor Counting (SLAC), Fixed Effect Likelihood (FEL), Internal Branch Fixed Effect Likelihood (IFEL), Random Effects Likelihood (REL), Mixed Effect Model of Episodic selection (MEME) and Fast Unbiased Bayesien Approximation (FUBAR)54.

Detection of putative recombination event

The sequences were analyzed for the identification of reliable putative breakpoints by different tools including SimPlot version 3.5.155, GARD (http://www.datamonkey.org/GARD), DAMBE version 5.2.3056 and RDP4 version 4.9557. However, owing to an enhanced accuracy, clarity and reliability of analysis, outcomes gained by RDP4 were considered conclusive for further interpretation. The RDP4 was preferred because it employs a combination of seven different algorithms named RDP, GENECONV, BootScan, MaxChi, Chimaera, SiScan and 3Seq to better unleash putative recombinant and parent isolates at p < 0.001. A putative recombination event was assumed to have occurred only when it was consistently identified by at least four of the above-mentioned algorithms at a probability threshold of 0.05.

Ethical approval and informed consent

This research did not involve human participants or animals. This article does not contain studies with animals or humans performed by any of the authors.