Blue-Winged Teals in Guatemala and Their Potential Role in the Ecology of H14 Subtype Influenza a Viruses

Wild aquatic birds are considered the natural hosts of 16 HA (H1–H16) and 9 NA (N1–N9) subtypes of influenza A viruses (FLUAV) found in different combinations. H14 FLUAVs are rarely detected in nature. Since 2011, H14 FLUAVs have been consistently detected in Guatemala, leading to the largest collection of this subtype from a single country. All H14 FLUAVs in Guatemala were detected from blue-winged teal samples. In this report, 17 new full-length H14 FLUAV genome sequences detected from 2014 until 2019 were analyzed and compared to all published H14 sequences, including Guatemala, North America, and Eurasia. The H14 FLUAVs identified in Guatemala were mostly associated with the N3 subtype (n = 25), whereas the rest were paired with either N4 (n = 7), N5 (n = 4), N6 (n = 1), and two mixed infections (N3/N5 n = 2, and N2/N3 n = 1). H14 FLUAVs in Guatemala belong to a distinct H14 lineage in the Americas that is evolving independently from the Eurasian H14 lineage. Of note, the ORF of the H14 HA segments showed three distinct motifs at the cleavage site, two of these containing arginine instead of lysine in the first and fourth positions, not previously described in other countries. The effects of these mutations on virus replication, virulence, and/or transmission remain unknown and warrant further studies.


Introduction
Influenza A viruses (FLUAV) infect a wide range of bird species and mammals, including humans [1]. The virus genome contains 8 segments of negative single-stranded RNA encoding 6 internal (PB2, PB1, PA, NP, M, and NS) and 2 surface (HA and NA) gene segments. Subtype classification is based on the antigenic properties of the surface glycoproteins. To date, 18 HA (H1-H18) and 11 NA (N1-N11) subtypes have been described [2][3][4]. Wild aquatic birds in the orders Anseriformes and Charadriiformes are considered the natural hosts for 16 HA (H1-H16) and 9 NA (N1-N9) subtypes, playing an important role in the perpetuation of FLUAVs in nature. H3, H4, and H6 subtypes are usually detected at high frequencies in the breeding grounds of North America, while others such as H8, H12, H14, and H15 are rarely found [5,6].
The H14 subtype was initially detected in 1982 in 4 virus isolates obtained from 3 mallard ducks and 1 herring gull in the former Soviet Union (Russia and Kazakhstan) [7].

Viral RNA Extraction
Viral RNA from combined tracheal and cloacal swabs was extracted from 250 µL of supernatant with Trizol LS reagent (Invitrogen, Carlsbad, CA, USA) using an organic extraction protocol described previously [16,17]. Extracted RNA was resuspended in 100 µL of DEPC-treated water. All extracted RNA was stored at −70 • C until use.

FLUAV Detection by RRT-PCR
Swabs were screened in duplicate for FLUAV using real-time reverse-transcriptase polymerase chain reaction (RRT-PCR) targeting the matrix gene, as described previously [13,18,19]. Briefly, a QuantiTect Probe RT-PCR (reverse transcription polymerase chain reaction) Kit (QIAGEN, Hilden, Germany) was used to perform RRT-PCR reactions in the ABI 7300 Real-Time PCR System (Applied Biosystems, Foster City, CA, USA). Each reaction contained 12.5 µL of kit-supplied 2X RT-PCR Master mix, 10 pmol of each primer, 0.3 µM probe, 0.25 µL of kit-supplied enzyme mix, 6.5 U RNase inhibitor, and 8 µL of RNA template, as described previously [13,15,16]. Thermal cycling conditions were as follows: one cycle of reverse transcription at 50 • C for 30 min and 94 • C for 15 min, followed by 45 cycles of denaturation at 94 • C for 1 s and combined annealing and extension at 60 • C for 27 s. Fluorescence signals were obtained at the end of each cycle after the annealing/extension step. Swabs showing Ct values < 40 were considered positives.

Virus Isolation
Influenza virus isolation was attempted in triplicate from samples that tested positive using RRT-PCR FLUAV but failed or yielded incomplete genomes after direct NGS sequencing. Each positive bird sample was tested in triplicate. For each bird sample, a pool of the tracheal/cloacal swab sample (200 µL) was inoculated into the allantoic cavity 9-day-old specific-pathogen-free (SPF) embryonated chicken eggs following the protocol described in the World Health Organization's Manual on Animal Influenza Diagnosis and Surveillance [21]. Collected allantoic fluids were tested using hemagglutination assay to assess for the presence of FLUAV. Allantoic fluids that tested negative were diluted 1:2 with sterile phosphate-buffered saline (pH 7.4) with 1X antibiotic and antimycotic solution (Corning, Corning, NY, USA) and inoculated in a new batch of three SPF eggs. The process was repeated three times as needed. Samples were considered negative for the presence of viable virus if no FLUAV growth was detected after three serial passages. Positive allantoic fluid was aliquoted and stored at −70 • C until later use.

Sequencing and Genome Assembly
MS-RTPCR products were sequenced using the Illumina platform as described previously [14] with minor modifications. Briefly, amplicons from MS-RTPCR reactions were cleaned using 0.45X of Agencourt AMPure XP Magnetic Beads (Beckman Coulter, Brea, CA, USA), according to the manufacturer's protocol, and eluted in 30 µL of HyClone molecular biology water (Genesee Scientific, San Diego, CA, USA). Amplicons were quantified using a Qubit buffer kit (ThermoFisher, Waltham, MA, USA) in a Qubit 3.0 fluorometer (Ther-moFisher, Waltham, MA, USA) and normalized to 0.2 ng/µL. Adaptors were added via tagmentation using the Nextera XT DNA library preparation kit (Illumina, San Diego, CA, USA). The reaction was set as 60% of the suggested final volume. Samples were purified using 0.7X of Agencourt AMPure XP Magnetic Beads and analyzed on a Bioanalyzer using a High Sensitivity DNA kit (Agilent, Santa Clara, CA, USA) to determine the distribution of fragment size. Libraries were pooled and normalized to 1-4 nM. After denaturation, the final loading concentration of the pooled libraries was 14 pM. Libraries were sequenced using a MiSeq Reagent Kit V2, 300 cycles (Illumina, San Diego, CA, USA). Genome assembly was performed using a customized pipeline developed at the Icahn School of Medicine at Mount Sinai [22].

Phylogenetic Analyses
Independent phylogenetic analyses were performed for the surface (HA and NA) and internal (PB2, PB1, PA, NS, NP, and M) gene segments. Sequences were aligned with MUSCLE 5.1 [23] and manually trimmed to keep the open reading frame (ORF) of the gene of interest. Background sequences were subsampled to remove identical sequences and those sequences with less than 80% of the total length of the ORF using a SeqKit 0.16.1 bioinformatic tool [24]. The best-fit model of nucleotide substitution was determined for each gene using the Bayesian information criterion (BIC) obtained using jModelTest 2.1.10 [25]. For N3, N4, and N5, as well as the internal gene segments, phylogenetic trees were constructed using maximum-likelihood (ML) inference methods, using RAxML 8.2.12 [26] under the general time reversible GTR + G nucleotide substitution model with 1000 bootstrap replicates. Trees were run at least twice to confirm topology. H14 phylogeny was inferred via time-scaled phylogenetic analysis using the Bayesian Markov chain Monte Carlo (MCMC) approach as implemented in BEAST 1.10.4 [27,28]. A relaxed uncorrelated lognormal (UCLN) molecular clock was used, with a constant population size, and a general-time reversible (GTR + G) model of nucleotide substitution with gamma-distributed rate variation among sites. Three independent analyses of 50 million generations were performed to ensure convergence, sampling every 5000 states. The burn-in percentage of each dataset was identified using Tracer v1.7.1 [29] and after its removal, results were combined using LogCombiner v1.10.4 [27,28]. The MCC trees were annotated using TreeAnnotator v1.10.4 [27,28]. Outlier sequences were identified using TemEst V1.5.3 [30]. When necessary, outlier sequences were removed, and the datasets were run again as described above. The resulting trees were visualized in FigTree 1.4.4 (http://tree.bio.ed.ac.uk/software/Figtree/) and aesthetically modified using Inkscape v0.48.1 (https://inkscape.org). H14 virus phylodynamic analysis was done in BEAST 1.10.4 [27,28] considering a discrete space of three locations (Eurasia, North America, and Guatemala), as well as discrete traits (HA cleavage site and NA), with an asymmetric substitution model. Significant transition rates were identified using Bayes factors (BF): ≥1000 was deemed as decisive support, 100 ≤ BF < 1000 as very strong support, 10 ≤ BF < 100 as strong support, and 3 ≤ BF < 10 as supported. The resulting MCC tree was visualized on SpreaD3 [31] and aesthetically modified using Inkscape v0.48.1 (https://inkscape.org). Gene lineages and genotype assignments were performed based on the phylogenetic trees as described previously [32]. Briefly, FLUAV genes with ≥95% nucleotide identity supported by bootstrap values >70% were classified in the same phylogenetic clade, whereas with ≥90% nucleotide identity within the same lineage.

Amino Acid Sequence Analysis
Nucleotide sequences were translated into amino acids using Geneious Prime 2021.2.2 (https://www.geneious.com). The resulting amino acid sequences were used to identify virulence markers and specific motifs in the sequences. from hunter-harvested wild aquatic birds. From those samples, 96.5% corresponded to BWT. Swab samples that tested positive (n = 234) for FLUAV's matrix (M) gene using RRT-PCR screening were used for MS-RT-PCR followed by NGS sequencing. Whole virus genome sequences were obtained from 89 out of 234 FLUAV-positive samples (38.2%) either from the swab material (n = 45) or virus isolate (n = 44), following the approach described by Ferreri et al. [14]. The HA and NA gene segments were successfully sequenced from an additional 12 swab samples that produced incomplete internal gene segment sequences. In addition, 15 additional swab samples produced at least one gene segment sequence; overall, 116 samples produced at least one complete gene segment.

Genetic and Phylogenetic Characterization of H14 FLUAVs from Guatemala
In this report, we focus the attention on the 17 H14 subtype strains identified (7 from swab material and 10 from viral isolate), all of which were found exclusively in BWT samples. The proportion of other identified subtypes can be found in Figure S1. In combination with previously reported H14 sequences for the 2011-2013 period, this new set of sequences brings the total number of Guatemalan H14 strains to n = 40 full H14 HA gene segment sequences (all of them with complete virus genome sequences). Subtype combinations of all H14 from Guatemala included H14N3 (n = 25), H14N4 (n = 7), H14N5 (n = 4), H14N6 (n = 1), and three mixed infections (2 H14N3/N5 and 1 H14N2/N3). To determine the genetic similarity among the Guatemalan H14 FLUAVs, mean pairwise percentage nucleotide identities were established using either concatenated or the major ORF sequences of each gene segment. Genetic comparisons across the concatenated gene segments (PB2, PB1, PA, HA, NP, M, and NS) from the 40 full genome sequences revealed 95.4% mean pairwise nucleotide identity ( Figure 1).

Genetic and Phylogenetic Characterization of H14 FLUAVs from Guatemala
During five wintering seasons (November-March) of FLUAV surveillance in avian species on the Pacific Coast of Guatemala (seasons 2014-2015, 2015-2016, 2016-2017, 2017-2018, and 2018-2019), 2294 paired tracheal/cloacal swab samples were obtained from hunter-harvested wild aquatic birds. From those samples, 96.5% corresponded to BWT. Swab samples that tested positive (n = 234) for FLUAV's matrix (M) gene using RRT-PCR screening were used for MS-RT-PCR followed by NGS sequencing. Whole virus genome sequences were obtained from 89 out of 234 FLUAV-positive samples (38.2%) either from the swab material (n = 45) or virus isolate (n = 44), following the approach described by Ferreri et al. [14]. The HA and NA gene segments were successfully sequenced from an additional 12 swab samples that produced incomplete internal gene segment sequences. In addition, 15 additional swab samples produced at least one gene segment sequence; overall, 116 samples produced at least one complete gene segment.
In this report, we focus the attention on the 17 H14 subtype strains identified (7 from swab material and 10 from viral isolate), all of which were found exclusively in BWT samples. The proportion of other identified subtypes can be found in Figure S1. In combination with previously reported H14 sequences for the 2011-2013 period, this new set of sequences brings the total number of Guatemalan H14 strains to n = 40 full H14 HA gene segment sequences (all of them with complete virus genome sequences). Subtype combinations of all H14 from Guatemala included H14N3 (n = 25), H14N4 (n = 7), H14N5 (n = 4), H14N6 (n = 1), and three mixed infections (2 H14N3/N5 and 1 H14N2/N3). To determine the genetic similarity among the Guatemalan H14 FLUAVs, mean pairwise percentage nucleotide identities were established using either concatenated or the major ORF sequences of each gene segment. Genetic comparisons across the concatenated gene segments (PB2, PB1, PA, HA, NP, M, and NS) from the 40 full genome sequences revealed 95.4% mean pairwise nucleotide identity ( Figure 1). At the individual gene level, the M1 ORF sequences presented the highest % mean pairwise identity (98.0%), followed by HA 97.7%, PB1 96.9%, NS1 96.6%, and NP 95.2%. The PA of H14 viruses from Guatemala belongs to the two North American lineages with 95.6% (n = 12) and 96.8% (n = 28) pairwise identity within each lineage. Identical H14 HA sequences were detected within the Guatemalan BWT samples (28 out 40), most of them identified from samples collected the same day. Detailed % mean pairwise identity can be found in Table S1-S11.
The 40 full-genome H14 virus sequences from Guatemala were compared to the 16 H14 full-length virus genomes from North America (n = 12, from Arkansas, Ohio, Texas, Wisconsin, California, Mississippi, and Missouri) and Eurasia (n = 4, from Russia, Ukraine, and Pakistan) available on BV-BRC (Bacterial and Viral Bioinformatics Resource Center, https://www.bv-brc.org/, accessed on 1 July 2022) and GISAID (https://gisaid.org/, accessed on 8 July 2022) databases. The similarity of the Guatemalan and North American H14 FLU-AVs is less than 95% based on concatenated internal gene segment sequences, suggesting active virus reassortment in Guatemala. Similar analyses revealed mean pairwise identity of 88.5% compared to H14 Eurasian sequences (Figure 1). Analyses using a single gene segment sequence showed mean pairwise percent identity between Guatemalan and North American sequences as follows:

H14 FLUAVs Were Introduced in Guatemala at Least Twice
Phylogenetic analyses were performed with all H14 HA sequences available to date (Figure 2A). The phylogenetic tree showed segregation of H14 FLUAVs into Eurasian and North American clades [11,13]. The Guatemalan H14 FLUAVs clustered together with viruses from North America with similar dates of collection (Figure 2A). These analyses further revealed one major introduction in 2011 [13] and then a second introduction in 2019. The closest ancestors of the HA of the H14 Guatemalan viruses were identified in the Mississippi flyway as previously observed [11,33]. After its first detection in 2011, the H14 viruses persisted in Guatemala until~2014.6 to 2014.9 (95% highest posterior density (HPD)). The most recent H14 virus identified in Guatemala (samples collected in 2019) clustered in a separate clade with other North American viruses, indicating a recent independent introduction of the H14 subtype in the country. Phylodynamic analysis considering a discrete space of three locations (Guatemala, North America, and Eurasia) with an asymmetric substitution model supports the Eurasian movement of H14 strains into the Americas (Bayes factor, BF = 11,050.9); however, due to the limited number of H14 sequences, it is not possible to establish significant movements between North America and Guatemala ( Figure 2B and Table S12).
It is well accepted that most FLUAVs rely on the host's trypsin-like proteases for cleavage of the HA into subunits, HA1 and HA2, to render the virus infectious. Among the H14 Guatemalan HA ORF sequences, three different predicted motifs at the HA cleavage site were identified. Sequences from samples collected in seasons 2010-2011 and 2012-2013 and samples (n = 3) collected in season 2018-2019 showed the motif PDKQTK'GLF, typical of H14 HA sequences previously identified in North American bird samples [10,11] and two other samples from Eurasia (Figures 2A and 3A). The cleavage site of these H14 HA sequences are unusual compared to other HA subtypes since they contained the amino acid lysine (K) instead of arginine (R) at the −4 and −1 positions on the HA1/HA2 cleavage site [7]. However, we have also detected additional alternative motifs:  amino acid lysine (K) instead of arginine (R) at the −4 and −1 positions on the HA1/HA2 cleavage site [7]. However, we have also detected additional alternative motifs: H14 HAs collected from late 2013 to mid-2016 (seasons 2013-2014, 2014-2015, and 2015-2016) contained either the motif PDKQTR'G (R, at the −1 position, n = 12 sequences) or PDRQTR'G (R, at the −1 and −4 positions, n = 2 sequences), the latter more typical of other HA subtypes.

Diverse Gene Constellations Detected in Blue-Winged Teals in Guatemala
All NAs and internal gene segments of the H14 FLUAVs from Guatemala grouped within the North American lineage (Figures 4 and S2-S5). Mean pairwise comparisons among the ORFs of distinct NA subtypes showed 97.9% sequence identity within the N3 (n = 27), 97.3% (n = 7) within N4, and 97.4% (n = 4) within N5. These analyses also revealed high nucleotide identity (>95%) with contemporary sequences of viruses circulating within the Pacific and Mississippi flyways, consistent with previous observations [13][14][15]. Of note, the NP gene of two H14N3 viruses collected in 2013 are recent introductions from the Eurasian lineage along with those from some contemporary North American viruses isolated from multiple flyways ( Figure 5) and >98% similar to a H3N8 virus from Alaska isolated in 2009 (Table S15) whose introduction from the Eurasian lineage has been previously documented [34]. A discrete traits analysis ( Figure 3A and Table S13) inferred that a PGKQAKG-to-PDKQTKG transition in the HA cleavage site occurred in Eurasia between the 1980s and early 2000s. The PDKQTKG motif persisted in Eurasian birds, where it was detected as recently as 2019. Two independent PDKQTKG-to-PDKQAKG transitions also occurred in Eurasia, but no PDKQAKG motifs have been observed since 2014. The first H14 viruses identified in North America retained the Eurasian PDKQTKG motif, which also persisted in North America until as recently as 2019. A PDKQTKG-to-PDKQTRG transition was observed in Guatemala in the early 2010s, followed by a PDKQTRG-to-PDRQTRG transition, although the latter transition was only observed in only two viruses. NA subtypes were swapped more frequently ( Figure 3B and Table S14). Decisive Bayes factor support (BF > 1000) was observed for the transition from N3 to N2 (BF= 18,833.4) and N3 to N5 (BF= 4703.6). There was no apparent association between the HA cleavage site motif and NA subtype.

Diverse Gene Constellations Detected in Blue-Winged Teals in Guatemala
All NAs and internal gene segments of the H14 FLUAVs from Guatemala grouped within the North American lineage (Figure 4 and Figure S2-S5). Mean pairwise comparisons among the ORFs of distinct NA subtypes showed 97.9% sequence identity within the N3 (n = 27), 97.3% (n = 7) within N4, and 97.4% (n = 4) within N5. These analyses also revealed high nucleotide identity (>95%) with contemporary sequences of viruses circulating within the Pacific and Mississippi flyways, consistent with previous observations [13][14][15]. Of note, the NP gene of two H14N3 viruses collected in 2013 are recent introductions from the Eurasian lineage along with those from some contemporary North American viruses isolated from multiple flyways ( Figure 5) and >98% similar to a H3N8 virus from Alaska isolated in 2009 (Table S15) whose introduction from the Eurasian lineage has been previously documented [34].  Phylogenetic inference via the maximum likelihood method revealed internal genes of H14 sequences grouped in subclades of reference sequences within a particular lineage ( Figures 5, 6 and S7-S11). Overall, the 40 H14 Guatemalan FLUAVs were associated with 20 clades with bootstrap support values of >70% of each segment (PB2 = 5, PB1 = 3, PA = 5, NP = 4, M = 1, and NS = 2) ( Figure 6). Interestingly, the internal gene segments of samples collected the same day and with identical H14 HA sequences were positioned in different clades on the phylogenetic tree, highlighting the high virus diversity found in these viruses in Guatemala. Phylogenetic inference via the maximum likelihood method revealed internal genes of H14 sequences grouped in subclades of reference sequences within a particular lineage ( Figures 5, 6 and S7-S11). Overall, the 40 H14 Guatemalan FLUAVs were associated with 20 clades with bootstrap support values of >70% of each segment (PB2 = 5, PB1 = 3, PA = 5, NP = 4, M = 1, and NS = 2) ( Figure 6). Interestingly, the internal gene segments of samples collected the same day and with identical H14 HA sequences were positioned in different clades on the phylogenetic tree, highlighting the high virus diversity found in these viruses in Guatemala.   Further analyses of the remaining predicted ORFs revealed neither mammalian-associated virulence markers in PB2 (E627K and/or D701N) [35,36] and PA (S409N) [37], nor NS1 protein (T92E) [38] or antiviral resistance markers (V27A and S31N on M2 [39] or H274Y on NA). The PB1-F2 N66 (asparagine) marker associated with increased virulence in mammals [40] was detected in 15 out 40 samples. The samples with the PB1-F2 N66 marker encoded an 87 aa long protein, whereas the rest carrying the PB1-F2 S66 (serine) signature encoded a 90 aa long protein. PB1-F2 proteins with either one of those lengths have been previously described in FLUAVs from ducks [39]. All Guatemalan H14 FLUAVs encoded a full-length PA-X protein of 252 amino acids in length. Further analyses of the remaining predicted ORFs revealed neither mammalianassociated virulence markers in PB2 (E627K and/or D701N) [35,36] and PA (S409N) [37], nor NS1 protein (T92E) [38] or antiviral resistance markers (V27A and S31N on M2 [39] or H274Y on NA). The PB1-F2 N66 (asparagine) marker associated with increased virulence in mammals [40] was detected in 15 out 40 samples. The samples with the PB1-F2 N66 marker encoded an 87 aa long protein, whereas the rest carrying the PB1-F2 S66 (serine) signature encoded a 90 aa long protein. PB1-F2 proteins with either one of those lengths have been previously described in FLUAVs from ducks [39]. All Guatemalan H14 FLUAVs encoded a full-length PA-X protein of 252 amino acids in length.

Discussion
During more than a decade of FLUAV surveillance in wild birds in Guatemala, an unusually large number of H14 subtype viruses were detected [13][14][15][16]41]. Prior to its initial isolation in North America in 2010, there was no evidence of significant circulation of this subtype in the Western hemisphere [9]. The H14 subtype FLUAVs appeared predominantly in Guatemala from the first detection in 2011 [13,14,16] through 2019. In this study, we analyzed a small data set of available sequences from H14 viruses worldwide (n = 62); however, our current observations are consistent with our initial findings [13]. Due to the limited information about this subtype, its ecology and how it is maintained in this region is unknown. One hypothesis for the high level of H14 detected in Guatemala is that it may be due to epizootic events originated from a single introduction, followed by local clonal expansion of the newly introduced subtype that later was maintained in local populations [42]. Studies of H3 subtype FLUAV epizootic events in mallards suggest that migrant birds play an important role as vectors of novel strains that are introduced and amplified by local bird populations [42]. It is feasible that H14 FLUAVs are maintained in susceptible hosts during non-migration seasons and later transmitted during migration season by a routinely sampled host, such as BWT. H14 viruses from Guatemala were exclusively detected from hunter-harvested BWT. Therefore, our sampling is biased towards game bird species that are in large numbers. BWT is one of the most abundant species of dabbling ducks in North America and has been found in breeding grounds from North America to wintering grounds in Mexico, Central America, and parts of South America [43]. Due to the wide distribution of BWT across the American continent, it is intriguing as to why the H14 subtype has become predominant only in Central America in BWT. We believe that yet to be defined local conditions in Guatemala may have provided the adequate environment for this subtype to thrive. We speculate that such conditions may explain in part the detection of H14 in Guatemala from 2011 to 2015, when all H14 HAs shared the same common ancestor. However, the distinct NA subtype combinations and internal gene segment constellations of Guatemalan H14 viruses within and between migratory seasons indicate active reassortment and independent introductions from unknown sources.
When the first H14 HA subtype was initially sequenced, one feature that stood out was the presence of K instead of R at the −1 position in the HA cleavage site [7,8,10,11,33]. Thus, previously described H14 HA from North America showed the PDKQTK'G motif, which is unusual among FLUAV subtypes. This rarer motif was detected initially during our surveillance, but it was also still present in Guatemala in 2019. Our continuous surveillance and sequencing approach directly from swabs revealed two additional cleavage site motifs in the H14 HAs, PDKQTR'G and PDRQTR'G, which place them closer to the more typical cleavage sites observed in other HA subtypes. It is commonly accepted that the HA cleavage is predicted to be a substrate for trypsin-like proteases found in the lumen of the intestinal tract of natural FLUAV hosts. Cleavage of HA is necessary for virus infection that modulates virulence. At this stage it is unclear how these different cleavage site motifs in the HA affect replication, virulence, and/or transmission of H14 viruses in natural and/or in potential accidental hosts (poultry, mammals) [44]. To our knowledge, it remains to be determined whether H14 FLAUV infections cause any disease in wild birds. Similarly, to our knowledge, there have been no reports of outbreaks of H14 FLUAVs in either poultry or mammals. It is also unclear whether our ability to detect different motifs is related to sequencing directly from swabs as opposed to sequencing from virus isolates obtained through passage in chicken eggs (an artificial substrate with potential bottleneck effects) [14]. Overall, continued detection of H14 FLUAVs in Guatemala appears dictated by potential multiple drivers in the avian reservoir that warrant further studies.

Supplementary Materials:
The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/v15020483/s1, Figure S1: Number of HA and NA identified in Guatemala; Figure S2: Detailed ML phylogenetic tree for N3; Figure S3: Detailed ML phylogenetic tree for N4; Figure S4: Detailed ML phylogenetic tree for N5; Figure S5: Detailed ML phylogenetic tree for N6; Figure S6: Detailed ML phylogenetic tree for PB2; Figure S7: Detailed ML phylogenetic tree for PB1; Figure S8: Detailed ML phylogenetic tree for PA; Figure S9: Detailed ML phylogenetic tree for NP; Figure S10: Detailed ML phylogenetic tree for M; Figure S11: Detailed ML phylogenetic tree for NS; Table S1: Detailed nucleotide pairwise identity of ORF sequences of concatenated genes; Table S2: Detailed nucleotide pairwise identity of ORF sequences of PB2; Table S3: Detailed nucleotide pairwise identity of ORF sequences of PB1; Table S4: Detailed nucleotide pairwise identity of ORF sequences of PA; Table S5: Detailed nucleotide pairwise identity of ORF sequences of HA; Table S6: Detailed nucleotide pairwise identity of ORF sequences of NP; Table S7: Detailed nucleotide pairwise identity of ORF sequences of M; Table S8: Detailed nucleotide pairwise identity of ORF sequences of NS; Table S9: Detailed nucleotide pairwise identity of ORF sequences of N3; Table S10: Detailed nucleotide pairwise identity of ORF sequences of N4; Table S11: Detailed nucleotide pairwise identity of ORF sequences of N5; Table S12: Bayes factor and posterior probability values using geographic locations; Table S13: Bayes factor and posterior probability values using HA cleavage site; Table S14: Bayes factor and posterior probability values using NA subtype; Table S15: Best ten BLAST hits of Guatemalan isolates.