Rabies is a zoonotic disease with an estimated global mortality of 59,000 people annually and a burden of more than 3.7 million disability-adjusted life years (DALYs) that is caused by a neurotropic lyssavirus^[1]. Dogs are the primary source of human rabies, as more than 95% of human cases can be traced to dogs^[2,3]. China faces a substantial burden of rabies, having endured three major human rabies epidemics, which occurred in the 1950s, 1981, and 2007^[4]. Implementation of various prevention and control measures has decreased the number of human rabies cases from 3,300 in 2007 to 167 in 2023. In China, Jiangsu is one of the regions most affected by rabies, although reported cases are dropping, from 201 cases in 2004 to 11 cases in 2023. However, it is worth noting that case numbers were slightly higher in 201, than in the previous year. According to statistics from the Jiangsu Health Commission, rabies ranks among the top five causes of death from notifiable infectious diseases in the province, and more than 90% of cases can be attributed to dog bites. Although the incidence of human rabies has decreased significantly, continued efforts are required to eliminate canine-transmitted rabies.

To date, the International Committee on Taxonomy of Viruses (ICTV) has recognized 17 Lyssavirus species, which are classified into three phylogroups based on serological and genetic distances. Rabies virus (RABV), in phylogroup I, is the primary pathogen responsible for dog-mediated rabies^[5]. RABV is a single-stranded negative-sense RNA virus with a ~12 kb genome that encodes five proteins: nucleoprotein (N), phosphoprotein (P), matrix protein (M), glycoprotein (G), and RNA polymerase (L)^[6]. Based on the global classification of RABV, all isolates from China are classified as China I–VII. China I, II, V, and VI are subclades of the Asian lineage; China III corresponds to the Cosmopolitan lineage; China IV is part of the Arctic/Arctic-like lineage; and China VII corresponds to the Indian subcontinent lineage^[7]. Understanding the genetic evolution of this pathogen is crucial for tracking variants and transmission pathways and can provide a molecular foundation for optimizing disease control strategies.

In this study, we analyzed 38 samples of dog brain tissue collected in Jiangsu Province between 2017 and 2020, two of which tested positive for RABV. Phylogenetic analysis of the whole-genome sequences of the two positive samples showed that both belonged to the China I lineage. We then conducted a genetic evolutionary analysis of all 37 RABV isolates collected in Jiangsu Province, including the isolates from our two positive samples. The results showed that there are three RABV lineages in Jiangsu Province: China I, II, and III. China I is the most prevalent lineage in Jiangsu Province, with transmission occurring mainly in dogs. This analysis offers valuable insights for the prevention and control of RABV and may aid in the development of more effective strategies. It aims to provide scientific evidence to optimize control policies in line with the goal of eliminating canine-transmitted rabies by 2030.

We collected brain tissue samples from 38 dogs in Jiangsu Province between 2017 and 2020 and used direct fluorescent antibody (DFA) and quantitative real-time PCR (qPCR) methods to detect RABV in the samples^[8].

We also collected 36 RABV whole genome sequences from GenBank, representing the seven RABV lineages prevalent in China (Supplementary Table S1). These 36 sequences, along with the sequences of the two new isolates were used as Database 1. All RABV isolates from Jiangsu Province (n=35) were also downloaded from GenBank (Supplementary Table S2). These 35 sequences and the two new sequences were used as Database 2. MAFFT software (v7.450) was used for multiple sequence alignment, and IQ-TREE (v1.6.12) was used to reconstruct maximum likelihood (ML) trees. The nucleotide and amino acid similarities of the N and G gene sequences were calculated using Geneious Prime software (v. 2025.0.2). ChiPlot (www.chiplot.online) was used to generate phylogenetic trees. A maximum clade credibility (MCC) tree was created using Beast software (v1.10.4) based on the G gene sequences of all RABV isolates from Jiangsu Province. FigTree (v1.4.4) was used to annotate and visualize the final tree. Protein structure was predicted using SWISS-MODEL (swissmodel.expasy.org), and key sites were visualized and annotated using PyMOL (v 3.1.6.1). Selection pressure was calculated using the online tools at datamonkey (www.datamonkey.org).

The DFA analysis showed clear apple green fluorescence signals in samples JSNJ1712 and JSNJ1801, indicating strong positivity for RABV (Figure 1A). The qPCR assay yielded Ct values of 12 for both JSNJ1712 and JSNJ1801. Sequencing analysis showed that the JSNJ1712 genome is 12,001 bp in length, containing the following gene-coding regions: N (1,353 bp; 85–1,437 bp), P (894 bp; 1,529–2,422 bp), M (609 bp, 2,518–3,126 bp), G (1,575 bp; 3,338–4,912 bp), and L (6,387 bp; 5,429–11,815 bp). The JSNJ1801 genome is 12,054 bp in length, and the coding regions for the N, P, M, G, and L genes are located at 157–1,509 bp, 1,601–2,494 bp, 2,582–3,190 bp, 3,402–4,976 bp, and 5,493–11,879 bp, respectively. We conducted a phylogenetic analysis of the two new isolates, which showed that JSNJ1712 and JSNJ1801 belonged to the China I group (Figure 1B). The nucleotide sequence similarity between the two isolates was 99.87%, and they shared the highest homology with isolate CSD0708D from Shandong Province.

Figure 1.  Detection of rabies virus (RABV) in dog brain tissue samples, and phylogenetic analysis of RABV from China. (A) Detection of RABV nucleoprotein antigen using direct fluorescent antibody (DFA) analysis. (B) Maximum likelihood (ML) tree of partial viral isolates using the full genome sequences from seven lineages in China (The green five stars represent the samples tested in this study; Bootstrap values = 1,000).

Because some isolates contained only N or G gene sequences, we conducted phylogenetic analyses using both the N and G genes. The results (Figures 2A, 2B) showed that the circulating RABV isolates in Jiangsu Province were primarily within China I, II, and III. Previous studies identified Jiangsu Province as a major source of China I lineage transmission to several eastern provinces, including Shandong, Fujian, Beijing, Anhui, Shanghai, and Zhejiang. The China II lineage in Jiangsu Province was transmitted from Henan^[9]. The China III lineage is mainly distributed in Inner Mongolia and Xinjiang, and it likely was introduced from Russia and Mongolia. However, its spread to other provinces has not been examined. The nucleotide and amino acid similarities of the RABV N and G genes in Jiangsu Province were analyzed, which showed that the N genes have nucleotide similarities between 86.29% and 100% and amino acid similarities from 95.77% to 100% (Figure 2C) and the G genes have nucleotide similarities ranging from 82.73% to 100% with amino acid similarities from 89.31% to 100% (Figure 2D).

Figure 2.  Phylogenetic and similarity analyses of RABV isolates in Jiangsu Province. (A and B) Maximum likelihood (ML) trees showing the genetic relationships based on the (A) RABV G gene and (B) RABV N gene, with detailed information on the isolates, including the host and year of isolation (Bootstrap values = 1,000). (C) Heatmaps of pairwise nucleotide and amino acid similarities in the G gene. (D) Heatmaps of pairwise nucleotide and amino acid similarity in the N gene.

Based on the GTR + G + I best-fit substitution and strict clock models, we constructed an MCC tree for the G gene sequences in RABV isolates from Jiangsu Province, with the time to the most recent common ancestor (tMRCA) shown at the major branch nodes (Figure 3A). Because there was only one sequence from the China III lineage, we focused on analyzing the tMRCA of the China I lineage. The tMRCA for China I is estimated to be 1971 (95% highest posterior density [HPD] interval: 1948.3 to 1987.2), with an average evolutionary rate of 3.729 × 10^-4 substitutions per site per year (95% HPD interval: 2.063 × 10^-4 to 5.835 × 10^-4). The Bayesian skyline plot (Figure 3B) revealed that the genetic diversity in this lineage remained relatively stable before 1960; however, between 1960 and 2010, it experienced a sharp decline, followed by a rebound and eventual stabilization, which aligns with the epidemic patterns observed in China. This can be attributed to the government’s implementation of strict dog management measures, including mandatory vaccination, culling of stray dogs, and the promotion of proactive vaccine use. However, the subsequent rise in genetic diversity corresponds with China's ongoing urbanization, which has resulted in more relaxed dog management practices in certain areas, particularly in remote rural regions, thereby facilitating the spread of the virus in these locations.

Figure 3.  Bayesian analysis and protein structure prediction. (A) Maximum clade confidence (MCC) trees (The red text indicates the new isolates in this article). (B) Bayesian skyline plot of RABV G genes. (C) Display of key amino acid sites on predictive protein models (Mutation sites are shown in red, the three main antigenic sites are shown in blue, and the two less-commonly recognized sites are shown in green).

To assess the selection pressure at each site of the G gene, we used the Mixed Effects Model of Evolution (MEME); Fixed Effects Likelihood (FEL); Fast, Unconstrained Bayesian AppRoximation (FUBAR); and Single-Likelihood Ancestor Counting (SLAC) methods on the Datamonkey website. Only the MEME method detected four positively selected amino acids (aa) sites: 54, 274, 305, and 490, whereas other sites were under negative selection. However, these mutations are not present at antigenic sites^[10], indicating that the existing vaccines can provide protection. Alignment of the amino acid sequences revealed that C54S was present in the 9811CHI isolate, D274V and A305V in the Jiangsu_Wx1 isolate, and Q490H in the JSNJ1712, JSNJ1801, CJS0841D, and SJS0842D isolates. Subsequently, we predicted the protein structure of the G protein based on the template (PDB: 7u9g) and selected the most reliable model based on the GMQE and QMEAN scores. Amino acid 305 is located in the α-helix, and amino acids 54 and 274 are found in the coil (Figure 3C). Because the amino acid at position 490 is located at the C-terminus, and it is difficult to accurately model the sequence at the C-terminus due to certain factors limiting model prediction, we did not include position 409 in the visualization. Since the coil and α-helix are involved in maintaining protein stability and function, mutations at these sites may help the virus adapt to different hosts or environmental conditions, thus continuously optimizing its functionality during evolution.