HTML
-
A total of 587 samples were collected from the lungs, lymph nodes, kidney, spleen, and brain of pigs suspected to be infected with PRV from intensive pig farms in eastern China including the Anhui, Fujian, Shandong, and Jiangsu provinces (Fig. 1) between 2017 and 2019 (Table 1). Sample homogenates were prepared by freeze–thaw cycles followed by centrifugation at 5000 ×g for 5 min. The supernatants were then harvested for DNA extraction. Positive samples confirmed by polymerase chain reaction (PCR) were used for virus isolation. For virus isolation, tissues were grinded, filtered using a 0.22-μm membrane, and inoculated into Vero cells. The cells were then incubated until obvious cytopathic effect (CPE) developed. Virus was further plaque purified.
Figure 1. Map indicating PRV-Central China positive farms in China. Different colors indicate the different provinces in China. The PRV reference sequences used in this study were from these colored provinces. The provinces marked with yellow stars are the provinces with positive samples in this study. Provinces with shadow coverage indicate significant association between virus and geographic location.
Strain name Origin Year Month No. GenBank gC gB gD gE ANHUI-1 China: Anhui 2018 Janurary MK922082 MK922098 MK610394 MK610413 ANHUI-2 China: Anhui 2018 Janurary MK922083 MK922099 MK610395 MK610414 ANHUI-3 China: Anhui 2018 March MK922084 MK922100 MK610396 – ANHUI-4 China: Anhui 2018 April MK922085 MK922101 MK610397 – ANHUI-5 China: Anhui 2018 April MK922086 MK922102 – – Anhui-CZ33 China: Anhui 2017 Janurary MK934521 MK922119 – MK934522 Anhui-BB11 China: Anhui 2017 Janurary MK922079 MK922120 – – Anhui-ZJ1 China: Anhui 2018 October – MK922121 – – FJ-1620 China: Fujian 2018 April – MK922103 MK610398 MK610415 FJ-2125 China: Fujian 2018 April – MK922104 MK610399 MK610416 FJ-5 China: Fujian 2018 May MK922087 MK922105 MK610400 MK610417 FJ-Z1 China: Fujian 2018 May MK922088 MK922106 MK610401 MK610418 FJ-Y21 China: Fujian 2018 May – MK922107 MK610402 MK610419 FJ-N1 China: Fujian 2018 May MK922089 MK922108 MK610403 MK610420 FJ-N2 China: Fujian 2018 May MK922081 MK922109 MK610404 MK610421 FJ-N3 China: Fujian 2018 May MK922090 MK922110 MK610405 MK610422 FJ-N4 China: Fujian 2018 May MK922091 MK922111 MK610406 MK610423 FJ-W2 China: Fujian 2018 May – MK922112 MK610407 MK610424 FJ-ZXF China: Fujian 2018 July MK922080 MK922113 MK610408 MK610425 FJ-YXJSJ1 China: Fujian 2019 Janurary MK922092 MK922116 MK610409 MK610428 FJ-YXJSJ2 China: Fujian 2019 Janurary MK922093 MK922117 MK610410 MK610429 FJ-YY China: Fujian 2019 Janurary MK922094 MK922118 MK610411 MK610430 FJ-QYQ3 China: Fujian 2019 Janurary MK922096 – MK610412 – FJ-QYQ4 China: Fujian 2019 Janurary MK922097 – – – FJ-QYQ2 China: Fujian 2019 Janurary MK922095 – – – FJ-SHS1 China: Fujian 2019 Janurary – MK922114 – MK610426 FJ-GSG5 China: Fujian 2019 Janurary – MK922115 – MK610427 Table 1. Sequence information of strains isolated in this study.
-
Virus DNA was extracted using the Virus Genomics DNA Isolation Kit (Tianlong Biotech, Suzhou, China) following the manufacturer's instructions. PRV was detected using a pair of specific PR-D-F/R primers (Supplementary Table S1) (Yue et al. 2009). PCR was performed in a 20-μL volume mixture comprising 10 μL of 2 × Taq Master mix (Vazyme Biotech, Shanghai, China), 7 μL of double distilled water (ddH2O), 1 μL of template DNA, 1 μL each of 10 pmol/μL forward and reverse primers. Thermocycler conditions used for PCR were 95 ℃ for 5 min, followed by 35 cycles of denaturation at 95 ℃ for 30 s, annealing at 59 ℃ for 30 s, and extension at 72 ℃ for 30 s, with a final extension at 72 ℃ for 10 min before storage at 4 ℃. The gB, gC, gD, and gE genes were also amplified (Supplementary Table S1) using the Phanta Max Super-Fidelity DNA polymerase (Vazyme Biotech, Shanghai, China) (Supplementary Table S1). Positive samples were sent to Tsingke (Nanjing, China) for DNA sequencing. Sequences were assembled using the BioEdit software (Hall 1999).
-
All PRV gB, gC, gD, and gE coding sequences from infected swine were collected from NCBI (https://www.ncbi.nlm.nih.gov/) (Supplementary Tables S2–S5), aligned using the MAFFT software (version 7.312) (Kazutaka et al. 2005) and manually adjusted using MEGA (version 7) (Kumar et al. 2016). The IQ tree software (version 1.6.5) was used to detect the best fit nucleotide substitution model according to the Bayesian information criterion (BIC) score (Lam-Tung et al. 2015). Maximum likelihood (ML) trees based on gB, gC, gD, and gE coding regions were reconstructed using RAxML (version 8.4.10) using the general time reversible (GTR) plus GAMMA distribution substitution model and 1000 bootstraps (Stamatakis 2014).
-
To analyze the correlation between each PRV sequence and geographical location after Markov Chain Monte Carlo (MCMC) analysis (Drummond and Rambaut 2007; He et al. 2019a), PRV sequences were first classified according to their country of isolation, while isolates from China were classified as eastern, northern, central, southern, and western China. The correlation was determined using the Bayesian Tip-Significance testing software (BaTS) (Parker et al. 2008). The parsimony score (PS) and association index (AI) statistics were calculated based on the gC gene MCMC analysis. When the P values of AI and PS were < 0.05, the correlation between PRV and geographical distribution was considered significant.
-
The gB, gC, gD, and gE ML trees were uploaded to Datamonkey (http://www.datamonkey.org/) to estimate sites and branched under selection (Delport et al. 2010; He et al. 2019b). The algorithms single-likelihood ancestor counting (SLAC), fast unconstrained Bayesian approximation (FUBAR), fixed effects likelihood (FEL), and mixed effects model of evolution (MEME) were used to identify selected sites. Positive selected branches were detected using an adaptive branch-site REL test for episodic diversification (aBSREL) algorithm (Kosakovsky Pond and Frost 2005; Murrell et al. 2012, 2013; Smith et al. 2015). A site was considered to be under positive selection if at least two algorithms were satisfied (P < 0.1 in SLAC, P < 0.05 in FEL and MEME, P > 0.9 in FUBAR).
Samples and Virus Isolation
DNA Extraction, PCR, and Sequencing
Sequence Alignment and Phylogenetic Analysis
Geographical Correlation
Selection Analysis
-
Samples are collected from four provinces in eastern China (Anhui, Fujian, Shandong, and Jiangsu), but PRV positive sample only detected in Anhui and Fujian Province (Table 1). Of 587 samples, 48 were positive for PRV (positive rate of 8.18%). Sequencing analysis of the main PRV glycoproteins gB, gC, gD, and gE revealed maximum nucleotide sequence divergences of 1.8%, 0.3%, 0.3%, and 0.6% within the isolates, and 2.6%, 5.6%, 2.0%, and 3.7% compared with the reference isolates from GenBank, respectively. Maximum amino acid sequence divergences were 4.0%, 0.7%, 0.3%, and 1.1% within the isolates for gB, gC, gD, and gE, and 4%, 9.5%, 3.5%, and 6.6% from other GenBank isolates, respectively. In addition, we observed some special amino acid site changes in the gC and gE proteins. For gC, more than half of the strains had a S to G change at position 41, and R162H and S345L changes were observed in four of the twenty strains. Compared with the Bartha vaccine strain, in addition to some specific amino acid mutations, there were some specific insertions from residues 59–65 on gC. A T386M change was observed in 8 gE strains, and 4 strains displayed a L575P change (Fig. 2). However, in case of gB, two strains (FJ-W2 and FJ-ZXF) were similar to the clade 1 strain NIA3.
Figure 2. Sequence alignment and amino acid differences between strains sequenced here and GenBank isolates. The strains sequenced here are indicated in red. Variant PRVs are indicated in orange, early clade 2 strains in blue, and clade 1 strains in green. The yellow bar indicates different site among these sequences.
-
We reconstructed phylogenetic trees of PRV based on the gB, gD, gE, and gC genes using all PRV sequences available from infected pigs in GenBank and the sequences generated here (Figs. 3, 4). Strains sequenced in this study were mainly similar with strains distributed in Guangdong, Guangxi, Henan, Hainan, Shandong, Shanxi, Fujian and Hubei Province. Based on the ML trees, we concluded that four trees had similar structures and PRV could be divided into two main clades, clade 1 and clade 2. Most of the sequences generated here belonged to clade 2, similarly to variant PRV (He et al. 2019a). This indicates that variant PRV is epidemic in eastern China. However, we also found that the sequences of different genes clustered with different clades. For example, gB of strains FJ-W2 and FJZXF located in clade 1 (containing the Bartha vaccine strain) and was genetically closer to strains from other countries, while gD, gE (Fig. 3B, 3C) and gC (Fig. 4) clustered in clade 2. This is indicative of inter-clade recombination events occurring in the PRV epidemic. In addition, sequences of different strains within clade 2 were grouped in different sub-clades for different genes, suggesting intra-clade recombination during PRV epidemic in eastern China.
Figure 3. Maximum likelihood trees reconstructed based on the gB, gD, and gE genes. Trees were reconstructed using RAxML (Version 8.4.10) with the general time reversible (GTR) plus GAMMA distribution substitution model, and 1000 bootstraps. The first layer of red dots indicates sequenced strains from this study. The second layer of colored circles indicates PRV host, and the third circle indicates country and regions of PRV. A gB gene. B gD gene. C gE gene. Red dots indicates PRV sequenced in this study.
Figure 4. Maximum likelihood trees reconstructed based on the gC gene. The tree was reconstructed using RAxML (Version 8.4.10) with the general time reversible (GTR) plus GAMMA distribution substitution model, and 1000 bootstraps. Colored lines indicate host and country of PRV. Red dots indicates PRV sequenced in this study.
-
Given that gC has the highest mutation rate and the largest number of reported sequences (Ye et al. 2015) among the glycoprotein-encoding genes, we studied the geographical correlation of PRV based on gC. We performed BaTS analysis and found that the P value of both AI and PS were < 0.01 (Table 2). Except for the P value of Malaysia, Belgium, Greece, United Kingdom, central China, and western China that equaled 1, the other 14 countries and regions showed a P value within 0.01. This is indicative of a strong correlation of PRV with geographical location, especially in some areas of China where the epidemic is more severe.
Statistic Observed mean Lower 95% CI Upper 95% CI Nulsl mean Lower 95% CI Upper 95% CI Significance AI 13.08 11.37 14.72 28.24 27.39 28.93 0.00 PS 102.28 98.00 107.00 201.63 198.27 204.85 0.00 MC Eastern China 11.06 11.00 12.00 1.74 1.48 2.09 0.01 MC USA 2.03 2.00 2.00 1.06 1.00 1.15 0.01 MC Brazil 3.52 2.00 6.00 1.37 1.14 1.94 0.01 MC Hungary 2.14 2.00 3.00 1.04 1.00 1.14 0.01 MC China 6.65 3.00 10.00 2.44 2.21 3.02 0.01 MC Malaysia 1 1.00 1.00 1.00 1.00 1.00 1.00 MC Spain 4.26 2.00 7.00 1.06 1.00 1.17 0.01 MC Germany 3.06 2.00 5.00 1.17 1.03 1.40 0.01 MC Slovakia 2.07 2.00 3.00 1.02 1.00 1.06 0.02 MC Italy 7.55 4.00 13.00 2.13 1.85 2.71 0.01 MC Argentina 2.9 1.00 5.00 1.09 1.00 1.31 0.01 MC Austria 1.8 1.00 2.00 1.01 1.00 1.03 0.01 MC Croatia 5.47 3.00 7.00 1.05 1.00 1.13 0.01 MC Belgium 1.47 1.00 3.00 1.09 1.00 1.29 1.00 MC Central China 1.45 1.00 2.00 1.20 1.06 1.39 1.00 MC North China 1.52 1.00 2.00 1.00 1.00 1.01 0.01 MC Western China 1.01 1.00 1.00 1.00 1.00 1.01 1.00 MC Greece 1.37 1.00 2.00 1.00 1.00 1.01 1.00 MC South China 3.03 3.00 4.00 1.56 1.32 2.02 0.01 MC United Kingdom 1 1.00 1.00 1.00 1.00 1.00 1.00 Bold number remind the difference in geographic associations significantly.
AI association index, PS parsimony score, MC monophyletic clade.Table 2. Geographical correlation analysis of PRV based on gC.
-
Using the FUBAR and MEME methods, we detected positive selection at gB residues 43, 834 and 908, gE residue 348 according to all methods, and two gC residues (75 detected with FEL, SLAC and FUBAR; and 194 detected with SLAC and FUBAR). No positive selection was detected on gD (Supplementary Table S6). As expected, some selected residues have been associated with relevant functions. For example, gB site 43 is near the B-cell epitopes within residues 59 to 126 (Zaripov et al. 1999), site 75 on gC is associated with virus adsorption (Karger and Mettenleiter 1993), and site 194 is on the 1/3 N terminal of the gC protein, which is associated with the HS receptor-binding domain (Flynn and Ryan 1995, 1996). In addition, the amino acid sites in which we found changes in relation to the reference strains had also important functions. In particular, change S41G on gC results in a change in the glycosylation site from NSS to NGS.