Three diarrheal samples were received by the laboratory from pigs exhibiting clinical signs of diarrhea at a pig farm in Zhejiang Province, China. Samples were collected using sterile sampling tools, immediately placed in 5 mL of viral transport medium, labeled, and transported on dry ice to the laboratory for further processing. Upon arrival at the laboratory, 800 μL of sterile 1 × PBS was added to the centrifuge tubes containing the fecal samples. Cotton swabs were agitated using a pipette tip, followed by vigorous vortexing for 2 min. The supernatant was transferred to a new centrifuge tube, followed by centrifugation at 4 °C, 8000 rpm for 2 min. The supernatant was transferred into new tubes, labeled, and stored at −80 °C. All sampling procedures were conducted with the consent of the farm owners and in accordance with institutional animal welfare guidelines.

The Viral RNA was extracted from processed fecal supernatants using the Trizol-based method. A 250 μL aliquot of the clarified supernatant was transferred to a 1.5 mL RNase-free microcentrifuge tube, followed by the addition of 1 mL Trizol reagent. The mixture was incubated on ice for 5 min, after which 200 μL chloroform was added. Following a 2-min incubation at room temperature, the sample was centrifuged at 12,000 rpm for 15 min at 4 °C. The aqueous phase (approximately 550 μL) was carefully transferred to a new tube and mixed with an equal volume of pre-chilled isopropanol. The mixture was incubated at −20 °C for 60 min and subsequently centrifuged at 12,000 rpm for 15 min at 4 °C. The resulting RNA pellet was washed once with 1 mL of ice-cold 75% ethanol and centrifuged at 12,000 rpm for 5 min at 4 °C. After removal of the supernatant, the pellet was air-dried and resuspended in 20 μL of DEPC-treated water. Purified RNA samples were stored at −80 °C until further use.

The viral RNA was extracted using the Trizol method and reverse transcribed into complementary DNA (cDNA) using the PrimeScript™ RT Reagent Kit (Takara, Japan), according to the manufacturer’s instructions. Quantitative real-time PCR (qPCR) was performed to detect the N genes of PDCoV, TGEV, PEDV, and SADS-CoV. Primer and probe sequences used for RT-qPCR are listed in Supplementary Table S1. Each 20 μL reaction contained TaqMan Fast Universal PCR Master Mix, specific primers and probes, and nuclease-free water. Thermal cycling was conducted under the following conditions: 50 °C for 30 min, 95 °C for 10 min, followed by 40 cycles of 95 °C for 15 s and 60 °C for 1 min. Samples with a cycle threshold (Ct) value ≤35 were considered positive.

The total RNA was extracted from fecal samples using Trizol reagent (Life Technologies, USA), and reverse transcription was performed using the PrimeScript™ RT Reagent Kit (Thermo Fisher, USA) to generate complementary DNA (cDNA). PCR amplification was conducted in a 20 μL reaction mixture containing 10 μL of 2 × Taq PCR Master Mix, 0.5 μL each of forward and reverse primers (10 μM), 2 μL of cDNA template, and RNase-free water. The thermal cycling protocol included an initial denaturation at 95 °C for 5 min, followed by 35 cycles of 95 °C for 30 s, 55 °C for 30 s, and 72 °C for 1 min, with a final extension at 72 °C for 5 min. PCR products were separated by electrophoresis on a 1.5% agarose gel and visualized under UV illumination.

The clinical fecal samples were collected from pig farms and subjected to standardized preprocessing procedures, including dehydration, homogenization, and filtration, to ensure sample purity and consistency. The Total RNA was extracted using commercially available RNA extraction kits, and RNA concentration and purity were assessed to confirm suitability for downstream applications. Depending on the experimental design, either mRNA-enriched or total RNA libraries were prepared. Library construction included reverse transcription to synthesize complementary DNA (cDNA), PCR amplification, and quality control of the final libraries. Sequencing was performed on the MGI high-throughput sequencing platform (e.g., MGISEQ or DNBSEQ series), generating whole-genome sequencing data. Raw reads were subjected to rigorous quality control procedures to remove low-quality sequences, and only high-quality reads were retained for subsequent analyses. Bioinformatic processing included sequence alignment and functional annotation using appropriate computational tools. These analyses enabled comprehensive microbial diversity profiling and the identification of potential functional genes and microbial biomarkers. A detailed analytical report was generated based on the results, and all data were securely archived in a dedicated database for future access and research applications (33).

A total of 193 complete PDCoV genome sequences were retrieved from the NCBI database. To assess the evolutionary relationships and recombination-associated phylogenetic incongruence, phylogenetic analyses were performed on seven datasets, including the complete genome, full-length S gene, ORF1b gene, and recombination breakpoint-defined regions of the genome and S gene. Multiple sequence alignments were conducted using MAFFT implemented in BioAider v1.727, and maximum likelihood phylogenetic trees were inferred using IQ-TREE v3.0.1 with 1,000 bootstrap replicates. The best-fit nucleotide substitution models were selected based on the Bayesian Information Criterion. All trees were visualized using iTOL (34).

The complete genome sequences of PDCoV, encompassing all open reading frames and coding regions, were retrieved from the GenBank database. Multiple sequence alignment was performed using MUSCLE (Windows 64) to ensure comparative analysis within a consistent reference framework. Codon usage analysis was subsequently conducted using CodonW software, and the Relative Synonymous Codon Usage (RSCU) values were calculated to evaluate codon preference across the viral genome. Codon usage patterns were visualized using an online tool (http://112.86.217.82:9929/), which generated bar charts and heatmaps illustrating codon usage frequencies in different genes or genomic regions. In addition to RSCU values, the Effective Number of Codons (ENC) was calculated to quantify overall codon usage bias, and the GC content at the third codon position (GC₃) was determined to provide further insight into compositional constraints influencing codon selection. To comprehensively assess codon usage bias, ENC-GC₃ plots (ENC plots) were generated using the Genepioneer online platform. These plots, which depict the relationship between ENC values and GC₃ content, enabled evaluation of the extent to which codon usage is influenced by mutational pressure versus natural selection. This integrated analytical framework provided a systematic understanding of codon usage characteristics in the PDCoV genome and their potential implications for viral gene expression and evolutionary adaptation.

The recombination events in the PDCoV genome were systematically analyzed using RDP4 (Recombination Detection Program, version 4) (35–37). Multiple sequence alignments of the full-length PDCoV genomes and spike (S) gene sequences were performed using the MAFFT algorithm implemented in BioAider v1.727. The aligned datasets were subsequently imported into RDP4 for recombination detection. Recombination analyses were conducted using multiple algorithms, including RDP, BootScan, MaxChi, and GENECONV, with default parameters applied. Recombination breakpoints and putative parental sequences were inferred by RDP4 based on consensus results from multiple methods. Only recombination events supported by supported by two or more independent algorithms and meeting the statistical significance threshold of p < 0.05 were considered reliable and retained for further analysis (Table 1).

This study performed comprehensive molecular and metagenomic analyses on clinical fecal samples collected from diarrheic pigs at a pig farm in Zhejiang Province, China, with the aim of identifying the principal causative agents and associated alterations in gut microbial composition during the outbreak. The qPCR revealed the presence of both PDCoV and SADS-CoV, while PEDV and TGEV tested negative (Figures 1A,B). The detection of SADS-CoV was further validated by agarose gel electrophoresis and sequencing of its N gene fragment (Figure 1C and Supplementary Figure S2A). High-throughput sequencing indicated that the microbial community was predominantly composed of bacteria, followed by viruses and a minor proportion of archaea (Figure 1D).

Within the virome, PDCoV accounted for the highest relative abundance (68%), followed by porcine circovirus (PCV) at 31%, with Posavirus 1 detected at low levels (Figure 1E). Although SADS-CoV was not identified by metagenomic sequencing, subsequent qPCR analysis detected partial genomic fragments, including ORF1a and the spike (S) gene regions. Due to the low viral load, de novo assembly of the complete SADS-CoV genome was not feasible (Supplementary Figure S1 and Supplementary Figure S2B), suggesting the virus may persist in a latent or low-level infection state, warranting continued surveillance. Further analysis of the sequencing data led to the successful identification and assembly of a novel PDCoV strain, designated PDCoV-ZJHZ2024. The complete genome sequence of PDCoV-ZJHZ2024 generated in this study has been deposited in GenBank under the accession number PX492333.

Metagenomic profiling also revealed pronounced dysbiosis in the gut microbiota of the affected pigs. The total bacterial abundance reached 21.6 million reads, with Firmicutes (18.7 M) markedly predominant over Bacteroidetes (1.23 M), resulting in a severely imbalanced Firmicutes/Bacteroidetes ratio. Members of the Bacilli (11.0 M) and Clostridia (4.63 M) classes were significantly enriched. Notably, the conditional pathogen Bacillus cereus was detected at a high abundance (7.47 M), indicating a potential contributory role in the pathogenesis of the outbreak. Although beneficial microbes such as Faecalibacterium prausnitzii, Lactobacillus spp., and Megasphaera elsdenii were present, their relative abundances were substantially reduced, suggesting a loss of ecological competitiveness and compromised intestinal barrier function. The dominance of PDCoV (Coronavirus HKU15) supports its role as the primary pathogen and further implies that viral infection may disrupt intestinal homeostasis, thereby promoting the proliferation of opportunistic bacteria and exacerbating disease severity (Figure 1D).

In summary, PDCoV was identified as the primary etiological agent responsible for this diarrheal outbreak, while SADS-CoV was present at low levels, indicative of latent infection. The pronounced dysbiosis of the gut microbiota—characterized by depletion of beneficial bacteria and expansion of conditional pathogens—was likely triggered by viral infection and may have contributed synergistically to disease progression.

Nucleotide BLAST analysis showed that the newly identified PDCoV strain PDCoV-ZJHZ2024 shared 98.38–99.63% nucleotide identity with previously reported PDCoV strains, with PDCoV-CHLNFX2022 being its closest relative. Comparative genomic analysis between PDCoV-ZJHZ2024 and PDCoV-CHLNFX2022 identified more than 60 nucleotide substitutions distributed across both coding and non-coding regions of the genome (Supplementary Table S2). Notably, a relatively high proportion of nonsynonymous substitutions was concentrated in the spike (S) gene, with approximately 40 amino acid variations detected, suggesting that this region may be subject to strong evolutionary selective pressure.

Several amino acid substitutions were identified in the S protein, including H19751Y, H19753Y, P20255L, and I20385T. These mutations may influence the structural properties or antigenic epitopes of the spike protein, thereby potentially affecting viral entry and host immune recognition.

Within the ORF1ab region, multiple nonsynonymous substitutions were distributed in functional domains associated with viral replication, which may affect viral replication efficiency. Overall, PDCoV-ZJHZ2024 exhibits substantial genetic diversity across multiple genomic regions, reflecting its dynamic evolutionary characteristics. The potential functional consequences of these mutations warrant further experimental investigation.

Using 193 reference sequences retrieved from the NCBI database, full-genome and spike (S) gene phylogenetic analyses were performed for PDCoV-ZJHZ2024 (Figures 2A,B). The full-genome phylogenetic tree revealed that PDCoV strains segregate into three major lineages: the Chinese lineage, the American lineage, and the Southeast Asian lineage. The Chinese lineage primarily comprises strains isolated within China, whereas the American lineage mainly includes strains from the United States, Peru, Japan, and South Korea. The Southeast Asian lineage consists of strains from Thailand, Vietnam, and Laos, exhibiting considerable genetic and geographic diversity.

Notably, several strains from South Korea and Japan clustered within the American lineage, suggesting possible introductions from North America. In addition, two strains from China (MN520191 and MN520198) were also grouped within the American lineage, indicating that American-lineage PDCoV strains may have been introduced into China through transboundary transmission.

In the full-genome phylogenetic analysis, PDCoV-ZJHZ2024 clustered within the Chinese lineage and showed the closest phylogenetic relationship to PDCoV-CHLNFX2022 (ON968724), suggesting a potential epidemiological linkage between the two strains. However, phylogenetic analysis based on the S gene (Figure 2B) revealed a partially incongruent clustering pattern compared with the full-genome tree. Specifically, the S gene of PDCoV-ZJHZ2024 showed the highest sequence similarity to PDCoV-CH-SDLY52-2021 (OM256446) rather than PDCoV-CHLNFX2022. This topological incongruence between the full-genome and S gene phylogenies suggests that a recombination event involving the spike gene may have occurred, a phenomenon commonly observed in coronaviruses, and this inference was further supported by subsequent recombination analyses.

Further comparison of the phylogenetic positions of PDCoV-ZJHZ2024 across different genomic regions indicated region-specific divergence. The full-genome phylogeny suggested a closer relationship to strains from Liaoning Province, whereas the S gene phylogeny showed greater similarity to strains from Shandong Province. This pattern supports the occurrence of inter-provincial genetic exchange among PDCoV strains circulating in China.

Codon usage bias in the newly sequenced PDCoV strain PDCoV-ZJHZ2024 was analyzed based on relative synonymous codon usage (RSCU) patterns and the distribution of the effective number of codons (ENC). RSCU analysis evaluates deviations between the observed frequencies of synonymous codons and their expected frequencies under equal usage assumptions, thereby reflecting codon usage preferences and potential evolutionary constraints within the viral genome.

The RSCU results showed that GCU, GGU, and CCA were the most frequently used codons across the complete genome (Figure 3A). In contrast, the S gene exhibited preferential usage of codons such as UAG, GCA, GAG, AGA, AGC, and GUU (Figure 3C), indicating distinct codon usage patterns among different genomic regions.

Furthermore, ENC–GC3 plot analysis revealed that data points from both the complete genome and the S gene of PDCoV-ZJHZ2024 deviated from the theoretical curve expected under mutation pressure alone (Figures 3B,D), suggesting that codon usage bias may be influenced by natural selection. These results indicate that, in addition to mutational pressure, evolutionary selection factors may play a role in shaping codon usage patterns in PDCoV.

Systematic recombination analyses revealed that PDCoV-ZJHZ2024 possesses a clear mosaic genomic architecture, indicating that its evolutionary trajectory has been strongly shaped by recombination events. At the whole-genome level, multiple algorithms implemented in the RDP package (GENECONV, BootScan, MaxChi, Chimaera, SiScan, and 3Seq) consistently detected a highly significant recombination event (Event No. 12; Table 2). This recombinant fragment was located at the ORF1b–S junction region, spanning alignment positions 18,958–20,798 (18,941–20,778 after gap removal), with both breakpoints supported by narrow 99% confidence intervals. The extremely low p-values obtained across methods (global KA, p = 1.39 × 10⁻¹⁹) further confirmed the robustness of this event. Consistently, MaxChi analysis revealed a clear parental lineage switch within the ORF1b–S region, with −log₁₀(p-values) exceeding the Bonferroni-corrected threshold (Figure 4A and Supplementary Figure S4A), while SimPlot analysis showed an abrupt shift in nucleotide similarity between PDCoV-ZJHZ2024 and the putative parental strains across the same region (Figure 4C and Supplementary Figure S4C), jointly supporting a reliable whole-genome-scale recombination event.

In addition to the large-scale recombination event, recombination signals were also detected within the spike (S) gene of PDCoV-ZJHZ2024 (Event No. 4). This event spanned alignment positions 130–1,828 of the S gene (130–1,823 after gap removal) and was identified by the MaxChi, Chimaera, and SiScan methods (Table 2). Although supported by fewer algorithms and showing lower statistical significance than the whole-genome event, MaxChi and SimPlot analyses still revealed local changes in sequence similarity between PDCoV-ZJHZ2024 and the putative parental strains within the corresponding region (Figures 4B,D and Supplementary Figure S4B). Consistently, region-specific phylogenetic analyses demonstrated a degree of topological incongruence among trees reconstructed from different S gene fragments (Supplementary Figure S3), suggesting a localized genetic rearrangement within the S gene, albeit with weaker supporting evidence.

To further validate the parental origins of the detected recombination events, phylogenetic trees were reconstructed separately based on recombinant and non-recombinant genomic regions. In trees inferred from the non-recombinant regions (alignment positions 1–18,957 and 20,799–25,422), PDCoV-ZJHZ2024 clustered with the inferred minor parental lineage, whereas reconstruction based on the recombinant fragment (positions 18,958–20,798) placed PDCoV-ZJHZ2024 within the major parental clade (Figures 5A,B). The pronounced topological incongruence between different genomic regions represents a hallmark of recombination-driven evolution and further supports the mosaic genomic structure of PDCoV-ZJHZ2024. Notably, phylogenetic analysis of the ORF1b region showed that PDCoV-ZJHZ2024 retained a conserved evolutionary position consistent with the non-recombinant portions of the genome (Supplementary Figure S4D), indicating that ORF1b was not involved in the major recombination event and may serve as a reliable reference for defining recombination boundaries. Collectively, these findings demonstrate that whole-genome-scale recombination constitutes the primary evolutionary mechanism shaping PDCoV-ZJHZ2024, whereas the recombination signal detected within the S gene likely represents a secondary or localized event.