Genetic Diversity and Phylodynamics of Avian Coronaviruses in Egyptian Wild Birds

Avian coronaviruses (ACoVs) are continuously evolving and causing serious economic consequences in the poultry industry and around the globe. Owing to their extensive genetic diversity and high mutation rates, controlling ACoVs has become a challenge. In this context, the potential contribution of wild birds in the disease dynamics, especially in domesticated birds, remains largely unknown. In the present study, five hundred fifty-seven (n = 557) cloacal/fecal swabs were collected from four different wild bird species from eight Egyptian governorates during 2016 and a total of fourteen positive isolates were used for phylodynamics and evolutionary analysis. Genetic relatedness based on spike (S1) gene demonstrated the clustering of majority of these isolates where nine isolates grouped within Egy/variant 2 (IS/885 genotype) and five isolates clustered within Egy/variant 1 (IS/1494/06 genotype). Interestingly, these isolates showed noticeable genetic diversity and were clustered distal to the previously characterized Egy/variant 1 and Egy/variant 2 in Egyptian commercial poultry. The S1 gene based comparison of nucleotide identity percentages revealed that all fourteen isolates reported in this study were genetically related to the variant GI-23 lineage with 92–100% identity. Taken together, our results demonstrate that ACoVs are circulating in Egyptian wild birds and highlight their possible contributions in the disease dynamics. The study also proposes that regular monitoring of the ACoVs in wild birds is required to effectively assess the role of wild birds in disease spread, and the emergence of ACoVs strains in the country.


Introduction
Coronaviruses are classified into three well-accepted genera; alpha-, beta-, and gamma-coronavirus whereas deltacoronavirus [1,2] is proposed as an additional group in the family. Alphacoronaviruses and betacoronaviruses have been isolated from mammals, whereas gammacoronaviruses, represented by the infectious bronchitis virus (IBV), are detected primarily in birds including poultry [2].

Polymerase (RdRp) and Spike (S1) Genes Amplification and Sequencing
Viral RNA was extracted from the allantoic fluids using Trizol LS®reagent kit (ThermoFisher, Life Technologies Corporation|Carlsbad, CA 92008, USA) according to the manufacturer's instructions. The eluted RNA in 30 µL of nuclease free water was stored at (−70 • C) until processing for the amplification. A single step RT-PCR was performed using the Verso One Step RT-PCR kit (ThermoFisher, Life Technologies Corporation|Carlsbad, California, USA) using previously reported primers for amplification of partial polymerase (RdRp) and full length S1 gene [6,21]. PCR products were gel excised and purified individually using QIAEX ® Gel Extraction Kit (Qiagen, Hilden, Germany). Gene sequencing was carried out using a BigDye Terminator v3.1 cycle sequencing kit (Qiagen, Hilden, Germany) in an ABI PRISM®300 using sequencing primers as described previously [22].

Sequence and Phylogenetic Analysis
The obtained sequences for RdRp (MF034358.1-MF034371.1) and S1 genes (MF034369.1-MF034385.1) were quality checked, assembled, edited using BioEdit software version 7.0.4.1 [23] and submitted to GenBank using BankIt tool of the GenBank (http://www.ncbi.nlm.nih.gov/WebSub/ ?tool=genbank). Analysis of full length S1 protein was conducted to determine the presence of potential N-linked glycosylation sites among the isolated Egyptian ACoV strains using BioEdit software version 7.0.4.1 [23]. To elucidate the phylogenetic relationships and high level clustering pattern based on full length S1 gene between ACoVs reported here and previously characterized from Middle East including Egypt and other parts of the world; we first compiled a dataset of all available S1 gene sequences. Then, this dataset was aligned in BioEdit version 7.0.1 [23] using Clustal W algorithm and was trimmed to equal lengths. To avoid sequence length biasness in evolutionary incursion, all sequences that aligned poorly or with incomplete information were excluded from the analysis. A phylogenetic tree was constructed using the nucleotide sequences of the S1 genes. For maximum-likelihood analysis of the phylogenetic relationship, a best-fit model was chosen on which further calculations and an ultrafast bootstrap equivalent analysis were based. IQ-tree software version 1.1.3 was used for these operations [24].

Selection Pressure Analysis
To explore the overall differences in selection pressure on the S1 gene, especially on the hypervariable regions (HVRs) sequences that define the cross-neutralization and escape mutant ACoV strains, we analyzed the occurrences of synonymous (dS) and non-synonymous (dN) substitutions using SNAP web tool (https://www.hiv.lanl.gov/content/sequence/SNAP/SNAP.html) [25], which plots the cumulative and per-codon occurrence of each substitutions. Six methods (RDP, GeneConv, BootScan, MaxChi, Chimaera and SiScan) integrated in the RDP v3 program [26] were applied on the avian CoVs sequences reported here to estimate any recombination event and to detect any putative recombination breakpoint using the nucleotide alignment of the S1 gene sequences. These methods were applied using the following parameters: window size = 20, highest acceptable p-value = 0.001 and Bonferroni correction. For reliable results, any putative recombination events detected by more than one method were considered.

Viral Sequences and Genomics Analysis
Comparison of nucleotide identity percentages, based up on S1 gene revealed that the fourteen isolates were 92-100% identity percentage (Table 3). In 2001, an Egyptian variant I strain (Egypt/Beni-Suef/01), closely related to the Israeli variant strain, was reported in different poultry farms [18]. Another variant was also identified in 2011 [27] and designated as Egyptian variant II. Deduced amino acids analysis was conducted to establish the genetic spectrum, origin, evolution, and the mutation trend analysis for Egyptian ACoVs. Hypervariable regions (HVRs) in the S1 gene showed distinct patterns among different viruses compared to the H120 and Ma5 strains (Table 4). In the analyzed sequence of the fourteen viruses in this study, amino acid mutations associated with virus tropism were identified at positions 38 (38N in only one virus while 38T in thirteen viruses) and 69 (69S in seven viruses and 69A in seven viruses) in the S1 protein of the variant strains. In addition, seventeen potential N-linked glycosylation sites in most of the Egyptian ACoV strains were found, with nine NXT and eight NXS sites in all the fourteen isolates reported in this study.

Phylogenetic Analysis
Based on the amplified CoV polymerase (RdRp) gene fragments, the phylogenetic analyses confirmed that the all detected viruses in this study belonged to gamma-coronaviruses ( Figure 2) along with other ACoVs. Phylogenomic analysis based on S1 gene of the ACoVs collected from Egypt and rest of the world, indicated that Egyptian viruses reported here were clustered with isolates from Libya, Oman and Kurdistan and shared ancestors (

Selective Pressure Sites and Recombination Analysis
A pairwise comparison bioinformatics approach (SNAP) was applied to determine the synonymous and non-synonymous substitution rates and selective evolutionary pressure for the S1 protein. The results presented in Figure 5 indicate the selection trends where numbers above zero indicate the positive selection, around zero shows the neutral selection, and below zero indicates the negative or purifying selection. The selection profiles of the amino acid sequence of all fourteen Egyptian ACoV strains showed different patterns within the S1 protein. Analysis of subtract synonymous and non-synonymous substitution rates in the isolates also demonstrated the presence of purifying (negative) selection. The cumulative difference between the non-synonymous substitution rate (dN) and the synonymous substitution rate (dS) (i.e., dN-dS) revealed that positive selection along the spike (S1) protein while negative and neutral selection was also detected ( Figure 5). When the S1 genes of ACoV strains were used as the query, a number of possible recombination sites were identified. Recombination analysis for the aligned sequences under study report the detection of possible recombination events; (1) AY091552.

Discussion
Continuous generations of natural mutant viruses or quasispecies subpopulations may result in the adaptation of such viruses to evade the host immune system and to gain higher transmissibility to new hosts [28]. ACoVs are continuously evolving, and the variant strains are one of the major emerging problems in the Egyptian poultry industry for the past five years and have resulted in significant economic losses. Furthermore, the infection predisposes the host for secondary bacterial infections resulting in an even higher morbidity and mortality rates [29]. Different serotypes have been reported worldwide and new variant serotypes are continuously being recognized [4]. However, there is no clear classification criterion for the circulating ACoVs, worldwide, which possesses complications in establishing epidemiological linking. Previous molecular studies revealed that the S1 subunit of S gene of ACoV is responsible for determining its serotypes [30] and new ACoV genotypes could arise as a result of substitutions in the amino acids of the S1 subunit [3]. S1 subunit has three HVRs, located within amino acids 38-67, 91-141, and 274-387 [3][4][5] and these HVRs were considered an essential determinant of coronavirus serotype specificity [31].

Discussion
Continuous generations of natural mutant viruses or quasispecies subpopulations may result in the adaptation of such viruses to evade the host immune system and to gain higher transmissibility to new hosts [28]. ACoVs are continuously evolving, and the variant strains are one of the major emerging problems in the Egyptian poultry industry for the past five years and have resulted in significant economic losses. Furthermore, the infection predisposes the host for secondary bacterial infections resulting in an even higher morbidity and mortality rates [29]. Different serotypes have been reported worldwide and new variant serotypes are continuously being recognized [4]. However, there is no clear classification criterion for the circulating ACoVs, worldwide, which possesses complications in establishing epidemiological linking. Previous molecular studies revealed that the S1 subunit of S gene of ACoV is responsible for determining its serotypes [30] and new ACoV genotypes could arise as a result of substitutions in the amino acids of the S1 subunit [3]. S1 subunit has three HVRs, located within amino acids 38-67, 91-141, and 274-387 [3][4][5] and these HVRs were considered an essential determinant of coronavirus serotype specificity [31].
Factors that contribute in the evolution and emergence of ACoV in commercial poultry are largely unknown. Evidences are that wild birds might constitute transmission vectors for ACoVs even in a huge geographic territory. It has also been proposed that when an ACoV infects the same wild bird host, this could lead to the emergence of new variants of ACoV [32] and the possibility of recombination between different ACoV strains. To assess the ACoV existence and carrying-potential of Egyptian wild birds for the first time in the country, 557 cloacal swabs were randomly collected from different wild bird species in different Egyptian governorates. Detection of gamma-coronaviruses (ACoVs) in wild bird samples from geographically separated populations highlights the potential of wild birds in the possible dissemination of ACoV. Since most of the ACoV positive wild birds were migratory and that Egypt has been proposed to be an important gateway in the Middle East and North of Africa for most of the emerging infectious diseases, it is likely that the mixing of birds and viruses from different hemispheres can lead to disease transmission. From this and other studies, it is reasonable to assume that a great variety of hitherto undetected ACoVs exist in wild bird species. Previous studies that examined the host range and genetic diversity of CoVs have revealed that ACoVs in wild birds are present mainly in wildfowl (Anseriformes) and waders (Charadriiformes) [22,33,34]. Our results corroborate these findings and indicate that ACoVs are common among different wild birds. Conclusively, the data show that there is circulation of genetically divergent ACoVs among wild bird populations in geographically distinct areas of Egypt. Although the number of samples analyzed in this study was limited, it can be assumed that the genetic variation of ACoVs among wild birds is much higher than previously thought. Isolates of ACoV that are sequenced and reported here not only show a high similarity with isolates from Israeli, Libya and Oman strains but are also highly similar (99-100%) to the strains previously reported from commercial poultry flocks, which implies that certain strains may have the potential to spread directly or indirectly to distant regions or to other countries.
Although no discrete and globally accepted classification system exists for ACoV, our results revealed that most of ACoVs available on GenBank are divided into variant and classic groups based on the sequence of the full length S1 genes. The classic group has only the Massachusetts serotype (mass-like/vaccine strains) (GI-1 lineage) according to Valsastro et al. [8]. In 2001, an Egyptian variant I strain (Egypt/Beni-Suef/01), closely related to the Israeli variant strain, was reported on different poultry farms [18]. Another variant was also identified in 2011 [27] and designated as Egyptian variant II. Variant genotype can be subdivided into at least eight subgenotypes according to their origin and evolution; 4/91 genotype, Taiwan and China like, Australia variant, QX variant, IS/1494/06 (Egy/Variant-1) and IS/885 (Egy/Variant-2). Egy/Variant I and II are belonging to GI-23 lineage, which represents a unique geographic wild-type cluster confined to the Middle East according to Valsastro et al. [8]. Some have become dominant in the majority of farms and are involved in respiratory and renal pathologies in the last five years.
However, previous studies reported that the variant group (GI-23 lineage) represents a unique geographic wild-type cluster, indigenous and predominate within the Middle East that includes IS/885/00 and IS/1494/06 [8,18,35], recent studies reported that the majority of the Egyptian ACoV strains (IBVs) belong to QX-like, 4/91, IS/1494/06 (Egy/Variant-1) and IS/885 (Egy/Variant-2) genotypes [27,36]. These genotypes were formerly present in the Egyptian commercial farms five years ago with increasing the mortality rates in broiler flocks but had not been described before in wild birds because of low surveillance and epidemiological studies directed towards wild birds worldwide. Our results indicated the prevalence of Egy/variant 1 (IS/1494/06 genotype) (5/14; 35.7%), in spite of relative abundance of Egy/variant 2 (IS/885 genotype) (9/14; 64.3%). The emergence of Egy/variant 1 (IS/1494/06-like) and Egy/variant 2 (IS/885-like) in Egypt has been introduced to the poultry population raises several speculations; it could be due to the introduction of carriers during importations of chicks from Middle East countries where this genotype was endemic or due to the role of migratory birds as a source of infection [19] that has been proposed in this study. Likewise, clustering of some ACoVs isolated from wild and commercial birds in Libya, Oman, and Kurdistan and allocated within the same Egyptian clusters either Egy/variant 1 and or Egy/variant 2 highlights the role of wild birds for transmission of such viruses. These results indicate the circulation of ACoVs in Egyptian wild birds, which could potentially be through the wild birds of neighboring countries and support the possible role of wild birds in the transmission of such viruses.
Recombination has been well reported and is thought to be a contributing factor in the emergence and evolution of different coronavirus genotypes as well as different species of coronavirus [37]. Thor and colleagues reported that the principal mechanisms for generating genetic and antigenic diversity within ACoV indicate a progressive evolutionary change that is the result of recombination events in the ACoV genome that play a major role in the origin and adaptation of the virus, leading to the emergence of new IB genotypes and serotypes [38]. Many recombination events have been reported in different ACoV strains, not only between field (wild-type) and vaccine viruses but also among field viruses either within the same genotype (intra-genotypic) or between different genotypes (inter-genotypic) [39,40], giving rise to new ACoV genotypes [41]. Recombination analysis for the aligned sequences under study report the detection of recombination events; 1) AY091552.2 IBV IS/720/99/S and JX173489.1 IBV/Eg/CLEVB-1/012 strains to produce a recombinant strain MF034384 ACoV/Teal/Kafr El Sheikh-Egypt/VRLCU-13/2016; 2) KC197206.1 IBV/ck/Egypt/12vir6109-79/2012 and JX174186.1 IBV Ck/Eg/BSU-3/2011 strains to produce a recombinant strain MF034381 ACoV/Cattle erget/ Kafr El Sheikh -Egypt/VRLCU-10/2016, which may be a precursor for further variants. Further studies are required to determine the pathobiological and clinical features of this virus in a chicken model. In addition, continuous sequence analysis for the currently circulating IBVs is highly recommended to study the possible spread of this virus and for better understanding of its phylogenetic relatedness to other viruses. These results indicate a progressive recombination of events that occur among wild birds as well as between the circulating variant field strains. Positively selected fragments of genes encoding viral proteins exposed on the surface of the capsid have been documented in other viruses [42,43]. There is an association between positively selected sites along the S1 subunit identified in this study and mapped neutralizing epitopes. It has been reported that mutations in the S1 protein often result in changes in antigenicity [5]. Likewise, parts of the hypervariable regions (HVR1, 2 and 3) defined in this study were shown to be under strong positive selection in the ACoV strains. Taken together, the strong positively selected motifs among the S1 protein may thus be associated with the immune response and receptor binding and would thus be important in future ACoV vaccine development.

Conclusions
This is the first report confirming the detection of ACoV in wild birds in different Egyptian governorates. Phylodynamic and evolutionary analysis indicate that wild bird ACoV strains are closely related to the infectious bronchitis viruses (IBVs) reported from commercial poultry, indicating their possible threat to commercial poultry. Intensified surveillance of wild birds is an important means of assessing the relative prevalence of ACoV variants, and this knowledge would aid risk assessments and risk management of these viruses. Wild bird surveillance that includes virus isolation may also be a tool for obtaining strains of ACoVs that can be used for vaccine development and diagnostics. Moreover, this study provides insights into the genetic diversity of ACoVs in wild birds' reservoirs. Further assessing the prevalence and effects of ACoV infections in wild birds will increase our knowledge on CoV interactions with their hosts and may suggest as yet unexploited avenues for combating CoV infections. There is a clear need for a better understanding of ACoV ecology through full genome sequencing that will provide detailed information on different subgroups, and this would require comprehensive data through better surveillance of wild birds.