The study was conducted across six major geopolitical regions in Nigeria: South-South (Akwa-Ibom), South-West (Ile-Ife), North (Plateau, Kano and Bauchi), and South-East (Enugu). Further analysis in Plateau state was dropped due to insecurity as further sampling could not be done. Within each region, sampling points were strategically selected to encompass major farms known to utilise animal waste manure for fertilisation and wastewater or polluted water sources for irrigation purposes. About 30 samples were collected from farms not using such manure, because there were cow droppings within and around the farmland (inadvertent manuring), which obviously were expected to have impacts on the crops.

A total of 208 fruit and vegetable samples were aseptically collected. Samples included Water Leaf, Cabbage, Carrot, Pumpkin, Lady Finger, Lettuce, Amaranth, Ugwu, Ora, Onugbu, Nchuanwu, Garden Egg, and Water Melon. For each produce type at a given location, a composite sample was created by pooling sub-samples from at least five different sampling points (e.g., five different vendors or farm plots) to obtain a representative average for that location. Each composite sample was placed in a sterile stomacher bag, stored in a cool box at 4°C, and transported to the laboratory for analysis within 6 hours of collection.

Quantification of Listeria spp. was performed using a standard plate count technique following a ten-fold serial dilution. Briefly, 25 g of each composite sample was homogenised in 225 mL of Buffered Peptone Water (BPW) for 2 minutes using a stomacher. The homogenate was then serially diluted from 10 ⁻¹ to 10 ⁻⁶ in 0.1% peptone water.

From appropriate dilutions, 1 ml aliquots were spread-plated in duplicate onto two selective agars: Oxford Listeria Selective Agar (Oxoid) and Microinstant ^® Listeria Selective Agar (Scharlab). All plates were incubated aerobically at 37°C for 24-48 hours. For Southwest analysis, pre-enriched incubation in serially diluted Half-Fraser Broth and Full-Fraser Broth (OXOID, CM0895 ISO, United Kingdom) at 37 ^oC for 24-48 h. Characteristic colonies on both media were counted, and the results were expressed as Colony Forming Units per 25 grams of sample (CFU/25g).

The mean count for each produce-type and region combination was calculated from the duplicate plate counts of the composite sample. This mean value, derived from the two-culture media, was used as the definitive concentration (Mean CFU/25g) for all subsequent risk assessment calculations.

Selected bacterial colonies were purified through subculturing onto fresh agar plates to obtain axenic cultures essential for subsequent phenotypic and genotypic characterisation. ¹⁶

Presumptive Listeria isolates were subjected to a series of biochemical assays for preliminary identification. Initial characterisation included Gram-staining, which revealed Gram-positive, rod-shaped morphology. This was followed by a catalase test, which yielded a positive result, and an oxidase test, which produced a negative reaction. Haemolytic reaction in 5% sheep red blood agar was examined with Listeria monocytogenes and L. ivanovii expected to show β haemolysis, while L. innocua is expected to show γ haemolysis. This specific combination of results is consistent with the typical biochemical profile of the genus Listeria. ¹⁶

2.5.1 DNA extraction, visualisation, and quantification

Genomic DNA was isolated from the purified bacterial cultures using the ZR Fungal/Bacterial DNA Miniprep kit (Zymo Research, USA), following the manufacturer’s instructions. To confirm the success of the extraction and assess DNA integrity, the eluate was analysed by agarose gel electrophoresis. A 1% agarose gel, prepared with Tris/Borate/EDTA (TBE) buffer and stained with ethidium bromide, was loaded with the DNA samples alongside a molecular weight marker. Electrophoresis was performed at 100 volts for one hour, after which the DNA fragments were visualised under ultraviolet light using a UV transilluminator (BIORAD, South Africa). The concentration and purity (A260/A280 ratio) of the extracted DNA were subsequently determined using a NanoDrop 2000 spectrophotometer (Thermo Scientific, UK). ¹⁷

2.5.2 Amplification and sequencing of genetic materials

The 16S rRNA gene was targeted for amplification to enable precise species identification. Each 25 μL polymerase chain reaction (PCR) mixture contained 12.5 μL of Taq 2X Master Mix (New England Biolabs, USA), 1 μL of each 10 μM forward (27F: 5'-AGAGTTTGATCMTGGCTCAG-3') and reverse (1525R: 5'-AAGGAGGTGWTCCARCCGCA-3') primers, 2 μL of DNA template, and 8.5 μL of molecular grade water. The thermocycling conditions comprised an initial denaturation at 94°C for 5 minutes; 36 cycles of denaturation (94°C for 30 seconds), annealing (55°C for 30 seconds), and extension (72°C for 45 seconds); followed by a final extension at 72°C for 7 minutes before holding at 10°C.

2.5.3 Visualisation of amplified fragments

To confirm successful amplification of the target gene, the PCR products were electrophoresed on a 1% agarose gel as described in the previous section, using a DNA ladder to estimate the amplicon size.

2.5.4 Gene sequencing and data analysis

Amplified PCR products were purified in preparation for sequencing. A sequencing reaction mixture was prepared by combining 1 μL of purified PCR product with 9 μL of Hi-Di Formamide. The samples were then sequenced bidirectionally using the Sanger method on a 3130xl Genetic Analyser (Applied Biosystems, USA).

2.5.5 Sequence submission to GenBank

The obtained 16S rRNA gene sequences were deposited into the GenBank database to obtain accession numbers, ensuring the data is publicly available for future research. ¹⁸

2.5.6 Sequence alignments and phylogenetic investigation

The resulting forward and reverse sequence reads were assembled and edited to generate a consensus sequence using BioEdit software (version 7.0.5.3). This consensus sequence was used to query the GenBank database via the BLASTN algorithm on the National Centre for Biotechnology Information (NCBI) platform to identify the closest known species matches. For phylogenetic inference, related sequences were downloaded and aligned using ClustalW. A phylogenetic tree was constructed using the Maximum Likelihood method based on the Tamura-Nei model in MEGA 7 software (version 7.0.63). ¹⁹

A quantitative microbial risk assessment (QMRA) was conducted to estimate the probability of illness from Listeria monocytogenes associated with the consumption of contaminated fresh produce. The assessment followed the internationally recognised four-step framework: hazard identification, exposure assessment, dose-response assessment, and risk characterisation. ^{14, 20}

2.6.1 Hazard identification

The pathogen of concern was identified as Listeria monocytogenes, a facultative anaerobic, Gram-positive bacterium known to cause the severe foodborne illness listeriosis. The hazard was characterised by its ability to grow at refrigeration temperatures and its high mortality rate, particularly among immunocompromised individuals, pregnant women, newborns, and the elderly.

2.6.2 Exposure assessment

The exposure assessment aimed to estimate the number of L. monocytogenes cells ingested per serving of a given produce type. The initial concentration of Listeria spp. C_listeria_spp, measured in colony-forming units per 25 grams (CFU/25 g) as detailed in Section 3.4, served as the starting point. This value was converted to a concentration per gram and adjusted to reflect the concentration of the pathogenic species, L. monocytogenes (C_path), using Equation 1. C _ path ( CFU / g ) = [ ( Mean _ CFU / 25 g ) / 25 ] × F _ path Equation 1

Where: ✓

Mean CFU/25g is the mean count derived from the two-enumeration media (Oxford and Microinstant ^® Listeria) for each produce-region combination.

The daily exposure dose (D) was then calculated using Equation 2. D ( CFU / serving ) = C _ path × R × P × S Equation 2

Where: ✓

R is the reduction factor accounting for preparation steps (e.g., washing, cooking). A value of 1 was applied, representing a worst-case scenario where no reduction occurs before consumption.

2.6.3 Dose-response assessment

The probability of infection per serving was estimated using the exponential dose-response model ( Equation 3), which is commonly applied for pathogenic bacteria including L. monocytogenes. P _ infection = 1 − exp ( − r × D ) Equation 3

Where: ✓

P_infection is the probability of infection per eating occasion.

The probability of illness (P_ill) was subsequently derived by accounting for the proportion of infections that lead to clinical symptoms, as shown in Equation 4. P _ ill = P _ infection × P _ ill | inf Equation 4

Where: ✓

P_ill|inf is the conditional probability of illness given that infection has occurred. A value of 0.5 (50%) was selected to represent a highly susceptible population subset. ²⁶

2.6.4 Risk characterization

The final annual risk of illness was calculated by integrating the per-serving risk with the frequency of consumption throughout the year. The probability of illness after f exposures was calculated using the exponential annual risk model ( Equation 5). P _ ill _ annual = 1 − ( 1 − P _ ill ) f Equation 5

Where: ✓

f is the annual consumption frequency, assumed to be 12 times per year (once per month) for each produce type.

For the very low values of P_ill calculated in this assessment, Equation 5 can be approximated by P_ill_annual ≈ P_ill × f.

The antibiotic susceptibility profile of the selected Listeria species isolates was determined using the standardised disk diffusion technique (Kirby-Bauer method) as outlined by the Clinical and Laboratory Standards Institute (CLSI, 2020). A panel of seven antibiotics was selected to represent major classes of antimicrobial agents. These included amoxicillin-clavulanic acid (AUG), cefoxitin (CFX), ceftriaxone (CRO), ofloxacin (OFX), imipenem (IMI 10), vancomycin (Va 30), and oxacillin (Ox 5). Before testing, each bacterial isolate was revived through two successive subcultures on Brain Heart Infusion agar at 37°C for 24 hours to ensure purity and metabolic activity. A standardised bacterial inoculum was prepared by suspending several well-isolated colonies in sterile saline to achieve a turbidity equivalent to 0. 0.5 McFarland standard, approximately 1–2 × 10 ⁸ CFU/mL.

This standardised suspension was then evenly swabbed onto the surface of Mueller-Hinton Agar plates. Antibiotic discs of specified potencies were aseptically placed onto the inoculated agar surface, ensuring sufficient spacing to prevent overlapping zones of inhibition. The plates were incubated at 35°C ± 2°C for 16–18 hours. Following incubation, the diameters of the zones of complete inhibition were measured to the nearest millimetre using a calibrated digital calliper. The results were interpreted based on the clinical breakpoints provided by CLSI (2020). It is important to note that, as specific interpretive criteria for Listeria species are not available for all antibiotics, the breakpoints for cefoxitin and oxacillin were applied based on standards for Staphylococcus spp., a common practice in studies of Listeria susceptibility. For clear and conservative data analysis in this study, any result interpreted as ‘Intermediate’ was categorised as ‘Resistant’ to provide a more stringent assessment of non-susceptibility.

The resulting susceptibility data were analysed and visualised using the Python programming language, version 3.12.11. The overall resistance prevalence was depicted using a simplified stacked bar chart, which clearly illustrated the percentage of isolates classified as resistant or sensitive for each antibiotic, providing an immediate macroscopic view of resistance trends. To explore the patterns of multidrug resistance and relationships between the isolates, a hierarchical clustering heatmap was constructed. This analysis employed the Jaccard distance metric and average linkage method to cluster both the bacterial isolates based on their antibiograms and the antibiotics based on their co-resistance patterns. This dual dendrogram heatmap effectively identified clusters of multidrug-resistant strains and groups of antibiotics to which resistance was commonly linked, offering a sophisticated visual representation of the complex susceptibility profile data.

The detection of antimicrobial resistance genes ( vanA, vanB, tetL, and tetM) and the virulence gene hly (encoding listeriolysin O) in the identified Listeria species was performed using multiplex polymerase chain reaction (PCR) assays. The protocol for amplifying the resistance genes was adapted from the multiplex PCR method established by Iweriebor et al., ²⁷ utilising the primer sequences and thermocycling conditions detailed in Table 1. Similarly, the detection of the hly virulence gene was carried out via a uniplex PCR assay following the methodology described by Kumar et al., ²⁸ which is a well-characterised and specific protocol for identifying this critical virulence factor in Listeria monocytogenes.

The polymerase chain reaction (PCR) was performed in a final volume of 25 μL. This reaction mixture contained 12.5 μL of a commercial master mix (OneTag Quick Load, New England BioLabs), 5 μL of DNA template, and primers at concentrations ranging from 0.5 μM to 1.25 μM, with the remaining volume adjusted using nuclease-free water. The thermocycling protocol began with an initial denaturation and enzyme activation step at 94°C for five minutes. This was followed by 35 amplification cycles, each involving denaturation at 94°C for 60 seconds, primer annealing at 55°C for 60 seconds, and strand extension at 72°C for 90 seconds. To ensure complete elongation of any nascent DNA strands, a final extension was carried out at 72°C for 15 minutes. The PCR products were then visualized by resolving them on a 2% agarose gel via electrophoresis at 110 V for 45 minutes. Amplicon sizes were determined by comparison to a 100 bp DNA ladder.

PCR was conducted at 60°C with a gradient of 2°C, utilizing Taq polymerase at a concentration of 0.5–3.0 units, MgCl ₂ at 1–3.0 mM, and primers at 25–100 ng. The annealing duration was 60 seconds, while the extension times were set at 30 seconds and 1 minute. The total number of cycles performed was 40, employing four sets of primers in the multiplex PCR assay.

The results of the microbiological analysis, which quantified Listeria species contamination across various fresh produce types from different geopolitical regions in Nigeria, are summarised in Table 2. The table presents the average bacterial counts, expressed in colony-forming units per 25 grams (CFU/25 g), obtained from two selective culture media: Oxford Listeria (Oxoid) and Microinstant ^® Listeria (Scharlab). The final composite average count for each produce-type and region combination was calculated from these values. A total of 208 samples were analysed, of which 80 were positive for Listeria species, resulting in an overall occurrence rate of 38.5%. The data reveal substantial variation in contamination levels and frequency, with produce from Northern regions such as Kano and Bauchi generally exhibiting the highest bacterial loads, notably in samples of Lady finger (8.33 × 10 ² CFU/25 g) and Amaranth (7.43 × 10 ² CFU/25 g). In contrast, produce including Pumpkin and Ora from various regions showed the lowest levels of contamination.

Table 3 (a) and (b) present the identities and accession numbers of Listeria strains characterised in this study, with emphasis on their genetic distinctiveness and relatedness to previously reported strains. Table 3(a) shows the de novo sequenced Listeria monocytogenes and Listeria innocua isolates that were successfully deposited in the NCBI database, thereby confirming their novelty and genomic uniqueness. In contrast, Table 3(b) highlights Listeria monocytogenes isolates with varying percentages of identity to established reference strains, reflecting the degree of genomic similarity or divergence among the isolates. Collectively, these results demonstrate both the presence of unique local strains and the close genetic relationship of others to globally reported Listeria strains, underscoring the diversity within the isolates recovered in this study.

To elucidate the genetic relatedness and evolutionary relationships among the Listeria isolates collected in this study, a phylogenetic tree was constructed based on partial 16S rRNA gene sequences. The analysis included a total of thirty-two isolates recovered from various locations and sources across Nigeria. The sequences were aligned, and a phylogenetic tree was inferred using the Maximum Likelihood method. As shown in Figure 1, the isolates clustered into distinct clades, confirming their identification as Listeria monocytogenes and revealing their phylogenetic grouping. Notably, the tree also illustrates the genetic diversity and potential evolutionary paths among the strains, with bootstrap values indicating the robustness of the major nodes. The inclusion of a Listeria innocua strain (PX282439) served as an outgroup to root the tree and provide a point of evolutionary reference.

Figure 2 depicts the antibiotic resistance profiles of 156 Listeria strains across seven different antibiotics. The stacked bar chart shows the percentage of strains classified as sensitive (S) and resistant (R) for each antibiotic. Notably, Vancomycin and Oxacillin exhibit high percentages of resistant strains (83.1% and 88.1% respectively), while Ofloxacin has the highest percentage of sensitive strains (81.4%). This figure offers a clear visual summary of the resistance patterns observed in the tested Listeria isolates. While Figure presents the antibiotic profile, other subsequent figures present the pattern and relationships with other factors

Figure 3 displays a hierarchical clustering heatmap that visualises the antibiotic resistance profiles and the presence of specific resistance genes across the Listeria strains analysed. The heatmap employs a colour gradient where blue represents sensitive or gene absence, red indicates resistance, and reddish-orange signifies the presence of a resistance gene. The dendrograms flanking the heatmap illustrate the hierarchical clustering of both the Listeria strains (rows) and the antibiotic/gene markers (columns), highlighting groups with similar patterns. The x-axis is labelled with the antibiotic names and italicised resistance gene names, while the y-axis shows the individual Listeria strains. This figure facilitates the identification of correlations between observed antibiotic resistance phenotypes and the presence of the investigated resistance genes, revealing potential relationships and clustering patterns within the Listeria isolates.

Figure 4 presents a hierarchical clustering heatmap illustrating the antibiotic resistance profiles and the presence of specific antibiotic resistance genes in a collection of Listeria strains. The heatmap is coloured to indicate sensitive/absent (blue, 0), resistant (red, 1), and gene presence (reddish-orange, 2) for each strain and antibiotic/gene. Hierarchical clustering was applied to both the Listeria strains (rows) and the antibiotic/gene markers (columns) to group those with similar patterns of resistance and gene presence, as depicted by the dendrograms on the left and top of the heatmap, respectively. The x-axis labels include the names of antibiotics and italicised names of resistance genes, while the y-axis labels represent the individual Listeria strains. This visualisation allows for the identification of clusters of strains with similar resistance phenotypes and genotypes, providing insights into the relationships between antibiotic resistance and the presence of these specific genes within the Listeria collection.

Figure 5 presents a heatmap visualising the pairwise correlations between the antibiotic resistance profiles and the presence of specific resistance genes. The heatmap utilises a colour scale to represent the strength and direction of the correlation coefficients, ranging from strong negative correlations (cool colours) to strong positive correlations (warm colours). Each cell in the heatmap displays the calculated Pearson correlation coefficient between the corresponding antibiotic/gene pairs. For example, AUG (augmentin) shows a positive correlation with CRO (ceftriaxone) with a coefficient of 0.64, suggesting that resistance to AUG is often associated with resistance to CRO in this dataset. Conversely, AUG exhibits negative correlations with the antibiotic resistance genes vanA (-0.30), vanB (-0.21), tetL (-0.14), and tetM (-0.14), indicating that resistance to AUG is generally less likely when these specific genes are present. This visualisation provides a clear overview of the interdependencies among the measured resistance phenotypes and genotypes, enabling the identification of highly correlated variables that may suggest potential associations between specific resistance mechanisms and observed resistance patterns.

The results of this QMRA are presented in Table 4. The table provides a comprehensive overview, detailing the produce type, region, sample size, original enumeration data from both culture media, the calculated mean count, and the resulting risk estimates. The data reveal significant variations in the estimated risk, differing by over two orders of magnitude between the highest and lowest risk-producing items. This stratification of risk directly correlates with the microbial load data, effectively identifying high-risk produce-region combinations that warrant targeted intervention strategies.