Rabies is one of the most dreadful zoonotic diseases, exhibiting almost 100% fatality rates in infected mammals if left untreated [1]. The etiological agent of rabies, Rabies Lyssavirus, otherwise called as Rabies Virus (RABV), is a nonsegmented negative-strand RNA virus, with the broadest host range among lyssaviruses, including humans, various domestic animals, livestock, and wildlife. In Georgia, rabies has been officially documented since 1930, where nearly one third of animal rabies cases were found to occur in livestock, including 83% of cases associated with cattle [1]. This pilot study provides some initial but important information on the population structure and mechanisms of evolutionary divergence of RABV strains recovered in Georgia, and on those of their genetically most closely related conspecifics from elsewhere.
We analyzed the DNA sequences of the nucleoprotein-encoding (N) gene loci (the 1350-bp region) of a total of 72 RABV isolates, recovered (during 2015-2016) from dogs, jackals, and cattle in Georgia. The information on both the above DNA sequences and the sources of these isolates (Table S1) was retrieved from the GenBank database of the National Center for Biotechnology Information (NCBI, https://www.ncbi.nlm.nih.gov/). The above sequences were submitted to the NCBI GenBank database by the Pox and Rabies branch, Centers for Disease Control and Prevention, Atlanta, GA, the United States of America. The Maximum Likelihood (ML) method, implemented in the Molecular Evolutionary Genetics Analysis software MEGA X (v. 10.1.8) [2], was applied to determine STs of the N gene loci of these RABV isolates, and to construct a phylogenetic tree. In the phylogenetic analysis, we also included their genetically most closely related conspecifics from elsewhere. The DNA sequences of these conspecifics were obtained using the Basic Local Alignment Search Tool (BLAST) with the megablast algorithm searching for the genetically closest homologs of the N gene loci of the determined RABV STs in the blastn database.
The Tajima's Neutrality Test [3] was applied to determine a number of segregating sites and an average number of nucleotide differences across the targeted N gene loci of the RABV isolates. Besides, the FUBAR (Fast, Unconstrained Bayesian AppRoximation) model [4] was applied to infer nonsynoymous (dN) versus synonymous (dS) substitution rates, as well as the dN/dS ratios (ա), on a per-site basis for the coding alignment of the above genetic loci. The FEL (Fixed Effects Likelihood) model [5] was applied to infer dN versus dS substitution rates on a per-site basis for the same coding DNA sequence alignment across the RABV phylogenetic clusters. Using the Tamura-Nei model (+G) [6], we also estimated a shape parameter for the discrete Gamma Distribution in these N gene loci. The Kimura 2-parameter model [7] was applied to determine a transition/transversion (ti/tv) bias (R) across the targeted nucleotide sequences by computing a tree topology with the maximum Log likelihood 3605.928. A maximum composite likelihood estimate of the pattern of nucleotide substitution was also determined as described previously [8].
The method of split decomposition [9], implemented in SplitsTree (v.4.14.4) [10], was used to reconstruct possible genetic recombination events of the N gene loci between the RABV strains. For the SplitsTree-generated splits graphs, the bootstrap values being ≥ 95 (from 1000 replicates) for each node of a parallelogram(s), and the fit values being ≥ 95 for each splits network, were considered to be statistically significant. SimPlot (v3.5.1) [43] [11] was applied to detect recombination hotspots across the above N gene loci.
The DNA sequence analysis discriminated the local 72 RABV isolates into 41 STs (Table S1). As shown, while certain STs of the N gene loci (STs: 9, 10, 12, 13, 15-20, 23, 25-28, 30-32, 35, 36, 38, 39, and 41) were uniquely associated with the RABV isolates from dogs, several STs (STs: 1, 7, 11, 29, and 37) were shared by the isolates recovered from both dogs and cattle. Two RABV isolates (MW055109.1 and MT079888.1), recovered from the jackals, were found to share the ST22 with the isolate (MT079902.1) from a dog. STs 4, 5, 6, 8, 33, and 34 appeared to be associated with the RABV isolates recovered from cattle in the neighbor regions across the Western Georgia (Imereti and Guria) or the north-western part of the country (Samegrelo-Zemo Svaneti), while STs 14, 21, 24, and 40 were uniquely associated with the South-East region (Kvemo Kartli) or the eastern part (Kakheti) of Georgia.
The BLAST-identified genetically most closely related conspecifics of the local RABV isolates were associated mainly with the strains from the neighbor countries, such as Russia, Azerbaijan and Turkey followed by Iran, Iraq, Hungary, and Tajikistan (Table S2), sharing the DNA-DNA identity in a range of 96.74-99.93%. The ancestral ST and the most recent ancestral ST were found to be associated respectively with the RABV strain (KY860610.1) isolated from a dog in Turkey, and the strain (MK598340.1) isolated from a red fox in Hungary, when the above conspecifics of the local RABV isolates were analyzed alone using ML. However, the ancestral ST versus the most recent ancestral ST appeared to be associated respectively with the strain (JQ944706.1) recovered from a dog in Russia, and two strains (MT079899.1 and MT079940.1 sharing the same ST37) from a dog and cattle in Georgia, when the entire RABV population, including the local isolates, were analyzed using the same algorithm. The phylogenetic analysis resulted in one large genetic clade with multiple clusters and sub-clusters incorporating a great majority of the STs, and one relatively distinct phylogenetic clade including only four STs 39, 41, 40 and 38 (Fig. 1). As shown, the former also included the RABV strains recovered in Turkey and Russia, while the latter encompassed the strains recovered in Azerbaijan, Hungary, and Tajikistan. For a better phylogenetic resolution, in this analysis, we also included multiple RABV isolates from Iran, because, based on the BLAST analysis, certain strains from this country were genetically closely related to the ones from Turkey, while others to a single strain from Iraq, constituting collectively another distinct phylogenetic clade (Fig. 1). The above phylogenetic clades may represent the RABV evolutionary lineages, although, a significantly larger strain collection from all these countries, including Georgia, need to be analyzed to support the above statement.
The analysis of molecular evolutionary patterns of the N gene loci of RABV isolates revealed that the ti substitution rates were notably higher than the tv substitution rates, especially the T (U)→C and C→T (U) rates (Table S3). The maximum composite likelihood estimates of the pattern of nucleotide substitutions resulted in the following nucleotide frequencies: 29.64% (A), 26.28% (T/U), 20.56% (C), and 23.52% (G). The ti/tv rate ratios respectively were k1 = 9.007 (purines) and k2 = 17.013 (pyrimidines), with the overall ti/tv bias being R = 6.214. The FEL-inferred dN/dS rates were in a range of 0.0035-0.05. The FUBAR could reveal pervasive negative/purifying selection at 147 sites across the N gene loci examined. No evidence for pervasive positive/diversifying selection in these gene loci could be found in the FUBAR analysis (Fig. S1). The Tajima’s neutrality test of the nucleotide variations of the N gene loci resulted in a negative value: D = - 1.237845 (Table 1).
SplitsTree generated the splits graphs with the strongest fit (100), exhibiting multiple parallelograms, some nodes of which were supported by the strong bootstrap values (94.9-99.9) (Fig. 2), when examining the above N gene loci. As shown, the genetic recombination events of the N gene loci, displayed in a capacity of these parallelograms, involved not only certain STs of the RABV isolates from Georgia (ST40, ST21, ST 15, ST17) (Fig. 2a), but also multiple strains from elsewhere, mainly from Iran (Fig. 2b-d). SimPlot could determine multiple genetic recombination hotspots across the N gene loci of some of these RABV strains (Fig. S2a-d). A great majority of the recombining RABV isolates were from dogs.
Based on the results of the above analysis having identified the specific STs shared by multiple RABV isolates from two different isolation sources (dogs and cattle), it is suggested that rabid dogs could contribute to the transmission of rabies in cattle in Georgia. Besides, as shown, while the distribution of RABV strains reflect some regional patterns (e.g., in Kvemo Kartli and Kakheti), certain strains from dogs seem to have spread across some regions of this country (such as Imereti, Guria, and Samegrelo-Zemo Svaneti). We assume that at least stray dogs could be the transmission vehicles of these strains across the above regions.
It must be noted that the ST22, shared by certain RABV strains recovered from the dog and jackals, belongs to one of the RABV cosmopolitan phylogenetic clades described previously [12]. While the immediate sources of rabies acquisition by the above animals remain unknown, jackals, as being generally a part of the dominant wildlife reservoirs for emerging rabies in domestic and agricultural animals [13], can be thought to contribute to this disease transborder transmission involving Georgia. As shown, the genetically most closely related conspecifics of the local RABV isolates appeared to be associated with Russia, Azerbaijan, Turkey, Tajikistan, Iran, Hungary or Iraq. The ML analysis results of the ancestral STs can be suggestive of a history of transborder transmission of certain RABV strains from Russia, Turkey, and Hungary, including Georgia that seems to have contributed to the recent events of the above process exhibiting this disease global transmission [14, 15].
The ti/ tv substitution rates, determined for the targeted N gene loci of the RABV strains in our analysis, are consistent with the results from the previous studies demonstrating a significantly higher frequency of C→U changes compared to that of their reverse (U→C) and other tis changes at site in many mammalian RNA viruses with <5% of heterogeneity [16]. Here, we show that the U→C and C→U substitutions were one of the major drivers of DNA sequence diversification of the N gene loci of the RABV strains examined. Besides, the ti/tv bias estimate, obtained from our analysis, was also very close to the ti/tv value (5.682) determined previously for the N gene loci of the RABV strains recovered elsewhere [17], thus suggesting some common evolutionary trends of this viral agent at least across certain regions globally.
The previous studies show that the evolution of the RABV N gene has been strongly subjected to purifying selection pressures [18], and that it has been highly constrained especially at nonsynonymous sites, with no clear evidence for positive selection [19]; however, it is noteworthy that some codons of this gene have been affected by nonsynonymous substitutions et relatively higher rates, which may signify the existence of localized selection pressures [17]. In our analysis, the FEL-inferred dN/dS rates on the per-site basis revealed no clear tendency towards amino acid changes, being further supported by the FUBAR inferences exhibiting pervasive negative/purifying selection at multiple sites across the targeted N gene loci. In contrast, no evidence for pervasive positive/diversifying selection could be found in the above gene loci. The Negative Tajima’s D negative value, obtained from the analysis of the RABV strains in this study, could be suggestive of negative selection and recent selective sweeps [20], as well as of a recent population expansion [21].
The transmission of RABV strains in non-flying mammals was shown to occur at a low rate at a global level, resulting in the phylogenetic structure that reflects a significant population subdivision rather than gene flow [22]. However, some studies provided initial evidence for genetic recombination events involving the polymerase gene [23] and the N gene [24] in some natural populations of RABV. Here, the results obtained from our SplitsTree and SimPlot-analyses, provide additional evidence for intra-species recombination, being in a strong agreement with the previous notion [24] suggesting that the N gene could be affected by homologous genetic recombination across the RABV natural populations. The above findings pinpoint collectively to the presence of rabies coinfections in certain rabid animals, such as e.g., rabid dogs. However, more in-depth studies are needed to determine whether HGT contributes to virulence and/or pathogenicity of this viral pathogen, as well as to its adaptation to various host environments.