Evolutionary dynamics of recent peste des petits ruminants virus epidemic in China during 2013-2014.

Peste des petits ruminants virus (PPRV) causes a highly contagious disease, peste des petits ruminants (PPR), in sheep and goats which has been considered as a serious threat to the local economy in Africa and Asia. However, the in-depth evolutionary dynamics of PPRV during an epidemic is not well understood. We conducted phylogenetic analysis on genomic sequences of 25 PPRV strains from China 2013-2014 outbreaks. All these strains clustered into a novel clade in lineage 4. An evolutionary rate of 2.61 × 10-6 nucleotide substitutions per site per day was estimated, dating the most recent common ancestor of PPRV China 2013-2014 strains to early August 2013. Transmission network analysis revealed that all the virus sequences could be grouped into five clusters of infection, suggesting long-distance animal transmission play an important role in the spread of PPRV in China. These results expanded our knowledge for PPRV evolution to achieve effective control measures.

Peste des petits ruminants (PPR) is one of the most important diseases of sheep and goats, which is caused by a morbillivirus, peste des petits ruminants virus (PPRV). Alongside sheep and goats, PPRV also affects wildlife species including captive, wild small ruminants and sub-clinically cattle and buffalo (Parida et al., 2015). PPR can cause a very high mortality reaching up to 100% in immunologically naïve populations. Thus, PPR is regarded as a significant hurdle to the development of sustainable agriculture. In March 2015, Office Internationale des Epizooties (OIE) and Food and Agriculture Organization (FAO) officially launched a global program to eradicate PPR by 2030. The genome of PPRV is 15,948 or 15,954 nucleotides (nt) in length and organized into six transcriptional units encoding six structural proteins, the nucleoprotein (N), the phosphoprotein (P), the matrix protein (M), the fusion protein (F), the hemagglutinin protein (H) and the polymerase or large protein (L) (Bailey et al., 2005, Bao et al., 2014). The P transcription unit also encodes two non-structural proteins, C and V.

After the first report of PPR in 1942 in Côte d'Ivoire, the disease has been reported in most parts of sub-Saharan Africa, North Africa, the Middle East, South Asia, Central Asia, and East Asia (Banyard et al., 2014, Dhar et al., 2002, Kwiatek et al., 2011). The first outbreak of PPR in sheep and goats in China was reported in Tibet, China in 2007 (Wang et al., 2009). Another outbreak of PPR in wild small ruminants in Tibet, China was reported in 2008 (Bao et al., 2011). No further spread of the disease has been reported in China until late 2013. On November 30th, 2013, an outbreak of PPR was reported by the Ministry of Agriculture (MOA) of China in a farm in Huocheng county of Yili region in Xinjiang province. Four more infected farms were identified in December 2013 in the same province. The disease has rapidly and widely spread to 21 provinces along with the movements of goats and sheep (Banyard et al., 2014, Wang et al., 2015, Wu et al., 2015). This epidemic of PPR in China continued to the end of June 2013, causing death of more than 16,000 sheep and goats. The spread of PPR in China has posed serious threat of this highly contagious disease on the large population of domestic small ruminants (more than 300 millions sheep and goats) and different species of wild small ruminants in China and neighboring countries.

Phylogenetic studies using partial N or partial F gene sequences have grouped PPRV strains into four lineages (Kwiatek et al., 2007, Shaila et al., 1996). The first three lineages were distributed in Africa. And lineage IV included strains from Asia and recently northern Africa (Kwiatek et al., 2011). Recent phylogenetic analysis revealed that the Chinese Tibetan strains and the Chinese 2013–2014 strains were separately grouped into two clusters in lineage IV (Wu et al., 2015). Historically, PPRV molecular epidemiology has focused on the partial N (255 nt region; nucleotide site 1360–1614 of PPRV genome) or partial F gene (322 nt region; nucleotide site 5779–6100 of PPRV genome) sequences, preventing the in-depth evolutionary analysis among strains (Cosseddu et al., 2013, Munir et al., 2011, Ozkul et al., 2002, Padhi and Ma, 2014, Soltan and Abd-Eldaim, 2014). Genome sequences for PPRV have recently been used for an evolutionary study on a global scale (Muniraju et al., 2014c). However, the evolutionary dynamics of PPRV in an epidemic outbreak has not been studied.

Analysis on the evolution of PPRV genome during the course of an epidemic may help the interpretation of field epidemiology data and the implementation of efficient control measures. Such in-depth analysis has been conducted for several other RNA viruses, such as human respiratory syncytial virus, hepatitis C virus, highly pathogenenic avian influenza virus and food-and-mouth disease virus (Agoti et al., 2015, Bataille et al., 2011, Lu et al., 2015, Otieno et al., 2016, Valdazo-Gonzalez et al., 2012). The outbreak of PPR in China during 2013–2014 provided a unique opportunity to study the in-depth epidemiological and evolutionary dynamics of PPRV during an epidemic. The diagnosis and the genetic characterization of the PPR viruses being involved in this epidemic have been reported (Wu et al., 2015). However, the molecular evolution of PPRV between infected farms during the outbreak remains unclear.

In this study, we present full-genome sequence data for the PPR epidemic in China from November 2013 to June 2014. Phylogenetic analysis was used to study the evolutionary dynamics of PPRV during an epidemic.

Material and methods

Field samples

Full-length PPRV genome sequences were achieved directly from clinical samples of infected animals taken from 25 infected farms in 21 provinces. The infected farms included in this study were identified from MOA's Animal Disease Surveillance PPR Data Archive (http://www.syj.moa.gov.cn/dwyqdt/). Clinical samples were screened with real-time RT-PCR with amplification targeting N gene for confirmation. Although 3–5 animals were sampled for diagnosis in each infected farms, samples from one PPR-confirmed animal from each infected farm were used for genome amplification. We used clinical samples of nasal swabs, mesenteric lymph node and spleen. All clinical samples used for genome amplification were detailed in Table 1.

Table 1

Details of the 25 PPRV strains analyzed in this study.

Strain	Collection date	Collection location	Species	Sexa	Sample	Accession no.	Clade
China/XJYL/2013	11/30/2013	Xinjiang	goat	M	mesenteric lymph node	KM091959	A
China/XJ2/2013	12/20/2013	Xinjiang	goat	M	spleen	KX421384	D
China/XJ3/2013	12/21/2013	Xinjiang	sheep	F	mesenteric lymph node	KX421385	C
China/XJ4/2013	12/22/2013	Xinjiang	goat	F	spleen	KX421386	D
China/XJ5/2013	12/29/2013	Xinjiang	goat	M	mesenteric lymph node	KX421387	E
ChinaGS2014	1/22/2014	Gansu	sheep	M	mesenteric lymph node	MF443351	E
ChinaNX2014	2/17/2014	Ningxia	sheep	M	mesenteric lymph nod	MF443340	A
ChinaLN2014	3/17/2014	Liaoning	goat	M	mesenteric lymph node	MF443341	E
ChinaCQ2014	3/30/2014	Chongqing	goat	M	mesenteric lymph node	MF443353	B
ChinaHLJ2014	3/31/2014	Heilongjiang	goat	F	nasal swab	MF443346	D
ChinaYN2014	4/1/2014	Yunnan	goat	M	spleen	MF443336	B
ChinaSaX2014	4/1/2014	Shaanxi	goat	F	mesenteric lymph node	MF443339	D
ChinaJX2014	4/1/2014	Jiangxi	goat	F	mesenteric lymph node	MF443342	E
ChinaJL2014	4/1/2014	Jilin	sheep	F	nasal swab	MF443344	D
ChinaJS2014	4/2/2014	Jiangsu	goat	F	mesenteric lymph node	MF443343	B
ChinaHeN2014	4/3/2014	Henan	goat	M	nasal swab	MF443347	D
ChinaHB2014	4/3/2014	Hubei	goat	M	mesenteric lymph node	MF443348	B
ChinaAH2014	4/3/2014	Anhui	goat	F	spleen	MF443354	E
ChinaSX2014	4/5/2014	Shanxi	goat	F	nasal swab	MF443337	E
ChinaGX2014	4/16/2014	Guangxi	goat	F	mesenteric lymph node	MF443350	E
ChinaGZ2014	4/21/2014	Guizhou	goat	M	nasal swab	MF443349	B
ChinaZJ2014	4/25/2014	Zhejiang	goat	F	nasal swab	MF443335	B
ChinaHN2014	4/25/2014	Hunan	goat	M	mesenteric lymph node	MF443345	D
ChinaGD2014	5/15/2014	Guangdong	goat	F	mesenteric lymph node	MF443352	E
ChinaSC2014	6/10/2014	Sichuan	goat	F	nasal swab	MF443338	B

F = Female; M = Male.

Genome sequencing

Viral RNA was extracted and used directly for viral genome determination. Fourteen pairs of oligonucleotide primers were used to amplify 14 overlapping fragments by reverse transcription-PCR as previously described (Bao et al., 2014). The PCR products were purified and sequenced with an ABI 3730XL genome sequencer (Applied Biosystems, USA). An addition of 15 PPRV genome sequences representing viral strains from different endemic country was obtained from GenBank. The GenBank accession number, the year and country of isolation were listed in Table 2. Vaccine strains were not included in this study.

Table 2

Details of the 15 PPRV reference strains used in this study.

Strain	Year	Country	Host	GenBank Acc. No.	Reference
China/Tibet/30/2007	2007	China	goat	FJ905304
China/33/2007	2007	China	goat	KX421388
China/Tibet/Bharal/2008	2008	China	bharal	JX217850	(Bao et al., 2012)
Côte d' Ivoire 89	1989	Cote d'Ivoire	goat	EU267273	(Chard et al., 2008)
Nigeria 76/1	1976	Nigeria	goat	EU267274	(Chard et al., 2008)
Turkey 2000	2000	Turkey	goat	NC-006383	(Bailey et al., 2005)
CIV 01 P 2009	2009	Cote d'Ivoire	goat	KR781451
Ethiopia 1994	1994	Ethiopia	goat	KJ867540	(Muniraju et al., 2014b)
Ethiopia 2010	2010	Ethiopia	goat	KJ867541	(Muniraju et al., 2014a)
Ghana NK1 2010	2010	Ghana	goat	KJ466104	(Dundon et al., 2014)
Oman 1983	1983	Oman	goat	KJ867544	(Muniraju et al., 2014b)
UAE 1986	1986	UAE	goat	KJ867545	(Muniraju et al., 2014b)
India TN Gingee 2014	2014	India	goat	KR261605	(Masdooq et al., 2015)
Uganda 2012	2012	Uganda	goat	KJ867543	(Muniraju et al., 2014b)
Morocco 2008	2008	Morocco	goat	KC594074	(Muniraju et al., 2013)

Phylogenetic analysis

All the genome sequences were aligned using MegAlign software in Lasergene package. A phylogenetic tree was inferred by maximum-likelihood method from the nucleotide alignment of genome sequences using MEGA version 4 software, assuming a TN93 model of base substitution (equal substitution rates among sites and between transitional and transversional substitutions) (Tamura et al., 2007). The statistical significance of phylogenies constructed was estimated by bootstrap analysis with 1000 repetitions.

The analysis of evolutionary rates and times to common ancestry

Temporal dynamics of PPRV were analyzed with time-resolved phylogenies using the Bayesian Markov Chain Monte Carlo (MCMC) method in the BEAST package 1.7.1 (Drummond et al., 2002). The JModelTest program was used to select the best-fitting nucleotide substitution model (Posada, 2009). The datasets were analyzed using the GTR+G substitution model under an uncorrelated exponential relaxed clock with an exponential growth model. The MCMC chains were run for 2 × 108 generations and sampled every 20,000 generation. Convergence was assessed from the effective sample size (ESS) with a 10% burn-in using Tracer v1.6. ESS value above 200 was accepted.

Bayesian Skyline Plot analysis

The Bayesian Skyline Plot, estimating the change in effective population size through time, was generated from the BEAST Bayesian MCMC output files using Tracer v1.6.

NETWORK analysis

The phylogenetic network was constructed using the Median Joining method in the program NETWORK (Hajimorad et al., 2003). This software uses a parsimony method to connect each sequence to its closest neighbor, and allows the creation of internal nodes, which could be interpreted as unsampled or extinct ancestral genotypes to link the existing genotypes in the most parsimonious way.

Nucleotide sequence accession numbers

The generated PPRV sequences were deposited in GenBank with the accession numbers given in Table 1.

Results

Genetic diversity in PPRV strains in China

A total of 25 complete genomes were obtained for the PPR viruses circulating in China from November 2013 to June 2014. These PPRV genomes were from 25 PPR infected farms in 21 provinces across China as shown in Fig. 1.

Fig. 1

Map indicating the locations of infected farms where the PPR viruses were collected and sequenced. Cases sequenced in this study are represented by black squares.

We firstly investigated the sequence differences in the partial N gene (nucleotide site 1360–1614 of these genomes) and partial F gene sequence (nucleotide site 5785–6106 of these genomes). Only one substitution site, C1457T, was found in the N partial region, showing two genotypes (T for strain ChinaGX2014 and C for other strains). One substitution site located in F partial region (C6090T), with a T in strain ChinaCQ2014 and C for other strains. We further studied the sequence differences along the whole genome. These 25 China PPRV genomes shared 99.86–99.99% nucleotide identities. A total of 96 substitution sites were found along the genome sequences of these 25 China PPRV strains ( Table 3). The number of pair-wise nucleotide differences across the entire genome varied from 1 (between strains ChinaJS2014 and ChinaHB2014) to 23 (between strains ChinaNX2014 and ChinaJL2014; ChinaJS2014 and ChinaGD2014). This increased resolution suggested the advantage for tracking PPRV short-term evolution using full-length genome sequence.

Table 3

Nucleotide and amino acid diversity of China 2013–2014 strains.

Region	Length (bp)	No. of variable sites among CN 2013–2014 strains		No. of variable sites between CN 2013–2014 and other strains
		Nucleotide	Amino acid	Nucleotide	Amino acid
Leader	55	1		0
N 3′ UTR	52	2		0
N cds	1578	4	1	11	5
N 5′ UTR	59	3		0
P 3′ UTR	59	1		0
P cds	1530	12	10	15	8
P 5′ UTR	66	0		3
M 3′ UTR	32	0		0
M cds	1008	5	1	9	0
M 5′ UTR	444	7		12
F 3′ UTR	639	12		21
F cds	1641	6	4	13	4
F 5′ UTR	136	5		5
H 3′ UTR	20	0		0
H cds	1830	7	4	14	10
H 5′ UTR	107	0		2
L 3′ UTR	22	0		0
L cds	6552	30	8	43	12
L 5′ UTR	69	1		2
Trailer	31	0		0
Total		96		150

A total of 32 substitution sites were found in the non-coding region. Sixty four substitutes were identified in the coding region, among which 28 substitutes were found non-synonymous (NS). The highest ratio of NS substitutes were found in P gene (10 NS sites out of 12 sites), followed by F gene (4/6 sites), H gene (4/7 sites), L gene (8/30 sites), N gene (1/4 sites) and M gene (1/5 sites). We investigated the distribution of substitution sites at different codon positions ( Table 4). A total of 17 nucleotide changes were found at the first codon position, among which 16 sites contributed to amino acid changes. At the second codon position, 16 nucleotide substitutions were found, with all of them contributing to amino acid changes. Although 37 nucleotide changes were found at the third codon position, only 2 sites lead to amino acid changes. The finding suggested that the nucleotide substitutions at the first and second codon position contributed mainly to the amino acid changes.

Table 4

Sequence variability at different codon position of China 2013–2014 strains.

CDS	No. of variable sites among CN 2013–2014 strains									No. of variable sites between CN 2013–2014 and other strains
	1st codon			2nd codon			3rd codon			1st codon			2nd codon			3rd codon
	nta	aab	%	nt	aa	%	nt	aa	%	nt	aa	%	nt	aa	%	nt	aa	%
N cds	0	0		1	1	100%	3	0	0%	3	3	100%	1	1	100%	7	1	14%
P cds	7	7	100%	2	2	100%	3	1	33%	4	3	75%	2	2	100%	9	3	33%
M cds	1	1	100%	0	0		4	0	0%	2	0	0%	0	0		7	0	0%
F cds	2	2	100%	1	1	100%	3	1	33%	2	1	50%	1	1	100%	10	2	20%
H cds	2	2	100%	2	2	100%	3	0	0%	4	4	100%	3	3	100%	7	3	43%
L cds	5	4	80%	4	4	100%	21	0	0%	6	3	50%	5	5	100%	32	4	13%
Total	17	16		10	10		37	2		21	14		12	12		72	13

nt = nucleotide.

aa = amino acid.

These 25 genomes obtained in this study were combined with the 15 genomes derived from GenBank to generate a set of 40 genomes from 11 different epidemic countries in Asia and Africa: India, Turkey, United Arab Emirates, Oman, China, Côte d'Ivoire, Nigeria, Ethiopia, Ghana, Uganda and Morocco. To identify the genome characteristics of the 25 China 2013–2014 PPRV genomes, we conducted multiple sequence alignments between these 25 sequences and the other 15 PPRV genomes. The genome length of all the China 2013–2014 PPRV strains was 15,954 nucleotides (nt) with a 6-nt insert (TCCCTC) at nucleotide 5215 within the 5’ UTR of F gene, comparing to the 15,948-nt-long genome of all the other PPRV strains. We identified a total of 150 nucleotide differences between the China 2013–2014 strains and all the other strains. Only 45 substitution sites were found in the noncoding region. In total, 105 substitutes were identified in the coding region, among which 39 substitutes were found non-synonymous (NS). The highest ratio of NS substitutes were found in H gene (10 NS sites out of 14 mutation sites), followed by P gene (8/15 sites), N gene (5/11 sites), F gene (4/13 sites), L gene (12/43 sites) and M gene (0/9 sites). The distribution of substitution sites at different codon position was found as following: 21 nucleotide differences (causing 14 amino acid changes) at the first base, 12 nucleotide (12 amino acid) differences at the second base, 72 nucleotide (13 amino acid) changes at the third base.

Phylogenetic analysis using maximum likelihood revealed that all the PPRV genomes from China between 2013 and 2014 could be grouped into a distinct clade in lineage IV. As shown in Fig. 2, lineage IV PPRV formed four clades: 4.1 (with strain Turkey 2000), 4.2 (with strains from Africa), 4.3 (with strains from India and Tibet of China), and clade 4.4 (with strains from China 2013–2014).

Fig. 2

Phylogenetic relationships of PPRV genome sequences. Shown is a maximum composition likelihood tree of 25 PPRV genome sequences generated in this study and 15 publicly available PPRV genomes. The tree is unrooted. The scale bar is given in numbers of substitutions per site. Bootstrap resampling (1000 replications) support values are shown at the nodes.

Evolution rate

To analyze the evolutionary rate and the time to the most recent common ancestors (TMRCA) of PPR viruses obtained in China, the 25 genome sequences obtained in this study were used to perform phylogenetic analysis by the Bayesian-based coalescent approach. The mean evolutionary rate of PPRV genome was estimated at 2.61 × 10−6 nucleotide substitutions per site per day (95% highest posterior density [95% HPD], 4.37 × 10−6; 7.74 × 10−7), corresponding to 9.54 × 10−4 nucleotide substitutions/site/year ([95% HPD], 1.59 × 10−3; 2.83 × 10−4). TMRCA estimation showed that the origin of the PPRV China 2013–2014 strains dated back to early August 2013 (95% HPD, July 2013; November 2013) ( Fig. 3). The latest divergence of PPRV from Xinjiang province to other provinces was dated between the end of October and the end of November in 2013. The earliest divergence of PPRV among the other provinces except Xinjiang was dated between the middle of November and the end of December in 2013. It is consistent with the one-month time difference between the outbreaks in Xinjiang province in December 2013 and the first outbreak outside Xinjiang province in the end of January 2014.

Fig. 3

Time-revolved phylogenetic tree based on PPRV genome using Baysian MCMC analysis. The tree was estimated using an uncorrelated exponential relaxed sclock model under an exponential growth model. The scale bar indicates time in day. The posterior probability is shown at each node.

Time-resolved phylogenetic analysis and transmission pathways

Time-resolved phylogenetic analysis was generated from all PPRV genomes from China, 2013–2014. Five clusters of sequences could be identified (Cluster A-E, Fig. 3). However, no distinct temporal and spatial distribution was found for samples in each cluster. As an example, Cluster D included two samples from Xinjiang province in December 2013 (China/XJ2/2013 and China/XJ4/2013), two samples from northeast area between March and April 2014 (ChinaHLJ2014 and ChinaJL2014), and three samples from middle area in April 2014 (ChinaSaX2014, ChinaHeN2014 and ChinaHN2014).

Median Joining phylogenetic analysis implemented by NETWORK also showed that all the 25 PPRV sequences derived from China, 2013–2014 were mainly grouped in five transmission clusters. All clusters were connected at the base of the network through four internal nodes (Mv1, Mv2, Mv3 and Mv4 in Fig. 4). Mv1 was identified as a putative common internal node of all the 25 sequences by comparing these 25 sequences with the most closely related sequence, China/Tibet/30/2007. China/XJYL/2013, which was sampled from the first reported outbreak during the 2013–2014 epidemic, was found to be most closely related to Mv1 by 3 nt substitutions. All sequences in Cluster D shared a same putative ancestor of Mv3, which was derived from Mv1 by 3 nt substitutions. All sequences in Cluster E were originated from a putative ancestor of Mv4, which was derived from Mv3 by 1 nt mutation. It was noteworthy that in each cluster the sequence collected from Xinjiang province was one of the most closely related to the putative ancestor.

Fig. 4

Median-joining phylogenetic network of PPRV during 2013–2014 China epidemic. The median-joining network was constructed from the genome sequences. Each connecting branch line represents a nucleotide substitution. White dots indicate putative ancestors. Each cluster is represented by a grey circle.

We identified two pairs of samples that were unambiguously connected in the phylogenetic network, which indicated the likely inter-farm transmission events. The first likely inter-farm transmission event involved a farm in Hubei province in middle China (ChinaHB2014) and the other farm in Jiangsu province along the east coast of China (ChinaJS2014), with a distance of 1000 km. The second likely inter-farm transmission event linked one farm in Xinjiang Province in western China (China/XJ5/2013) and the other farm in Shanxi province in the middle of China (ChinaSX2014), with a distance of 3000 km.

Population dynamics of PPRV

We estimated the population dynamics of the genome of PPRV using Bayesian skyline plot (BSP) analysis. The plot depicted a consistent pattern of changes in the effective number of infections (N t) through time. A sudden increase of viral population was found during January 2014, followed by a stable platform until June 2014 ( Fig. 5).

Fig. 5

Bayesian skyline plot analysis of the PPRV genomes during China 2013–2014 epidemic. The Y-axis shows the effective population size and the x-axis shows generation time in day. Mean effective population size is indicated by a black line. The 95% HPD of the effective population size are indicated by grey lines.

Discussion

Here we present the first intra-epidemic analysis of PPRV genome by studying viruses from 25 farms infected during the 2013–2014 PPR outbreaks in China. This has demonstrated the potential of such studies to provide in-depth insights into the evolution of PPRV in endemic country to achieve effective control measures.

The 25 PPRV genomes obtained during November 2013 and June 2014 in China showed considerable sequence variety. A total of 96 substitution sites were found along the genome sequences, among which only two sites were found in N gene or F gene partial region. Five clusters were confirmed for the phylogeneitc analysis using genome sequences of these strains. However, only two genotypes were found for these strains in N gene or F gene partial region. It is suggested that although N gene or F gene partial sequence was widely used for PPRV phylogenetic analysis, whole-genome sequences were more useful for detailed transmission study in a short period of time.

Sequence comparison showed that the distribution of nucleotide substitution along the PPRV genome during this endemic was random. The genetic diversity of China 2013–2014 PPR viruses was mainly due to synonymous mutation. It is noteworthy that the highest ratio of NS substitutes was found in P gene. Previous studies revealed that P protein of Rinderpest virus played an important role in transcription and replication of the genome (Kaushik and Shaila, 2004, Saikia et al., 2008). The bottleneck effects occurring during virus transmission may result in the fixation of mutation in P gene. However, the limitation of data scale (only 25 sequences sampled during a short period of 7 months) deduced the failure of selection pressure analysis in this study (data not shown). Further research should be conducted to cover longer time scale to reveal the selection pressure on the evolution of PPRV in China.

Our estimation of the evolutionary rate of PPRV genome in China during 2013–2014 was 2.61 × 10−6 nucleotide substitutions per site per day (9.54 × 10−4 nucleotide substitutions/site/year), which is consistent with previous prediction for the evolutionary rate of other paramyxoviruses (10−3–10−4 nucleotide substitutions per site per year) (Chong et al., 2010, Furuse et al., 2010). Muniraju et al. estimated that the evolutionary rate of PPRV genome was 9.09 × 10−4 nucleotide substitutions/site/year according to the analysis based on 12 PPRV genomes (Muniraju et al., 2014c). The TMRCA of PPR viruses collected in China during 2013–2014 is estimated to be in early August 2013, which is three months earlier than the confirmation of the first outbreak in Yili region of Xinjiang province. It is possible that the virus was transmitted into Xinjiang, China and kept undetected for several months. The presence of a single putative common ancestor around August 2013 for all PPRV China 2013–2014 genome sequences provided clear evidence for a single introduction of the virus. The same conclusion was drawn by previous phylogenetic analysis on N gene or F gene partial sequence (Wu et al., 2015). Phylogenetic analysis based on partial sequence of N gene revealed that China/XJYL/2013 was most closely related to PPR viruses circulating in middle-Asian countries (Kock et al., 2015, Wu et al., 2015). Since no genome sequence of PPRV Central-Asian strain was available, it is not possible to investigate the phylogenomic relationship between PPRV China 2013–2014 strains and central-Asia strains in this study. More genome sequences of PPRV strains circulating in neighboring countries especially in central Asia will help locate the origin of PPRV China 2013–2014 epidemic.

The virus caused two waves of outbreaks after introduction. The first wave of outbreaks was limited in Xinjing province between November and December 2013, involving 5 sporadic outbreaks. The second wave of outbreaks started from the end of January 2013 to the beginning of June 2013. During this wave of outbreaks, more than 200 sporadic outbreaks were reported in 21 provinces. It is of question how the virus spread out of Xinjiang province and then across most part of China. BSP analysis conducted in this study revealed that a sudden increase of viral population was found during January 2014, which was consistent with the occurrence of the second wave of outbreaks. The timing of this wave coincides with the Spring Festival (31st January 2014) which is associated with rising of lamb-consuming and therefore increased risk of long-distance animal transmission. Field investigation revealed that the infected farms in Xinjiang province had intended to fatten the purchased goats and sell to other provinces. In this study, we constructed the transmission pathway of the spread of virus basing on the full-length genome sequences. It revealed that viruses collected from the second wave of outbreaks shared three putative ancestor viruses with those from the first wave of outbreaks. It indicated that the second wave of PPR in China was caused by multi-origin transmission from Xinjiang Province.

Recently, much has been conducted to track the high-resolution transmission history of viruses by using full-genome sequencing. Transmission network of the Middle East respiratory syndrome coronavirus in Saudi Arabia was reconstructed using full genome sequences (Assiri et al., 2013, Cotten et al., 2013). Emergence and evolution of Norovirus in Sydney was investigated by using genetic sequences (Eden et al., 2014). Evolution of Rift Valley fever virus during an outbreak in Kenya during 2006–2007 was described using genome sequences (Bird et al., 2008). In this study, we constructed the transmission pathway of the spread of PPR virus basing on the full-length genome sequences and identified five clusters of infections. Farms in each clusters showed dispersed geographical location, with farm-to-farm distance ranging from 360 km to 3400 km. It indicated that such transmission between farms was caused by long-distance animal movement. Putative-ancestor-virus-centered radial shape of each cluster suggested that samples collected from different province might come from a same virus origin. Previous epidemiologic analysis suggested that infected animals in different provinces could be traced back to some national livestock markets (Fan et al., 2015). These animals were found bought from the market and transferred thousands of kilometers away to the farms by trucks. The incubation period of PPR is typically 4–6 days, which may range between 3 and 14 days. According to the convenient highway network in China, it is suggested that the infected animals which were during incubation period were transferred from Xinjiang province to some national livestock market and then mixed with sheep and goats in the market. It is of most possibility that the suspected animals in the market were infected with PPR and were transferred to the farms in other provinces by long-distance translocation during incubation period.

We identified two pairs of samples that were unambiguously connected in the phylogenetic network, which indicated the likely inter-farm transmission events. The virus from ChinaJS2014 was one nucleotide different from the virus from ChinaHB2014. It indicated a single farm-to-farm transmission despite the long distance between these two infected farms. The network analysis revealed a number of long branches (more than 10 nt substitutions) between the putative ancestor virus and viruses collected from infected farms, which may indicate the presence of intermediate undetected infected farms. Although more than 200 infected farms were reported during 2013–2014 PPRV epidemic in China, only 25 of them were full-genome sequenced and analyzed in this study. Another possibility is that long branches were a result of intra-farm evolution of the virus within an infected farm because we analyzed only one sample from each farm. In further research, analysis with samples collected from more infected farms and more samples from one infected farm will provide deeper information on transmission pathway.

The possibility of further spread of PPR to neighboring countries is of great concern. Among the 14 neighboring countries of China, PPR was endemic in five countries (India, Bhutan, Nepal, Pakistan and Afghanistan) in South Asia and has been reported in two countries (Kazakhstan and Tajikistan) in central-Asia (Banyard et al., 2010, Lundervold et al., 2004). Before 2013, no occurrence of PPR was reported in East Asia and Southeast Asia. In September 2016, two outbreaks of PPR were reported to the Office Internationale des Epizooties in western part of Mongolia (http://www.oie.int). Although no evidence was available to show that the disease in Mongolia was transmitted from China, the widespread of PPR in China has put high risk of this disease to the neighboring PPR-free countries.

This study showed that full genome sequence analysis can make an important contribution to understanding the intra-epidemic transmission of PPRV. Analysis of these data provided useful information for effectively control and finally eradication of the disease.

Funding information

This work was financially supported by The National Key Research and Development Program of China (2017YFD0502300), the Animal Diseases Surveillance Fund of the Ministry of Agriculture, China (08-52) and a scholarship from the China Scholarship Council to J.B (20143012).