Fecal samples were obtained from diseased and healthy cats and dogs in Sichuan province. Ascites samples were obtained from cats developing peritonitis in Sichuan province. Viral RNA genomes were extracted using Sangon Viral RNA Kit (Sangon, Shanghai, China) and stored at −80°C. The cDNA was synthesized using the Vazyme HiScript III RT SuperMix Kit (Vazyme, Nanjing, China) or the RT Easy™ I Master Premix Kit (Foregene, Chengdu, China). FCoV and CCoV were detected using the previously reported primers (Herrewegh et al., 1995; Naylor et al., 2001). The S1 gene, N gene, and M gene of the positive samples were amplified with the PrimeSTAR ^® Max DNA Polymerase (Takara, Japan). Purified RT-PCR products were cloned into the pMD-19T vector. The positive recombinant plasmids were sent for sequencing (Sangon, Shanghai, China). Multiple sequence alignment was performed using MACSE (Ranwez et al., 2011) in PhyloSuite (version 1.2.2) (Zhang D. et al., 2020) and then manually adjusted in MEGAX (Kumar et al., 2018).

A nucleotide identity matrix was generated by Bio-Aider (Zhou et al., 2020). Nucleotide diversity (π) and haplotype diversity (Hd) of FCoV were calculated using DNAsp5 (Librado and Rozas, 2009). F_ST, a population genetic differentiation value, was calculated using Arlequin 3.5 (Excoffier and Lischer, 2010). K_ST and Snn values were calculated using DNAsp5 (Librado and Rozas, 2009). The hypothesis of deviation from zero population differentiation was implemented by 1,000 permutations of the initial data.

The discriminant analysis of principal components (DAPC) for FCoV was performed with the R package and adegenet library by geographic location (Jombart et al., 2010). Due to insufficient CCoV samples in different geographic regions, only FCoV DAPC analysis was completed. The dataset was first imported and transformed based on PCA. K-means algorithm was performed for clustering, and LD analysis was performed to evaluate the clustering condition under different K values.

Gene recombination analysis was performed using RDP5 (Martin et al., 2021) for different S1, N, and M genes. According to previous reports, TGEV is also involved in the recombination; hence, the corresponding gene sequences of TGEV were added. We used seven methods such as RDP (Martin and Rybicki, 2000), 3Seq (Lam et al., 2018), GENECONV (Sawyer, 1999), Chimera (Posada and Crandall, 2001), SiScan (Gibbs et al., 2000), MaxChi (Smith, 1992), and BootScan (Salminen et al., 1995) in RDP5. Credible recombination events could be detected by at least three of the seven methods with a P-value cutoff of 10^–6. Simplot was applied to further show the breakpoint position of the recombination event.

The maximum likelihood phylogeny for S1, M, and N genes was inferred using IQTREE version 2.1.3 (Nguyen et al., 2015) in PhyloSuite (version 1.2.2) (Zhang D. et al., 2020) with SeACoV reference sequence as outgroup. The substitution model calculated by ModelFinder (Kalyaanamoorthy et al., 2017) in PhyloSuite (version 1.2.2) (Zhang D. et al., 2020) was GTR + F + R10 for S1 gene, GTR + F + R5 for N gene, and GTR + F + R6 for M gene. All of the maximum likelihood trees were with 10,000 ultrafast bootstrap replicates to estimate branch support. Finally, the ML trees were visualized and annotated in iTOL v6¹.

To reveal the natural selection between the different lineages of FCoV and CCoV, the branch model was performed in EasyCodeML (version 1.4) (Gao et al., 2019) by a non-temporal ML tree and a PML format non-recombinant sequence file. In addition, to find the site that was positively selected on the S1 gene, we implemented a site model for each subgenotype sequence set. Only when the probability of identified positive sites is greater than 0.95, it could be credible.

We performed the phylodynamic analysis using the S1 gene of FCoV-I and CCoV-IIa. Firstly, we used TempEst (version 1.5.3) to preliminarily assess a temporal signal of S1 gene dataset and then remove some sequences with an abnormal temporal signal. Finally, we used Bayesian evaluation of temporal signal (BETS) (Duchene et al., 2020) method to evaluate the final time structure. The generalized stepping-stone sampling method was performed in BEAST (version 1.10.4) to obtain the marginal likelihood estimation of heterochronous model (Mhet) and isochronous model (Miso) and infer Bayes factor (BF). If log BF > 5, it meant that the dataset had an excellent temporal structure.

Path-sampling and stepping-stone sampling methods are used to calculate the best molecular clock model and tree prior model (Suchard et al., 2018). The result showed that the strict clock model with the Bayesian Skygrid coalescent model (Gill et al., 2013) is the best model for FCoV S1. Also, the uncorrelated lognormal relaxed clock model with the Bayesian skyline coalescent model (Minin et al., 2008) is the best model for CCoV S1. BEAST runs were performed using four independent chains with a chain length of 1 × 10⁸ for FCoV and one chain with a chain length of 5 × 10⁷ for CCoV with GTR + F + I + G4 substitution model. After data convergence (ESS > 200) in Tracer (version 1.7.1), the log and tree file were mixed by LogCombiner (version 1.10.4) with a burn-in period of 10% total chain length and all trees. The MCC tree was annotated by TreeAnnotator (version 1.10.4) and visualized by Figtree (version 1.4.4). In addition, we used Skygrid analysis and Bayesian skyline analysis by Tracer (version 1.7.1) to reconstruct the effective FCoV and CCoV population history over time, respectively.

To allow better bioinformatics analysis, 28 S1 genes, 34 N genes, and 34 M genes (coding-complete sequences) were obtained from FCoV-positive samples, while 10 S1 genes, 12 N genes, and 12 M genes were obtained from CCoV-positive samples (Supplementary Table 1). The nucleotide sequences of the entire genome, S gene, M gene, and N gene of selected FCoV and CCoV strains were retrieved from GenBank (Supplementary Table 1).

The average nucleotide identity of FCoV S1, N, and M genes from Sichuan, China, was 79.71% (41.6–99.96%), 92.9% (90.04–99.91%), and 88.67% (79.82–99.87%) (Supplementary Tables 2–4), while that of CCoV was 81.18% (38.29–98.64%), 91.13% (77.11–99.91%), and 91.61% (80.88–99.49%), respectively (Supplementary Tables 5–7). In addition, the average nucleotide identity between FCoV and CCoV was 48.52% (39.66–99.05%), 80.59% (76.85–94.8%), and 80.1% (79.25–99.75%), respectively. We observed that the S1 gene of either FCoV or CCoV showed very low similarity (about 40% at the lowest), even within one species. Those findings suggest that there might be a genetic relationship between them or a gene exchange and co-evolution.

We identified 22 haplotypes among the FCoV S1 genes in our study, with a haplotype diversity of 0.958 and a nucleotide diversity of 0.171 (Supplementary Table 8). The nucleotide diversity of FCoV S1 in Sichuan was lower than that of other isolates in China, Japan, and the United States, but higher than that of isolates in the United Kingdom and Netherlands (Supplementary Table 8). A total of 25 and 26 haplotypes were identified among the N and M genes in our study with a respective haplotype diversity of 0.975 and 0.979 and a nucleotide diversity of 0.063 and 0.108, respectively (Supplementary Table 8). By contrast, the nucleotide diversity of FCoV N and M genes in Sichuan was higher than that of other isolates in China. K_ST, S_nn, and F_ST, all three population genetic differentiation values, were statistically significant (p < 0.05) (Supplementary Table 9). Compared with other strains in China, FCoV N genes in Sichuan were highly differentiated, followed by S genes (Supplementary Table 9). In general, Sichuan strains are rigorously differentiated in comparison with the strains in the United States and the United Kingdom, but were less differentiated than the Dutch strains.

Discriminant analysis of principal components showed that the optimal cluster K value of S1, N, and M genes was seven, four, and six, respectively (Supplementary Figures 1–3). It was obvious for all viral genes that Sichuan FCoV groups shared many genetic components with the Dutch FCoV groups, but only a few with the Japanese and American FCoV groups in terms of S, M, or N genes (Supplementary Figures 1–3). The DAPC scatter plot also indicated that Sichuan FCoV groups and other FCoV populations in China were genetically related to populations in Netherlands, Japan, and the United States (Supplementary Figure 4).

We found multiple reorganization events, including those reported in previous studies and newly discovered. A total of 12 S1, 4 N, and 8 M gene recombination events were detected, respectively, by RDP5 (Supplementary Table 10). CCoV NC0604 had the same recombination pattern as the previously reported CCoV-IIb (Figure 1A). This was the first time that a recombination event had been reported in mainland China. The C-terminal domain of CCoV NC0604 S1 gene exhibited high similarity with CCoV-IIa, while the N-terminal domain exhibited high similarity with TGEV (Figure 1A). We found a putative recombination receptor-binding sequence “TTACTACAG” in CCoV NC0604 strain, which shared high similarity to the putative donor-binding sequence “TTACGCAA” (Decaro et al., 2009). In our study, recombination events also occurred between CCoV NC0521 and FIPV MY0628 in the N-terminal domain of the S1 gene (Figure 1B). This was similar to the recombination pattern of CCoV CB/05, because of their roughly same breakpoints (Supplementary Table 10). We also found that the recombination in the C-terminal domain, such as FIPV CD0402, was similar to CCoV HLJ/HRB/2016/13 (Figure 1C). Obviously, these recombination events did not occur among inter-lineages. In our study, we detected two N gene recombination strains (FIPV DY0523 and FIPV DY0612) through RDP5 (Figure 1D). The recombination breakpoint of FIPV DY0523 was 741 bp and, which was similar to FCoV UU30 (HQ392472). The recombination breakpoint of FIPV DY0612 was 581 bp, which was similar to PC/M477/06 strain (GU017103) (Figure 1D). It was obvious that the N gene recombination events occurred more often in the C-terminal domain. Similarly, we found that the M genes of FIPV LS0526 and FIPV MY0628 had recombination events at 271 and 341 bp, respectively, which were similar to other strains (FCoV KUK-H/L and FCoV 79-1683) (Figures 1E,F). These popular recombination events could be related to virus assembly, invasion, and pathogenicity.

Natural selection pressure of inter-lineage was calculated by a branch model in EasyCodeML with FCoV-I ω = 0.18993 (Figure 2A), FCoV-II ω = 0.00496 (Figure 2A), CCoV-I ω = 0.0001 (Figure 2B), CCoV-IIa ω = 0.22195 (Figure 2B), CCoV-IIb ω = 0.52488 (Figure 2B), and CCoV-IIv ω = 0.00148 (Figure 2B). The statistical P-value of all analyses was significant (p < 0.05). Two FCoV subtypes were subject to negative selection, but FCoV-II was subject to stronger constraints than FCoV-I. All subtypes of CCoV were constrained, and the order of the strength from strong to weak was CCoV-I, CCoV-IIv, CCoV-IIa, and CCoV-IIb. Interestingly, when analyzing CCoV-IIb from different hosts, we found that the canine-derived CCoV-IIb group was in negative selection (ω = 0.05935) (Figure 2C), while the human-derived CCoV-IIb group was in strong positive selection (ω = 20.89251) (Figure 2C). This result might be related to the sample size.

Then, we applied the site model in EasyCodeML to analyze positive selection in each intra-lineage with a sufficient sample size. The P-value of all positive selection analyses was also significant (p < 0.05) (Supplementary Table 11). Nine positive sites were detected in FCoV-I subtype (reference: FCoV black, EU186072) (Table 1). Three positive sites were detected in FCoV-II subtype (reference: FCoV 79-1146, AY994055) (Table 1). We detected nine and ten positive sites in CCoV-IIa (reference: CCoV strain 171, KC175339) and CCoV-IIb (reference: CCoV UCD-1 AF116248), respectively (Table 1). The positive sites of FCoV-I and CCoV-IIb were more evenly distributed in the S1 gene, and those of FCoV-II and in CCoV-IIa group were mainly located in the N-terminal domain, the receptor-binding domain (Figure 3). The positive sites 786, 787, and 799 of FCoV-I were located on both sides of the S1/S2 recognition sequence “RRXRR,” which could be related to the cleavage of S protein during the cell invasion process of the virus. Each positive selection site had high variability and the properties of variable amino acid residues were different, which could help the virus to face different living environments (Table 1).

We identified 37 S1 genes from our samples, of which 27 strains were clustered into FCoV-I subtype, eight strains were clustered into CCoV-IIa subtype, and only one FCoV-II strain and one CCoV-IIb strain were recognized (Figure 4). All FCoV-I strains and CCoV-I strains in our study were clustered in the same branch as the reference strains, such as UU strain from Netherlands and other strains from China (Figure 4). The FCoV-II and CCoV-IIa strains of our samples were clustered together with the Vietnamese strain CCoV/dog/HCM27/2014, the Tanzanian strain SH36_2004, and the Japanese strain Tokyo/cat/130627 (Figure 4). The only CCoV-IIb strain NC0604 was clustered together with the British CCoV strain 2020/7 and the Australian CCoV strain CCoV/7/2020/AUS, but it diverged from the human-derived CCoV (Figure 4). The FCoV strains from Taiwan, China were clustered together with the Japanese strains and showed certain genetic differences from the strains in the Chinese mainland (Figure 4).

According to the phylogeny results of N genes, most FCoV-I strains were clustered together with the Dutch strains, while most CCoV strains and Vietnamese strains were clustered together (Supplementary Figure 5). Their certain species characteristics could assist in species identification. Unlike the S1 genes and N genes, the M genes of FCoV and CCoV have an obvious connection (Supplementary Figure 6). In general, part of the M genes identified from our FCoV-positive samples was clustered with strains from the United States and Brazil. The M genes of seven strains were still in the same branch as the Dutch strains. Interestingly, there were five samples in a separate clade clustered with the Japanese strain FCoV Tokyo/cat/130627. Most of the CCoV M genes were clustered with the Chinese strains.

Due to the limitations of the number of sequences, we only performed phylodynamic analysis on the S1 gene of the classical strains from FCoV-I and CCoV-IIa subtypes. The dataset had an optimal temporal signal (FCoV R² = 0.5816, CCoV R² = 0.706) using TempEst (Supplementary Figure 7). Further BETS analysis results showed that log BF(FCoV) = 2400.29 > 5 and log BF(CCoV) = 19.439 > 5, so Mhet was better than Miso and the dataset had an excellent temporal signal (Supplementary Table 12). The MCC trees showed the temporal historical origin of FCoV and CCoV based on the S1 gene. The mean evolutionary rate of FCoV-I S1 genes was estimated at 1.244 × 10^–3 (95% HPD: 1.13 × 10^–3−1.34 × 10^–3) subs/site/year. The time of the most recent common ancestor (tMRCA) of FCoV-I was estimated at 1822 (95% HPD: 1801–1841), which was probably the origin time of the British FCoV population and time of differentiation with the Netherlandish population (Figure 5A). The estimated tMRCA of the Netherlandish FCoV-I groups was 1862 (95% HPD: 1855–1876), and that of the United States population was 1887 (95% HPD: 1876–1899) (Figure 5A). We also estimated that the earliest origin of the Chinese FCoV-I population was 1876 (95% HPD: 1866–1884) and the tMRCA of the Chinese population was 1886 (95% HPD: 1876–1892) (Figure 5A). Then, the mean evolutionary rate of CCoV-IIa S1 gene was estimated at 1.281 × 10^–3 (95% HPD: 9.40 × 10^–4−1.63 × 10^–3) subs/site/year. The tMRCA of CCoV-IIa was 1901 (95% HPD: 1870–1929) derived from the German CCoV population (Figure 5B). The origin time of CCoV-IIa of Tanzania was 1917 (95% HPD: 1886–1946) (Figure 5B). The CCoV-IIa of a few Chinese strains was originated in 1979 (95% HPD: 1968–1991), but we could clearly find that the origin time of most Chinese strains in the dataset was 1983 (95% HPD: 1975–1993) (Figure 5B). By combining the estimation of the S1 genes of FCoV-I and CCoV-II, we concluded that the tMRCA of FCoV and CCoV was 4,598 years ago, which is 2578 BC (95% HPD: 4808 BC–500 BC) (Supplementary Figure 8).

The analysis of population historical dynamics shows that the FCoV population was in a stable state before 1920 (Figure 6A). After that, the population size began to expand rapidly. However, it began to shrink around 1975 (Figure 6). Although the population size fluctuated slightly after 2000, the overall trend was shrinking at a historically low level (Figure 6A). The CCoV population size had grown slowly before 1970. After that, it has begun to shrink slowly (Figure 6B). However, it began to shrink rapidly in 2000 and expand sharply in 2005 (Figure 6B). Although there was a slight decline thereafter, the overall population size was at a historically high level.