Spatiotemporal Dynamics, Evolutionary History and Zoonotic Potential of Moroccan H9N2 Avian Influenza Viruses from 2016 to 2021

The H9N2 virus continues to spread in wild birds and poultry worldwide. At the beginning of 2016, the H9N2 Avian influenza virus (AIV) was detected in Morocco for the first time; despite the implementation of vaccination strategies to control the disease, the virus has become endemic in poultry in the country. The present study was carried out to investigate the origins, zoonotic potential, as well as the impact of vaccination on the molecular evolution of Moroccan H9N2 viruses. Twenty-eight (28) H9N2 viruses collected from 2016 to 2021 in Moroccan poultry flocks were isolated and their whole genomes sequenced. Phylogenetic and evolutionary analyses showed that Moroccan H9N2 viruses belong to the G1-like lineage and are closely related to viruses isolated in Africa and the Middle East. A high similarity among all the 2016–2017 hemagglutinin sequences was observed, while the viruses identified in 2018–2019 and 2020–2021 were separated from their 2016–2017 ancestors by long branches. Mutations in the HA protein associated with antigenic drift and increased zoonotic potential were also found. The Bayesian phylogeographic analyses revealed the Middle East as being the region where the Moroccan H9N2 virus may have originated, before spreading to the other African countries. Our study is the first comprehensive analysis of the evolutionary history of the H9N2 viruses in the country, highlighting their zoonotic potential and pointing out the importance of implementing effective monitoring systems.


Introduction
The H9N2 AIV was initially isolated from turkeys in Wisconsin in the United States of America in 1966 (A/turkey/Wisconsin/1/1966(H9N2)) [1]. In Asia, the first cases of H9N2 infection in poultry were recorded in Hong Kong in the 1970s [2]. In the 1990s, the H9N2 spread became more extensive with several outbreaks affecting industrial poultry in China, Germany, South Korea, Italy, Iran, Pakistan, the USA, South Africa, and Ireland [3][4][5]. This lead to a persistent viral circulation in many countries in Asia, the Middle East, and North Africa [6], where the virus has attained endemic status in poultry populations [7][8][9]. Today, H9N2 is considered the most prevalent, and damaging, low-pathogenic avian influenza (LPAI) subtype in poultry worldwide [7], presenting a threat to both animal and human health [9,10]. H9N2 infections usually cause minor clinical disease in poultry, which makes its identification in the field rather challenging [11]. However, despite their low pathogenicity, H9N2 viruses may cause variable pathogenicity and mortality in domestic poultry [5,6,12], especially when associated with immunosuppression or coinfection with other viruses or bacteria [3,11,13,14], resulting in significant economic losses [5,6,12].
The H9N2 virus has a zoonotic potential with interspecies transmission between poultry and mammals, including swine and humans [15][16][17][18]. At least 59 human cases have been reported to date [9]. Furthermore, H9N2 viruses have proved able to provide partial or even complete sets of internal gene segments to zoonotic avian influenza viruses (AIVs), leading to the emergence of human-lethal reassortants, as in the case of China in 2013 with the appearance of H7N9 and H10N8 viruses, and, a year later of the H5N6 subtype [19][20][21][22][23][24][25][26] posing serious risks to public health. Moreover, there are multiple reports of H9N2 natural infection in other mammals, such as ferrets [27], mice [28], minks, foxes, raccoon dogs [29], cats, and dogs [30].
Since the early 2000s, the emergence of H9N2 viruses has been reported in several North African countries, including Libya in 2006 [33], Egypt in 2006 [34], and Tunisia in 2009 [35]. Additionally, following the first isolation of the LPAI H9N2 virus in Morocco in January 2016 [36], H9N2 viruses have been extensively circulating in North Africa, spreading to West Africa, as well as to East Africa. Since 2017, H9N2 viruses have been detected in Algeria [37] and in several Sub-Saharan African countries, including Burkina Faso [38], Ghana [32], Uganda, Kenya, Senegal [39], Benin Republic, and Togo [40]. All H9N2 viruses detected in poultry in Africa belong to the G1 lineage and are related to those isolated in the United Arab Emirates.
In Morocco, following its first identification in 2016, the H9N2 virus has spread rapidly across the country, causing significant economic losses in poultry [36]. A previous study revealed that the Moroccan isolate might be related to the Middle East isolate A/chicken/Dubai/D2506.A/2015 [41]. The virus showed high pathogenicity under field conditions with aggravation of the disease in case of coinfection with other respiratory diseases [41,42]. Despite the swift implementation of a mass poultry vaccination program following the onset of the disease, the virus became endemic in all poultry sectors of the country [43]. Knowledge of the evolutionary history and the circulation of the Moroccan H9N2 viruses is crucial for public and animal health in Morocco. A better understanding of the epidemiological and geographic origins of the Moroccan H9N2 viruses might help to improve the control and management of future outbreaks. In this study, we characterized the whole genome of H9N2 AIVs collected in Morocco from 2016 to 2021, exploring the spatiotemporal dynamics and the acquisition of mutations of particular concern for animal and human health.

Sample Collection
From 2016 to 2021, thanks to the establishment of a program to monitor the presence of influenza viruses in poultry in Morocco, samples were collected from 81 poultry farms located in seven different regions in the country, including layers, breeders, turkeys, and broiler flocks of varying ages. Most flocks (72.84%, 59/81) were vaccinated against AIV H9N2. Oropharyngeal, cloacal swabs and organ samples were collected from chickens and turkeys, stored at 2-8 • C, and sent to the laboratory within 24 h for further analysis.

Real Time RT-PCR for H9 Gene Detection
Swabs, as well as tissue samples, were re-suspended in 1 mL of phosphate-buffered saline (pH 7.4 +/− 0.2, supplemented with 10,000 IU/mL penicillin, 10 mg/mL streptomycin, 0.25 mg/mL gentamycin, and 5000 IU/mL nystatin), thoroughly homogenized by using Bead Blaster 24 Homogenizer (Benchmark Scientific, 2600, Sayreville, NJ, USA). After centrifugation at 1500× g for 10 min at 4 • C, total RNA was purified with the NucleoSpin ® RNA Virus Kit (Macherey-Nagel, Düren, Germany) following the manufacturer's instructions. RNA was then eluted in 50 µL RNase-free water. All samples were tested by RT PCR for the influenza A matrix gene as described by Spackman, E. et al. (2002) [44]. The positive samples were then tested for the H9 gene detection by RT-PCR according to Monne, I. et al. (2008) [45] by using the AgPath-IDTM One-Step RT-PCR Kit (Ambion, Applied Biosystems, Grand Island, NY, USA). The rRT-PCR was performed with Aria Mx (Agilent Technologies, Santa Clara, CA, USA) and ABI 7500 Fast (Applied Biosystems, Foster City, CA, USA) 96-well plate real-time PCR devices.

Whole Genome Sequencing
The complete genomes were amplified using SuperScript™ III One-Step RT-PCR System with Platinum™ Taq High Fidelity DNA Polymerase (Invitrogen, Carlsbard, CA, USA), as previously described by Zhou B. et al. (2009) [46]. Amplicons were purified by using Agencourt AMPure XP (Beckman Coulter Inc., Brea, CA, USA) and quantified with Qubit™ DNA HS Assay (Thermo Fisher Scientific, Waltham, MA, USA). Sequencing libraries were prepared with the Illumina Nextera XT DNA Sample Preparation Kit (Illumina, San Diego, CA, USA) and sequenced on the Illumina MiSeq platform (2 × 250 bp Paired-End; Illumina, San Diego, CA, USA).
Phylogenetic network analyses using the Median Joining method implemented in the program NETWORK 10 (https://www.fluxus-engineering.com/, accessed on 4 June 2021) were performed on the HA gene sequences of the H9N2 influenza viruses collected in Morocco in the period 2016-2021 [55]. This method uses a parsimony approach to reconstruct the relationships between highly-similar sequences, and allows the creation of "median vectors," which represent un-sampled sequences used to connect the existing genotypes in the most parsimonious way. The evolutionary rate and the time to the most common recent ancestor (tMRCA) were estimated for all the gene segments using BEAST v1.10.4 [56]. Markov chain Monte Carlo (MCMC) sampling analyses employing an uncorrelated lognormal relaxed molecular clock, that allows for rate variation across lineages, and the SRD06 substitution model (HKY85 + Γ4 with two partitions-1st + 2nd positions vs. 3rd position -, base frequencies, and Γ-rate heterogeneity unlinked across all codon positions) were implemented [57]. Tree Annotator v1.10.4 (http://beast.bio.ed.ac.uk/TreeAnnotator/, accessed on 11 August 2021) was used to summarize the Maximum Clade Credibility (MCC) trees, which were then visualized in FigTree v1.4.4 (http://tree.bio.ed.ac.uk/software/figtree/, accessed on 11 August 2021).

Analysis of Selection Pressures
Site-specific selection pressures for the HA gene segment of the Moroccan H9N2 viruses were measured as the ratio of nonsynonymous (dN) to synonymous (dS) nucleotide substitutions per site. Specifically, the Mixed Effects Model of Evolution (MEME) and the fixed-effects likelihood (FEL) methods were used to infer sites subjected to episodic or pervasive positive selection [58,59]. Analyses were performed using the Datamonkey online version of the Hy-Phy package [60]. All analyses used the GTR model of nucleotide substitution and ML phylogenetic trees.

Phylogenetic and Evolutionary Analyses
In the framework of the surveillance activity implemented in Morocco between 2016 and 2021, field samples collected from 81 poultry farms were submitted to the laboratory to be tested for the presence of AIV. Overall, 76 out of 81 (76/81) tested positive for AIV and specifically for the H9 subtype ( Figure 1 shows the location of the sampled and the positive farms). In total, 28 representative positive samples were selected based on the year and area of collection for whole genome sequencing. Specifically, we obtained the complete genome of 22 viruses and the partial genome of 6 viruses (Table 1).  Figures S1-S7). The tree topologies indicate that the H9N2 viruses from Morocco had evolved into multiple genetic groups within the country. However, the lack of a consistency of the groups among the different phylogenies is suggestive of the occurrence of multiple inter-subtype reassortment events.  Virus evolution into multiple genetic clusters is evident also from the HA medianjoining network analysis. Although all the 2016-2017 sequences (yellow in Figure 3) show a high genetic relationship (from 1 to 16 nucleotide differences) and cluster together on the centre of the network, the viruses collected in 2018-2019 (green in Figure 3) and 2020-2021 (blue in Figure 3) separate from the 2016-2017 progenitor viruses by long branches (from 15 to 54 nucleotide differences) that spread from the centre in multiple directions, forming separate clusters. The high number of median vectors, representing the lost ancestral sequences, clearly indicate a relevant data gap. No particular grouping was observed according to the geographical origin, which suggests a wide virus dispersal within the country.  To investigate the role of Morocco in the geographic spread of the H9N2 viruses of the G1 lineage, we performed a Bayesian phylogeographic analysis. Our results showed a virus spread from the Middle East to Morocco, and subsequently from Morocco to Algeria and from Morocco to West Africa. All of these transitions were well supported by Bayes Factors higher than 20 (BF > 20, Figure 4). This finding suggests that Morocco might have played a central role in the spread of the H9N2 virus in the African continent, although this interpretation may have been biased by the higher number of sequences from Morocco compared to the ones from other African countries.

Molecular Analysis and Positively Selected Sites
The molecular analysis revealed that the H9N2 viruses from Morocco present the same cleavage site (RSSR/GLF) as the previously identified North African H9N2 viruses of the G1 lineage. In addition, most of the Moroccan viruses possess the amino acid Leucine at position 226 and Threonine at position 155 (H3 numbering) of the HA receptor binding site; these mutations imply preferential binding to human-like α2-6-linked sialic acid (SA α2-6) receptors [61,62]. Of note, site 226 has been identified as under positive selection by both MEME and FEL methods (p-value < 0.05) ( Table 2). The H9N2 Moroccan viruses show three different amino acids in this position: Leucine (L), Glutamine (Q), and Arginine (R). These amino acid substitutions may affect not only the virus receptor binding properties but also influence its antigenic specificity [63]. Although amino acids L and Q have frequently been identified in viruses belonging to the G1 lineage, R has rarely been observed in the H9N2 viruses (one virus from Algeria and one from Morocco). Previous studies demonstrated that this mutation is associated with change of the pandemic 2009 (H1N1) viral HA binding preference, from the human-type receptor, α2,6-linked sialic acid, to the avian-type receptor, α2,3-linked sialic acid [64,65].  In addition to 226, in the HA gene we identified three further sites (sites 190, 160, and 186) under positive selection: one (site 190) using both methods, while the other two were recognized only by FEL. Of note, all these sites fall within the receptor-binding domain and/or antigenic sites [63,66,67].

Discussion
In Morocco, the first outbreak of LPAI H9N2 was described in broiler flocks in January 2016 [36]. The virus was later isolated from layers and breeders in different regions of the country. As part of a surveillance program, the national veterinary services implemented mass vaccination against H9N2 using homologous inactivated vaccines in all sectors of poultry production, first for breeders and layers, and a few days later for turkey and broiler flocks. The outbreak was controlled in 3-4 months, after which the vaccination was authorized for all types of poultry production using autogenous and heterologous inactivated vaccines. Over the past 5 years, 13 H9N2 vaccines originating from different H9 strains, i.e., from Morocco, UAE, Korea, Turkey, Italy, Egypt, were available and commercialized in Morocco (see Supplementary Materials Table S2). Despite this vaccine availability, a recent study on the effectiveness of the four most used vaccines in Moroccan poultry farms has confirmed the need to continuously adapt the vaccine strains to the virus circulating in the field [69]. In the neighboring North African countries, vaccination against AIV H9N2 is also implemented in poultry production sectors using different vaccine types, such as inactivated homologous vaccines and vectored vaccines in Algeria. In Tunisia, vaccination using inactivated homologous vaccines for chicken and turkey breeders and layers has been performed since 2015, but to the present day, no vaccination has been used in broiler flocks. In Egypt, vaccination is currently implemented using local and non-local strains, with variable efficacy. In contrast, vaccination in the field has never been reported in Libya [70]. Despite such vaccination approaches, LPAI H9N2 has become enzootic in North Africa and continues to cause economic losses even in the vaccinated flocks. Thus, circulating H9N2 viruses in North Africa need to be further investigated, and at the same time, revision of employed vaccines and of regional surveillance activities should be implemented.
Our results indicate that all Moroccan H9N2 viruses from 2016 to 2021 belong to the G1 like-lineage and are related to the viruses from the Middle East and North and West Africa, as described in previous reports [43,71]. Results obtained from the spatial analysis demonstrated that AIV H9N2 had been introduced into Morocco from the Middle East, probably through the import of falcons for hunting, and has been circulating in the country for the last six years. The study also showed that the H9N2 G1 like-lineage is the predominant H9 genotype in Moroccan poultry farms. From Morocco, the virus spread to Algeria and West Africa, suggesting the key role of Morocco in the H9N2 virus dispersal on the continent, as reported in previous studies [32,[37][38][39][40]. Such virus-spread might have occurred through the movement of poultry professionals [38] and/or the trade of poultry products, since the export market of day-old-chicks and hatching eggs from Morocco to West Africa is well documented [72]. However, the gap in epidemiological data, as well as the limited number of sequences available from the other African countries, might have affected our results and biased our interpretation. Evolutionary and network analyses of the H9N2 viruses from Morocco indicate a high similarity among all the 2016-2017 HA sequences, which is in accordance with previous studies [43,71]. Long branches separate the 2018-2019 and 2020-2021 viruses from their ancestors (2016)(2017). This might be explained by the implementation of vaccination from 2016 onwards. The high vaccination coverage and the use of multiple types of vaccines may have driven the virus evolution in multiple directions, with the consequent generation of multiple separate genetic groups in the country. This was further confirmed by the identification of several sites in the HA protein under positive selection located in the RBS or antigenic sites, which may facilitate vaccine escapes, considering the fact that most (71.05%, 54/76) of the H9N2 positive flocks were vaccinated.
As previously observed in the H9N2 African viruses [32,38,40,73], most of Moroccan H9N2 HA protein harbored mammalian-adaptive mutations, including the hemagglutinin human-like marker 216L (226 in H3 numbering), which promotes preferential binding to the α2-6-linked sialic acid receptors of human epithelial cells [60] and was described as an immune escape mechanism [62]. The evidence for positive selection of the 226L marker during serial passage of the virus in chickens was provided by Jegede et al., 2018 [74], suggesting that the H9N2 viruses acquired the mammalian adaptation while circulating in poultry. However, specific receptor-binding studies will be instrumental to assess the zoonotic potential of the H9N2 Moroccan viruses, as other mutations in the HA and internal proteins may be necessary for a successful avian to human adaptation. The present work is the first comprehensive study that provides an overview of the historical origins and evolutionary features of Moroccan H9N2 AIV, based on the Bayesian phylodynamic analysis and whole genome characterization. Over the past five years, the H9N2 AIV in Morocco has shown a rapid genetic diversification, likely directed by selective pressure exerted by vaccine-induced antibodies. The observed mutations and reassortments might extend the range of host species and could lead to the emergence of new strains [75,76]. In addition, our study demonstrates that viruses with a marker well-described in literature associated with zoonotic potential (226L, H3 numbering) may represent a serious risk for public health in Morocco. Hence, the circulation of H9N2 AIV should be investigated in other species, such as wild birds and mammals, including humans, to better understand virus ecology in Morocco. All the more so as a recent molecular study on the presence of the AIV in wild birds staying in the corridor of wetlands during their winter crossing in Morocco, showed a high positivity-rate of the influenza virus in different birds species [77]. Continuous systematic surveillance of molecular evolution of H9N2 AIV is strongly needed to anticipate the emergence of novel H9N2 strains escaping vaccine protection or with an increased human pandemic potential.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/v14030509/s1, Figure S1: Maximum Likelihood phylogenetic tree of the NA gene segment (IQ-TREE v.1.6.8). The H9N2 viruses from Morocco analyzed in this report are marked in blue. Ultra-fast bootstrap supports equal to or higher than 80% are indicated next to the nodes; Figure S2: Maximum Likelihood phylogenetic tree of the PB2 gene segment (IQ-TREE v.1.6.8). The H9N2 viruses from Morocco analyzed in this report are marked in blue. Ultra-fast bootstrap supports equal to or higher than 80% are indicated next to the nodes; Figure S3: Maximum Likelihood phylogenetic tree of the PB1 gene segment (IQ-TREE v.1.6.8). The H9N2 viruses from Morocco analyzed in this report are marked in blue. Ultra-fast bootstrap supports equal to or higher than 80% are indicated next to the nodes; Figure S4: Maximum Likelihood phylogenetic tree of the PA gene segment (IQ-TREE v.1.6.8). The H9N2 viruses from Morocco analyzed in this report are marked in blue. Ultra-fast bootstrap supports equal to or higher than 80% are indicated next to the nodes; Figure S5: Maximum Likelihood phylogenetic tree of the NP gene segment (IQ-TREE v.1.6.8). The H9N2 viruses from Morocco analyzed in this report are marked in blue. Ultra-fast bootstrap supports equal to or higher than 80% are indicated next to the nodes; Figure S6: Maximum Likelihood phylogenetic tree of the M gene segment (IQ-TREE v.1.6.8). The H9N2 viruses from Morocco analyzed in this report are marked in blue. Ultra-fast bootstrap supports equal to or higher than 80% are indicated next to the nodes); Figure S7: Maximum Likelihood phylogenetic tree of the NS gene segment (IQ-TREE v.1.6.8). The H9N2 viruses from Morocco analyzed in this report are marked in blue. Ultra-fast bootstrap supports equal to or higher than 80% are indicated next to the nodes; Figure S8: the HA median-joining network for Moroccan regions; Table S1: Moroccan H9N2 viruses from 2016 to 2021 included in the study; Table S2: List of H9N2 vaccines available in Morocco.

Data Availability Statement:
The primary data used to support the findings of this study are available from the corresponding author upon request.