Emergence and spread of two SARS-CoV-2 variants of interest in Nigeria

Lineage dynamics of SARS-CoV-2 in Nigeria over the first three epidemic waves

We generated a total of 1577 genomes from samples obtained between March 2020 and September 2021, across the first three waves of the pandemic (Fig. 1A and B). We collected samples from 25 of the 36 states and the Federal Capital Territory of Nigeria (Fig. 1C). The northern states were highly undersampled relative to the southern states, with sampling highly uneven across time (Fig. 1A and B).

• Download figure
• Open in new tab

Fig. 1:

Genomic epidemiology of SARS-CoV-2 in Nigeria.A) Reported daily news Epicurve of SARS-CoV-2 cases in Nigeria (blue bars) in Nigeria and genomes produced assembled during the first three waves (red line) with the y-axis transformed to a log-scale. B) Time-varying sampling fraction of genome produced in this study per new casesC) Geographic distribution of sequences generated in the current study. The map of Nigeria shows the number of genomes from states across the country as region colour and marker size display lineages per state. D) Lineage frequency profile for the study period. The frequency of PANGO lineages across the country over a period of 80 weeks. Lineages that are not VOC or VOIs, or that do not appear more than three times over the course of the pandemic, are grouped a “others”. E) Human mobility data from Google for retail shops, pharmacies, parks, public transport, workplaces, and residential activity places from January 2020 to November 2021, as compared to the pre-pandemic baseline (January to February, 2021).

To investigate the lineage dynamics across the different outbreak waves, we used the pangolin nomenclature tool to assign all sequences to corresponding lineages. We detected more than 35 lineages across the first three waves (Fig. 1D). We found that a number of ancestral lineages circulated during the first wave, including A, B, B.1, and B.1.1 (Fig. 1D). The onset of community transmission in Nigeria was initially delayed, with the first wave initiated two months after the first case detection in March 2020. A nationwide lockdown and ban on local and international air travel were enforced, with the first wave ending in August 2020 after four months (Fig. 1B). During the lockdown, we found that there was up to a 50% decline in activity at workplaces, retail, supermarkets, and public transport stations, and up to a 25% increase in residential activity (Fig. 1E).

After this strictly adhered-to lockdown, Nigeria opened its airspace to international traffic in October 2020. This was shortly followed by an increase in new cases in November 2020, marking the start of the second wave (Fig. 1B), which was characterised by B.1.1.7, B.1.525, B.1.1.318 (VOC and VOI), and other variants such as L.3 and B.1.1.487 (Fig. 1D).

Our secondary analysis focused on the lineages (B.1.525 and B.1.1.318) that contributed to the surge during the second wave and have also been shown to possess genotypic traits for increased transmissibility and also increased virulence in human and animal cell lines, such as the E484K spike protein RBD mutation. We excluded B.1.1.7 from our analysis as it has already been well studied^22,23, and it was not the dominant variant in Nigeria during the second wave, unlike in other places where it was reported, such as the UK and the USA^24,25,26. Moreso, the second wave began after travel restrictions were lifted and human mobility was back to normal, thus the need to investigate the emergence of B.1.1.318 and B.1.525 in Nigeria. The third wave was initiated in June 2021, with the Delta variant (B.1.617.2) and its sub-lineages sweeping to dominance in sampling (Fig. 1B and D). In the secondary analysis described in this work, we also sought to identify the timings of emergence and number of introductions of these VOI.

Regional connectivity drove the introduction of B.1.1.318 into Nigeria

Using Bayesian phylogeographic reconstructions, we investigated the timing and origin of the B.1.1.318 emergence to better understand Nigeria’s connectivity in the global SARS-CoV-2 transmission networks (see Methods). B.1.1.318 was first detected in Nigeria in Lagos State in December 2020. The lineage was detected in multiple countries within and outside Africa, notably resulting in a large outbreak in Mauritius¹⁹.

In our phylogeographic reconstruction, we found that B.1.1.318 emerged in Africa (root state posterior support = 0.95-6) in early August [mean tMRCA = 5 August, 95% HPD 25 June to 20 September] across two replicates (Fig. 2A). We estimated that B.1.1.318 was introduced to Nigeria on at least 53 independent occasions [mean introductions, 95% HPD 50-59] beginning in November 2020, after travel restrictions were lifted in October (Fig. 2B). We found that the majority of introductions originated from other African nations (Fig. 2B), indicating Nigeria’s strong connectivity to the region. However, the number and origin of introductions and exports estimated from genomic data is sample dependent, with bidirectional transmission to and from undersampled regions obscured by uneven sampling globally. To overcome this, we quantified air travel patterns to and from Nigeria over time to better understand which countries were most likely to act as transmission sources or sinks based on their connectivity to Nigeria.

• Download figure
• Open in new tab

Fig. 2:

A) Time-resolved B.1.1.318 phylogeny. Branches are coloured by region-level geographic state reconstruction. Internal nodes annotated with black points represent posterior support > 0.75. B) Mean number of introductions (Markov jumps) of B1.1.318 into Nigeria from all source regions, binned by month. The mean number of introductions overall (inset). C) Volume of monthly air passengers inbound to Nigeria across 2020-2021 by source region (bar, left axis) and fold increase in travel volume compared to May 2020 baseline. D) Estimated Introduction Intensity Index for Nigeria from 2020-2021. E) Percentage of inbound air travel to Nigeria by source region across 2020-2021. See Supplementary Figure 2 for country-level data.

Before travel restrictions were lifted in October 2020, we found that the low levels of incoming air travel volume predominantly originated from other African nations, followed by the Middle East (driven by the UAE) and Europe (driven by the UK) (Fig. 2C-D, SFigure 2). This suggests that the introduction risk was predominantly driven by regional connectivity during the period when travel restrictions were in place. We sought to investigate this further by quantifying the introduction risk for all countries connected to Nigeria by air travel. We estimated the introduction intensity index (III), which accounts for the number of cases in the source countries as well as the travel volume to Nigeria (see Methods). Overall, we found that the introduction risk from air travel was very low during the period of travel restrictions (May 2020 to October 2020), as incoming travellers originated largely from other African nations with low incidence (Fig. 2F).

We found that the number of incoming air passengers peaked in December 2020 - January 2021 (Fig. 2C-D). We also observed that the number of inbound passengers from Europe, the Middle East, and Northern America increased relative to other African nations after restrictions were lifted, resulting in comparable volumes from Africa, Europe, and the Middle East (Fig. 2C-D). The 30-fold increase in incoming passengers by December 2020 was reflected in the corresponding peak in the III, attributable to increased introduction risk from Europe and North America from November 2020 into early 2021 (Fig. 2F). The surge was predominantly driven by the high-incidence second waves of the UK and USA, respectively (SFigure 1).

We found that the III peak but not the estimated sources of risk were consistent with the number and origin of introductions of B.1.1.318 in our phylogeographic reconstructions (Fig. 2B). We estimated the highest introduction risk for Europe after travel restrictions were lifted, whereas the majority of introductions estimated from genomic data originated in other African countries. The III from other African nations was relatively low, despite comparable incoming air travel volume. This is likely both a factor of the relatively low incidence in most African countries during this period (Fig. 2F, SFigure 3) as well as the fact that our travel data does not account for land-based travel, which will expectedly severely underestimate the introduction risk from surrounding countries. The III from the African region was largely driven by South Africa, which was experiencing a peaked incidence from its second wave (SFigure 3). In our III analyses, we found negligible introduction risk from outside of Europe and North America.

• Download figure
• Open in new tab

Fig. 3:

A) Time-resolved B.1.525 phylogeny. Branches are coloured by region-level geographic state reconstruction. Internal nodes annotated with black points represent posterior support > 0.75. B) Mean number of exports (Markov jumps) of B.1.525 from Nigeria to destination regions, binned by month. Mean number of exports overall (inset). C) Volume of monthly air passengers outbound from Nigeria across 2020 - 2021 by destination region (bar, left axis) and fold increase in travel volume compared to baseline (May 2020). D) Percentage of outbound air travel from Nigeria by destination region across 2020-2021. See Supplementary Figure 1 for country-level data. E) Mean number of introductions (Markov jumps) of B1.525 into Nigeria from all origin regions, binned by month. The mean number of introductions overall (inset). F) Estimated Exportation Intensity Index (EII) for Nigeria from 2020-2021.

The Eta variant of interest (B.1.525) likely emerged in Nigeria

The B.1.525 (Eta) lineage drove the second wave from January to March 2021. B.1.525 was initially reported to have originated in Nigeria or the UK^16,17,18. We investigate the timing, origin, and transmission dynamics of B.1.525 with Bayesian phylogeographic reconstructions (see Methods). We estimated that the B.1.525 lineage emerged in Nigeria (root state posterior support = 0.998) in late July [mean tMRCA 23 July, 95% HPD 7 June to 4 September], during the final weeks of Nigeria’s first wave (Fig. 3A). This suggests that this lineage circulated cryptically for several months before its first detection in Nigeria on 12 December 2020, likely owing to low-incidence-associated sparse sampling (Fig. 1A).

We sought to better understand Nigeria’s connectivity in global SARS-CoV-2 transmission networks by identifying the most likely destinations of B.1.525’s exportation, as the lineage 1) emerged in Nigeria and 2) is Nigeria’s most sampled lineage, excluding the Delta variant sublineages. To investigate the spread of B.1.525 from Nigeria, we reconstructed the timing and pattern of geographic transitions out of and into Nigeria across the full posterior of our Bayesian phylogeographic reconstructions. We estimated that B.1.525 was exported from Nigeria a lower bound of 295 times [mean Markov jumps, 95% HPD 259-335] (Fig. 3B). The majority of sampled exports were destined for Europe, followed by the rest of Africa and Northern America from December 2020 to March 2021 (encompassing Nigeria’s second wave) (Fig. 3B). During this period, we also found that B.1.525 was re-introduced from Europe a lower bound of 20 [95% HPD 6-37] times (Fig. 3E). These estimated source-sink dynamics support Nigeria’s strong connectivity, particularly to Europe, but will underestimate exports to undersampled regions such as neighbouring African countries. To mitigate the effect of global sampling biases on genomic estimates of transmission dynamics, we again analysed changes in air travel to and from Nigeria over time to understand the sources and seeds of Nigeria’s bidirectional transmission dynamics.

African nations dominated the reduced level of outbound travel during the period of travel restrictions (May to October 2020) (Fig. 3C-D), suggesting that regional connectivity drove Nigeria’s export risk during this period. We integrated the air travel data with Nigeria’s epidemic incidence data to quantify an exportation intensity index (EII)²⁷. The EII quantifies the temporal trend in the daily estimated number of viral exports from Nigeria to sink countries (see Methods). Given the data available, Nigeria’s estimated export intensity was low overall compared to countries with high air traffic data to Nigeria, as the country’s SARS-CoV-2 incidence remained comparatively low on a global-scale for the entire period under investigation (Fig. 3F). Nigeria’s first wave had relatively low incidence, peaking at about 4000 cases in a week in late June (week 18) (SFigure 4B). Combined with reduced travel volumes, the period from May 2020 to September 2020 was characterised by negligible export intensity (Fig. 3F).

Outbound air travel recovered gradually over 2020, peaking in December 2020 - January 2021 (Fig. 3C-D). We found that this peak corresponded with the peak of the larger second wave in Nigeria and therefore also the EII, as well as the distribution of exports of B.1.1.525 from Nigeria (Fig. 3B, F). The proportion of outbound air travel destined for countries in Europe, the Middle East, and North America increased relative to travel destined for African countries after travel restrictions were lifted (Fig. 3C-D). We estimated the highest EII for Europe (driven predominantly by travel to the UK) and Northern America (driven by travel to the USA) from December 2020 to February 2021 (SFigure 4). The export intensity to Europe is consistent with the high number of estimated exports in the genomic data, though North America was disproportionately underestimated in the genomic data (Fig. 3B). Overall, air travel volume destined for Europe was only moderately higher (10%) than Africa, North America, and the Middle East from December (after first detection and at the start of the second wave) to April 2021 (end of the second wave) (Fig. 3C-D). However, Europe had 4-fold more sampled introductions of B.1.525 than African nations or the USA. Notably, the comparatively high EII for the Middle East region (representing 20% of outbound volume overall) was not represented in our phylogeographic reconstructions (Fig. 3B). This highlights how unevenly distributed surveillance capacity can lead to an underestimation of transmission events (SFigure 3A).

Sample preparation

RNA was extracted from nasopharyngeal swabs in viral transport media and saliva/sputum in PBS using the QiAmp viral RNA mini kit (Qiagen) according to the manufacturer’s instructions and MagMax pathogen RNA/DNA kit (Applied Biosystems) using a Kingfisher Flex purification system (ThermoFisher Scientific) according to the manufacturer’s instructions.

Probe-based SARS-CoV-2 real-time RT-PCR

RT-qPCR screening of suspected samples was carried out targeting N, ORF1ab and RdRP genes of the virus using commercially available kits: genesig® (Primerdesign Ltd, UK), DaAnGene (Daan Gene Co., Ltd, China), Liferiver (Shanghai ZJ Bio-Tech Co., Ltd, China), Genefinder (OSANG Healthcare Co., Ltd, Korea), and Sansure (Sansure Biotech Inc., China).

SARS-CoV-2 whole-genome sequencing

Samples with moderate to high viral load detected with RT-qPCR (Ct value < 30) were selected for sequencing, which were samples collected from genomic surveillance and hospitalised individuals spatially and temporally. The sampling gap was the inability to sequence samples from every state during all different waves of the pandemic in Nigeria. The Illumina COVIDseq protocol was used for sequencing preparation (https://www.illumina.com/products/by-type/clinical-research-products/covidseq-assay.html). RNA was converted to cDNA, followed by tiling amplification (400 bp) with Artic V3 primer pools, which covers the whole genome of the virus. Nextera DNA flex libraries were made from the amplicons. Libraries were pooled and sequenced on the Illumina MiSeq, NextSeq 2000, and NovaSeq 6000 platforms at the African Center of Excellence for Genomics of Infectious Diseases (ACEGID), Redeemer’s University, Ede, Nigeria.

Genome Assembly

We used the viral-ngs pipeline v2.1.19 (https://github.com/broadinstitute/viral-ngs) for demultiplexing, quality control, and genome assembly, consistent with prior work^33,34,35. The demux_plus applet of viral-ngs was used for demultiplexing basecalls to BAM files. Genome assembly was carried out using the assemble_refbased workflow on the viral-ngs pipeline which maps each unmapped BAM file to the SARS-CoV-2 reference genome (NC_045512.2) to generate coverage plots and FASTA files of successfully assembled genomes.

Lineage/clade assignment

We used Phylogenetic Assignment of Named Global Outbreak LINeages (PANGOLIN)³⁶ to identify the lineages of SARS-CoV-2 circulating within Nigeria and to identify the lineages circulating in specific states of Nigeria. Nextclade (https://clades.nextstrain.org/) was used to assign the sequences to globally circulating viral clades and to investigate mutations in the genome that could affect primer amplification during PCR.

Maximum likelihood phylogenetic analysis

1,577 genomes with sequence length ≥ 20,930 (covering 70% or more of the reference - NC_045512.2) were assembled and aligned using MAFFT³⁷. We used IQTREE³⁸ and the GTR (generalised time-reversible) model with a bootstrap value of 1000 to construct the phylogenetic tree. Treetime³⁹ was used for quantifying the molecular clock and identifying ancestral phylogeny in the augur pipeline of Nextstrain⁴⁰.

Bayesian phylogenetic analysis

All B.1.1.318 (n = 3,858 with major locations: Europe, USA, and Africa accounting for 94%) and B.1.525 (n = 8,278 with major locations: Europe, USA, and Africa accounting for 86%) sequences were downloaded from GISAID⁴¹ on 18 August 2021. Sequences with > 5% ambiguous nucleotides, a length < 95%, or incomplete dates were discarded. Sequences were aligned to the reference (NC_045512.2) using minimap2, with the 5’ and 3’ UTRs and known problematic sites (GitHub - W-L/ProblematicSites_SARS-CoV2 3) masked⁴². Both lineage-specific datasets were downsampled for representativeness and computational tractability with a phylogenetic-informed downsampling scheme. A maximum-likelihood phylogeny was reconstructed for each lineage with automatic model selection in IQTREE⁴³. The phylogeny was downsampled by root-to-tip tree-traversal, with all internal nodes subject to two rules: 1) if 95% of the leaves subtended by the internal node represented a single country, the earliest representative was retained alongside a random sequence; 2) if leaves from the same location were separated by a zero-branch length, the earliest representative was retained alongside a random sequence. All countries with more than the median number of sequences across countries were then downsampled randomly to the median, this yielded a total number of 1,118 and 1,746 genomes for B.1.1.318 and B.1.525 respectively. All Nigerian sequences were retained (n = 73 for B.1.1.318; n = 256 for B.1.525).

We reconstructed Bayesian time-scaled phylogenies for each lineage using BEAST 1.10.5⁴⁴. For both lineages, the time-scaled phylogenies were reconstructed under an HKY substitution model with a gamma-distributed rate for variation among sites⁴⁵, a relaxed molecular clock with a log-normal prior⁴⁶, and an exponential growth coalescent tree prior⁴⁷. For each dataset, we combined two independent Markov Chain Monte Carlo (MCMC) chains of 200 million states, run with the BEAGLE package⁴⁸ to improve run time. Parameters and trees were sampled every 20,000 steps, with the first 20% of steps discarded as burn-in. Convergence and mixing of the MCMC chains were assessed in Tracer v1.7, to ensure the effective sample size (ESS) of all estimated parameters was > 200⁴⁹.

We performed an asymmetric discrete trait analysis using BEAST version 1.10.5 to reconstruct the location-transition history across an empirical distribution of 4000 time-calibrated trees (sampled from each of the posterior tree distributions estimated above). We aggregated the country of sampling on a regional level for computational tractability. We used Bayesian stochastic search variable selection (BSSVS) to infer non-zero migration rates and identify the statistically supported transition routes into and out of Nigeria by a Bayes factor test⁵⁰. In addition to the discrete trait analysis, we used a Markov jump counting approach to estimate the timing and origin of geographic transitions into Nigeria to account for uncertainty in phylogeographic reconstruction associated with sparse sampling and low sequence variability⁵¹. We used the TreeMarkovJumpHistoryAnalyzer from the pre-release version of BEAST 1.10.5 to extract the Markov jumps from posterior tree distributions⁵². We used TreeAnnotator 1.10 to construct Maximum clade credibility (MCC) trees for all datasets. Trees were visualised using baltic (https://github.com/evogytis/baltic).

Air traffic and human movement data

In order to associate international air travel and local human movement with variant movement across borders, we used air traffic data from the International Air Transportation Association (IATA) obtained from BlueDot (https://bluedot.global/) to quantify the volume of international travellers to and from Nigeria. The dataset included the number of travellers by origin and destination country, aggregated by month, across all international airports in Nigeria. Data was obtained for the period May 2020 to April 2021, encompassing the emergence and spread of B.1.525 and B.1.1.318. We also curated human movement data within Nigeria from Google Mobility (https://www.google.com/covid19/mobility/) from January 2020 to November 2021, which reflects a baseline before lockdown (January - March 2020) and the study’s period of interest. We used R to display movement patterns after grouping human movement into six divisions: retail & recreation, grocery & pharmacy, parks, transport stations, workplaces, and residential areas.

Estimated introduction and exportation intensity index

From genomic data, estimates of the number and origin/destination of bidirectional transmissions with Nigeria are dependent on the sample’s representativeness of the population and are therefore limited by global and local sampling biases. To limit these sampling biases, we supplemented our phylogeographic analyses by estimating the introduction (III) and export intensity index (EII) according to Du Plesis et al²⁷. The III estimates the daily risk of introductions into Nigeria from each country as a product of the number of asymptomatically infected individuals in each source country on that day (estimated from the time series of deaths) and their likelihood of traveling to Nigeria (based on the volume of inbound air travel from the source country). The EII is calculated similarly, based on the number of asymptomatically infected individuals in Nigeria and their likelihood of travel to each destination by air. Notably, we could not obtain data for land-based travel. Connectivity to regional and neighbouring countries and the associated introduction and export intensity are therefore likely severely underestimated. A time series of reported deaths for each country were collected from the outbreak.info R package⁵³. As air travel data was aggregated by month, we conservatively assumed that travel was uniform across all days of the month. We quantified the III and EII for the period May 2020 to April 2021, encompassing the emergence and spread of B.1.525 and B.1.1.318.