Genomic Evolution and Selective Pressure Analysis of a Novel Porcine Sapovirus in Shanghai, China

1. Introduction

Sapovirus (SaV), a member of the Caliciviridae family alongside Norovirus, Lagovirus, and Vesivirus, stands out as a leading cause of acute gastroenteritis in humans, especially among younger children[1]. SaVs, with a non-enveloped structure, possess a linear genome comprising polyadenylated, single-stranded, positive-sense RNA. The SaV genome features two primary open reading frames (ORFs): ORF1 encodes nonstructural proteins VP1, while ORF2 encodes the minor protein VP2 with an unknown function[2]. Genetically diverse, SaVs fall into 19 genogroups[3]. SaVs have been identified globally[4,5,6] and also in various species such as swine[7], mink[8], cows[9], dogs[10], and bats[11]. Therefore, it is essential to focus attention on this agent, given its potential for cross-species dissemination[12].

Porcine SaV (PoSaV) was first observed in fecal specimens from piglets in the USA using electron microscopy in 1980[13]. The initial gastroenteritis outbreak in Chinese piglets caused by PoSaV was reported in Shanghai in 2008. Despite its global presence, genomic characterization in China remains limited. Moreover, the elusive viral isolation poses a significant obstacle to understanding the mechanism and developing vaccines. A comprehensive investigation into the genetic evolution of PoSaV is pivotal for pathogenesis studies. Considering the potential for zoonosis, the continuous evolution of PoSaV variants poses a major concern for global porcine production. Phylogenetic and molecular evolutionary analyses can uncover codons under negative/purifying selection in a species. Thus, this study aims to genetically characterize prevalent PoSaV strains, determining their genetic divergence and selective pressure. These findings serve as a guide for effective defense and containment of PoSaV infections in the country.

2. Materials and Methods

2.1. Metagenomic Sequencing and Data Analysis

Twenty piglet diarrhea fecal samples were collected from three herds in Shanghai, divided into two groups (PCM1 and PCM2), and subjected to metagenomic sequencing. Genomic DNA/RNA was extracted from the pretreatment samples with MagPure Viral DNA/RNA Mino LQ Kit (Magen kitR6662-02). The viral RNA underwent reverse transcription by SuperScript III reverse transcriptase (Invitrogen). Subsequently, the SMARTer Ultra Low Input RNA kit was used to synthesize double-strand cDNA, and the quality of viral nucleic acids was assessed using the NanoDrop spectrophotometer (Thermo Fisher Scientific) and 1.5% agarose electrophoresis. The Next® UltraTM DNA Library Prep Kit for Illumina® (New England Biolabs) was employed for generating sequencing libraries. High-throughput sequencing occurred on an Illumina Novaseq 6,000, producing 150 bp paired-end reads by the Magigene Company (Guangzhou, China).

The raw sequencing reads were processed to obtain clean data using Soapnuke (v2.0.5) for further analysis. Quality-filtered reads were then de novo assembled to generate the metagenome. Different virus families were classified according to the annotation information from the NCBI taxonomy database.

2.2. Complete Genome Sequencing of PoSaV and Genomic Evolutionary Analysis

To analyze piglet-prevalent PoSaV’s genomic characteristics, the full-length sequence of a GIII genotype isolate was amplified. The PCM1 and PCM2 samples were then mixed, followed by extracting the total RNA using TRIzol reagent (TAKARA, Japan) followed by cDNA synthesis using the downstream primer (5’-CGGTACGCGTAACCAGGGAAAGA-3’). PCR, employing 2 × Phanta Max Master Mix (Vazyme, China), utilized prime pairs designed based on reference sequences KT922089. The amplification conditions involved initial denaturation at 98 ℃ for 2 min, followed by 30 cycles of 98 ℃ for 20 s, 60 ℃ for 20 s, and 72 ℃ for 1 min, with an additional extension step at 72 ℃ for 5 min. Each fragment underwent cloning into pEASY Blunt Zero Vector and sequencing by Sangon Biotech (Shanghai) Co., Ltd. (China). Furthermore, the full-length sequence was subsequently compared and spliced using Clustal and Lasergene 7.0 software.

For phylogenetic analysis, PoSaV ORF2 genes were obtained from GenBank, forming the basis for a neighbor-joining phylogenetic tree. This tree, annotated via the interactive Tree Of Life (iTOL) software (http://itol.embl.de/)[14], aimed at illustrating the evolutionary relationships. To assess PoSaV evolution rates, RdRp sequences from GenBank-derived PoSaVs underwent analysis. This involved constructing both a phylogenetic tree and a time tree, using the maximum likelihood method within MEGA 11 with 1000 bootstrap replicates. A 50% cut-off value was applied to condense the tree. The RelTime with Dated Tips (RTDT) method estimated divergence times for all tree branching points. This analysis referred to 50 PoSaV strains identified between 1998 and 2022.

2.3. Recombination and Selective Pressure Analysis

Recombination analysis utilized Recombination Detection Program (RDP) software version 4, examining the entire genome of PoSaVs[15].

The amino acid homology of PoSaV ORF1 proteins was assessed using Larsergene 7.1 and MEGA 11 for selective pressure analysis. Synonymous nonsynonymous analysis was conducted through SNAP services (http://hcv.lanl.gov/content/sequence/SNAP/SNAP.html) to determine the nonsynonymous to synonymous substitutions ratio (dN/dS) at each ORF1 coding region amino acid site, indicating evidence of positive or negative selection[16].

3. Results

3.1. Sapovirus Is Highly Prevalent in the Piglet Diarrhea Feces

To perform a comprehensive search for piglet gut viruses, twenty diarrhea feces were collected and divided into two groups: PCM1 and PCM2. The sapovirus abundance in piglet diarrhea feces was estimated using metagenomics. Results showed that RNA viruses, excluding phages, constituted 55.47%, with Caliciviridae and Astroviridae accounting for 27.54% and 21.78%, respectively (Figure 1a). Sapovirus emerged as the sole Caliciviridae member in this metagenomic sequencing. Analyzing the sapovirus ratio among total RNA viruses confirmed infection rates of 55.36% and 29.89% in PCM1 and PCM2 samples, respectively (Figure 1a). These findings indicate sapovirus as one of the predominant viruses in piglet diarrhea guts.

3.2. Phylogenetic and Time Evolutionary Analysis

To understand the genomic evolution of sapovirus in piglet diarrhea, a full-length sequence of sapovirus was amplified and spliced, named as SHCM/Mega2023. The phylogenetic tree, based on the ORF2 genes of SHCM/Mega2023 and corresponding GenBank strains, identified its classification as GIII genotype. SHCM/Mega2023 demonstrated the closest relationship with the 2016 USA-isolated JJ259 strain (KT922089) (Figure 2).

Subsequently, employing the RTDT method based on the RdRp genes, a time tree was calculated to examine sapovirus divergence times (Figure 3). The GII genotype showed the earliest origin among eighteen genotypes, dating back to the 1980s. GIV emerged in 1987. GXVII, GVII, and GV shared a common divergence time in 2005. GVIII distinguished itself between 1998 and 2009, while GXI emerged between 2001 and 2014. Notably, the GIII genotype exhibited the most extended span of divergence times, ranging from 1998 to 2014. SHCM/Mega2023, along with Orf1 (MH490911), diverged around 2014.

3.3. Recombinant Analysis

A complete genome recombination analysis utilized the RDP5 software to examine SHCM/Meg2023 and 33 reference strains. Results indicated that SHCM/Meg2023 potentially resulted from recombination between MK962338 (major parent) and KT922089 (minor parent), supported by six detection algorithms (RDP 6.821×10^-20, BootScan 6.750×10^-18, MaxChi 1.314×10^-14, Chimaera 6.32×10^-14, SiScan 1.743×10^-24, 3Seq 6.574×10^-11) (Figure 4). SAV/Pig-wt/ESP/P452/2017 (MK962338) and JJ259 (KT922089) are GIII PoSaV isolates from Spain (2017) and the USA (2016), respectively. Further analysis identifies position 2969-5132 in the SHCM/Meg2023 genome as the predicted recombinant fragment, with the breakpoint within ORF1 (Figure 4).

3.4. Evolutionary Selective Pressures on the ORF1 Codons

To further evaluate the extent of selective pressure on the ORF1 gene, the standardized differences in dN/dS were calculated for each position. Scores > 0 suggest heightened diversifying selection, highlighting the propensity for ORF1 protein regions to diversify. Purifying selection (dN/dS < 1) was indicated by BUSTED models for SHCM/Mega2023 (Figure 5a). In the ORF1 protein, seven amino acids (2 aa, 110 aa, 529 aa, 620 aa, 1668 aa, 1977 aa, and 2232 aa) were subject to diversifying positive selection (ω3 = 10.8), while 1314 amino acids experienced purifying selection (ω1 = 0.0396) (Figure 5b). Three positively selected sites (1668 aa, 1977 aa, and 2232 aa) were situated within ORF1 at capsid protein VP1.

4. Discussion

Pigs in all growth stages can contract PoSaVs, with a higher infection rate observed in post-weaning pigs[17,18]. No significant differences in SaV prevalence were found between age-matched groups of pigs with and without diarrhea in the field. However, PoSaVs frequently co-infect with other enteric pathogens, such as porcine epidemic diarrhea virus, porcine astrovirus, porcine kobuvirus, porcine deltacornavirus, etc.[19]. While most PoSaV genogroups remain unknown[20], a genomic evolutionary investigation of PoSaVs is necessary to provide effective clues for further research on cellular adaptability, pathogenic mechanisms, and prevention and control strategies.

To date, SaV has been categorized into 19 genogroups and 52 subtypes based on VP1 gene sequences, revealing four genogroups (GI, GII, GIV, and GV) and eight genogroups (GIII and GV-GXI) in humans and pigs, respectively[21]. GIII predominantly circulates in global swine herds worldwide[14], consistent with our study's findings. Our metagenomic data confirmed porcine sapovirus as the primary pathogen in piglet diarrhea feces, warranting attention. Phylogenetic analysis identified the prevalent PoSaV within the GIII genogroup. Given sapovirus's current division into 19 genogroups, our study pioneered the calculation and comparison of their evolutionary rates. Results affirmed GIII's widest span of divergence among genogroups, spanning from 1998 to 2014. We hypothesize that GIII exhibits the highest rates of transmission and variation.

SaV exhibits a broad host range with the potential for zoonotic transmission[22]. Recombination, acknowledged as an evolutionary force, can yield new viruses with distinct pathogenesis and virulence[23]. Intra- and inter-genogroup recombinant strains have been reported[7,24,25]. Here, we present the novel PoSaV isolate SHCM/Mega2023, characterized by recombination with MK962338 (major parent) and KT922089 (minor parent), with the breakpoint identified at 2969–5132. This differs from the intra-recombinant PoSaV (p2 strain) reported in Shanghai in 2015[26], which featured a breakpoint in the RaRp-capsid junction region, aligning with most SaV recombination events[27]. Although both SHCM/Mega2023 and p2 strains were reported in Shanghai, they exhibited distinct recombinant events, with the event in SHCM/Mega2023 being novel. The relationship between the novel recombinant event and pathogenesis remains unknown and requires further exploration.

Investigation into the selective pressures acting on a gene is crucial for comprehending its evolutionary dynamics and may unveil potential drug-targeting regions, as amino acid substitutions at negatively selected sites are likely intolerable[28,29]. Positively selected sites aid in identifying epitope-contacting residues[30]. Meanwhile, the presence of positive selection signatures in the genome is a characteristic indication of adaptation, revealing ongoing, recent, or ancient responses to environmental changes throughout a population's evolution[31]. In SaV, the polyprotein encoded by ORF1 undergoes cleavage to form various viral proteins, including 2C-like NTPase, 3C-like protease, VPg, 3D-like RNA-dependent RNA polymerase, and the major capsid protein (VP1). In this study, the selective pressure force of the ORF1 protein of PoSaVs was initially calculated, revealing seven positively selected sites in the ORF1 protein of SHCM/Mega2023. Notably, three positively selected sites (1668 aa, 1977 aa, and 2232 aa) within the capsid protein VP1 of ORF1 are potentially associated with antigenicity. VP1, the main structural protein of PoSaV, safeguards the viral nucleic acid and features multiple immunogenic epitopes that elicit both humoral and cellular immunity.