Emergence of Different Recombinant Porcine Reproductive and Respiratory Syndrome Viruses, China

Epidemiological investigations were conducted on recently emerging porcine reproductive and respiratory syndrome virus (PRRSV) strains in Shandong province in 2014–2015. The proportion of the NADC30 strain identified by ORF7 sequence alignment has been gradually increasing. Three emerging PRRSV strains were successfully isolated, and the complete genomic sequences were determined. Our results indicate the importance of recombinant strains in Shandong province, China. There was a varied degree of recombination of two or three strains (classical, HP-PRRSV and/or NADC30). Moreover, the recombination strains affected the pathogenicity of newly emerged strains.


Results
Epidemiology of PRRSV and variability of ORF7. PRRSVs of recent outbreaks (since 2014) are characterized by severe sow abortion and respiratory symptoms of conservation piglets. PRRSV strains were present in 13 different regions of Shandong province. Within a total of 347 clinical samples of diseased pigs (182 lung and lymph node samples and 165 serum samples), 243 samples were PRRSV-positive in RT-PCR using primers that targeted the PRRSV ORF7 (N) gene. The percentage of PRRSV positive samples was 40%-50% from June 2014 to April 2017 ( Supplementary Fig. 1). Comparative analyses of ORF7 sequences showed that ORF7 primers were able to distinguish classical PRRSV (represented by the virus VR2332), HP-PRRSV (represented by the virus JXA1) and NADC30. However, the proportion of NADC30 strains in the PRRSVs of recent outbreaks (since 2014) was 50.62%, which was higher than that of classical PRRSV (6.67%) or HP-PRRSV (42.80%) ( Table 1).
To better understand the genetic relationship and evolution of PRRSV in Shandong province, China, an ORF7-based phylogenetic tree of 243 positive PRRSVs and 54 reference viruses was constructed. According to the phylogenetic tree ( Fig. 1), these PRRSVs were divided into two distinct genotypes: genotype I and genotype II. All the PRRSVs of Shandong province belonged to genotype II, which was divided into three monophyletic lineages: the first sub-group contained 16 isolates belonging to classical PRRSV (represented by the virus VR2332), the second sub-group contained 104 isolates related to HP-PRRSV (represented by the virus JXA1), and the third sub-group contained 123 isolates (76 isolates from 2017) that were closely related to a group represented by NADC30.
Identification of PRRSV isolation. Three PRRSV strains (SDhz1512, SDlz1601 and SDYG1606) that caused piglets to have clinical respiratory distress and interstitial pneumonia were recovered by PAMs. PAMs inoculated with four PRRSV strains (SDhz1512, SDlz1601, SDYG1606 and SX-1) showed green fluorescence, and the positive control cells infected with SX-1 had died. No fluorescence signal was observed in the control cells (Fig. 2).

Genomic characterization of 3 PRRSV isolates.
The complete genome sequences from three PRRSV isolates were determined and deposited in GenBank (accession no. KX980392, KX980393 and KY053458). They were classified as genotype II isolates but differed in length: 15069 nt for SDlz1601, 15367 nt for SDhz1512, and 15074 nt for SDYG1606. Full-length sequence analyses of SDlz1601, SDhz1512 and SDYG1606 were performed together with reference PRRSV strains including two North American isolates (VR2332 and NADC30)  Table 2). The nucleotide and amino acid homology of the 5′UTR, 3′UTR, ORF1a, ORF1b, and ORF2~7 genes in the three strains and the representative strains, including NADC30, VR2332 and JXA1, were compared. The results showed that the 5′UTR, ORF2~7 and   Amino acid analysis of NSP2, secondary structure and 3D-structure prediction. Non-structural protein 2 (NSP2) had the highest genetic diversity in the PRRSV genome. To further characterize the deletion regions, the NSP2 of the three isolated strains (SDhz1512, SDlz1601 and SDYG1606) was compared with 8 representative genotype II isolates. Multiple sequence alignment revealed that SDlz1601 and SDYG1606 had three discontinuous deletions (a total of 131 amino acids) in NSP2 that resembled previous NADC30-like PRRSVs, which can be used as molecular markers to distinguish them from other PRRSVs 15 . These deletions were not present in the Chinese HP-PRRSV strains (Fig. 3). Furthermore, the SDYG1606 isolate had the same 9 amino acid  deletions (793-801) in NSP2 and was closely related to NADC30 (Fig. 3), but SDhz1512 had 30 amino acids deletion and was more similar to the HP-PRRSV strains (HUN4 and JXA1). Surprisingly, the 9 amino acid deletions (814~822) of NSP2 in SDlz1601 differed from those of previously sequenced NADC30-like isolates from China. NSP2 secondary structures of the three isolates (SDhz1512, SDlz1601 and SDYG1606) and three reference strains (VR2332, JXA1 and NADC30) were predicted using protein software of DNAStar version 7.1, Madison WI ( Supplementary Fig. 2). SDhz1512 was consistent with the highly pathogenic strain JXA1, in which the α-helix was deleted and the β-strand was increased. Compared with the classical strain VR2332, SDlz1601 and SDYG1606 were more consistent with the new NADC30-like strain, and the deletion of a total of 111 amino acids (323~433) resulted in the lack of a large α-helix, β-strand and random coil. In addition, the Antigenic Index in this region (320-360) was reduced to a low level using Protein software (DNAStar, version 7.1, Madison WI).
In conclusion, the absence of amino acid groups in the discontinuous large fragments in NSP2 of NADC30-like strains resulted in significant changes in the secondary structure of the protein.
Recombination analyses in SDhz1512, SDlz1601 and SDYG1606. In this study, a comprehensive and novel approach was employed to characterize recombination in the complete genome of PRRSV. Using seven algorithms, RDP, GENECONV, MaxChi, Chimaera, Bootscan, SiScanand, and 3Seq methods in RDP4.70 software provided strong evidence of recombination events in most of the sequences analyzed. Seventeen recombination events were detected in the three isolates (Table 3). SDYG1606 had 3 potential recombination signals, whereas SDhz1512 had 8 potential recombination signals and SDlz1601 had 6 potential recombination signals (Table 3). Of all 16 potential recombination events, 14 recombination events were detected in SDYG1606. SDlz1601 had a remarkably high degree of certainty, with a P-value in at least four algorithms of <1 × 10 −6 ( Table 3). The other 14 recombination events with a high degree of certainty due to a recombinant score >0.  (Table 3). In addition, the remaining 3 recombination events had a fair likelihood since their recombinant score was between 0.429 and 0.549 (Table 3). Furthermore, there were similar recombinant events with the same locations and parents in different isolates including a recombination located at 1944-2021 with parents NADC30/JXA1 in SDhz1512 (event 7) and SDlz1601 (event 17) (Table 3).
To further validate and confirm the putative recombination events and their breakpoint positions in isolates detected with RDP, the SimPlot method in the SimPlot v3.5.1 program was employed. This method records the coherence of the sequence relationship over the entire length of the isolate and its potential parents 23 . This study revealed that SDhz1512, SDlz1601 and SDYG1606 were the result of recombination between the NADC30-like viruses and the classical strains and HP-PRRSV strains present in outbreaks in China. From the similarity plot, we identified seven recombination breakpoints in SDhz1512: one was located in NSP1α (nucleotides [nt] 352), two points were located in NSP7 (nt6349 and 6938), two points were located in NSP9 (nt8009 and 8830), one point was located in NSP10 (nt10000), and the last point was located in ORF2a (nt12125) (Fig. 5Aa). The analysis revealed that SDlz1601 was the result of recombination between the NADC30 viruses and HP-PRRSV strains circulating in China. There were seven recombination breakpoints detected in the SDlz1601 isolate: two points were located in NSP2 (nt 2082 and 3919), one was located in NSP9 (nt8774), two points were located in ORF4 (nt 13398 and 14050), one was located in ORF5 (nt14344), and one was located in ORF6 (nt15015) (Fig. 5Ba). The analysis revealed that SDYG1606 was the result of recombination between the NADC30 viruses and the classic HP-PRRSV strains circulating in China. From the similarity plot, we identified three recombination breakpoints: one was located in NSP2 (nt1942) and the others were located in NSP4 (nt5759) and NSP9 (nt9296) (Fig. 5Ca). The breakpoints of SDhz1512 separate the genome into eight regions: three are closely related to the HP-PRRSV strains in China (represented by JXA1) (Fig. 5Ab), four are closely related to classical PRRSV strains from China (represented by China1a) (Fig. 5Ac), and remaining region is closely related to the PRRSV strains from North America (represented by NADC30) (Fig. 5Ad). The breakpoints of SDlz1601 separate the genome into eight regions: four are closely related to the classic HP-PRRSV strains in China (represented by JXA1) (Fig. 5Bb) and the other four are closely related to PRRSV strains from North America (represented by NADC30) (Fig. 5Bc). The breakpoints of SDYG1606 separate the genome into four regions: two are closely related to the classic HP-PRRSV strains in China (represented by JXA1) (Fig. 5Cb) and the other two are closely related to the PRRSV strains from North America (represented by NADC30) (Fig. 5Cc).
Interestingly, eight of the above 6 reported and 3 isolated strains had different parents of recombination with other PRRSV strains such as classical strains (VR2332 and/or CH1a), HP-PRRSV (JXA1 and/or 09NEN1), and NADC30 ( Table 4). The recombination breakpoints analysis of the above recombinants demonstrated that they . Tertiary structure models for NSP2 of VR2332, SDlz1601, SDYG1606, SDhz1512, JXA1 and NADC30 using I-TASSER Server. TM-score values range from [0, 1]. A higher score indicates a better structural match. Statistically, a TM-score <0.17 means a randomly selected protein pair with the gapless alignment taken from PDB; a TM-score >0.5 corresponds approximately to two structures of similar topology. C-score values are in the confidence interval [−5, 2], with higher values indicating better models. The TM-score and RMSD are the indexes for evaluating the molecular structure. The correlation of the C-score with RMSD and the TM-score is extra high.
were similar to the phenomenon of fracture recombination, and the recombinant positions mainly occurred in NSP2, NSP9, ORF2, ORF4 and ORF7 (Table 4). These proteins may have been the hot regions for recombination event breakpoints.

Phylogenetic analysis.
To determine the genetic relationship of the HP-PRRSV strains that have recently emerged in China, a phylogenetic tree based on the full genome sequences of SDhz1512, SDlz1601, SDYG1606 and another 73 published PRRSV isolates that were identified between 2013 and 2016 was generated. The phylogenetic analysis of the entire PRRSV genome was performed using RaxML and the BI tree, and its confidence was evaluated with 1000 bootstraps. As shown in Fig. 6, the results showed that genotype II isolates in China could be divided into three sub-genotypes: those closely related to NADC30-like isolates, China classical isolates, and HP-PRRSV isolates. The phylogenetic tree further revealed that SDhz1512, SDlz1601, and SDYG1606 were recombinant strains and had lower identity to other isolates. The complete genome of the three isolates shared only 89-92% identify, which differed from the published NADC30-like strains. This suggests that the virus gained genetic diversity by recombining with local HP-PRRSV, classical strains and North America isolates (represented by NADC30) in China, and these recombinants are becoming more common. The MCC tree revealed, with high resolution, the evolutionary history of NADC30-like viruses and related viruses in China, from which we could infer that NADC30-like PRRSVS were more closely related to NADC30 and were clustered into a separate branch, thus distinguishing them from the HP-PRRSV cluster represented by JXA1 and HuN4. The Bayesian skyline plot revealed two periods in which there was a marked increase in relative genetic diversity (Fig. 7) and whose timing and scale matched well with the epidemiological record of the HP-PRRSV and NADC30 PRRSVS outbreaks. PRRSV recombination strains exhibited pathogenicity for piglets. The pigs challenged by PRRSV SDhz1512, SDlz1601 and SDYG1606 strains developed similar clinical signs, such as fever (≥40 °C) and anorexia, and the clinical signs became more severe over the next few days. Three pigs in the SDlz1601 or SDYG1607 strain challenge group had respiratory distress, characterized by dyspnea, tachypnea and coughing, whereas two pigs in the SDhz1512 strain challenge group displayed similar symptoms. The piglets in the control group remained clinically normal throughout the study. The viral titers of the three isolates in sera reached peak levels at 14 days post-inoculation (dpi) (Supplementary Fig. 3), whereas all the pigs in the mock-inoculated control group were negative for PRRSV throughout the entire experimental period. The viruses were recovered from the sera of pigs in the three infected groups, and ORF7 genes were sequenced to confirm they were the original virus.
The major gross pathology finding was different levels of interstitial pneumonia in the PRRSV-challenged pigs (Fig. 8). Interstitial pneumonia was characterized by thickening of the alveolar septa and infiltration of infiltrating lymphocytes and macrophages. Necrotic bronchial epithelial cells could also be found in the trachea. No pathological lesion was identified in the control pigs. The morbidity was determined by clinical symptoms and pathological changes (

Discussion
Since 2013, the phenomenon of sow abortion has occurred frequently, resulting in significant losses to the global pig industry. The genetic diversity of HP-PRRSV has greatly increased by rapid evolution and recombination events 14,19,24 . Recently, several outbreaks of NADC30-like PRRSVs were reported in different provinces of China 18 . Strikingly, all NADC30-like strains carry the genetic marker of PRRSV MN184 strains, namely, the exact discontinuous deletions in the NSP2-coding region 20 . Some workers claimed that geographic separation is a factor influencing PRRSV evolution, based on ORF5 and genome sequences 25,26 . Yoon et al. 27 have also reported that there was no immediate relationship between the date or place of collection and the topological distribution of ORF7 in the different PRRSVs. Conversely, Hao X et al. 28 claimed that PRRSV evolution was also based on the ORF7 sequence. However, regarding the phylogenetic relationships among the PRRSVs, our results based on the whole genome sequence data of 73 clinical samples obtained between 2013 and 2016 using RaxML and a BI tree revealed that the Chinese and North American genotypes isolates could be divided into three sub-genotypes using ORF7 primers: those closely related to NADC30 isolates, classical isolates and HP-PRRSV isolates. Moreover, three PRRSV strains (SDlz1601, SDhz1512, and SDYG1606), which had more than 96% homology with NADC30     Recombination is not uncommon in PRRSV, and there are more and more reports about the recombination of PRRSVs 5,19 . Recently, Li Y et al. 16 reported PRRSV HNyc15 isolation characterized by NADC30-like recombination with classical PRRSV strain VR2332 and CH-1a between ORF2 and ORF4. Zhao K et al. 19 also indicated that the NADC30-like PRRSV strain from North America has spread to several provinces in China and recombined with local HP-PRRSV strains. Zhao H et al. 24 reported that nine of 28 isolates and one isolate from another laboratory were potential complicated recombinants between the vaccine JXA1-R strains and predominant circulating strains, and the PRRSV recombination rate is increased under the current vaccination pressure. However, our study showed that SDlz1601 and SDYG1606 exhibited a 111 aa deletion at positions 323 to 433 and a 19 aa deletion at positions 506 to 524, which can be considered to be NADC30-like strains. Conversely, the SDhz1512 strain was a recombination strain of NADC30, classical and HP-PRRSV. A 30 aa sequence in the NSP2 gene was detected that was consistent with HP-PRRSV. Accordingly, not only NADC30-like strains but also recombinant strains of NADC30 with classical PRRSV and/or HP-PRRSV have become the current popular trend.
A recent publication indicated that piglets infected with NADC30-like viruses, such as HNjz15, CHsx1401, and FJ1402, showed fever and respiratory disorders, but no death 29,30 . Obviously, the virus NADC30-like pathogenicity is different from the Chinese HP-PRRSV, which is fatal for piglets 31 . In this study, animal challenge was used to determine the relationship between recombination and the virulence of three isolated PRRSV strains. The SDYG1606 strain led to more severe interstitial pneumonia than SDlz1601 although the sequence characteristics of it were similar to NADC30-like strains. The reason for the difference in virulence may result from varied recombination breakpoints. However, the piglet deaths were caused by secondary infection of bacteria. The new variant strain, SDhz1512, which was a triple recombination strain including NADC30, classical and HP-PRRSV, caused more minor clinical symptoms than the NADC30-like strains. Our present data indicated that the virulence of different PRRSV strains with different recombination breakpoints was varied, but the virulence of the recombinant strain with classical PRRSV was slightly weaker.
Bayesian phylodynamic models showed one remarkable improvement compared to traditional methods and could make use of associated epidemiological information to infer genetic relations. These inferences could be used to identify viral dispersion routes that correspond with transportation patterns involving high PRRSV risk 32 . In this study, 75 NSP2 sequences of genotype II PRRSVs from 2000 to 2016 of China obtained from GenBank were used to estimate the evolutionary history of PRRSV in China. The Bayesian skyline plot revealed two periods in which there was a marked increase in relative genetic diversity. The first growth phase corresponded to the initial HP-PRRSV outbreak in 2005~2010, whereas the second growth phase was mainly associated with the NADC30-like outbreak in 2013~2017. It matched well with the epidemiological record of the HP-PRRSV and NADC30 PRRSV outbreaks.  Currently, there are a variety of commercial PRRSV vaccines in the Chinese market. The current modified live vaccines do not provide complete cross-protection against heterologous PRRSV strains 20 . Recent studies suggest that several newly identified virulent PRRSV isolates have been introduced into swine populations through the inoculation of PRRSV-derived inactivated vaccines 33,34 . Recently, Shi M et al. 13 indicated that recombinants appeared to be highly pathogenic, such that the recombination events that generated them either preserved or increased the pathogenicity of the parental strains. Zhao K et al. 19 also confirmed the high pathogenicity of JL580 recombinants between NADC30-like and HP-PRRSV strains in piglet challenge. Several studies reported the outbreaks of NADC30-like PRRSVs in vaccinated pig herds with 30-50% fatality 15 . Thus, the current prevalence of recombinant strains may be the result of PRRSVs escaping vaccine immunization. However, whether the existing vaccine can provide protection for recombinant strains requires further study.

Material and Methods
Clinical samples. During  were isolated from serum using porcine alveolar macrophages (PAMs) as previously described 35 . The inoculated cells were maintained at 37 °C in a 5% CO 2 atmosphere for 72 h. The PRRSV strains were confirmed with an  immunofluorescence assay using a monoclonal antibody directed against the PRRSV N protein, which was produced from hybridoma cells derived from Sp2/0 myeloma cells and spleen cells of BALB/c mice immunized with the N protein of PRRSV strain SX-1 (accession no. GQ857656.1).
Primers designed to determine the complete genome sequence. To determine the complete genomes of SDhz1512, SDlz1601 and SDYG1606 strains, 12 primers were designed on the basis of the genomic sequence of NADC30 available in GenBank (accession No. JN654459). Two Chinese isolates (JXA1 and SD1100) available from the National Center of Biotechnology Information (NCBI) were used as references. The primers used for full genome sequencing of the three isolates are shown in Table 6.
Sequence assembly and alignment. Consensus 41 . Default settings were used and the threshold p-value set at 0.05 using the Bonferroni correction. To avoid false positive results, recombination events supported by at least six different methods were considered. The putative recombination events and their breakpoint positions were further validated and confirmed via the SimPlot program (version 3.5.1), which performs similarity comparisons within a 500-bp window sliding along the genome (10 bp step size) using the similarity score and Bootscan plot 42 based on 26 representative genome sequences. The parental sequences detected were further checked for the presence of recombination. A cluster analysis maximizing the value of x2 was then used to select breakpoints among the clusters 38 . These breakpoints were used to divide the alignment into segments for phylogenetic tree construction.
Pathogenicity study of the three PRRSV isolated strains. To determine the pathogenicity differences among SDhz1512, SDlz1601 and SDYG1606, an animal experiment was designed. Twenty four-week-old piglets were randomly divided into four groups and maintained in individual rooms. All the piglets were negative for antibodies of porcine circovirus type 2 (PCV2) and PRRSV before the experiment. Each piglet in the PRRSV strain group (SDhz1512, SDlz1601 or SDYG1606) was intranasally administered 2 ml of the virus containing 2 × 10 5 TCID 50 . Each animal in the control group was given the same dosage of PAM culture supernatant. The viremia of PRRSV-inoculated animals was detected at 0, 3, 7, 14, and 21 dpi by an IFA-microtitration infectivity assay. The animals were euthanized and lung samples were collected for histopathology at 14 dpi.
Phylogenetic classification of PRRSV. The ORF7 nucleotide phylogenetic and molecular evolutionary analyses of 243 positive PRRSVs were conducted using MEGA version 7, along with 54 sequences of representative PRRSVs available in GenBank from various countries and areas. Seventy-five complete genomes representing the main PRRSV strains were obtained from GenBank. Multiple sequence alignments were performed in MAFFT version 7.263, and the nucleotide substitution model was evaluated by j-ModelTest v2.1.10 using the Akaike Information Criterion (AIC) 43 . A general-time reversible (GTR) model of nucleotide substitution was used with a gamma-distributed (Γ) rate variation among sites, with a proportion of invariant sites (I). A maximum likelihood (ML) tree was reconstructed using RaxML v8.2.10 44 , and its confidence was evaluated by 1000 bootstraps.
Time dynamic evolution analysis of PRRSV. According to the description of Shi M et al. 13 and Sook Hee Yoon 26 , 75 NSP2 gene sequences of genotype II PRRSVs from 2000 to 2016 of China obtained from GenBank were used to estimate the evolutionary history of PRRSV in China. Alignments with recombinant regions removed were utilized for phylogenetic analyses using the Bayesian Markov chain Monte Carlo (MCMC) approach available in BEAST v1.8.2 45 . A Bayesian skyline non-parametric coalescent model and uncorrelated lognormal (UCLN) relaxed molecular clock were selected. Bayesian MCMC chains were run for 500 million generations, 10% of which were removed as burn-in, and sampled every 50,000 steps. Convergence and uncertainty in parameter estimates were evaluated by calculating the active sample size 46 in Tracer v1.6.
Comments. By submitting a comment you agree to abide by our Terms and Community Guidelines. If you find something abusive or that does not comply with our terms or guidelines, please flag it as inappropriate.