Anal swabs (cotton swabs soaked in physiological saline and inserted 2–3-cm deep into the anus of the sampling dog; these swabs were gently rotated, applied inside the anus, and then inserted into a clean test tube) from 276 dogs (182 with diarrhea and 94 healthy) were obtained from pet hospitals in five provinces (Henan, Inner Mongolia, Jiangsu, Hubei, and Anhui) of China between 2020 and 2023.

Each swab was washed with 1 mL of phosphate-buffered saline solution (0.01 mol/L), and 200 μL of suspensions were separated for viral DNA/RNA extraction using Simply Viral DNA/RNA Coextraction Kit (Bioer Biotechnology, Inc., Hangzhou, China) according to the manufacturer’s instructions. The extracted DNA and RNA were stored at −80°C for future use.

As previously reported, CachaV detection was performed using nested polymerase chain reaction (nt-PCR) as follows: the first PCR used the outer primer set IF (5′- CAACTAGCCGAATGCAGGGA-3′) and IR (5′-CGATAACATCCCCGGACTGG-3′) and the nested PCR used the inner primer set IF (5′-AGCTCAGTTTGGCCCAGATC-3′) and IR (5′-AGAGGGATCGCTGGATCTGT-3′) (10). In addition, canine parvovirus (CPV-2) (21), canine distemper virus (CDV) (22), canine coronavirus (CCoV) (23), and canine bufavirus (CBuV) (24) were screened in these collected samples using PCR/RT–PCR as previously described. The infection and co-infection statuses of the samples were summarized and displayed using UpSet plot packages (25). The prevalence of CachaV was compared between healthy and diarrheal dogs using Fisher’s exact test via GraphPad Prism 9.5 (San Diego, California, United States). Statistics were deemed significant at p < 0.05.

The nearly complete genome sequences of the nine CachaVs were amplified using specific primers as described previously (19). Overlapped amplification of CachaV genome segments was performed using PCR in a 20-μL reaction system containing a template (>100 ng/L), 10 pmol forward/reverse primer sets, and Ex-Taq polymerase (TaKaRa Biotechnology Co., Ltd., Dalian, China). The reaction procedure was as follows: predenaturation at 94°C for 5 min; followed by 35 cycles of combined denaturation at 94°C for 45 s, annealing at 56°C for 45 s, extension of 72°C for 70 s; and final extension at 72°C for 10 min. After ligation of the PCR amplicons using ClonExpress Ultra One Step Cloning Kit (Vazyme Biotechnology Co. Ltd., Nanjing China), positive clones were sent for sequencing to Generalbiol -Biotechnology, Chuzhou, China.

To analyze the genome sequence of CachaV and demonstrate its phylogenetic relationship, the sequence segments were assembled using DNAStar7.0 software (DNASTAR Inc., Madison, WI, United States) and aligned with all chaphamaparvoviruses. We used BioAider kit for the differential analysis of the 9 studied CachaV strains and 24 CachaV strains. A phylogenetic evolutionary tree was constructed using the genomes of the 9 studied CachaV strains; 24 CachaV strains detected in dogs, wolves, or cats in China or other countries (Supplementary Table S1). According to the MODELS program in MEGA 11 software, the HKY + G + I model, JTT + G model, and JTT + G model were used for the whole genome phylogenetic tree, NS1 phylogenetic tree, and VP1 phylogenetic tree, respectively. Evolutionary trees were further constructed based on the amino acid (aa) sequences of VP1 and NS1 using 1,000 guided replicates and the maximum-likelihood method via MEGA 11 software (26). Furthermore, recombination was predicted in the strains evaluated in this study using RDP4. Potential recombination events detected via three or more programs along with the identity analysis of the parents were considered potential events, with the highest acceptable p-value cutoff of 0.05 (27). RDP4 is a computer program typically used for recombination prediction; it includes guided scan, MAXCHI, mosaic, 3SEQ, gene cloning, LARD, and SISCAN in addition to the traditional RDP method that can predict recombination (28).

Based on the Meg-Align results, the aa sequences of NS1 and VP1 of the nine studied strains were compared with the reference CachaV strains identified in dogs or cats. Antigenic epitope prediction was performed using DNAMAN 5.2.2 for the resulting and reference strains (LynnonBiosoft, America) (Supplementary Table S2). To understand the molecular characteristics of the coding proteins, IDEXX1 (accession number: MH893826) and two variant CachaV strains (CHN230216 and CHN220916) were subjected to structural homology modeling for NS1 and VP1 using SWISS-MODEL2 and were visualized using PyMOL Molecular Graphics System 2.3 (DeLano Scientific LLC, America).

Among the 276 samples, 9 (3.3%) tested positive for CachaV after viral screening, including 8 (8/182, 4.4%) from dogs with diarrhea in Henan and Jiangsu provinces and 1 (1/94, 1.1%) from a healthy dog in Anhui province. Notably, more than half of the total positive cases were detected in Henan and Jiangsu provinces. Table 1 shows the clinical information of the CachaV-positive dogs. Statistical analysis revealed no association between CachaV infection and clinical symptoms (p > 0.05). Meanwhile, two dogs with diarrhea were quadruply infected (CachaV + CPV-2 + CCoV + CDV), one dog with diarrhea was triply infected (CachaV + CDV + CBuV), and one dog with diarrhea was dually infected (CachaV + CPV-2). The infection and co-infection statuses of each pathogen screened are shown via UpSet and Venn plots in Figure 1.

The nine identified strains had the same genome structure as IDEXX1, each with two main codes: NS1 (663 aa) and VP1 (504 aa). Notably, the aa sequence identities of NS1 and VP1 in the nine CachaV strains with IDEXX1 were 97.89–98.94% and 97.42–99.8%, respectively. Compared with the 24 previously reported CachaV strains, the aa sequence identities of NS1 and VP1 in the 9 obtained CachaV strains were 97.59–99.55% and 96.23–99.8%, respectively. Moreover, the NS1 of CHN230521 showed the highest identity (99.85%) with NWT-W88 and CP-T015, whereas CHN220318 showed the lowest identity (96.23%) with OM640109. For the VP1 protein, CHN230521 exhibited the highest identity (100%) with NWT-W88 (accession no.: OK546101, originated from wolf, Canada, 2009), CP-R107C (accession no.: OP225937), CP-T015 (accession no.: OP225942), and CP-T046 (accession no.: OP225944), whereas CHN220916 showed the lowest identity (95.83%) with CY56 (accession no.: OM640109).

Phylogenetic analysis was performed using the genome sequences and aa sequences of NS1 and VP1 in the 9 identified strains and 24 CachaV strains (Figure 2). As shown in the evolution trees, CHN230521were closely related to some CachaV strains collected in Thailand, Canada and the United States, and the remain eight obtained strains mainly clustered with Chinese CachaV strains and displayed different branching, respectively. Most strains obtained in this study belonged to the same branch with the two cat-originated strains previously reported in China. Notably, no obvious recombination signal was observed in our comparative analysis with all chaphamaparvoviruses, despite the complex whole genome evolution of these nine strains.

Through alignment of the derived aa sequences of the NS1 protein between the obtained CachaV and 24 CachaV strains, the following mutation sites were detected only in the nine obtained strains: Met104Leu, Asn207Asp., His220Arg, Gln469Arg, and Gly519Ser (Supplementary Table S3). The CHN210713 and CHN220916 strains showed unique site mutations at 99 (Phe → Ser), whereas the remaining seven strains had undergone aa changes at 456 (Leu → Thr). By comparing the aa sequences of CachaV-NS1, continuous mutations of SerGlyGlyTyr252-255CysLeuThrPhe were identified, which were common in other CachaV reference strains detected in China (Supplementary Table S4). For NS1, the CHN210713 strain identified in this study had mutations at loci 131 (Phe → Ser), 238 (Asp→Gly), and 247 (Val → Ala). Notably, the mutations generated at N9D in VP1 were present only in the CHN220916 and CHN230216 strains. To better understanding these mutation sites, mutations related to the tertiary structure change in NS1 were predicted and are displayed in Figure 3. The NS1 tertiary structure is predicted to vary at S252C which caused the change ofα-helix and random coil region, whereas these sites in VP1 have not been predicted to undergo structural changes.