Hybrid Genome Sequencing and Comparative Analysis of Three Novel Listeria monocytogenes Strains: Insights into Lineage Diversity, Virulence, Antibiotic Resistance, and Defense Systems

TL;DR

This study sequenced and analyzed three new Listeria monocytogenes strains to understand their genetic diversity, virulence, and resistance traits.

Contribution

The study introduces three novel L. monocytogenes strains with complete genomes and identifies unique virulence and resistance features.

Findings

Phylogenomic analysis showed a conserved core genome but diverse accessory genome among the strains.

Strain A2D10 contains the LIPI-3 operon, suggesting a hypervirulent phenotype linked to a fatal clinical case.

All strains showed similar antimicrobial resistance profiles and diverse defense systems.

Abstract

Listeria monocytogenes is a major foodborne pathogen, responsible for severe listeriosis outbreaks associated with contaminated foods. This study reports the comparative genomic analysis of three novel L. monocytogenes strains C5, A2D9 and A2D10, obtained from dairy and clinical sources. Hybrid genome sequencing with Oxford-Nanopore and Illumina technologies provided high-quality complete chromosomes. Phylogenomic analysis revealed a highly conserved core genome alongside accessory genome diversity. Strain C5 belonged to sequence type ST2, while A2D9 and A2D10 were assigned to ST155 and ST1, respectively. All strains exhibited close genomic relatedness to isolates from dairy animals and/or the dairy environment. Functional analysis identified conserved metabolic functions across all genomes. A total of 40 virulence genes were detected, including the LIPI-1 island in all strains and the LIPI-3 operon exclusively in A2D10, indicating a potential hypervirulent phenotype consistent with its ST1 background and the associated fatal clinical outcome. All strains exhibited similar antimicrobial resistance profiles typical of L. monocytogenes and diverse defense systems. The newly sequenced strains provide a valuable resource for functional analyses of the mechanisms underlying adaptation of L. monocytogenes to diverse environments.

Figures

Full text

1. Introduction

Listeria monocytogenes is a Gram-positive, facultative intracellular pathogen that causes listeriosis, a severe foodborne disease with high mortality rates among vulnerable populations such as neonates, pregnant women, elderly, and immunocompromised individuals [1,2,3]. This bacterium thrives in different environments, and it is primarily transmitted to humans through the consumption of contaminated foods, such as soft cheeses, smoked fish, and ready-to-eat (RTE) products. It has the ability to survive in a variety of food environments due to its ability to tolerate cold, high salinity, and disinfectants [1,4]. While cases of listeriosis remain relatively rare, it is of major public health concern due to the severity of clinical outcomes and the challenges associated with its control in the food chain [1].

The genomic architecture of L. monocytogenes reflects its ecological adaptability and pathogenic versatility. Its genome consists of a circular chromosome approximately of 2.9 Mb in size, with a high coding density and relatively conserved gene order across strains [5]. However, extensive comparative genomic analyses have revealed substantial interstrain variability, largely attributed to mobile genetic elements such as prophages, genomic islands, and plasmids. These elements contribute to genome plasticity, influencing traits such as stress tolerance, antimicrobial resistance, and virulence [5].

L. monocytogenes strains are classified into four major evolutionary lineages (I–IV), four molecular serogroups (IIa, IIb, IIc and IVb) and 12 serotypes (1/2a, 1/2b, 1/2c, 3a, 3b, 3c, 4a, 4b, 4c, 4d, 4e and 7), which are associated with distinct, but overlapping ecological niches and clinical outcomes [6,7,8]. Lineage I includes strains with the highest virulence potential [1]. Most strains belong to lineages I and II. Serotypes 1/2b and 4b (lineage I) and 1/2a (lineage II) are frequently associated with human listeriosis cases [6,8]. Lineage II strains, particularly serotype 1/2a, are also commonly isolated from food and farm environments and are responsible for cases of listeriosis in both humans and animals [8,9]. In contrast, lineages III and IV are predominantly associated with animal hosts and are not common [7,8,10].

Comparative studies have demonstrated lineage-specific differences in genome content, including the presence or absence of key virulence factors, metabolic genes, and regulatory systems. For example, the master regulator prfA, listeriolysin O (hly) and internalin genes (inlA, inlB) are central to host invasion and intracellular survival, but their sequence and expression can vary significantly between strains, influencing pathogenic potential [11,12]. L. monocytogenes also exhibits an intrinsic antimicrobial resistance (AMR) background characterized by the presence of core genes such as fosX, norB, mprF, lin, and mdrL, suggesting resistance to multiple antibiotic classes [13,14]. L. monocytogenes can harbor different classical antiviral defense systems like restriction-modification (R-M) and CRISPR-Cas [5,15,16,17]. Multiple defense systems can also be combined in a strain-specific manner [18].

Despite the availability of many sequenced genomes, questions remain regarding the evolutionary forces shaping genome reduction, the role of horizontal gene transfer, and the extent of genomic conservation across environmental and clinical isolates. While some L. monocytogenes strains exhibit streamlined genomes with loss of accessory genes, others maintain large, mosaic structures enriched with horizontally acquired elements. Understanding these differences is essential to unravel the molecular mechanisms underlying adaptation and virulence [19].

In the present study, we conducted a comparative genomic analysis of three novel L. monocytogenes strains, comprising one isolate obtained from a dairy farm environment and two clinical isolates most likely originating from foodborne infections. Our objective was to generate high-quality, complete genome assemblies with minimal gaps and assembly errors using a hybrid sequencing approach that combined Oxford Nanopore Technologies (ONT) long reads with Illumina short reads. We also aimed to investigate in depth key genomic features, including the distribution of virulence-associated genes, antimicrobial resistance determinants, and other adaptive traits. These newly characterized strains expand the genomic landscape of L. monocytogenes and may support future comparative and functional studies.

2. Materials and Methods

2.1. Isolation and Antibiotic Resistance Testing

A total of three L. monocytogenes strains were used in this study. Specifically, L. monocytogenes C5 was previously isolated from the dairy farm environment as already described [20], while two additional clinical strains A2D9 and A2D10 were isolated from patients with listeriosis, most likely of foodborne origin (please see below). Strain A2D10 was associated with a fatal clinical outcome.

Strains were cultivated in Brain Heart Infusion (BHI) medium (Oxoid, Basingstoke, UK) at 37 ± 0.2 °C overnight, and genomic DNA was extracted with the NucleoSpin DNA RapidLyse kit (Macherey-Nagel, Düren, Germany). DNA integrity and approximate fragment size were assessed on agarose gels. Their antimicrobial susceptibility testing was assessed based on the European Committee on Antimicrobial Susceptibility Testing (EUCAST) using the standardized disk diffusion method [21].

2.2. Complete Genome Sequencing

The complete genomes of L. monocytogenes strains were obtained through a hybrid assembly approach using both ONT long-read and Illumina short-read sequencing technologies.

Long-read sequencing was performed on a GridION X5 Mk1 device (ONT, Oxford, UK) using R9.4.1 flow cells and chemistry (FLO-MIN106, SQK-LSK109, EXP-NBD104, EXP-NBD114). DNA isolates were measured by fluorometry (Qubit4, Thermo Fisher Scientific, Waltham, MA, USA) and spectrophotometry (NanoDrop One, Thermo Fisher) to obtain precise quantification and purity data. 600 ng of each bacterial sample was used as initial input for the library preparation workflow. Nanopore sequencing preparation followed standard ONT protocols including multiple cleanups. Appropriate bead ratios and long fragment buffer were employed to enrich for longer fragments. In brief, fragments were first end-prepped and then barcoded. Barcoded samples were pooled using equal DNA mass (100 ng) and ligated to sequencing adapters. Approximately 50 fmol of final library, based on fragment length estimation, was loaded onto the flow cell. Strains were sequenced to 166× (A2D10), 221× (A2D9) and 251× (C5) coverage with Nanopore, allowing for selection of the longest, highest quality reads for assembly.

Short-read sequencing was conducted by Novogene (Cambridge, UK) using the Illumina NovaSeq 6000 platform (Illumina, San Diego, CA, USA) with 2 × 150 bp paired-end reads, yielding greater than 300x coverage for all bacterial genomes. Library preparation followed standard workflows including adapter ligation, sequencing and quality control.

2.3. Bioinformatics Analysis

A dedicated snakemake [22] workflow (https://github.com/simakro/CobraHy) was implemented for efficient and reproducible data processing of each sample. The workflow takes long Nanopore and short Illumina reads as input to generate high quality hybrid assemblies and includes taxonomic classification, plasmid identification, genomic annotation and quality assessment. In detail Nanopore long reads were subjected to quality control and adapter trimming with NanoPlot [23] and Porechop [24], respectively. Illumina reads also underwent adapter trimming and quality control with FastQC [25]. Long reads were assembled using Flye [26], to generate the initial draft assembly, which was subsequently polished first with nanopore data using Medaka [27], then with Illumina short reads aligned with the bwa-mem2 2.2.1 software package [28], in three iterations of Pilon [29]. All contigs were checked for circularity and oriented using the Circlator tool [30]. Plasmid identification was performed with PlasClass [31]. Assemblies were taxonomically classified with GTDBtk [32], and the closest reference genome was automatically retrieved. The reference was used to guide annotation with Prokka [33], as well as assessment of completeness with BUSCO [34]. Two assemblies were generated for each strain using the described snakemake workflow.

To maximize contiguity and precision of the primary assembly, an optimized read selection algorithm was employed to choose nanopore reads with the highest length and quality until a coverage of 60x was reached, yielding minimum read lengths and qualities of 15 kb and Q15 (A2D9 and C5), and 12 kb and Q12 (A2D10), respectively. A secondary assembly, employing a lower size cutoff of 700 bp, was performed to facilitate detection of smaller plasmids potentially missed in the primary assembly due to selection for very long reads.

The generated primary assemblies and associated output of the snakemake workflow were used for further and in-depth analysis. Completeness and contamination of assemblies were assessed with CheckM2 v 1.1.0. with the UniRef100 v1. reference database [35]. Clustering of strains was conducted with the Type (Strain) Genome Server (TYGS) [36], including the L. monocytogenes EGD-e as reference strain (accession number NZ_CP023861.1). The genomes of strains were aligned and visualized using progressiveMauve [37] and Circoletto [38] with default BLAST parameters. Comparative and functional genomics were performed using EDGAR 3.5 to compute the pan-genome, core-genome and average nucleotide identity (ANI) [39,40]. Typing of strains was performed using the BIGSdb-Lm [2]. JSpeciesWS was used to identify the closest related genomes to each of the three L. monocytogenes strains [41] with tetra correlation search (TCS). Additionally, multilocus sequence typing (MLST) and core-genome multilocus sequence typing (cgMLST) were carried out on assembled genomes using the CLC genomics workbench v24.0.2 (Qiagen, Hilden, Germany), employing L. monocytogenes MLST [42] and cgMLST schemes. Reconstruction of metabolic pathways was achieved through BlastKOALA using the KEGG v116.0 database [43]. Clusters of orthologous groups (COG) functional categories were obtained with eggNOG-mapper v2 [44] against the eggNOG 5.0 database. Finally, genes related to virulence and antimicrobial resistance (AMR) were also identified with ABRicate v1.0.1 [45] against the ARG-ANNOT [46], CARD [47], MEGARes [48], and VFDB databases [49]. Only genes meeting the default cut-off criteria of ABRicate (sequence identity ≥80% and alignment coverage ≥80%) were considered present in the genome. Defense mechanisms were explored with the DefenseFinder [50], prophage regions using PHASTER tool [51] and genomic islands (GI) were predicted by IslandViewer 4 [52]. All tools were run with default parameters unless otherwise specified.

3. Results and Discussion

3.1. General Characteristics of the Novel L. monocytogenes Complete Genomes

According to CheckM2, all three assemblies displayed excellent completeness (>99.9%) and minimal contamination (<1%) (Table 1) with a GC content of 38%, which is consistent with previous reports for L. monocytogenes [53].

Further genome annotation using Prokka revealed highly consistent genomic features across the three L. monocytogenes strains (Table 2).

The strains contained between 2837 and 2971 predicted genes, of which 2769 to 2903 were identified as protein-coding sequences (CDS). In all strains, 67 tRNA genes and 18 rRNA genes were identified, indicating a conserved core genomic organization of RNA-coding genes. These results are in accordance with the previous studies of L. monocytogenes genomes [5]. Finally, no plasmids were identified in any of the strains.

3.2. Phylogenomic and Whole-Genome Comparative Analysis

Phylogenomic analysis of L. monocytogenes revealed the evolutionary relationships among the L. monocytogenes strains and the type strains from the TYGS database (Figure 1A). Genome alignments with Mauve suggested a uniform coverage across most genomic regions (Figure 1B). To complement these observations, Circoletto was employed to visualize pairwise genomic alignments between the reference and each of the other strains (Figure 1C). This analysis clearly indicated the close relatedness of EGD-e and A2D9, as well as C5 and A2D10. The ANI matrix results indicated that all pairwise comparisons exceeded 98%, confirming the pairs of strains showing closer relations as suggested by the Circoletto analysis (Figure 1D).

3.3. Typing of L. monocytogenes Strains

Initial typing with BIGSdb-Lm predicted the classification of strains in two serogroups. Strains C5 and A2D10 were most closely related to serogroup IVb, while A2D9 to serogroup IIa. L. monocytogenes strains belonging to serogroups IVb and IIa have been previously isolated from ready-to-eat foods of animal origin, raw meat and food production environments [6,54].

Further investigation using JSpeciesWS revealed that genomes of each of L. monocytogenes strains showed high similarity (Z-score ≥ 0.999) to L. monocytogenes originating from specific environments. Strain C5 displayed the strongest correlation with L. monocytogenes Scott A isolated from milk during a listeriosis outbreak [55] and ATCC 19117 isolated from sheep [56]. Strain A2D9 showed the closest genomic relationship to foodborne isolates L. monocytogenes Finland 1998 [57] and QOC2 [58], both linked to outbreaks associated with the consumption of dairy products. Finally, strain A2D10 showed the highest similarity to L. monocytogenes strains NTSN and F2365, which were originally isolated from sheep and cheese, respectively [59,60]. The close genomic relationship between the studied strains and those associated with the dairy environment suggests that they may have originated from this environment as well.

The strains were further compared using the MLST L. monocytogenes scheme with seven housekeeping genes (abcZ, bglA, cat, dapE, dat, ldh, and lhkA) [61] (Figure 2). L. monocytogenes C5 strain was assigned to sequence type (ST) 2 of clonal complex (CC) 2 [62]. Isolates belonging to ST2 have been detected in dairy, meat, and fish products [63,64,65], and clinical cases [66]. Strain A2D9 was assigned to ST155 of CC155 [61], which has been identified in fish, meat products and clinical samples [67,68]. L. monocytogenes A2D10 strain was assigned to ST1 of hypervirulent CC1 [59,65,69,70]. ST1 strains have been verified to be associated with clinical cases [61], as well as with food-related environments [71,72]. L. monocytogenes EGD-e was assigned to ST35 of CC9 as previously described [73]. Members of ST35 have been associated with food products of animal origin [69,74]. Finally, cgMLST analysis assigned L. monocytogenes C5 to cgST9209, A2D9 to cgST21252, and A2D10 to cgST32884, although these assignments were inconclusive. The inconclusiveness arises from the presence of several nearly identical cgMLST profiles in the database, each exhibiting minimal allelic differences and comparable k-mer similarity scores. As more L. monocytogenes genomes are deposited and assigned to cgMLST profiles in future updates, the increased database coverage may help refine these assignments.

The comparative genomic and MLST analyses revealed that all three newly sequenced L. monocytogenes strains share close genomic relationships with isolates previously recovered from dairy and other food-related environments. Although strains A2D9 and A2D10 were obtained from clinical cases, their clustering with foodborne and dairy-associated reference strains suggests a possible food-related origin or transmission pathway, consistent with the well-documented ability of L. monocytogenes to cross the food chain barrier and cause human infection.

3.4. Pan-Genome and Core-Genome Analysis

Comparative genomic analysis using EDGAR 3.5 was performed on the three strains and the reference genome. The results showed that the genomes were highly conserved, sharing a large proportion of common genes in the core-genome, while only a small fraction of genes were unique to individual isolates (the accessory genome) (Figure 3). The analysis revealed that the core-genome accounted for approximately 80% of all genes, indicating a very close genetic relationship between the isolates.

Additionally, the comparison between the core-genome remained relatively stable showing a slight decline with each added genome (Figure 4). Pan-genome continued to increase with the addition of more strains. These results indicate considerable genomic diversity among the isolates, supporting the concept of L. monocytogenes as a species with an evolving and expanding gene repertoire.

The Venn diagram displaying the distribution of orthologous genes indicated that 2594 genes were shared by all strains (Figure 5).

Surrounding regions denoted strain-specific or partially shared genes. For instance, strain A2D10 possessed 48 unique genes and A2D9 had 47, while L. monocytogenes C5 harbored 93 genes absent from the other isolates. The reference genome also contributed 87 unique coding sequences. Intermediate overlaps, reflected genes conserved in subsets of the strains, indicating functional divergence and variable genome content.

3.5. Functional Annotation

Functional analysis revealed that proteins from all three L. monocytogenes strains and the reference strain were assigned to similar KEGG pathways. Most genes were found related to metabolic functions, including carbohydrate, amino acid, energy, nucleotide and lipid metabolism, as well as genetic and environmental information processing and cellular processes (Figure 6A). The core- and pan-genome exhibited comparable distribution of functions to those of the individual strains. In contrast, the total functional distribution of singletons varied significantly. Similar patterns between strains, the core genome and the pan-genome were also observed in COG functional categories, while singletons again showed distinct profile (Figure 6B). Overall, both KEGG and COG annotations indicated that the main metabolic pathways and cellular processes are largely conserved across all four L. monocytogenes strains, while the presence of singletons likely reflects strain-specific adaptations to distinct environmental niches.

3.6. Virulence Potential

A total of 40 virulence-associated genes were identified in the L. monocytogenes strains (Table 3). Certain genes were detected across all four genomes, including key pathogenicity determinants such as prfA, plcA, plcB, hly, actA, and mpl, which constitute the L. monocytogenes pathogenicity island 1 (LIPI-1) essential for intracellular survival and spread within the host [11,12,75,76,77,78]. The genomes also carried multiple internalin genes (inlA, inlB, inlC, inlF, and inlK) mediating adhesion, invasion, and immune evasion [12,79], together with the surface adhesin vip gene whose product interacts with host Gp96 and promotes epithelial internalization [80,81]. In addition, the lntA gene, encoding a nuclear-targeted protein that modulates host immune responses through interference with the BAHD1 chromatin repressor [82], was also present in all strains. Further virulence factors involved in host interaction and colonization, including oatA, bsh, prsA2, iap, lap, lapB, and fbpA [83,84,85,86] were found to be present.

Strains A2D9 and EGD-e both carried the additional internalin gene inlJ, which promotes bacterial adherence [87]. Strain A2D9 also harbored pdgA, involved in immune evasion [88,89], as well as the aut and ami genes, which contribute to bacterial adhesion, invasion and binding to host cells [90,91]. The presence of these genes suggests an extended repertoire of surface- and cell-wall–associated pathogenicity traits, potentially enhancing the capacity of A2D9 to adhere, invade, and persist within the host.

The genome of A2D10 also harbored the complete lls gene cluster (llsA, llsB, llsD, llsG, llsH, llsP, llsX, and llsY), forming the L. monocytogenes pathogenicity island 3 (LIPI-3), which encodes the thiazole/oxazole–modified microcin listeriolysin S (LLS) [77,92,93]. This peptide, produced by hypervirulent L. monocytogenes strains, has been associated with severe listeriosis outbreaks and enables modification of the host gut microbiota, promoting efficient intestinal colonization and invasion of internal organs [93]. These findings are consistent with the fatal clinical outcome of this strain. Strain C5 did not harbor any unique virulence genes compared to other strains.

3.7. Antimicrobial Resistance

In vitro antimicrobial testing of all L. monocytogenes strains revealed identical susceptibility profiles. All isolates were resistant to cefoxitin, oxacillin, daptomycin, and moxifloxacin, exhibited intermediate susceptibility to ciprofloxacin and remained sensitive to ampicillin, benzylpenicillin, gentamicin, erythromycin, clindamycin, quinupristin/dalfopristin, linezolid, teicoplanin, vancomycin, tigecycline, and trimethoprim/sulfamethoxazole. Resistance to cefoxitin and oxacillin is recognized as intrinsic in L. monocytogenes, reflecting the low affinity of its penicillin-binding proteins for cephalosporins and certain β-lactams antibiotics [13]. Reduced susceptibility to daptomycin has also been reported in L. monocytogenes [94]. Similarly, intermediate or reduced responses to fluoroquinolones such as ciprofloxacin and moxifloxacin are consistent with previous observations [95,96], while the overall susceptibility to ampicillin, benzylpenicillin, gentamicin, erythromycin, clindamycin, linezolid, and vancomycin reflects the typical antimicrobial profile of this species [13,97].

Genome screening revealed the presence of intrinsic AMR genes fosX, norB, mprF, and lin (Table 4). These genes confer resistance to fosfomycin, fluoroquinolones, cationic antimicrobial peptides, and lincosamides, respectively [14,98]. The presence of fosX, norB, and mprF is consistent with the observed resistance to β-lactams, daptomycin, and fluoroquinolones [14,99].

Furthermore, all strains carried the mdrL gene, which encodes a multidrug efflux pump responsible for the export of macrolides, cefotaxime, heavy metals, ethidium bromide and benzalkonium chloride [13]. The combined presence of mdrL, fosX, norB, mprF, and lin highlights the conserved multidrug resistance profile of L. monocytogenes strains supporting the findings of our in vitro antimicrobial tests [13,100,101].

3.8. Defense Mechanisms and Miscellaneous Traits

A variety of defense mechanisms were identified across all L. monocytogenes strains (Table 5). Strain C5 possessed type I, II and IIG R–M systems, reflecting an enhanced ability to recognize and degrade invading foreign DNA [102,103,104]. Strain A2D9 carried a CRISPR–Cas Type IB system together with an Iet system composed of the ATPase and protease components [105], as well as a retron I-C defense module [50,106], suggesting an expanded defense repertoire. A type II R-M system was also predicted for strain A2D10 alongside an RloC protein, shown to function through tRNA cleavage [107]. Additionally, this strain harbored abortive infection system (Abi) [108,109] and defense-associated reverse transcriptase (DRT) systems, which are known to restrict phage infection and inhibit protein translation [109,110]. Finally, the R-M type II and Abi defense systems were also predicted for EGD-e strain.

All strains contained one questionable and one incomplete prophage region, except for EGD-e, which harbored two questionable regions. Additionally, strain C5 contained one intact prophage region. Furthermore, strains C5, A2D9 and EGD-e each carried four genomic islands, while strain A2D10 harbored two. The observed differences in defense systems, prophage regions, and genomic islands highlight the genomic plasticity of the L. monocytogenes strains and their ongoing adaptation through horizontal gene transfer. It should be noted that all findings of this study are based on bioinformatics annotation and have not undergone any experimental functional validation.

4. Conclusions

The comparative genomic analysis of three novel L. monocytogenes strains provided comprehensive insights into their genetic composition, functional potential, virulence traits, antimicrobial resistance, defense mechanisms and genomic plasticity. All genomes exhibited high completeness and a large core-genome, while the presence of unique genes highlighted genomic diversity. MLST analysis confirmed that the strains represented distinct sequence types. Functional annotation revealed that most genes were involved in fundamental metabolic processes, while the distribution of singletons suggested strain-specific functions that may contribute to environmental persistence or host adaptation.

Variation in virulence genes suggested strain-specific differences in pathogenic potential, most notably the presence of the LIPI-3 island in A2D10, which is associated with enhanced virulence and may help explain the severe clinical outcome linked to this strain. The conserved resistance and susceptibility patterns across all genomes indicated the typical antimicrobial profile of L. monocytogenes, supporting the use of ampicillin or benzylpenicillin as first-line therapy in clinical practice. The identification of resistance determinants such as fosX, norB, mprF and lin further confirms the current reliability of first-line therapies but underscores the necessity of ongoing genomic surveillance to detect any future emergence of high-risk acquired traits. Furthermore, the presence of multiple bacterial defense systems, prophage regions and genomic islands suggests mechanisms that may facilitate adaptation to diverse ecological niches, including food-processing environments and the human host.

The identified genomic features reveal key traits relevant both to clinical management strategies and the survival of L. monocytogenes in food-processing environments. Further investigation of these strains is warranted to elucidate the functional significance of the identified genomic features and their contribution to L. monocytogenes virulence and environmental adaptation.