Introduction

The rapid and widespread global circulation of influenza A(H1N1)pdm09 virus in 20091,2,3,4 and severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) in 20205,6,7,8,9 demonstrate how respiratory viruses can quickly spread globally following their emergence. Influenza A viruses (IAV) and novel coronaviruses remain a major public health threat with potential to cause pandemics resulting in economic losses, social lockdowns, and millions of deaths10,11,12,13. Therefore, global surveillance efforts to monitor these pathogens through various initiatives, for example, the Global Influenza Surveillance and Response System (GISRS)14 remain justified to improve understanding of their global transmission networks and dynamics for appropriate interventions.

The global surveillance of influenza through GISRS has resulted in the generation of geographically and temporally extensive virus sequence data, which has provided a unique opportunity to investigate the global spread of influenza viruses15,16,17,18,19,20. Three key hypotheses have been previously proposed to explain how influenza viruses spread globally. First, the “source-sink” model, whereby East and Southeast (E-SE) Asia represent a global source population of novel influenza virus strains, while temperate regions represent ecological sinks15,16,18. Second, influenza viruses occur as temporally migrating metapopulations, whereby new virus strains can emerge in any geographical region, with the location of the source population changing from season-to-season17. Third, the global patterns of spread of seasonal influenza viruses have been shown to vary with the rate of genetic and antigenic evolution of different influenza virus types and subtypes19. However, despite documented high disease burden21,22,23, the role of sub-Saharan African countries in the global spread of influenza viruses remains unclear due to insufficient surveillance and sequence data from these regions22,24.

To improve understanding of IAV dispersal pathways in sub-Saharan Africa and the context of locally observed virus diversity, we studied the spread of A(H1N1)pdm09 and A(H3N2) viruses in the region using phylogeographical methods based on codon-complete sequences from samples obtained from a multi-country surveillance involving five sub-Saharan African countries. We then compared these data with sequences that were available from other countries in sub-Saharan Africa and around the globe.

Materials and methods

Sample sources and molecular screening

The samples analysed here were collected under three distinct studies: (a) the Pneumonia Etiology Research for Child Health (PERCH) study, which enrolled children aged between 1 and 59 months admitted to hospital with severe or very severe pneumonia in surveillance sites in Kenya, Zambia, Mali, Gambia, and South Africa25,26,27 from August 2011 through January 2014; (b) a Kenya-wide surveillance of influenza viruses in Kenya among patients of any age hospitalized with severe acute respiratory illness (SARI)28 from January 2011 through December 2013; and (c) inpatient paediatric (1–59 months) pneumonia surveillance of influenza viruses in Kilifi County and Referral Hospital (KCH)29 from January 2011 through December 2013. Detailed descriptions of the study sites and population, sample collection, transportation, storage, and processing have been described elsewhere; PERCH study25,27,30, the countrywide surveillance study28, and the inpatient paediatric pneumonia surveillance study29. For the PERCH study, a total of 138 IAV positive specimens were analysed from five surveillance sites from five countries (Fig. S1A): Kenya, 27; South Africa, 40; Zambia, 27; The Gambia, 22; and Mali, 22. A further 106 A(H1N1)pdm09 and 16 A(H3N2) virus sequences were available from the countrywide surveillance study and the paediatric KCH study, respectively.

RNA extraction and multi-segment real-time PCR (M-RTPCR) for IAV

Viral RNA extraction and M-RTPCR was undertaken as previously described29. Briefly, we performed viral nucleic acid extraction from IAV positive samples (Cycle threshold (Ct) < 35.0) using the QIAamp Viral RNA Mini Kit (Qiagen). We reverse transcribed ribonucleic acid (RNA) and amplified the codon-complete region of influenza in a single M-RTPCR using the Uni/Inf primer set31 in 25 μL PCR reactions. We evaluated successful amplification by running the M-RTPCR amplicons on 2% agarose gel and visualized the gels looking for expected bands on a UV transilluminator after staining with RedSafe Nucleic Acid Staining solution (iNtRON Biotechnology Inc.,).

IAV codon-complete sequencing and genome assembly

Following M-RTPCR, we purified, quantitated and normalized amplicons to 0.2 ng/μL as described previously29. Briefly, the amplicons were purified with 1X AMPure XP beads (Beckman Coulter Inc.,), quantified with Quant-iT dsDNA High Sensitivity Assay (Invitrogen), and normalized to 0.2 ng/μL. Indexed paired end libraries were generated from 2.5 μL of 0.2 ng/μL amplicon pool using Nextera XT Sample Preparation Kit (Illumina) following the manufacturer’s protocol. Amplified libraries were purified using 0.8X AMPure XP beads, quantitated using Quant-iT dsDNA High Sensitivity Assay (Invitrogen), and evaluated for fragment size in the Agilent 2100 BioAnalyzer System using the Agilent High Sensitivity DNA Kit (Agilent Technologies). Libraries were then diluted to 2 nM in preparation for pooling and denaturation for running on the Illumina MiSeq (Illumina). Pooled libraries were NaOH denatured, diluted to 12.5 pM and sequenced on the Illumina MiSeq using 2 × 250 bp paired end reads with the MiSeq v2 500 cycle kit (Illumina). Five percent Phi-X (Illumina) spike-in was added to the libraries to increase library diversity by creating a more diverse set of library clusters. We carried out contiguous nucleotide sequence assembly from the sequence data using the FLU module of the Iterative Refinement Meta-Assembler (IRMA)32 using IRMA default settings.

Collation of contemporaneous global sequence dataset

Global comparison datasets for influenza A(H1N1)pdm09 and A(H3N2) viruses were retrieved from the GISAID EpiFlu™ database (https://platform.gisaid.org/epi3/cfrontend; accessed 20 March 2022). The datasets were prepared to determine the relatedness of the viruses in this report to those circulating around the world thus understand their global context. Only codon-complete full-length genome sequences sampled between January 2010 and December 2013 were included to improve the phylogenetic resolution of the analyses and investigate reassortment patterns. For A(H1N1)pdm09 virus, the final dataset of 460 global sequences sampled from January 2010 to December 2013 was available (numbers in parenthesis indicate number of sequences): Africa (13); Asia (148); Europe (82); North America (170); South America (2); Oceania (39). For A(H3N2) virus, the final dataset of 749 global sequences sampled from January 2010 to December 2013 was available (numbers in parenthesis indicate number of sequences): Africa (7); Asia (196); Europe (38); North America (178); South America (199); Oceania (138). The accession numbers for the global genome sequences are available in the study’s GitHub repository, (https://github.com/DCollinsOwuor/Phylogeography-and-reassortment-patterns-of-human-influenza-A-viruses-in-sub-Saharan-Africa).

Phylogenetic analysis

We aligned and translated consensus nucleotide sequences in AliView v1.2633 and concatenated individual gene segments into codon-complete genomes using SequenceMatrix v16.0.134. We constructed codon-complete gene segment and concatenated genome phylogenetic trees of A(H1N1)pdm09 and A(H3N2) viruses with maximum-likelihood (ML) and bootstrap analysis of 1,000 replicates. We inferred the best-fit nucleotide substitution models using IQ-TREE v1.6.1135,36 and implemented those chosen by the Bayesian Information Criterion for each concatenated virus genome. We visualized and annotated the phylogenetic trees using ggtree R package v2.3.537 and used the full-length HA sequences of all virus sequences generated from the PERCH, countrywide surveillance, and inpatient paediatric pneumonia surveillance study to characterize A(H1N1)pdm09 and A(H3N2) viruses into clades using PhyCLIP v2.038 and IAV clade representative strains.

Reassortment analysis

We used tanglegrams of the time-resolved trees from TreeTime v0.8.539 to visualize the phylogenetic relationships between individual codon-complete gene segments of A(H1N1)pdm09 and A(H3N2) viruses using dendextend R package v1.15.240 and explored the pairwise phylogenetic congruence between IAV gene segments to identify intra-subtype reassortment. Tanglegrams are sets of two trees representing different segments of IAVs, with tips from matching viruses connected by lines. We compared topologies of tanglegram trees that matched the phylogeny of HA segment with that of other segments, coloring tanglegram twines based on observed reassortment patterns. We then estimated the reassortment rates (the number of reassortment events per lineage per year (events/lineage/year)) among the PERCH study and sub-Saharan African IAVs using a coalescent reassortant constant population model (CoalRe)41 in BEAST2 v.2.6.6 (https://www.beast2.org/) with parameters: GTR + I + G4 substitution model, strict clock, prior for reassortment rate = exponential with mean 0.25, chain length of 500 million, and 10 per cent burn-in. We then used the Graph-incompatibility-based Reassortment Finder (GiRaF) tool42 to computationally assess reassortment among IAVs. Our GiRaF runs, that computationally assessed reassortments, converged with an average standard deviation of split frequency < 0.01, tree length estimated sample size (ESS) value > 100, tree length potential scale reduction factors (PSRF) of 1.0, and a maximum split frequency PSRF of 1.00–1.002 (Tables 1 and 2). We manually discarded the first 25 percent of trees from each of the two runs (burn-in), then combined and processed the remaining trees. Only reassortment events with a high confidence value (≥ 0.70) and reassortant sets predicted in at least three of twenty-eight gene pairwise comparisons, as recommended by GiRaF developers were reported. We ran two simultaneous runs of MrBayes v3.2.743 on the datasets (A(H1N1)pdm09 viruses, 200 million generations; A(H3N2) viruses, 100 million generations) under a GTR + I + Γ substitution model and sampled trees at every 200,000th generation.

Table 1 MrBayes convergence diagnostics for A(H1N1)pdm09 viruses.
Full size table
Table 2 MrBayes convergence diagnostics for A(H3N2) viruses.
Full size table

Spatial dynamics of spread of IAVs

We conducted a phylogeographical analysis to assess the global spread of A(H1N1)pdm09 and A(H3N2) viruses using methods implemented in BEAST v1.10.4 package44, with an asymmetric discrete trait approach that applied the Bayesian stochastic search variable selection (BSSVS) model. To reduce the complexity of the maximum clade credibility (MCC) inference, we categorized location states into geographical regions (“Africa”, “Asia”, “E-SE Asia”, “Europe”, “North America”, “South America”, or “Oceania”) and visualized phylogeographic inferences with the SPREAD3 software v0.9.7.1c package45. Due to the important role of E-SE Asia in the global spread of influenza viruses, we partitioned Asia into E-SE Asia and rest of Asia to elucidate the geographical transition states associated with E-SE. To visualize the geographic spread of the virus over time, we generated a D3 file using SPREAD3 v0.9.7.1c package and used a global geo.json file for visualization. We then used the resulting log files to calculate the estimated virus migration rates between regions and Bayes factor (BF) values for significant migration rates between discrete locations: we deemed ≥ 1,000 as decisive support, 100 ≤ BF < 1000 as very strong support, 10 ≤ BF < 100 as strong support, and 3 ≤ BF < 10 as supported rates.

We also used the ML tree topologies for A(H1N1)pdm09 and A(H3N2) viruses from our mugration models to estimate the number of virus transmission events between sub-Saharan Africa and the rest of the world. Using the date and location-annotated tree topologies, we counted the number of transitions within and between sub-Saharan Africa and the rest of the world and plotted using ggplot2 v3.3.3 R package37.

Ethics

Scientific and ethical clearance for the study was obtained from institutional ethics review boards within each study country in Africa, the UK, and USA26. Additional ethical approval was sought and received from KEMRI Scientific Steering Committee (SSC# 1055 and 1433) and Oxford Tropical Research Ethics Committee (OxTREC# Ethics ID: 60-90). Informed consent was sought and received from the study participants for the study.

Results

IAV sequencing and codon-complete genome assembly of IAV genomes

Of the 138 IAV positive samples that were available from the PERCH study, 100 (73%) were successfully sequenced. The recovered genome sequences were classified into A(H1N1)pdm09 (n = 31) and A(H3N2) (n = 69). These were combined with the additional sequences from Kenya (106 for A(H1N1)pdm09 and 16 for A(H3N2)). Among A(H1N1)pdm09 viruses from the PERCH study, all viruses from 2011 to 2013 fell into six clades namely, clade 3 (n = 5), clade 7 (n = 1), clade 8 (n = 2), clade 6A (n = 1), clade 6B (n = 2), and clade 6C (n = 20) (Fig. 1 and Fig. S1B). Among A(H1N1)pdm09 viruses from the countrywide Kenya study, all viruses from 2011 to 2013 fell into clade 6 (n = 96) and clade 6C (n = 10) (Fig. 1 and Fig. S1B). For A(H3N2) viruses, all PERCH Africa study viruses from 2011 to 2013 fell into five clades namely, clade 3A (n = 6), clade 3B (n = 2), clade 3C.2 (n = 1), clade 3C.3 (n = 39), and clade 7 (n = 21) (Fig. 2 and Fig. S1C). Among A(H3N2) viruses from the KCH paediatric study, all viruses from 2011 to 2013 fell into clade 3B (n = 4) and clade 7 (n = 12) (Fig. 2 and Fig. S1C).

Figure 1
figure 1

Phylogenetic trees of the eight individual gene segments of influenza A(H1N1)pdm09 viruses from sub-Saharan Africa from the 2011–2013 influenza seasons. The tree branches and tips are coloured by virus clade, as shown on the colour-coded key. The positions of intra-subtype reassortant viruses in each phylogeny are indicated by filled stars for PERCH study viruses while viruses from other sub-Saharan African countries are indicated by filled diamonds. PERCH pneumonia etiology research for child health.

Full size image
Figure 2
figure 2

Phylogenetic trees of the eight individual gene segments of influenza A(H3N2) viruses from sub-Saharan Africa from the 2011–2013 influenza seasons. The tree branches and tips are coloured by virus clade, as shown on the colour-coded key. The positions of intra-subtype reassortant viruses in each phylogeny are indicated by filled stars for PERCH study viruses while viruses from other sub-Saharan African countries are indicated by filled diamonds. PERCH pneumonia etiology research for child health.

Full size image

Patterns of reassortment of IAVs from sub-Saharan Africa

We investigated the reassortment patterns by comparing the position of viruses on the eight segment-specific phylogenies of A(H1N1)pdm09 and A(H3N2) viruses to identify inconsistencies arising from intra-subtype reassortment. For most of the virus strains, those belonging to the same clade, as inferred from HA phylogeny, clustered together on phylogenies generated from the seven remaining gene segments for A(H1N1)pdm09 (Fig. 1) and A(H3N2) viruses (Fig. 2), respectively. However, in 44 (29%) A(H1N1)pdm09 virus (Fig. 1) and eight (9%) A(H3N2) virus (Fig. 2), intra-subtype reassortants were suspected based on tip locations of viruses in the gene segment phylogenies. Of these reassortants, 31 A(H1N1)pdm09 and four A(H3N2) viruses were identified among our PERCH study samples. GiRaF correctly identified all the 52 intra-subtype reassortant viruses. Interestingly, ten of the 44 A(H1N1)pdm09 and two of the eight A(H3N2) reassortant viruses were recently identified in a continentwide reassortment analysis of IAVs in Africa using GiRaF46.

The relatedness of these reassortant viruses were visualized using tanglegrams based on pairwise phylogenetic congruence between codon-complete gene segments for A(H1N1)pdm09 (Fig. 3) and A(H3N2) (Fig. 4) viruses, with tanglegram twines colored by virus clade to highlight reassortant strains. Additionally, for the viruses detected in sub-Saharan Africa, coalescent model analysis estimated a mean reassortment rate of 0.0225 [95% HPD, 0.1272–0.3231] events/lineage/year among A(H1N1)pdm09 viruses and a mean reassortment rate of 0.0577 [95% HPD, 0.0146–0.1027] events/lineage/year among A(H3N2) viruses. We further observed that the reassortant A(H1N1)pdm09 and A(H3N2) viruses persisted in circulation for 1–2 consecutive years in sub-Saharan Africa based on their date-annotated tree topologies.

Figure 3
figure 3

Evolutionary relationships of each gene segment of influenza A(H1N1)pdm09 viruses from sub-Saharan Africa, displaying intra-subtype reassortment with the subtype, highlighting the reassortant strains by clade from the 2011–2013 influenza seasons.

Full size image
Figure 4
figure 4

Evolutionary relationships of each gene segment of influenza A(H3N2) viruses from sub-Saharan Africa, displaying intra-subtype reassortment with the subtype, highlighting the reassortant strains by clade from the 2011–2013 influenza seasons.

Full size image

Viral imports and exports in sub-Saharan Africa

We assessed how 150 A(H1N1)pdm09 virus and 94 A(H3N2) virus genomes from sub-Saharan Africa compared to 441 A(H1N1)pdm09 and 749 A(H3N2) contemporaneous virus genomes, sampled globally by inferring their phylogenies. For A(H1N1)pdm09 virus, we inferred 15 importations originating from outside sub-Saharan Africa (seven from North America, six from Asia, and two from Europe), which represent 15 independent introductions into sub-Saharan Africa from geographical locations outside the region (Fig. 5A). We also inferred four virus location transition events among the sub-Saharan African countries. Additionally, we captured six virus export events from sub-Saharan Africa to North America, Asia, and Europe, all occurring from Kenya (Fig. 5A). For A(H3N2) viruses, we inferred 10 importations originating from outside sub-Saharan Africa (three each from North America and Asia and two each from Oceania and South America), which represent 10 independent introductions into sub-Saharan Africa from geographical locations outside the region (Fig. 5B). Furthermore, we inferred 13 virus location transition events among the sub-Saharan African countries. Additionally, six virus export events from sub-Saharan Africa to North America and Oceania were inferred (Fig. 5B).

Figure 5
figure 5

Number of viral imports and exports from sub-Saharan Africa shown as an alluvium plot using virus sequences from sub-Saharan African countries, Asia, Europe, North America, South America, and Oceania for (A) influenza A(H1N1)pdm09 viruses and (B) influenza A(H3N2) viruses.

Full size image

Global migration dynamics of IAVs

For A(H1N1)pdm09 virus, we observed significant migration pathways from E-SE Asia (0.81–3) into multiple geographical regions including sub-Saharan Africa, Asia, Europe, North America, and Oceania. Additionally, we observed significant migration pathways from North America (0.52–0.73) to E-SE Asia, Europe, and Oceania (Table 3). The observed migration pathways corroborate these results (Fig. S2). For A(H3N2) virus, we observed significant migration pathways from E-SE Asia (0.62–3.34) into all other geographical regions including sub-Saharan Africa, and from North America (0.6–1.63) to Europe, Oceania, and South America (Table 4). The global migration pathways for A(H3N2) corroborate these results (Fig. S3).

Table 3 Asymmetrical migration rates between global location states inferred using the BSSVS model for influenza A(H1N1)pdm09 virus.
Full size table
Table 4 Asymmetrical migration rates between global location states inferred using the BSSVS model for influenza A(H3N2) virus.
Full size table

Discussion

We observed that IAV strains from 2011 to 2013 from our study evolved over consecutive influenza seasons and fell into distinct A(H1N1)pdm09 and A(H3N2) virus clades, some of which persisted in circulation as intra-subtype reassortants for consecutive influenza seasons. Furthermore, we observed multiple introductions of IAVs into sub-Saharan Africa over consecutive influenza seasons, with viral importations originating from multiple global geographical regions, for example, North America, Asia, Europe, and Oceania. Additionally, we inferred virus transfer between the sub-Saharan African countries and virus exports from sub-Saharan Africa to other geographical regions, for example, North America, Asia, and Europe. On a global scale, IAVs spread from multiple geographical regions to multiple geographical destinations including sub-Saharan Africa.

Co-circulation of different IAV subtypes and clades is common during epidemics of seasonal influenza, which may facilitate intra-subtype reassortment events among circulating viruses47,48,49,50,51. Overall, the reassortment rates among the sub-Saharan African A(H1N1)pdm09 viruses (0.1272–0.3231) were comparable to the 0.15–0.8 events/lineage/year estimated among global A(H1N1)pdm09 viruses41. However, sub-Saharan African A(H3N2) viruses reassorted at a lower mean rate of 0.0577 [95% HPD, 0.0146–0.1027] than the global estimate of 0.35–0.65 events/lineage/year41. Although A(H3N2) viruses undergo frequent reassortment events globally, the low reassortment rate among some A(H3N2) viruses might be due to the resulting reassortant viruses being unfit and negatively selected, thus not detected in appreciable frequencies52,53. Circulation of intra-subtype reassortants in our study for 1–2 consecutive years demonstrates that IAVs can persist in the population and spread geographically. A recent continentwide study of the rates and patterns of reassortment of IAVs in Africa revealed that novel reassortant viruses emerged every year, with some reassortant viruses persisting in different countries and regions for 1–5 consecutive years46. Persistence of intra-subtype reassortants has been demonstrated previously52,54, although the factors associated with such persistence patterns at a population level require further investigation. For example, a reassortant A(H3N2) virus subclade, designated 3C.3A2/re, emerged in 2016–2017 influenza season and resulted in more hospitalizations and deaths in the 2017–2018 North American influenza season52. Increased whole genome sequencing over consecutive influenza seasons, particularly in understudied sub-Saharan Africa regions, would allow for an improved understanding of the frequency and timing of such intra-subtype reassortments and the contribution to the evolutionary and transmission dynamics of seasonal influenza.

Globally, we observed significant migration pathways for IAVs. Taken together, these results suggest that the global spread of A(H1N1)pdm09 and A(H3N2) virus epidemics are driven by different geographical regions, which also includes sub-Saharan Africa, in which E-SE Asia and North America are major transmission sources17,19,20,55. Our findings support the notion that influenza viruses persist as temporally migrating metapopulations in which new virus strains can emerge in any geographical region, with the location of the source population changing regularly17. This underscores the need for improved influenza surveillance particularly in understudied sub-Saharan Africa regions for complete understanding of the global patterns of spread of influenza.

The paucity of sequence data from other African countries limited our analysis of the patterns of spread and persistence of IAVs in the continent, which might have been useful to demonstrate intra-continental spread of influenza viruses in detail, since persistence may be facilitated by climatic variability that generates temporally overlapping epidemics in neighboring countries. Despite this paucity, our study leveraged on existing national surveillance studies, for example, from the countrywide surveillance of influenza viruses in Kenya where A(H1N1)pdm09 virus sequences were included. It is possible that Africa plays a bigger role in the spread of influenza, which we might not have captured. Furthermore, the analysis in this report only involved the coding regions of the A(H1N1)pdm09 and A(H3N2) virus gene segments. Although noncoding regions are conserved, mutations that affect viral replication may occur, and this information may not have been captured in this study.

In conclusion, despite sparse data from influenza surveillance in sub-Saharan Africa, our findings support the notion that influenza viruses persist as temporally structured migrating metapopulations in which new virus strains can emerge in any geographical region, including in sub-Saharan Africa; these lineages may have been capable of dissemination to other continents through a globally migrating virus population. Further knowledge of the viral lineages that circulate within understudied sub-Saharan Africa regions is required to inform vaccination strategies in those regions.