Molecular Detection and Characterization of Porcine Epidemic Diarrhea Virus and Porcine Aichivirus C Coinfection in México

Swine enteric viral infections are responsible for substantial economic losses in the pork industry worldwide. Porcine epidemic diarrhea (PEDV) is one of the main causative agents of diarrhea in lactating pigs, and reports of PEDV coinfection with other enteric viruses highlight the importance of viral interactions for disease presentation and outcomes. Using next-generation sequencing (NGS) and sequence analyses from samples taken from piglets with acute diarrhea, we explored the possible interactions between PEDV and other less reported pathogens. PEDV coinfection with porcine kobuvirus (PKV) was detected in 36.4% (27/74) of samples. Full genomes from porcine coronavirus and kobuvirus were obtained, as was a partial porcine sapovirus genome (PSaV). The phylogenetic results show the clustering of these strains corresponding to the geographical relationship. To our knowledge, this is the first full genome and isolation report for porcine kobuvirus in México, as well as the first phylogenetic analysis for porcine sapovirus in the country. The NGS approach provides a better perspective of circulating viruses and other pathogens in affected production units.


Introduction
Diarrheic diseases in swine have strong economic impacts on production units around the world. In recent years, porcine coronaviruses have represented some of the principal causes of these diseases; however, reports of naturally occurring coinfection of coronaviruses with other viral agents highlight the importance of less characterized viruses on disease severity and outcome [1,2].
Porcine epidemic diarrhea virus (PEDV) induces an enteric disease with high mortality rates in suckling piglets characterized by acute diarrhea and severe dehydration. It was first identified in Europe in 1978 and in Asia in 1982. Since then, PEDV has had a great impact on the Asian swine industry. In May 2013, the presence of PEDV was reported in México along with the United States and Canada, resulting in heavy economic losses for the swine production units [3,4].
PEDV is a member of the Coronaviridae family, subfamily Orthocoronavirinae, and genus Alphacoronavirus and is an enveloped virus approximately 120 nm in diameter with a

Sample Collection
Samples were collected from 10 different farms distributed in four Mexican states, and all tested animals presented clinical manifestations compatible with gastrointestinal disease at the time of sampling. The commercial farrow-to-finisher units tested were located in the states of Nuevo León, Guanajuato, and Estado de México, and samples from the state of Chiapas belonged to backyard swine herds. The production located in Estado de México had a history of previous diarrheic outbreaks caused by PEDV. For the rest of the tested farms, no previous outbreaks caused by PEDV were reported.
A total of 68 rectal swabs from adult swine and six intestinal tissue samples from piglets in litters presenting clinical signs were taken following euthanasia according to the American Veterinary Medical Association (AVMA) Guidelines for the Euthanasia of Animals: 2013 Edition on the farm. The samples were aliquoted and stored at −80 • C in sterile plastic tubes with Dulbecco's modified Eagle's medium (DMEM) until processing. tion, and quality assessment, the extractions were divided for use in molecular testing and sequencing.

Molecular Testing
cDNA was synthesized using Superscript VILO (Invitrogen, Hercules, CA, USA) and PCR for PEDV detection was performed using Platinum™ Taq DNA Polymerase High Fidelity (Cat no: 11,304,011 Invitrogen™, Hercules, CA, USA) following the protocol described by Li et al. [27] The reaction was conducted under the following conditions: 94 • C for 2 min, (94 • C for 1 min, 55 • C for 1 min, and 72 • C for 90 min) for 40 cycles, with a final period of 10 min at 72 • C using the following primers targeting the PEDV M gene: 5 -AGTCTTACATGCGAATTGACC-3 and 5 -AGCTGACAGAAGCCATAAAGT-3 .
For the detection of PKV, the protocol described by Liu et al. [28] was employed using the primers PK-F10 5 -GGTGGACTCATTGAGTAC-3 and PK-R10 5 -CCTCCCTGGGTGCA GCTTC-3. The samples were also tested for the presence of Deltacoronavirus using the protocol used by Hu et al. with the primers forward 5 -CGCGTAATCGTGTGATCTATGT-3 and reverse 5 -CCGGCCTTTGAAGTGGTTAT-3 [29]. Two PEDV-and PKV-positive intestinal tissue samples of piglets with acute clinical signs were selected for whole genome sequencing and sent to the next-generation sequencing core facility at the Biotechnology Institute, UNAM in cold chain transport.

Sequencing and Bioinformatics Analyses
RNA sequencing of positive samples was performed at the next-generation sequencing core facility at the Biotechnology Institute, UNAM. RNA concentrations were measured (Genova Nano micro-volume spectrophotometer 737,501, Jenway, Staffordshire, UK) and sequencing libraries were constructed using the TruSeq Stranded mRNA Sample Preparation Kit (Illumina, San Diego, CA, USA) using 5 uL of RNA. The samples were then deep sequenced on the Illumina NextSeq500 system, generating paired-end reads of 75 bases. After quality assessments by FASTQC (Babraham Bioinformatics, UK) [30] and adapter elimination using Cutadapt v1.11, the sequence reads were first filtered using Bowtie2 Alignment software against the Sus scrofa reference genome (GCF_000003025.6).
The remaining data were assembled with the SPAdes genome assembler v3.12.0 (Center for Algorithmic Biotechnology, St. Petersburg, Russia) and SSPACE v3.0 tool (BaseClear B.V., The Netherlands) [31][32][33]. The assembled contigs were subject to taxonomic classification using the Kraken metagenomic classification pipeline. To identify contigs corresponding to viral genomes, Basic Local Alignment Search Tool (BLASTn) was used against a viral genome database obtained from the NCBI website (available from: https: //ftp.ncbi.nlm.nih.gov) (accessed on 20 February 2020) [34,35]. Finally, the original reads were mapped to the viral contigs using the BBMap Alignment software to verify the coverage [25]. Prediction of homologous recombination events was performed using RDP4 (Recombination Detection Program) [26]. The analyses were executed at the Center for Genomic Sciences, UNAM facilities.
For the phylogenetic analyses, complete genome sequences of PKV, as well as sequences corresponding to PEDV and the capsid coding region of PSaV, were downloaded from the NCBI database. Then, representative sequences of clusters were selected by CD-HIT-EST [22], and sequences from México were included if available. The representative sequences and obtained genomes were aligned using MAFFT v7.45 (Osaka University, Japan) [39]. For PKV, only the coding sequence (CDS) was used for the phylogenetic reconstruction due to the variability in the lengths of the available sequences. For the PSaV phylogeny, the capsid protein coding region was employed.

Viral Culture
The two original intestinal samples that tested positive for PKV and PEDV by RT-PCR and confirmed by sequencing were propagated separately for virus isolation. Intestinal samples were prepared as 10% homogenates in DMEM and passed through 0.45 µm filters (GE Healthcare, Chicago, IL, USA). The resultant supernatants were inoculated with Vero cells (African green monkey kidney) cultured at 37 • C with 5% CO2 in Dulbecco's Modified Eagle's Media (DMEM; GIBCO, Gaithersburg, MD, USA). Cytopathic effects (CPEs) were evaluated on a daily basis and analyzed for viral replication by RT-PCR after every passage.

Prevalence of Porcine Kobuvirus and Porcine Coronavirus in Diarrheic Pigs
Coinfection of PKV and PEDV was detected in 80% (8/10) of the analyzed farms located in Nuevo León, Guanajuato, and Estado de México states. All the farms presented clinical manifestations at the moment of sample collection. Of the samples used in this study, 51.4% (38/74) were positive for PKV and 67.6% (50/74) were positive for PEDV. Overall, PKV and PEDV coinfection was detected in 36.4% (27/74) of samples. The Deltacoronavirus genome presence was only detected in one farm located in Guanajuato with no PEDV or PKV coinfection.

Viral Genome Detection
Positive PEDV and PKV intestinal samples belonging to two piglets with acute clinical signs were selected for next-generation sequencing. After the elimination of low-quality reads, an average of 7,029,014 reads per sample were obtained. These reads were assembled, and contigs were taxonomically annotated using Kraken and BLASTn. The Kraken results show the identification of known swine pathogens classified in families Coronaviridae, Picornaviridae, and Caliciviridae; other viral agents of families Tospoviridae, Bromoviridae, and Marseilleviridae were also detected. The assembled contigs of swine viruses and their length are shown in Table 1. At the species level, sample 205H showed coinfection by PEDV, PKV, and PSaV, while sample 103H showed evidence of PEDV and PKV coinfection, but no presence of PSaV. Two complete PEDV genomes (28,029 and 28,034 nt) were assembled in addition to one full-length sequence corresponding to porcine kobuvirus Aichivirus C (8116 nt), a partial porcine kobuvirus Aichivirus C sequence (580 nt), and a partial porcine sapovirus sequence (2977 nt) ( Table 2). BLASTn analyses using an NCBI viral genome database showed high coverage and identity between both PEDV genomes and a USA strain (KF272920.1).

Genome Analysis of Porcine Kobuvirus
PEDV and Aichivirus C coinfections were identified in the two analyzed samples. A complete sequence of the porcine kobuvirus genome (GenBank Accession no. MT211964) was obtained in sample 205H, and a contig of 10% of the length of the genome was also generated for sample 103H (Table 2). Phylogenetic analysis based on the complete polyprotein gene revealed that the strain detected in this study formed a cluster with isolate OH/RV11/2011 from the USA and Chgz/2015 (Figure 1).
Sequence analyses of the polyprotein gene revealed a nucleotide identity range of 89.70-90.66% and an amino acid identity between 96.91 and 97.84% to reference strains from the USA and China (GenBank: MF506730.1, KC424640.1), respectively. The previously reported 90 nt motif duplex present in strains from Hungary, Spain, and China and the insertion present at position 135-137 nt of US/OH/RV50/2011 were not present in the assembled sequences [14,21]. No evidence of recombination was found.

Genome Analysis of the Porcine Epidemic Diarrhea Virus
For both samples, full-length coronavirus genomes were assembled. These sequences presented a ≈99.6% identity to local Mexican sequences. Full genome phylogenetic analysis showed that the assembled genomes were located in a cluster with other previously reported no-indel Mexican strains ( Figure 2). ORF3 of the two genomes consisted of 224 aa, with a high nucleotide identity (≈99.8%) to local strains, and, as expected from the phylogeny results, did not have the deletion present in attenuated strains [40,41]. However, the two ORF3 sequences showed mutations at the amino acid level: V → A at amino acid position (aa) 21, V → I at aa 54, F → V at aa 80, and A → T at aa 100.
The analysis of the S gene of both samples showed a nucleotide similarity with local strains of 99.47-99.66% and 99.42-99.78% at the amino acid level (MH004420.1, KJ645708.1) and a 99.66% nt/96.64% aa identity to the reference strain CV777 (NC_003436.1), with most amino acid changes located in the S1 domain. In addition, multiple mutations were located among regions containing neutralizing epitopes at residues 592-607: G → S at aa 594, A → E at aa 605, and in the epitope region 636-789: Y → S at aa 766, N → D 707, N → S at aa 719, N → S at aa 724, and Y → S at aa 766 of the reference genome CV777 S protein. The obtained sequences were uploaded to the GenBank database: MT490315 (103H), MT490316 (205H).

Sequence Analysis of PSaV
A sequence of 2977 nt corresponding to PSaV partial ORF 1 and ORF 2 was assembled from one of the samples, displaying high nucleotide homology (91.13%) to an available Sapovirus GIII sequence (KX688107.1). For VP2, the similarity was 93.36% at the nucleotide

Genome Analysis of the Porcine Epidemic Diarrhea Virus
For both samples, full-length coronavirus genomes were assembled. These sequences presented a ≈99.6% identity to local Mexican sequences. Full genome phylogenetic analysis showed that the assembled genomes were located in a cluster with other previously reported no-indel Mexican strains (Figure 2). ORF3 of the two genomes consisted of 224 aa, with a high nucleotide identity (≈99.8%) to local strains, and, as expected from the phylogeny results, did not have the deletion present in attenuated strains [40,41]. However, the two ORF3 sequences showed mutations at the amino acid level: V → A at amino acid position (aa) 21, V → I at aa 54, F → V at aa 80, and A → T at aa 100.
The analysis of the S gene of both samples showed a nucleotide similarity with local strains of 99.47-99.66% and 99. .78% at the amino acid level (MH004420.1, KJ645708.1) and a 99.66% nt/96.64% aa identity to the reference strain CV777 Figure 1. Phylogenetic tree using the full CDS of PKV. The tree was inferred using the maximum likelihood method and GTR model with gamma distribution and SH-like branch support using PhyML software v3.3.21 (Université de Montpellier, France). The sequence obtained in this study is marked in color.
Viruses 2021, 13, 738 7 of 12 (NC_003436.1), with most amino acid changes located in the S1 domain. In addition, multiple mutations were located among regions containing neutralizing epitopes at residues 592-607: G → S at aa 594, A → E at aa 605, and in the epitope region 636-789: Y → S at aa 766, N → D 707, N → S at aa 719, N → S at aa 724, and Y → S at aa 766 of the reference genome CV777 S protein. The obtained sequences were uploaded to the GenBank database: MT490315 (103H), MT490316 (205H).

Sequence Analysis of PSaV
A sequence of 2977 nt corresponding to PSaV partial ORF 1 and ORF 2 was assembled from one of the samples, displaying high nucleotide homology (91.13%) to an available Sapovirus GIII sequence (KX688107.1). For VP2, the similarity was 93.36% at the nucleotide level and 97.08% at the amino acid level to the USA 2015 sequence (MK965901.1). The capsid gene presented a homology of 97% at the amino acid level with the GIII 2019 USA strain (MK965900.1). Phylogenetic analysis based on the capsid gene region of PSaV showed a clustering with samples from the USA and China (Figure 3).

Viral Culture
Following inoculation and passage with PKV-positive samples, a cytopathic effect was visible 3 to 4 days post-infection (p.i.). After harvesting, PKV RNA was consistently detected in two isolates after several passages. Porcine coronavirus RNA was not detected by RT-PCR in the cells or supernatants in passages in cell culture, and no cell fusion or syncytia characteristic of porcine coronavirus replication was observed (Figure 4).

Viral Culture
Following inoculation and passage with PKV-positive samples, a cytopathic effect was visible 3 to 4 days post-infection (p.i.). After harvesting, PKV RNA was consistently detected in two isolates after several passages. Porcine coronavirus RNA was not detected by RT-PCR in the cells or supernatants in passages in cell culture, and no cell fusion or syncytia characteristic of porcine coronavirus replication was observed ( Figure 4).

Viral Culture
Following inoculation and passage with PKV-positive samples, a cytopathic effect was visible 3 to 4 days post-infection (p.i.). After harvesting, PKV RNA was consistently detected in two isolates after several passages. Porcine coronavirus RNA was not detected by RT-PCR in the cells or supernatants in passages in cell culture, and no cell fusion or syncytia characteristic of porcine coronavirus replication was observed ( Figure 4).

Discussion
Since the introduction of PEDV in México in 2013, it has become one of the primary causative agents for diarrheic manifestations in swine productions. Additionally, PKV was detected in pigs from the majority of analyzed herds. Although PEDV was identified as being responsible for the majority of the outbreaks in the analyzed production units, we were able to identify the presence of PKV in most of these farms. PKV was unknown to be circulating in México and had not been considered in the control measures at the affected farms due to the lack of previous reports in the country.
The presence of both viruses in the same herds with similar clinical signs could indicate a high association of kobuvirus infection with diarrhea but not necessarily as the causative agent since no affected farm was solely positive for PKV. In addition to coronavirus infection by PEDV, diarrheic outbreaks in Mexican porcine farms have been attributed to Deltacoronavirus infection; however, the only farm positive for Deltacoronavirus in this study was negative for coinfection with PEDV and PKV [42].
PKV has been detected in pigs with and without clinical symptoms [43], with the highest prevalence occurring in suckling pigs and 40% positivity in non-diarrheic piglets, as reported by Park, 2010 [13]. It is possible that PKV is an opportunistic pathogen of gastroenteritis in pigs when primary infection by PEDV or other pathogens occurs; however, the possible effects in disease outcome, recovery times, and daily weight gain in surviving animals are unknown. Although the PEDV and PKV coinfection pathogenesis remains unclear, Zhao et al. reported a higher prevalence of this coinfection over the prevalence of infections by PEDV alone in China, which could indicate a possible greater impact and dissemination for coinfection cases [18].
Based on the whole genome phylogenetic analyses, both PEDV and PKV showed a close relation to strains reported in México, China, and the USA. No local sequences for PKV were available at the time of this study. Another agent detected in the NGS analyses was PSaV. The PSaV capsid analysis reflected the same geographical association as PKV and PEDV. The clustering of these viruses reflects the close commercial relation of the swine industry between the USA and México. To our knowledge, this is the first full genome report of porcine kobuvirus in México, as well as the first phylogenetic analysis for porcine PSaV in the country.
The taxonomic classification of de novo assembled contigs by BLASTn reported the presence of several viral agents in the samples, mainly members of Coronaviridae, Picornaviridae, and Caliciviridae families in sample 205H known to cause similar clinical presentations. However, the assembly of full-length genomes for PKV in sample 103H was not possible, and PSaV was not detected. The difference in the detected pathogens may derive from the disparity in the sequencing efficiency, as both animals were housed in close proximity in the same area. Additional analyses based on contigs and filtered reads detected the presence of known commensal and environmental viruses coherent with the type of sample used-likely originating from the digested intestinal content.
During an active outbreak in a production unit, diagnosis is restricted to a small array of agents known to be present in the geographical region; thus, the most prevalent causes such as PEDV are likely to be investigated and controlled, leaving possible secondary infections undiagnosed, which could lead to complications, particularly for feedback-based control. While the implementation of regular viral metagenomics as a regular diagnostic tool could be prohibitive for most farms, an NGS analysis of representative samples from farms could provide a better perspective and understanding of the circulating viruses and other pathogens in the production unit and could be available from a commercial provider. Further studies must be conducted to estimate the prevalence of PKV and PSaV in swine production units and their clinical implications.