From January to March 2021, 460 feces samples of black-headed gulls were collected at seven sampling points near Erhai Lake in Dali, Yunnan Province, and Dianchi Lake in Kunming, Yunnan Province. All samples were placed in virus transport medium (VTM) and kept in dry ice for transportation to the laboratory and stored at −80°C until use.

Viral RNA was extracted from the fecal samples using MagaBio plus Virus DNA/RNA Purification Kit II (BioFlux, China) following the manufacturer’s instructions. RNA extracted from fecal was washed off with 70-μl RNase-free water and stored in −80°C as a template for reverse-transcriptional PCR (RT-PCR). Initial CoV screening was performed by amplifying a 440-bp fragment of the rdrp gene of CoVs using primers DCoV-F (5′-GTGGVTGTMTTAATGCACAGTC-3′) and -R (5′-TACTGYCTGTTRGTCATRGTG-3′) (Lau et al., 2018). Primescript™ one-step RT-PCR kit ver. 2 (Takara) was used for PCR amplification. About 20 μl PCR mixture included 0.8 μl Primescript 1 step enzyme, 10 μl 2 × 1 step buffer, 5.2 μl ddH₂O, and 2 μl extracted RNA. The mixtures were amplified for 40 cycles as 94°C for 30 s, 48°C for 30s, and 72°C for 40 s and a final extension at 72°C for 10 min in an automated thermal cycler (Applied Biosystems). To confirm the bird species, the mitochondrial cytochrome b (cytb) gene was amplified and sequenced of the samples.

The PCR products were gel-purified and sequenced with an ABI Prism 3,700 DNA analyzer (Applied Biosystems), using the PCR primers. The sequences of the PCR products were compared with known sequences of the rdrp genes of CoVs in the GenBank database.

Three complete genomes of BHG-DCoVs were amplified and sequenced using RNA extracted from black-headed gull feces as templates. RNA was amplified with degenerate primers, which were designed by multiple alignments of other coronavirus genomes with the complete genome, using Primescript one-step RT-PCR kit version 2. Additional primers were designed according to the results of the first and subsequent rounds of sequencing. The 5′ and 3′ genome end sequences were obtained by 5′ and 3′ race (Roche), respectively. The expected size of PCR products was purified by gel and directly sequenced. The sequence was assembled to obtain the full-length genome sequence.

Multiple sequence alignment was performed by MAFFT v7.149 in BioAider v1.334 (Katoh and Standley, 2013; Zhou et al., 2020). The maximum likelihood phylogenetic tree of CoVs’ RdRp was constructed by IQ-tree v1.6.10 in PhyloSuite program with 10,000 ultrafast bootstraps, and the most appropriate substitution model of aa was calculated using ModelFinder according to the Bayesian information criterion (BIC) method (Nguyen et al., 2015; Zhang et al., 2019). The phylogenetic tree of ORF1ab, NIRAN, 3CL-pro, HEL1, and S protein was performed using MEGA7 with 1,000 bootstraps in neighbor-joining method (Kumar et al., 2016).

The rdrp gene was aligned using MAFFT program with codon method in BioAider v1.334. Then, we detected the temporal structure in these rdrp gene sequences by TreeTime program (Sagulenko et al., 2018). The correlation coefficient of R² was 0.01, and there was very weak signal between sampling time and genetic distance in the data (Supplementary Figure 1). Therefore, we used a uniform distribution priori value (from 8 × 10^–5 to 2 × 10^–4 subs/site/year) according to the latest report of evolution rate of rdrp gene in DCoV (Wang et al., 2021). Then, we ran a Markov chain of 10 million steps with sampling every 1,000 steps in BEAST v1.10.4 (Drummond et al., 2012). The mean evolution rate and time of the most recent common ancestor (tMRCA) were calculated under uncorrelated lognormal relaxed clock. The most appropriate substitute model of nucleotide was GTR + F + G4 calculated by ModelFinder according to BIC method. We checked the effective sample size (ESS) of parameters in Tracer v1.7 program and ensured that they all reached convergence (ESS > 200). Finally, the maximum clade credibility (MCC) tree was obtained by discarding the first 10% of states in Tree Annotator package.

Paired evolutionary trees of DCoVs and their corresponding hosts were used to examine the evolutionary relationship. Some DCoV strains without host records were excluded in this analysis. In addition to our samples, only one or two strains (from different countries) were reserved as the representative for these highly similar sequences with the same host. The host evolutionary tree and topologies were obtained from the TimeTree website¹ by entering the Latin text of all species to be analyzed (Kumar et al., 2017). Then, we constructed the viral phylogenetic tree based on the rdrp gene using Mrbayes v3.2.6 program (Ronquist et al., 2012). We run three hot and one cold Markov Chain Monte Carlo (MCMC) chains in Mrbayes, and the appropriate substitution and site-variation model of nucleotides was set to GTR + F + I + G4 which was tested and estimated by ModelFinder software according to BIC method. We performed two independent runs with total generations of 2,000,000 for each run and sampled every 100 generations. The first 10% burned samples were discarded, and the 50% majority rule tree was generated using the remaining samples. We checked the ESS value of each parameter in Tracer v1.6 program and ensured that they all reach convergence (ESS > 200).

According to previous research on the evolutionary history of RNA virus (Shi et al., 2018), we examined the co-divergence (coevolution) and cross-host transfer events among DCoVs and their hosts using eMPRess program (Santichaivekin et al., 2021). The consistency test of paired host–virus tree in eMPRess was conducted based on the duplication-transfer-loss (DTL) model, using a maximum parsimony approach (MPR) to find a “best” mapping of the virus tree onto the host tree which minimized the total “cost.” The null hypothesis is that the reconciliation score of the original dataset was due to change in eMPRess. The “cost” of every event was set as follows: co-divergence = 0 (fixed), duplication = 1, transfer (aka host switch) = 1.2, and loss = 1. The statistical test of co-divergence was conducted by comparing the estimated cost to null distributions calculated from 100 random permutations of the host tip mappings.

A total of 460 fecal samples were collected. Reverse-transcriptional PCR detection showed that three samples were positive of DCoV, with a total positive rate of 0.65%. These positive samples were collected in Cuihu Park, Kunming (Table 1). The three PCR products were sequenced, and based on the sequences, we constructed a phylogenetic tree together with other coronaviruses, preliminarily confirming that the sequences we detected could be clustered into the genus of deltacoronavirus (Figure 1). Especially, the RdRp sequences were highly similar to those of DCoVs HKU27, HKU28, and HKU29 detected in Dubai in 2018 (Lau et al., 2018), with identities of 99.290–99.645%.

In order to further reveal the genetic and evolutionary characteristics of the potential DCoVs detected in black-headed gulls, we amplified and sequenced the whole genome of the three DCoVs positives. As a result, three whole-genome sequences of 26,159 nt were obtained. Then, we annotated the genomes by referring to HKU29 genome (LC364344.1) and identified 10 typical open reading frames (ORFs) of coronavirus which were aligned as 5′ untranslated region (5′ UTR)-replicase ORF1ab-spike (S)-envelope (E)-membrane (M)-Nucleocapsid (N)-3′ UTR. Moreover, additional ORFs coding non-structural (NS) proteins (NS6, NS7a, NS7b, NS7c, and NS7d) were identified (Figure 2). The identities among the three genomes ranged from 99.382 to 99.957%, and they showed similarity of 96.606–96.911% to DCoVs HKU27, HKU28, and HKU29 (Table 2). Hence, we named the three DCoV strains in black-headed gulls as BHG-DCoV HNU4-1, HNU4-2, and HNU4-3.

The 3CL-pro, NiRAN, RdRp, ZBD, and HEL1 are the hallmark domains for classification of coronaviruses (International Committee on Taxonomy of Viruses Executive Committee et al., 2020). In order to confirm the taxonomic status of HNU4-1, HNU4-2, and HNU4-3, we compared the aa sequences of these domains with the known DCoVs White-eye coronavirus HKU16 (WECoV-HKU16), HKU27, HKU28, and HKU29. The result showed that all the eight viruses could be classified into the same species under the genus of deltacoronavirus (Table 2). Although the eight DCoVs were classified into the same species, the structural proteins of the BHG-DCoV, HKU27, HKU28, and HKU29 shared high similarity while WECoV-HKU16’s were quite distinct from the others (Table 2).

In order to further compare the differences among BHG-DCoV and HKU27, HKU28, and HKU29, we compared the amino acids encoded in their coding regions. The results showed that there were amino acid changes in ORF1ab. We used the cleavage sites of 16 NSPs of ORF1ab analyzed previously to predict the cleavage sites of ORF1ab of BHG-DCoV (Lau et al., 2018). The results showed that the mutations of ORF1ab mainly occurred in NSP3 and NSP6 (Figure 3). BHG-DCoV and HKU27, HKU28, and HKU29 also had different degrees of amino acid mutations in M, N, NS7a, NS7b, and NS7c. Among them, the most amino acid mutation was N protein, with a total of five amino acids mutated, and the least was NS7c, with only one amino acid mutation (Figure 3).

The difference of S protein among BHG-DCoV, HKU27, HKU28, HKU29, mammalian DCoV, and sparrow coronavirus HKU17-USA (close to PDCoV) was compared (Figure 4 and Supplementary Figure 2). The S protein was quite different between avian DCoV (except HKU17-USA) and mammalian DCoV, and the BHG-DCoV only shared about 45% aa identity with mammalian DCoV (Supplementary Table 1). Besides, there were multiple deletions/insertions and mutations in the NTD domain of S1 region, which may be a possible region for host receptor utilization restriction (Figure 4). However, the difference in the S2 region was relatively small (Supplementary Figure 2). We further analyzed the variation in BHG-DCoV compared with HKU27, HKU28, and HKU29 in S1 region. The result showed that there were nine sites with completely different amino acids between BHG-DCoV, HKU27, HKU28, and HKU29 in S1-NTD (D59A, N/S65Q, Q162L, G/S163I, S179E, T180K, L190F, Q207K, and A253R), and two sites in S1-CTD (Y355H and N419T). Besides, HKU28 was the closest BHG-DCoV on S protein and showed about 95.7% aa identity.

To assess the impact of mutations on S protein structure and function, we compared the differences in primary and tertiary structures of HNU4-1 and HKU28. In primary sequence, there were a total of 50 aa mutations and one aa insertion in S protein of HNU4-1 compared to HKU28 (Supplementary Table 2). Of these, 16 mutations were in S1-NTD and five mutations were in S1-CTD. We subsequently predicted the tertiary structures of S protein of HNU4-1 and HKU28 by homology modeling method, respectively (Figure 5). The tertiary structure of HKU28 covers the region from aa 51 to 1061, and the tertiary structure of HNU4-1 covers the region from aa 46 to 1062, the chain A as shown in Figure 5A. In general, the tertiary structure of the S protein of HNU4-1 and HKU28 was very similar. However, some mutation sites caused significant changes in residues conception, such as D59A, Q162L, S163I, S179E, T180K, Q207K, N256T, and P258N in S1-NTD, and Y355H, N419T, and K455Q in S1-CTD (Figure 5B). Notably, Q162L, S163I, S179E, T180K, and Q207K are at the top of the NTD, and Y355H, N419T, and K455Q are also at the top of the CTD; we speculate that it may be related to host receptor utilization. However, avian receptor utilization for DCoV is currently unclear, which greatly limits the in-depth exploration of the impact of mutations on function.

Based on the amino acid sequences of ORF1ab, S, HEL1, 3CL-pro, and NIRAN, we constructed the phylogenetic trees of BHG-DCoVs and the other 30 DCoVs. The results showed that HNU4-1, HNU4-2, HNU4-3, HKU27, HKU28, and HKU29 were clustered as a sub-branch within a branch including WECoV-HKU16 (Figure 6). Compared with HKU27, HKU28, and HKU29, the three BHG-DCoV are more closely related. It is worth noting that in the phylogenetic tree of ORF1ab, HEL1, 3CL-pro, and NIRAN, they were also clustered with GPCoV found in Antarctica. In the phylogenetic tree constructed by S fragment, except HKU27, HKU28, and HKU29, they also clustered with a Arenaria interpres DCoV (ARDCoV) found in Canada (Figure 6; Wille et al., 2021).

The estimated mean evolutionary rate of DCoV rdrp gene was about 1.49 × 10^–4 subs/site/year (9.10 × 10^–5–2.0 × 10^–4, 95% HPD). The time of the most recent common ancestor (tMRCA) differentiation time of DCoV was about 2977 B.C (6617 B.C–1035 B.C, 95% HPD). The tMRCA of White-eye coronavirus HKU16 species was about 1200 A.D (594 A.D–1508 A.D, 95% HPD). The tMRCA of HNU4-1, HNU4-2, and HNU4-3 was about 1964 A.D (1914 A.D–1993 A.D, 95% HPD). However, the tMRCA of HKU27, HKU28, and HKU29 was about 1952 A.D (1896 A.D–1983 A.D), and it suggests that they may differentiate earlier than BHG-DCoV (Figure 7).

The distribution of virus taxa on the host tree showed the host range of DCoV, including 18 species of the Aves class and two species of the Mammalia class (Figure 8 and Supplementary Table 3). The analysis using eMPRess of virus–host mapping detected a total of 23 events, including five co-divergence events, four duplication events, and 14 host-switching events (Figure 8). The global co-divergence test indicated that the co-divergence signal was weak in the host–virus tree, because the “cost” of the original host–virus tree was not significantly less than the randomness (p-value > 0.01) (Supplementary Figure 3). Frequent host-switching events indicate that DCoV was actively jumping among different hosts. Especially, during the history of evolution, HKU28 underwent several reliable host-switching events (with supporting values of 100) among which one led to the generation of HNU4 (Figure 8). This indicated that HNU4 might originate from HKU28 via cross-species transmission among different avian orders.

Besides, the DCoV in houbara (Gruiformesor) looks “restless” and spread to other hosts, such as falcon (Falco peregrinus), pigeon (Columba livia), and white-eye (Zosterops japonicus). Consistent with the previous speculation, the host transfer between houbara and falcons may be caused by the food chain (Lau et al., 2018). It is worth noting that the DCoV in white-rumped snowfinch (Montifringilla taczanowskii), sparrow, and munia shared co-divergence with three hosts. Then, the sparrow DCoV spread to common-magpie (Pica pica), quail (Coturnix coturnix), and mammal (porcine and Asian leopard cat). The munia DCoV may spread to magpie-robin (Copsychus saularis). In general, many host-switching events were predicted in DCoV, which imply that DCoV was very “active” in their hosts. However, the current sample of bird DCoV is still limited, and only a single host species has been detected in some orders. Therefore, more convincing evolutionary relationships require more samples to support in future. The large number of potential host-switching events indicates that the prevalence of DCoV in birds is wide and complex, and it is necessary to strengthen monitoring.