Detection of Coronaviruses and Genomic Characterization of Gammacoronaviruses from Overwintering Black-Headed Gulls (Chroicocephalus ridibundus) in Yunnan Province, China

1. Introduction

Since the 21st century, coronaviruses (CoVs) have caused three pandemics, namely severe acute respiratory syndrome (SARS) caused by SARS-CoV in 2002, Middle East respiratory syndrome (MERS) caused by MERS-CoV in 2012, and coronavirus disease 2019 (COVID-19) caused by SARS-CoV-2 in 2019, which have caused huge threats and losses to the lives and economies of people worldwide. These viruses are believed to originate from zoonotic viruses harbored in wild animals and then infect humans through accidental cross-species transmission[1,2,3,4,5,6,7,8,9].

Coronavirus (CoV) is a positive-sense, single-stranded, non-segmented RNA virus with typical spherical enveloped virions. Its viral glycoproteins (spikes) surround a helical nucleocapsid, and the whole genome is approximately 27-32 kb. The 2/3 region at the 5’ of the genome encodes two large open reading frames (ORFs): ORF1a and ORF1b, which are translated into two large polyproteins (PPs), PP1a and PP1b. The remaining 1/3 region of the genome encodes four major structural proteins: spike (S), envelope (E), membrane (M), and nucleocapsid (N). Some accessory proteins are interspersed among the structural proteins[10,11]. CoVs belong to the order Nidovirales, family Coronaviridae, subfamily Orthocoronavirinae[12]. The International Committee on Taxonomy of Viruses (ICTV) divides the subfamily Orthocoronavirinae into four genera, Alphacoronavirus (α-CoV), Betacoronavirus (β-CoV), Deltacoronavirus (δ-CoV)and Gammacoronavirus (γ-CoV), based on the five conserved protein domains used for classification in the order Nidovirales: 3C-like protease (3CLpro), RNA-dependent RNA polymerase (RdRp), Nidovirus RdRp-Associated Nucleotidyltransferase (NiRAN), zinc-binding domain (ZBD), and superfamily 1 helicase (HEL1)[13]. The genera α-CoV and β-CoV mainly infect mammals, while the genera δ-CoV and γ-CoV mainly infect poultry and birds, and a small number also infect mammals[14].

Infectious bronchitis virus (IBV), a representative CoV of the genus Gammacoronavirus, was the first reported in the United States in 1931. IBV is highly contagious and mainly infects the respiratory tract, kidneys, and reproductive system of poultry, leading to respiratory distress, kidney damage, and decreased egg production, which causes huge economic losses to the poultry industry. Since the discovery of IBV, it had been the only known γ-CoV for more than 50 years. As time passed, the viral diversity and host diversity of the genus Gammacoronavirus have gradually been discovered[15,16,17,18,19,20,21,22,23].

Black-headed gull (Chroicocephalus ridibundus) belongs to the family Laridae of the order Charadriiformes, is a migratory bird, widely distributed in North America, Eurasia, northern, central and southern Africa, the Indian subcontinent, the Indo-China Peninsula, the southwest and southeast coastal areas of China, and also on the islands in the Pacific Ocean[24]. In China, Xinjiang, Inner Mongolia and Heilongjiang are its main breeding areas. In winter, the black-headed gull migrates from Siberia to the south of China for wintering, and Yunnan province is a typical wintering habitat for it. Black-headed gulls are hosts for some zoonotic viruses, and as migratory birds they have the ability to fly long distances, and they are abundant, gregarious and in very close contact with humans and other animals. It is due to these characteristics that black-headed gull play an important role in the transmission of pathogens by influencing their dynamics and thus the ecology and evolution of various viruses and bacteria[14,19,25].

Previous studies have found that the black-headed gull is a natural host for δ-CoV and γ-CoV [26,27]. We also conducted an investigation on the CoVs of the black-headed gull and analyzed the characteristics of the three whole genomes (HNU4-1, HNU4-2, HNU4-3) of black-headed gull deltacoronaviruses (BHG-DCoVs)[28]. In this study, we collected fecal samples from black-headed gulls at three sampling sites in Kunming City, Yunnan Province, China, and detected the presence of δ-CoV and γ-CoV using RT-PCR. For the positive samples of black-headed gull gammacoronaviruses (BHG-GCoVs), we performed whole-genome amplification to obtain the whole genomes of the two strains of BHG-GCoVs and were bioinformatically analyzed, suggesting that they may have the potential for cross-species transmission.

2. Results

2.1. Detection of CoVs in the Feces of Black-Headed Gulls

A total of 59 CoVs were detected from the 509 collected fecal samples of black-headed gulls in three sites (Cuihu Park, Daguanlou Park and Haigeng Dam) , with a prevalence of 11.59% (59/509, 95%CI: 8.80% - 14.38%). Among them, the prevalence of BHG-DCoVs was 3.54% (18/509, 95%CI: 1.90% - 5.10%), and that of BHG-GCoVs was 8.06% (41/509, 95%CI: 5.68% - 10.43%). The prevalence of BHG-GCoVs was higher than that of BHG-DCoVs (χ² = 9.518, P < 0.01) (Table 1).

We constructed a phylogenetic tree for the 59 detected BHG-CoVs based on the partial RdRp (δ-CoVs about 435bp, γ-CoV about 393bp) sequencing results, using the neighbor-joining method (Figure 1). The results showed that 41 strains of BHG-GCoVs (GenBank accession numbers: PQ676674, PQ676677- PQ676682, PQ676685- PQ676687, PQ676689,PQ676690, PQ676692, PQ676694, PQ676695, PQ676697- PQ676701, PQ676704, PQ676706, PQ676710-PQ676715, PQ676717- PQ676720, PQ676723- PQ676729, PQ676731,PQ676732) had a relatively close genetic relationship with the ruddy turnstone CoV (MT993597) from Australia. The partial RdRp nucleotide similarity among all the detected BHG-GCoVs in this study ranged from 97.96% to 100%. In addition, BHG-DCoVs formed two subclades. Among them, 10 strains (GenBank accession numbers: PQ676675, PQ676676, PQ676683, PQ676684, PQ676688, PQ676691, PQ676693, PQ676702, PQ676707, PQ676716) formed a separate branch, and they seemed to be the ancestors of HKU27 (LC364342), HKU28 (LC364343), and HKU29 (LC364344). The other 8 strains (GenBank accession numbers: PQ676696, PQ676703, PQ676705, PQ676708, PQ676709, PQ676721, PQ676722, PQ676730) had a relatively close genetic relationship with HNU4-1 (OL311150), HNU4-2 (OL311151), and HNU4-3 (OL311152). The partial RdRp nucleotide similarity among all the detected BHG-DCoVs in this study ranged from 93.79% to 100%. Host molecular identification was carried out on all the positive samples through cytb gene amplification, and it was determined that the host were both the black-headed gulls.

2.2. Whole-Genome Structural Analysis of BHG-GCoVs

To further reveal the genetic and evolutionary characteristics of the GCoVs detected in the feces of black-headed gulls, we performed whole-genome amplification and sequencing on two GCoVs-positive samples, which resulted in the acquisition of two whole-genome sequences of BHG-GCoVs with 27358 nt and 37355 nt respectively, and the nucleotide similarity of the two whole genomes was 98.75%, which were designated as DLU1 (GenBank accession number PQ676672) and DLU2 (GenBank accession number PQ676673).

Genome annotation of 2 strains of BHG-GCoVs with reference to AvXc-GCoV (PP845453) from the terek sandpiper (Xenus cinereus) in China identified the genomic order of DLU1 and DLU2 as 5’UTR-ORF1ab-spike (S) -envelope (E) -membrane (M) -ORF5a-ORF5b- 5a-5b-Nucleocapsid (N) -3’UTR. Comparison with AvXc-GCoV, Shelduck_GCoV, and Duck_GCoV 2714 showed that both DLU1 and DLU2 contained ORF1a/b, S, E, M, and N genes typical of the orthocoronavirus subfamily. In addition, all five strains contained the 5a and 5b accessory proteins, but the BHG-GCoVs lacked the ORF9a and ORF9b putative proteins relative to AvXc-GCoV, and the 3a and 3b accessory proteins relative to Shelduck_GCoV and Duck_GCoV 2714 (Figure 2).

2.3. Genome Similarity Analysis of BHG-GCoVs.

We performed BLAST on the NCBI based on the whole genomes of BHG-GCoVs and selected 21 strains of GCoVs (Supplementary Table 4) that were significant for comparison with BHG-GCoVs (https://blast.ncbi.nlm.nih.gov/Blast.cgi. accessed on 19 December 2024). Among them, 13 strains were the whole genomes recently uploaded by Chinese scientists. They came from 13 different hosts, including 3 families of 2 orders, mainly the Scolopacidae and Anatidae. Another 3 strains were GCoVs detected in Australia and Poland, and their whole genomes or relatively long genomes had been obtained. The other 5 strains were representative strains of the genus Gammacoronavirus.

The whole genomes of BHG-GCoVs have a high nucleotide similarity with the GCoV (AvXc-GCoV) from Xenus cinereus and the GCoV (AvNp-GCoV) from whimbrels (Numenius phaeopus), ranging from 94.29% to 94.70%. However, there is a significant difference in nucleotide similarity with Duck_GCoV 2714, the representative strain of Gammacoronavirus anatis, which is only 74.38% (Supplementary Table 1).

Five conserved replicase domains in non-structural proteins(3CLpro, NiRAN, RdRp, ZBD, and HEL1) are used for the classification of coronaviruses[29]. To confirm the classification, we compared the nucleotide (nt) and amino acid (aa) similarities of these five domains between BHG-GCoVs and the selected 21 strains of GCoVs. The results showed that the average aa similarity of the five domains of BHG-GCoVs was relatively high with AvXc-GCoV, AvNp-GCoV, common gull GCoV, and ruddy turnstone GCoV, ranging from 99.14% to 99.53%. The aa similarity with Duck_GCoV 2714 was 90.75%. However, the aa similarity with the other four representative strains was only 67.68% - 83.42% (Supplementary Table 1).

Among the structural proteins, the gene of the S protein is the most prone to mutations in CoVs. The S protein of BHG-GCoVs contains 1,123 aa and has a unique S1/S2 cleavage site (520-FTPEKQVLDIKRFTSNYTEYAPMVFSD-546)[27], forming two subunits, S1 (546 aa) and S2 (577 aa). The results of the similarity analysis of the S protein showed that the aa similarity between BHG-GCoVs and the GCoV of ruddy turnstone (Arenaria interpres) (AvAi-GCoV) was relatively high, ranging from 96.26% to 96.44%. In the S1 subunit, there were 37 aa differences between BHG-GCoVs and AvAi-GCoV, with a similarity of 93.22%. In the S2 subunit, there were 2 aa differences between BHG-GCoVs and AvAi-GCoV, with a similarity of 99.65%. The tertiary structure of the S protein of DLU1 and AvAi-GCoV is shown in Figure 3(A). Interestingly, the aa of the S protein of BHG-GCoVs had a relatively high similarity with GCoVs from Charadriiformes (AvXc-GCoV, AvNp-GCoV, AvAi-GCoV, AvCt-GCoV, AvCa-GCoV, excepting for the common gull GCoV), ranging from 90.21% to 96.44%, while the similarity with GCoVs from Anseriformes (AvAc-GCoV, AvAa-GCoV, AvAz-GCoV, AvMs-GCoV, Shelduck_CoV, Duck_GCoV 2714) was only 31.15% - 54.81%. It is worth noting that although both the common gull and the black-headed gull belong to the family Laridae, their aa similarity of the S protein of GCoVs was only 54.51% - 54.59% (Supplementary Table 2).

In the S1 subunit, the aa similarity between BHG-GCoVs and GCoVs from Charadriiformes (AvXc-GCoV, AvNp-GCoV, AvAi-GCoV, AvCt-GCoV, AvCa-GCoV, excepting for the common gull GCoV) ranged from 85.74% to 93.05%. The aa similarity between BHG-GCoVs and that from Anseriformes (AvAc-GCoV, AvAa-GCoV, AvAz-GCoV, AvMs-GCoV, Shelduck_CoV, Duck_GCoV 2714) was only 16.15% - 39.65%, and the similarity that with the common gull GCoV was 37.63% - 37.81%.In the S2 subunit, the aa similarity between BHG-GCoVs and that from Charadriiformes (AvXc-GCoV, AvNp-GCoV, AvAi-GCoV, AvCt-GCoV, AvCa-GCoV, excepting for the common gull GCoV) ranged from 94.45% to 99.83%. The aa similarity that with Anseriformes (AvAc-GCoV, AvAa-GCoV, AvAz-GCoV, AvMs-GCoV, Shelduck_CoV, Duck_GCoV 2714) was only 47.19% - 69.26%, and the similarity that with the common gull GCoV was 70.12% (results not shown).

Figure 3. (A) Tertiary structure of the S protein of DLU1 and AvAi-GCoV. Green represents the S1 subunit, yellow represents the S2 subunit, and red represents the amino acid differences in the S proteins of DLU1 and AvAi-GCoV. (B) Multiple alignment of the amino acids of the S protein of DLU1, AvAi-GCoV and common gull GCoV. The red box represents the amino acid differences between DLU1 and AvAi-GCoV.

Figure 3. (A) Tertiary structure of the S protein of DLU1 and AvAi-GCoV. Green represents the S1 subunit, yellow represents the S2 subunit, and red represents the amino acid differences in the S proteins of DLU1 and AvAi-GCoV. (B) Multiple alignment of the amino acids of the S protein of DLU1, AvAi-GCoV and common gull GCoV. The red box represents the amino acid differences between DLU1 and AvAi-GCoV.

2.4. Phylogenetic Analysis of the Whole Genome BHG-GCoVs

Based on the whole-genome nucleotide sequences, a phylogenetic tree was constructed using the maximum likelihood method, with GTR+G+I as the best substitution model, and with 1,000 bootstrap replicates with other CoVs . The phylogenetic tree indicated that the two BHG-GCoVs clustered into the species Gammacoronavirus anatis in the subgenus Igacovirus of the genus Gammacoronavirus, which is represented by the strain Duck CoV_2714 (NC_048214). BHG-GCoVs clustered into a sub-branch with AvXc-GCoV and AvNp-GCoV that were recently detected in Shanghai, China (Figure 4).

2.5. Phylogenetic Evolution of the Domains of BHG-GCoVs

We constructed phylogenetic trees with other GCoVs based on the aa sequences of the five conserved replicase domains (3CLpro, NiRAN, RdRp, ZBD, HEL1) in non-structural proteins and the S protein. The results showed that BHG-GCoVs clustered in the same sub-branch with AvXc-GCoV and AvNp-GCoV in the 3CLpro, ZBD, and HEL1 domains, and clustered in the same sub-branch with common gull GCoV in the NiRAN and RdRp domains. Additionally, it is worth noting that in terms of the S protein, BHG-GCoVs clustered into the same sub-branch with AvAi-GCoV (Figure 5).

2.6. Recombination Analysis of BHG-GCoVs

We found that 2 strains of BHG-GCoVs changed their positions in the phylogenetic trees based on different domains, indicating that recombination events might occur frequently during the process of virus evolution. Screening for potential recombination events was carried out on the whole-genome sequences of the two BHG-GCoVs in this study. The best possible recombination events in DLU1 and DLU2 were determined by the GENECONV (p-value: 8.66×10⁻²⁴³) and RDP (p-value: 2.035×10⁻²⁵¹) methods provided by the RDP v4.100 program. These recombination events with strong p-values were further confirmed in the SimPlot v3.5.1 program. The results showed that DLU1 and DLU2 were most likely to be recombinant strains of common gull GCoV and AvAi-GCoV, and there were breakpoints at genomic positions 20,156 nt and 24,432 nt, which were close to the start and end regions of the S protein (Figure 6).

2.7. Evolutionary Rate and tMRCA

The average evolutionary rate of CoVs in the RdRp gene was estimated to be 1.3 × 10^-4 nucleotide substitutions per site per year using Bayesian Skyline under a relaxed molecular clock model with an uncorrelated lognormal distribution [14]. Molecular clock analysis using the RdRp gene showed that the most recent common ancestor (tMRCA) of avian γ-CoV and beluga whale γ-CoV diverged at approximately 3934 B.C (95% highest posterior density (HPD): 9437 B.C-542 B.C), and that the subgenus Brangacovirus and Igacovirus tMRCA differentiation time is approximately 246B.C (95% HPD: 2244 B.C-887 AD), for subgenus Igacovirus it is approximately 788 AD (95% HPD : 153B.C-1323AD), for AvXc-GCoV the tMRCA differentiation time was approximately 1984 AD (95% HPD: 1950AD-2006AD), and that of BHG-GCoVs was approximately 1994AD (95% HPD: 1966AD-2011AD) (Figure 7).

3. Discussion

Since the CoVs that can infect humans are mainly α-CoV and β-CoV, the current surveillance of CoVs mainly focuses on α-CoV and β-CoV as well as their host animals, such as bats and rodents[30,31,32]. However, on a global scale, birds are also important hosts for some zoonotic viruses. The wide variety and abundance of bird species, the ability of some birds to fly long distances, and their close relationships with humans and domestic animals also play an important role in the transmission of pathogens, making it essential to monitor pathogens carried by birds. Although there have been no reports so far on the cases of CoVs directly infecting humans from poultry and wild birds, the case of human infection by porcine deltacoronavirus (PDCoV) in Haiti reported by US scientists in 2021, indicating that δ-CoVs and γ-CoVs , which mainly infect poultry and birds, also have the possibility of cross-species transmission to livestock and then infecting humans[33]. This has aroused our vigilance.

Black-headed gulls are the natural hosts of δ-CoVs and γ-CoVs. The present study was carried out on CoVs in black-headed gull, which seems to have a higher prevalence than in other areas, for example, the prevalence of BHG-CoVs was 8.20% (5/61) in the Bering Strait and the prevalence of BHG-CoVs was 4.53% (35/773) in Poland. The prevalence of γ-CoVs was higher than that of δ-CoVs in this study, as in other studies [22,26].

In this study, we obtained CoV-positive samples by detection and selected two BHG-GCoVs derived from black-headed gulls for whole-genome sequencing, obtained the whole genomes of the two BHG-GCoVs. As far as we know, the whole genome of BHG-GCoV has not been reported, DLU1 and DLU2 are the two whole genomes of BHG-GCoVs obtained so far.

Based on the BLAST of the whole genomes of BHG-GCoVs on NCBI, we found that BHG-GCoVs belong to the subgenus Igacovirus of the genus Gammacoronavirus. The subgenus Igacovirus is divided into three species, namely Gammacoronavirus anatis, Gammacoronavirus galli, and Gammacoronavirus pulli (https://ictv.global/taxonomy, accessed on 19 December 2024). The first two species both belong to infectious bronchitis virus (IBV), while BHG-GCoVs have a higher similarity with the species Gammacoronavirus anatis. Previously, the number of samples with whole-genome sequences or relatively long sequences obtained in the species Gammacoronavirus anatis was very limited, mainly coming from Poland, Australia, and China[27,34,35,36,37]. Recently, Chinese scientists have uploaded some whole-genome sequences belonging to the species Gammacoronavirus anatis from 13 different hosts of 3 families of 2 orders, which have provided an important reference basis for our research on the cross-species transmission of BHG-GCoVs. Based on the whole-genome phylogenetic analysis, the results showed that BHG-GCoVs belong to the species Gammacoronavirus anatis of the subgenus Igacovirus, and BHG-GCoVs clustered into the same sub-branch with AvXc-GCoV and AvNp-GCoV, which indicated that they have a relatively close genetic relationship. However, they have a relatively distant genetic relationship with the representative strain of Duck_GCoV 2714.

Based on the whole genome and five domain sequences of BHG-GCoVs, similarity was compared with other GCoVs. From the analysis results, the whole genomes of BHG-GCoVs were similar in length to those of other GCoVs from Charadriiformes, which was about 27,300 nt in length, while the whole genomes of GCoVs from Anseriformes were longer than 28,000 nt in length (Supplementary Table 1). Moreover, we also found that GCoVs carried by Charadriformes were more similar to BHG-GCoVs than those carried by Anseriformes. Therefore, it can be seen that GCoVs have a higher similarities within hosts from the same order. Similarly, the similarity of the S gene of BHG-GCoVs to that of other GCoVs from Charadriiformes was significantly higher than that of GCoVs from Anseriformes. Additionally, it is worth noting that although both the common gull and the black-headed gull belong to the family Laridae, the differences in the S protein of the GCoVs they carry are relatively large, only 54.51% - 54.59%, whether this is due to species differences or regional differences needs to be confirmed by more investigations.

The S protein of coronavirus plays an important role in binding to receptors on host cells and the fusion of the viral membrane with the host cell membrane. During the maturation of CoVs, the extracellular region of the S protein is recognized and cleaved into two subunits, S1 and S2, by host proteases[38]. The S1 subunit is mainly responsible for binding to cell receptors for virus adsorption, and the S2 subunit mainly mediates membrane fusion with the host cell and the internalization of the virus[39].Further analysis of the S protein of BHG-GCoVs found that the mutated regions of the S protein of BHG-GCoVs are mainly in the S1 subunit, while the S2 subunit is relatively more conserved, which is consistent with the previous research results of Polish scientists[27]. We speculate that it is due to the wide variety of cell types and receptors in different hosts that CoVs need to constantly mutate to attach to new receptors, and the virus internalizes in order to replicate and survive. Overall, the tertiary structures of the S proteins of DLU1 and AvAi-GCoV are very similar. However, some mutation sites cause significant residue variations.

Interestingly, the positions of the two strains of BHG-GCoVs changed in the phylogenetic tree based on different genes, especially the S protein, and so we hypothesized that a recombination event might have occurred. We obtained confirmation by further analysis that DLU1 and DLU2 were recombinant strains of common gull GCoV and AvAi-GCoV. In addition, estimates based on temporal divergence indicated that the tMRCA of GCoVs from Scolopacidae preceded that of BHG-GCoVs, and therefore, it is likely that BHG-GCoVs evolved through GCoVs from Scolopacidae and common gull GCoVs.

Health threats from CoVs are persistent and long-term, and γ-CoVs have more genetic diversity and host diversity, but our knowledge of them is still very limited. Long-term surveillance of CoVs carried by birds and poultry is necessary to prevent future emerging infectious diseases.

4. Materials and Methods

4.1. Sample Collection

In January 2023, a total of 509 fecal samples of black-headed gulls were collected from three sampling sites (Cuihu Park, Daguanlou Park, and Haigeng Dam) in Kunming City, Yunnan Province, China. All the samples were placed in cryopreservation tubes filled with 1 ml of virus transport medium (VTM). After collection, the cryopreservation tubes were immediately put into the car refrigerator for storage and then transported to the laboratory and stored at -80 °C until use.

4.2. RNA Extraction and CoV Screening by RT- PCR

Approximately 70 μl of viral RNA was extracted from the fecal samples using the MagaBio plus Viral DNA/RNA Purification Kit III (Bioer, Hangzhou, China) through an automatic nucleic acid extractor (Bioer, Hangzhou, China) according to the manufacturer’s instructions, and then stored at −80 °C as a template for reverse transcription-polymerase chain reaction (RT-PCR). CoVs were initially screened by nested RT-PCR. For the first round, the forward primer AvCoV-F1 (5’-GGKTGGGAYTAYCCKAARTG-3’) and the reverse primer AvCoV-R1 (5’-TGYTGTSWRCARAAYTCRTG-3’) were used[23]. The first round of PCR amplification was carried out using the HiScript II One Step RT-PCR Kit (Vazyme, Nanjing, China). The 25 μl PCR mixture included 12.5 μl of 2 × One Step Mix (Dye Plus), 1.25 μl of One Step Enzyme Mix, 1 μl each of the forward and reverse primers, 6.25 μl of ddH₂O, and 3 μl of the extracted RNA template. The mixture was subjected to reverse transcription at 50 °C for 30 min, pre-denaturation at 94 °C for 3 min, denaturation at 94 °C for 30s, annealing at 48.0 °C for 30s, and extension at 72 °C for 36s in an automatic thermal cycler (Applied Bio-systems, Shanghai, China). The above three steps of denaturation, annealing, and extension were repeated for 35 cycles, and finally, an extension was carried out at 72 °C for 7 min. This was performed to amplify a 602 bp fragment of the RdRp gene of CoV. For the second round, the forward primer AvCoV-F2 (5’-GGTTGGGACTATCCTAAGTGTGA-3’) and the reverse primer AvCoV-R2 (5’-CCATCATCAGATAGAATCATCAT-3’) were used[23]. The second round of PCR amplification was carried out using 2 × Rapid Taq Master Mix (Vazyme, Nanjing, China). The 25 μl PCR mixture included 12.5 μl of 2 × Rapid Taq Master Mix, 1 μl each of the forward and reverse primers, 9.5 μl of ddH₂O, and 1 μl of the product from the first round of PCR. The mixture was subjected to pre-denaturation at 95 °C for 3 min, denaturation at 95 °C for 15s, annealing at 48.5 °C for 15s, and extension at 72 °C for 26s. The above three steps of denaturation, annealing, and extension were repeated for 35 cycles, and finally, an extension was carried out at 72 °C for 5 min. This was done to amplify a 440 bp fragment of the RdRp gene of CoV.

The PCR products were gel-purified (OMEGA Bio-tek, Norcross, GA, USA) and sent to Sangon Biotech for bi-directional sequence determination[32]. The sequences of the PCR products were compared with known sequences of the RdRp genes of CoVs in the GenBank database.

To confirm the host species, the forward primer L0 (5’-GGACAAATATCATTCTGAGG-3) and the reverse primer H0 (5-GGGTGTTCTACTGGTTGGCTTCC-3’) were used to amplify and sequence the mitochondrial cytb gene of the samples[40].

4.3. Complete Genome Sequencing

Two complete genomes of BHG-GCoVs 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 onestep 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, Basel, Switzerland), respectively. The expected size of PCR products was purified by gel and directly sequenced. The sequence was assembled to obtain the whole genome sequence[28].

4.4. Genomic and Phylogenetic Analysis

The whole genome was annotated through Geneious Prime software. Multiple sequence alignments with other CoVs were conducted using MAFFT (v7.490)[41]. The most appropriate aa substitution model was calculated by the Find Best DNA/Protein Models tool in the MEGA X software, and the model with the lowest Bayesian Information Criterion (BIC) score was taken as the best substitution model[42].

The phylogenetic tree was constructed based on the whole genomes of CoVs using the maximum likelihood method. The phylogenetic trees for other non-structural and structural protein sequences were constructed using the neighbor-joining method, and the bootstrap replicates were all set to 1000 times. The phylogenetic trees were visualized using the iTOL v6 online tool[43].

The similarity comparison of the genomes was carried out using BioAider v1.334 software[44]. The recombination events of BHG-GCoVs were detected using the RDP package v.4 and Simplot v3.5.1[45].

4.5. Protein Tertiary Structure Analysis

Use AlphaFold3 to predict and model the tertiary structure of the S protein (https://alphafold3.org/), and use ChimeraX v1.8 software to visualize the model results.

4.6. Estimation of Divergence Dates

The RdRp gene was aligned using MAFFT program with codon method in BioAider v1.334. We detected the temporal structure in these RdRp gene sequences by TreeTime program. The correlation coefficient of R² was 0.01, and there was very weak signal between sampling time and genetic distance in the data. 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[28]. Then, we ran a Markov chain of 10 million steps with sampling every 1,000 steps in BEAST v1.10.4[46]. 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.

4.7. Statistical Analysis

Use IBM SPSS Statistics 25 software for statistical analysis. Use the chi-square (χ²) test to calculate the significant differences in prevalence.