Molecular characteristics and pathogenicity of porcine epidemic diarrhea virus in some areas of China from 2015 to 2018

Since The S gene of PEDV in collected samples were amplied by RT-PCR. A phylogenetic tree was constructed using MEGA6.0 software with the neighbor-joining method to analyze the evolutionary relationship. Nucleotide and deduced amino acid (AA) sequences of S were aligned using the MegAlign program of DNASTAR7.1 software to determine sequence homology. PEDV GDgh16 strain isolation, IFA identication and titer detection were performed in Vero cells. Six 4-day-old healthy colostrum deprived suckling piglets were used for challenging experiment of PEDV GDgh16 strain. Virus copies from the small intestine were detected by RT-qPCR. The other section was stained with the anti-N protein McAb at 1:1000 dilution for immunohistochemical (IHC) examinations. The G2-b strains became dominant in 2018. Compared with previous strains, these strains had multiple variations in the SP and S1-NTD domain and in the neutralizing epitope of the S protein. Furthermore, we successfully isolated and identied a new virulent G2-b strain, GDgh16, which was well adapted to Vero cells and had a high mortality rate in piglets. Our study provides full insight into the genetic characteristics of prevalent PEDV strains in parts of China, which suggests that the development of novel effective vaccines is necessary and pressing.

that the disease was caused by a highly pathogenic PEDV variant. Compared with the PEDV variant strains, the S genes of classical CV777 and new OH851 had the same insertions and deletions (S-INDEL strains) [11,12].
The genome of PEDV is approximately 28 kb in length and consists of seven open reading frames (ORFs), which encode four structural proteins and three non-structural proteins [13]. S is the largest structural protein, which contains neutralizing antibody epitopes and a speci c receptor binding site for virus entry [14]. At present, four antigenic epitopes has been characterized in the S protein, including the CO equivalent (COE) domain (aa positions 499-638), the epitopes SS2 (aa positions 748-755) and SS6 (aa positions 764-771), and the epitope 2C10 (1368 GPRLQPY 1374) [15,16]. Because of the vital role of the S protein and extensive mutation of the S gene, it is often used as a target gene for the analysis of virus genetic variation. Based on the S gene,PEDV strains can be classi ed into genogroup 2 (G2) and genogroup 1 (G1). At present, most isolated recovered in China belong to G2 [17]. Recently, a new mutation in the S gene of PEDV has been reported [18]. Different recombinant PEDV strains have also been reported in different areas of China [19,20]. Studies have shown that one province in particular has the co-existence of different genotyped PEDV strains. These results indicate that PEDV has been continuing to spread widely to most areas of China and has caused serious economic loss to the pig industry, thereby manifesting the complex evolution of the virus. Therefore, extensive research of the evolutionary pathogenic mechanism is essential in China.
To control the PEDV spread, the classical CV777-derived vaccine has been widely used in many areas of China; however, it cannot provide adequate protection against PEDV invasion. In contrast, the wide-scale use of the vaccines has increased the environment stress and led to PEDV variation to escape immune protest. To further and fully understand the prevalence and evolution of PEDV in South China, in this study, the diarrhea samples of piglets were collected, and the variation of the S genes of PEDV positive samples were analyzed by sequence alignment and a phylogenetic tree.

Sample collection
A total of 362 diarrheic samples from small intestine tissues or diarrhea feces of suckling piglets in pig farms in ve provinces of China (Guangdong, Guangxi, Jiangxi, Hu'nan and Hainan) were collected from June 2015 to October 2018. The piglets had severe watery-diarrhea and dehydration. The diarrhea feces were re-suspended in 1 mL phosphate buffer saline (PBS) solution in 1.5 mL Eppendorf tubes. After centrifugation at 10,000 × g for 5 min, 200 µL supernatants were transferred into new tubes for RNA extraction and virus isolation.

RNA extraction and sequencing
The total RNA of collected supernatants was extracted using TRIzol reagent (TaKaRa) according to the manufacturer's instructions. To extract RNA, reverse transcription PCR (RT-PCR) was performed using three pairs of newly designed primers for PEDV S gene ampli cation and detection ( Table 1). The overlapping three PCR products were identi ed by 1.5% agarose gel electrophoresis. The positive PCR products were sequenced by Sangon Biological Engineering Co. Ltd., and the entire sequence of the S gene was obtained by using DNAStar 7.1 software. The complete S gene sequence was submitted to GenBank, and the accession no. is listed in Table 2.  S gene sequence analyses The representative strains of complete genome sequences that were available in GenBank were collected and used for phylogenetic analyses ( Table 3). The phylogenetic tree of all the S genes of the representative strains and isolates was constructed by the neighbor joining method with 1000 bootstrap replicates using Molecular Evolutionary Genetics Analysis (MEGA) software (version 6.0) (http://www.megasoftware.net/). Titer detection for the virus growth curve Vero cells cultured in 24-well cell culture plates were infected with PEDV at a MOI 0.01. The cells and supernatants were collected at 12, 24, 36, 48, 60, 72 and 96 hours postinfection (hpi). Then, the cells were frozen and thawed three times. After centrifugation at 10,000 × g for 5 min at 4 ℃, the supernatants were collected and TCID 50 was determined using a microtitration infection assay.

Piglet challenge experiment
To determine the virulence of isolated strain GDgh16, six 4-day-old healthy colostrum deprived suckling piglets were obtained and were arti cially fed with bovine milk from birth. All the piglets without colostrum were randomly divided into two groups with three piglets in each group. One group was challenged orally with 0.5 mL

Statistical analysis
The numerical data were expressed as the mean ± SD, and the data were analyzed using GraphPad Prism software (version 5.02 for Windows; GraphPad Software Inc.). Differences between groups were assessed using ANOVA. Differences were considered statistically signi cant at a P value of < 0.05 and were extremely signi cant at a value of P < 0.01 or P < 0.001. Sixty-two S genes of the tested strains and the representative strains that were published in GenBank were analyzed by a phylogenetic tree. As shown in Fig. 1, the phylogenetic analysis showed that these strains could be divided into two groups, namely, G1 and G2. G1 included the classical strains (CV777 and SM98) and some isolates from China, USA and Japan after 2010. Thus, G1 was further divided into three subgroups: G1-a, G1-b and G1-c. G1-a and G1-b were classical S-INDEL strains, and G1-c was a new S-INDEL strain. G2 was a non-S-INDEL strain and was also divided into two subgroups, G2-a and G2-b, which consisted of a number of severely virulent strains from all over the world since 2010. The strains in our study belonged to G1-b, G1-c, G2-a and G2b. One strain, GDjm18-2, was categorized as subtype G1-b, which included classical vaccine strains CV777/attenuated and JS2008. One strain GDjm17-1 belonged to the G1-c cluster. The other strains that were identi ed in our study formed eight clusters. In these strains, 25 isolates from Guangdong, 3 isolates from Fujian and 1 isolate from Jiangxi formed three clusters and belonged to G2-b, having high similarity with GD-A and CH-GXNN-2012. The other 34 isolates formed ve clusters and belonged to G2-a. In the 34 strains, JXyc15 has a close relationship with the C4 cluster (North American strains). The other strains had a closer identity with CH-ZMDZ-11, CH-HNAY-2015 and CH-HNCDE-2016L. As shown in Table 4, all the isolated strains from 2015 belonged to G2-a (100%). In 2016 and 2017, there were 46.15% and 43.75% isolated strains belonging to G2-a, respectively. Compared with G2-a, the rate slightly increased, and there were 53.84% and 50% isolated strains belonging to G2-b, respectively. However, in 2018, there were 72.22% isolated strains that belonged to G2-b, which was much higher than that of G2-a in 2017 (22.22%).  To further analyze the AA mutations in different domains of the S protein in the isolates, the different domains of the S protein were aligned with CV777, and the average number of AA mutations each year were computed.
SS6, HR1, HR2 and TM in the S2 protein did not obviously change from 2015 to 2018.

Discussion
At present, PEDV has become a vital viral causing diarrhea and has caused much damage to pig farms worldwide. Because there is no effective vaccine against the emerging prevalent strain in China, the variant PEDV strain has been prevalent in many farms of different areas [21]. Considering the viral variant and limited protection of commercial vaccines, it is necessary to fully understand the genetic variation and epidemiology of PEDV for next-generation vaccine development.
In our study, the genetic variation of PEDV in parts of China was analyzed from 2015 to 2018.
The S gene encodes the largest structural proteins and could stimulate the body to produce neutralizing antibodies. Because of its extensive variant, the S gene has been commonly used as a target gene in studies on the genome characteristics of PEDV [22]. Phylogenetic analysis showed that the strains of four subgroups existed from 2015 to 2018, and G2-a and G2-b are the two most prevalent subgroups in China. From 2015 to 2018, eight strains that belonged to four subgroups (G1-b, G1-c, G1-a and G1-b) were epidemic in Jiangmen of Guangdong, which suggested that PEDV had mutated widely and the epidemic of PEDV was becoming more complex. These results agree with Wen et al.'s report [23]. In 2015, all the isolated strains belonged to G2-a, but in 2018, there were 72.22% strains belonging to G2-b, and only 22.22% strains belonging to G2-a. Interestingly, different from G2-a, which included the strains of other countries, such as America, South Korea, and Japan, the G2-b subgroup only contained Chinese-isolated strains. Combined with previous studies, these results suggest that G2-b strains might be the dominant strains in the future in China [24,25].
The S protein is much more variable, and many studies showed that the AA changes in the S protein might affect the virulence and pathogenicity. Our study showed that the number of AA mutations in the SP1 and S1-NTD domains increased in 2017 and 2018. It had been reported that S1-NTD might be a vital domain related to viral virulence [26,27]. The conformation change of S1-NTD might be related to the high pathogenicity of the PEDV strain FJzz1 [24,28]. Recently, increasingly more virulent PEDV strains emerged [24,28,29]. Whether these mutations changed the major conformations and altered the pathogenicity of these strains will be further explored in the future. Our data showed that the positive rate increased from 2015 to 2016, but it decreased from 2016 to 2018, which might be due to the improvement of disease prevention and control strategies. Many pig farms used the mode of "feed-back" to give sows immunity for protecting piglets against PEDV invasion. This was an effective measure to prevent PED, but a risk of virus dispersal also existed, which was why many recombinant PEDV strains were reported [30][31][32][33][34].
At present, four neutralizing epitopes of PEDV S protein have been determined, which were the COE domain (499-638), epitope SS2 (748-755), epitope SS6 (764-771) and epitope 2C10 (1368-1374) [15,16]. In our study, there were AA changes in as many as 35 positions in the COE domain. Moreover, one strain, GDhz16, had four continuous AA mutations in epitope SS6. Epitopes SS2 and 2C10 also had AA substitutions. Because of the mutation, especially some insertions and deletions in the S protein, antigenicity, pathogenicity and neutralization properties of isolated strains have changed [35,36]. That was why the prototype strain CV777-derived vaccine could protect against the disease induced by classical strains but not prevent the disease induced by variant strains [37,38]. Whether these AA changes affect the antigenicity and neutralization properties of the four neutralizing epitopes needs be explored in the future.
Based on previous epidemiological and clinical observations of eld strains since 2010, the emerging G2 strains were highly pathogenic [39]. To investigate the pathogenicity of isolated variant strains, three piglets were infected orally with GDgh16. The results showed that the piglets in the infected group began to have clinical signs of diarrhea at 12 h, and the piglets developed a typical symptom of PED at 16 h. The morbidity reached 100%. The piglets began to die at 24 hpi, and all had died by 48 hpi. Moreover, the small intestine had high viral copies and many viral antigens, which indicated that GDgh16 was a highly pathogenic strain. Other researchers demonstrated that different types of pigs infected with variant PEDV strains shared consistent outcomes [40][41][42][43]. These results indicated that the variant strains were a large threat to the pig industry, and how to control the PED spread has become a critical issue.

Conclusion
The  Amino acid sequence analysis of the neutralizing epitope in the S protein The amino acid sequence alignment of S protein neutralizing epitopes of isolated strains and reference strains using the Clustal W method.