Study period and location
Samples were collected once a month from March 2021 to October 2021 (active surveillance), apart from August 2021 when sampling was not possible for security reasons. Samples were collected from five households, one livestock market and one transhumance location in both Bassa and Jos South local government areas (LGAs) in Plateau State in northern Nigeria (Fig. 1a). These LGAs had been identified as being at high risk of FMD based on serological testing of samples from small ruminants [5] and FMD outbreaks reported in 2020 [21]. To be eligible for recruitment, households had to raise both cattle and small ruminants (sheep and/or goats) and agree participate in the study. Eligible households were identified and selected with input from local contacts in each LGA. Participation was voluntary and no incentives were given to take part in the study. Transhumance locations were defined as a location where herdsmen settle for up to two weeks before continuing their journey. There is only one livestock market and transhumance site in Jos South. Livestock markets and transhumance sites in Bassa were selected based on location, access, and agreement from people in charge to take part in the study. For transhumance sites, animal availability at the time of the sampling was also a factor. In addition, samples were collected, by the same field teams, from households reporting clinical outbreaks in the study area and in neighbouring LGAs during the same study period (outbreak investigation).
In both Bassa and Jos South, the households that were sampled changed between the March and April visits, after which samples were taken from the same households on all subsequent visits (i.e. May-October). Results are presented for all households, but statistical analysis was only carried out using the results from samples collected between April and October.
Animal sampling
During the first household visit nine animals (3 sheep, 3 goats and 3 cattle) were selected systematically for sampling. If fewer than three animals of a species were present, the total number was completed by sampling other species. At subsequent visits, the same animals were sampled, if possible, but this was often not the case as some animals were sold or slaughtered during the study period. At livestock markets and transhumance sites, animals were selected at each visit from various locations within the market or site. Each animal sampled was given a unique ID number which was linked to the site where the animal was held during the month of the visit. In addition, the age and sex of animals sampled were collected. Clinical examination was conducted on the animals sampled by a qualified veterinarian in the field team, and animals showing FMD-like lesions on the day of the visit were recorded.
Five millilitres of blood were collected from the jugular vein of each selected animal using pre-labelled vacutainer tubes (Becton Dickinson, USA). A sterile swab stick (SkyHealth, China) was used to swab the oropharyngeal/oral cavity of an animal and immediately put into PBS, as the transport medium. Once an animal had been sampled it was identified with a line in the ear using a non-toxic colour marker to avoid double sampling. All samples were kept at 4°C and sent to the National Veterinary Research Institute (NVRI) in Vom as soon as possible. On arrival at NVRI serum was separated into Eppendorf tubes clearly labelled with the unique animal ID number. Individual serum samples and oral swabs were stored at -20°C and shipped on dry ice to The Pirbright Institute, UK for testing.
Environmental sampling
At each sampling site, electrostatic dust cloths were used to swab areas of the environment where contact with secretions and excretions from infected animals was deemed likely (e.g. food troughs, hard floor surfaces, boots and tether ropes, transport vehicles and herder’s sticks). Up to 10 environmental samples per site per visit were collected. Each environmental sample was given a unique ID number which was linked to the site, place from which the sample was collected and month of the visit. The environmental samples were processed in the field by eluting the swabs in PBS and then adding the samples directly into lysis buffer (MagMAX Core or RLT buffer, Thermo Fisher Scientific, UK) at a ratio of 1:1. All samples were stored at 4°C and shipped on dry ice to The Pirbright Institute for testing.
Sample processing
Environmental samples, serum samples and oral swabs were tested for the presence of FMDV RNA by rRT-PCR. Viral RNA was extracted from samples using the KingFisher Flex automated extraction platform (Thermo Fisher Scientific, UK) with the MagMAX™ CORE Nucleic Acid Purification Kit (Thermo Fisher Scientific, UK). FMDV RNA was detected by rRT-PCR on the ABI 7500 system (Applied Biosystems, UK) using an assay that targets the 3D region of the FMDV genome [22]. Serum samples were also tested for the presence of antibodies against FMDV non-structural proteins (NSP). Samples were heat inactivated (at 56°C for 30 minutes, using a heat block) before testing with a PrioCHECK FMDV NS Antibody ELISA kit (Thermo Fisher Scientific Prionics AG, Waltham, MA, USA). Kits were used as per the manufacturer’s instructions.
Sequencing
Selected FMDV-positive samples were subjected to probe-enriched Illumina-based next generation sequencing. The selected samples had a mean RT-PCR CT value of 27.2 (range: 20.9 to 32.9; Additional file 1). First and second strand synthesis of total nucleic acid was performed as described previously using the SuperScript™ double-stranded cDNA synthesis kit (ThermoFisher) [23]. Libraries were prepared following the Illumina DNA prep with enrichment protocol. The first stage libraries were pooled and subsequently enriched using a library comprising 26,275 unique biotinylated oligos designed using 622 complete genomes available in GenBank [24]. Final libraries were diluted and run on the Illumina MiSeq using a V2 300 nano cartridge generating 2 × 150 paired-end reads.
Reads were initially assembled into contigs using SPADES, and each contig in turn was queried against a database of FMDV sequences to identify FMDV-specific contigs. Each FMDV contig was further subjected to BLAST online to identify the closet related sequence available. The 1D region of each most closely related virus was selected and used as a reference sequence for the reference assembly using BWA-MEM of each sample. Consensus sequences were extracted following reference assembly using VCF tools. The evolutionary history was inferred by using the Maximum Likelihood method based on the Tamura 3-parameter model [25]. The tree with the highest log likelihood (-2172.95) is shown. The percentage of trees in which the associated taxa clustered together is shown next to the branches. Initial tree(s) for the heuristic search were obtained automatically by applying Neighbor-Join and BioNJ algorithms to a matrix of pairwise distances estimated using the Maximum Composite Likelihood (MCL) approach, and then selecting the topology with superior log likelihood value. A discrete Gamma distribution was used to model evolutionary rate differences among sites (5 categories (+ G, parameter = 0.5492)). The tree was drawn to scale, with branch lengths measured in the number of substitutions per site. The analysis involved 20 nucleotide sequences. All positions with less than 95% site coverage were eliminated (i.e. fewer than 5% alignment gaps, missing data, and ambiguous bases were allowed at any position). Evolutionary analyses were conducted in MEGA7 [26].
Statistical analysis
McNemar's test for paired data and the kappa statistic were used to assess whether the results for oral swabs and serum samples tested by rRT-PCR from the same animals differed significantly from one another.
The proportion of positive samples by rRT-PCR or by ELISA was analysed using a binomial generalised linear mixed model with a logit link function. Because of the marked difference in the number of positive rRT-PCR results between LGAs, data from Bassa and Jos South were analysed in separate models (Additional file 2). By contrast, the number of positive ELISA results were similar in both LGAs and, accordingly, the ELISA results were analysed in a single model (Additional file 2). In each model the response variable was whether or not a sample was positive. Explanatory variables considered in the models included the LGA, month the sample was taken, sample type and species as fixed effects (see Additional file 2 for more details), and sampling location as a random effect. Models were implemented in a Bayesian framework using OpenBUGS (version 3.2.3; https://www.mrc-bsu.cam.ac.uk/software/bugs/openbugs/). Diffuse priors were used for model parameters: normal with mean 0 and standard deviation 10 for regression coefficients; and exponential with mean 100 for the random effect variance. Two chains of 60,000 iterations were run with the first 10,000 samples discarded to allow for burn-in of the chains. Chains were subsequently thinned by selecting every fifth iteration to reduce autocorrelation amongst the samples. Convergence was monitored visually and using the Gelman-Rubin statistic in OpenBUGS. Models including different explanatory variables were compared using the deviance information criterion (DIC) [27].
The force of infection (i.e. the rate at which susceptible individuals acquire infection) was estimated from age-seroprevalence data for each household. Specifically, the age and serological status of each animal (i.e. negative or positive by NSP ELISA) at its first sampling was extracted from the data-set (Additional file 1 and Additional file 2). To simplify the analysis, only home-bred animals (n = 306) were included in the analysis; those bought-in to the household (n = 18) were excluded. The force of infection was estimated using a catalytic model [28], so that
\(p(a)=1 - \exp ( - \lambda a),\)
where p is the seroprevalence at age a and λ is the force of infection. Heterogeneities in the force of infection were incorporated to allow it to vary amongst species and households/LGAs. In this case, the force of infection for an individual is given by
\(\log \lambda =\alpha +{\beta _{SPP}}+{\gamma _{LGA,HOUSE}},\)
where α is the baseline, β is the effect of species and γ is the effect of LGA and household. Here γ is treated as a random effect and drawn from a normal distribution with mean 0 and standard deviation σγ. Parameters in the model were estimated in a Bayesian framework using OpenBUGS (version 2.3.2). A Bernoulli likelihood was used with diffuse priors for model parameters (normal with mean 0 and standard deviation 10 for α and β; exponential with mean 100 for σγ). Two chains of 30,000 iterations were run with the first 5,000 samples discarded to allow for burn-in of the chains. Chains were subsequently thinned by selecting every fifth iteration to reduce autocorrelation amongst the samples. Convergence was monitored visually and using the Gelman-Rubin statistic in OpenBUGS. Models including different explanatory were compared using the DIC (Additional file 5). The final model was checked by comparing the observed data to the posterior predictive distribution [29].