Pressure treatment and reduction of viable cell counts
Reduction of viable cell counts (colony forming units, cfu) after pressure treatment at 400 and 600 MPa for 1 minute showed that pressure resistance differs between strains. The variance in cfu of the treated samples was significantly larger compared to control, indicating large variance in the level of resistance of L. monocytogenes. At 400 MPa for 1 minute, the 10 strains exhibited an average log10 reduction of 0.57 cfu/ml, ranging from 0.05 log10 cfu/ml for the most barotolerant strain (RO15) to 2.07 log10 cfu/ml for the most pressure sensitive strain (EGD-e). Similar variance was observed at 600 MPa. Here, the 10 strains exhibited an average log10 cell number reduction of 7.06 cfu/ml, with a range of 5.42 to 8.27 log10 cfu/ml (Table 1). Strain RO15 was also the most barotolerant strain based on an initial screening including several other L. monocytogenes strains (Supplementary text 1, Table S1, Figure S1).
Table 1
Reduction in viable cell counts after 1 minute 400 and 600 MPa pressure treatment. Table shows the average reduction in log10(cfu/ml) of selected L. monocytogenes strains after 1 minute at 400 or 600 MPa pressure compared to untreated control samples (n = 6 to 9 treated/untreated samples of per strain).
| Log10 reduction compared to control -Δlog10 (cfu/ml) |
Strain | 400 MPa | 600 MPa |
RO15 | 0.05 | 5.42 |
2HF33 | 0.15 | 6.97 |
MB5 | 0.22 | 8.27 |
AB199 | 0.26 | 6.95 |
AB120 | 0.35 | 7.30 |
C7 | 0.37 | 7.52 |
F2365 | 0.47 | 7.10 |
ScottA | 0.68 | 5.86 |
RO4 | 1.00 | 7.70 |
EGD-e | 2.07 | 7.49 |
One-way analysis of variance (ANOVA), comparing the pressure (400 MPa) treated samples of each strain to the control samples, showed that the mean cfu/ml for treated samples was significantly lower than the control for all strains (p < 0.02), except for strain RO15 (p = 0.15) (Fig. 1a). The same statistical analysis for HPP at 600 MPa indicated that all strains including strain RO15 had significantly lower cfu/ml in treated samples compared to controls (Fig. 1b).
Genome sequencing, general features and RNA-Seq
We sequenced and assembled seven L. monocytogenes strains (RO15, 2HF33, MB5, AB199, AB120, C7, and RO4). Sequenced seven strains and three additional strains (ScottA, F2365 and EGD-e, genome sequences of which are available in public databases) were used for comparative genome analysis. All genomes had a similar size, GC content and number of coding sequence (CDS) (Table 2). Strain 2HF33 had the largest genome with the highest number of CDS (3.11 Mbps; 3092 CDS) whereas strain RO4 had the smallest genome size with the least number of CDS (2.88 Mbps; 2806 CDS). As RO15 was defined as barotolerant strain, we sequenced the genome of L. monocytogenes strain RO15 using PacBio RSII long read technology. Obtained data was assembled using HGAP3 resulting in one continuous 3,042,507 bp sized chromosome (Fig. 2) at average 308X sequencing coverage. In addition to the chromosome, a contig with a complete, circular prophage sequence of 38,811 bp was obtained, which was also found as part of the chromosome (2729417–2770759 bp).
Table 2
General genome summary of selected L. monocytogenes strains. The table shows the genome size, number of predicted CDS and serotypes of selected strains.
Strain | Size, Mbp | #CDS | Serotype | MLST | Lineage | Genome Source | Source | Reference |
RO15 | 3.04 | 3022 | 1/2a | 155 | II | This study (Complete) | Herring with spices | [16] |
2HF33 | 3.11 | 3092 | 1/2a | 121 | II | This study (Contigs) | Salmon filleting section | [20, 21] |
MB5 | 2.97 | 2975 | 1/2a | 7 | II | This study (Contigs) | Salmon gutting machine | [20] |
AB199 | 2.99 | 2914 | 1/2a | 204 | II | This study (Contigs) | Drain/processing room | [16] |
AB120 | 3.00 | 2922 | 1/2a | 204 | II | This study (Contigs) | Sausage filler machine | [16] |
C7 | 3.09 | 3050 | 1/2a | 8 | II | This study (Contigs) | Salmon gutting machine | [20] |
F2365 | 2.91 | 2808 | 4b | 1 | I | GenBank: AE017262.2 | Mexican-style cheese | [17] |
ScottA | 3.02 | 2966 | 4b | 290 | I | GenBank: CM001159.1 | Clinical isolate | [18] |
RO4 | 2.88 | 2806 | 1/2a | 20 | II | This study (Contigs) | Dry cured salami | [16] |
EGD-e | 2.94 | 2875 | 1/2a | 35 | II | GenBank: AL591824.1 | Rabbit tissue | [19, 22] |
The serotype of the strains was reported in previous studies [16, 17, 18, 19, 20] (Table 2). Although RO15 was reported as a serovar 4b strain [16], our genome-based prediction suggests that it belongs to serovar 1/2a. Sequences for ORF2110 and ORF2819 primers used for identification of serovar 4b strains were not found in the genome of strain RO15, whereas sequences for primers targeting lmo0737 and prs, which are indicators of serovar 1/2a strains were present. Multilocus sequence typing (MLST) of the strains were also assigned based on the genome sequences (Table 2). Lineage assignment showed that only ScottA and F2365 belong to lineage I and all other strains were members of lineage II (Table 2, Figure S2).
In addition to genome sequencing, transcriptome analysis using RNA-seq for strain RO15 and ScottA provided a basic view of transcriptional activity of the genomes (Table S3).
Methylated DNA Motifs
High PacBio sequencing coverage allowed us to analyze DNA methylation modifications motifs in strain RO15 in addition to assembling the genome. Using the long read modification data eleven methylated sequence motifs were detected in the genome (Table 3). None of the detected motifs had a partner motif, i.e. a reverse-complementary sequence, but all detected motifs were only partially modified in the genome with less than 50% methylated motifs. While most of the motifs have already been deposited in the REBASE database [23] DADGYATYA, WNNTVVGCNTWNH, AHNBAACA, AGNNARNWW were novel, i.e. they have not been described as potential methylation sites previously. None of the detected motifs have been reported as recognition sequence motif for a restriction enzyme in the REBASE database.
Table 3
A summary of motifs. The table shows the detected methylated motifs and total number of motifs in the genome.
Motifs | Type | % Motifs Detected | # of Motifs Detected | # of Motifs in Genome | mean Coverage |
ADGYACYTV | m6A | 44.03% | 420 | 954 | 149.9 |
ADDTGGCA | m6A | 30.97% | 455 | 1469 | 148.5 |
TVVARARG | unknown | 22.20% | 1835 | 8265 | 148.8 |
ANNYASYA | m6A | 22.10% | 3289 | 14879 | 149.0 |
DADGYATYA | m6A | 21.02% | 326 | 1551 | 147.7 |
WNNTVVGCNTWNH | unknown | 18.14% | 582 | 3208 | 149.2 |
AHNBAACA | m6A | 13.36% | 839 | 6282 | 150.2 |
AGNNARNWW | m6A | 9.85% | 2225 | 22596 | 148.6 |
TNNNDNNH | unknown | 9.22% | 114413 | 1240313 | 148.9 |
TNNNCRVHNH | unknown | 7.59% | 6652 | 87598 | 149.0 |
TVNNNNNG | unknown | 3.26% | 7834 | 240556 | 149.6 |
We predicted one type II cytosine-5 DNA methyltransferase gene (OCPFDLNE_00657), and three type II N4-cytosine or N6-adenine DNA methyltransferase genes (OCPFDLNE_02168, OCPFDLNE_02626, OCPFDLNE_02808) in strain RO15 genome. In addition, genes for a type IV methyl-directed restriction enzyme (OCPFDLNE_00324) and type II restriction enzymes (OCPFDLNE_00658, OCPFDLNE_02625, OCPFDLNE_02807) were predicted. These predicted methyltransferases and restriction enzymes had high e-value (> 1E-50) BLASTP hits in the REBASE protein sequences database [23]. However, their recognition sequences are unknown. OCPFDLNE_02625 and OCPFDLNE_02807 type II restriction enzyme genes in prophage regions were expressed based on RNA-seq data.
Genomic comparison of the selected strains.
In the genome assembly we also identified a circular phage as an independent contig. While looking after phage originating genome parts prophage prediction revealed five prophage regions (10.7 kb to 47.9 kb sized) in the RO15 genome (Fig. 2). Prophage (region) 5 had the same sequence as separate circular phage. RNA-Seq count data suggested that most of the genes of prophage 1,3 and 5 in RO15 were transcribed. Following, the prophage prediction was performed for all strains, which showed that all strains had at least one or more prophage region in the genome with a maximum of six prophage regions in strain C7 (Table 4). One of the predicted prophage regions in ScottA was also transcribed based on the RNA-seq data.
Table 4
Comparison of CRISPR/Cas systems, anti-CRISPR genes, and prophage regions amongst the selected L. monocytogenes strains.
Strain | Average reduction of log10(cfu/ml) at 400 MPa | CRISPR/Cas systems | Number of spacers | Self-targeting spacer predicted | anti-CRISPR gene predicted | Number of prophage predicted |
RO15 | 0.05 | RliB-CRISPR, CRISPR I | 64 | yes | AcrIIA1, A2 | 5 |
2HF33 | 0.16 | RliB-CRISPR, CRISPR I, CRISPR II | 62 | yes | AcrIIA1, A2, A3, A4 | 5 |
MB5 | 0.22 | RliB-CRISPR, CRISPR II | 39 | yes | AcrIIA1, A2, A3 | 5 |
AB199 | 0.26 | RliB-CRISPR | 3 | no | no | 1 |
AB120 | 0.34 | RliB-CRISPR | 3 | no | no | 1 |
C7 | 0.37 | RliB-CRISPR | 7 | no | no | 6 |
F2365 | 0.48 | RliB-CRISPR | 3 | no | no | 1 |
ScottA | 0.68 | RliB-CRISPR | 3 | no | no | 3 |
RO4 | 1.00 | RliB-CRISPR | 5 | no | AcrIIA1, A2, A3 | 1 |
EGD-e | 2.07 | RLIB-CRISPR | 4 | no | AcrIIA1, A2, A3 | 2 |
Additional files |
Additional File 1: Supplementary Text 1. L. monocytogenes strain selection, including Table S1, Figure S1. |
Additional File 2: Figure S2. Maximum-likelihood tree shows lineages of strains. Concatenated sequences of the seven MLST gene fragments were used to create approximately maximum-likelihood tree to assign lineage of strains. Bold written strains were used in this study. The lineage information of other strains (not bold written) in this figure were shown in a previous study [64]. L. innocua strain Clip11262 was used as outgroup. |
Additional File 3: Table S2. Table shows Average Nucleotide Identity (ANI) (based on BLAST) calculation results of each strain. For strain ScottA and F2365 the ANI score with other strains were lower than 95. ANI score between ScottA and F2365 were higher than 95. For the other strains ANI scores between each other were more than 98. |
Additional File 4: Table S3. Gene count data for both strain RO15 and ScottA based on RNA-Seq. |
Additional File 5: Table S4 and Figure S3. Table S4 shows identified antibiotic resistance genes and percentage Identity of matching region in all strains using CARD database. Figure S3 shows average distance trees from alignments of protein sequences of identified antibiotic resistance genes. |
Additional File 6: Table S5. T-test against RO15 log10 reduction. |
Additional File 7: Table S6. Barotolerant strains specific genes based on pan-genome wide association analysis. Table shows genes and its annotation that have p < 0.01 based on pan-genome wide association analysis for barotolerant strains. Gene ID of strain RO15 was used for this table. Orthologs of these genes with 95% identity was also seen in strain MB5 and 2HF33. |
Genome alignment of the selected nine strains against the strain RO15 (Fig. 3) showed large sequence gaps between strain RO15 and the other strains. The gaps were mainly related to prophage regions with prophage region 3 being specific for strain RO15. Prophage region 2 was only seen in the strain 2HF33. Prophage region 4 was partially seen in strains 2HF33, EGD-e, MB5 and ScottA. Similarly, prophage region 5 was partially seen in strains 2HF33, C7, and RO4. In addition, strains ScottA and F2365 had less aligned regions compared to the other strains. Average nucleotide identity (ANI) results showed that ANI scores of strain ScottA and F2365 with other strains were lower than 95 (Table S2).
Annotation of CRISPR systems revealed that RliB-CRISPR system was seen in all the strains, which is in line with previous RliB-CRISPR system studies [24, 25]. Interestingly, CRISPR I or CRISPR II system genes were present only in strain RO15, 2HF33, and MB5 (Fig. 4), which exhibited the lowest reduction in log10 cfu/ml with the 400 MPa pressure treatment (Table 1). The alignment of the spacer sequences of the CRISPR systems back to the genome itself revealed that only strain RO15, 2HF33, and MB5 contained self-targeting spacers with 100% identity. As expected, spacer sequences were aligned to the prophage regions, except one spacer sequence in RO15. It aligned to the gene ATP-dependent helicase addB (OCPFDLNE_02427) located on the chromosome.
We predicted that strain RO15, 2HF33, MB5, RO4, and EGD-e had anti-CRISPR genes in the prophage regions (Table 4). Homologues of all four previously annotated Listeria anti-CRISPR genes, i.e. acrIIA1, acrIIA2, acrIIA3, and acrIIA4 [26], were seen in strain 2HF33. Gene acrIIA4 was not seen in strains MB5, RO4, and EGD-e, but the rest of the anti-CRISPR genes were present. RO15 contains two copies of acrIIA1 (OCPFDLNE_02770, OCPFDLNE_02583) and acrIIA2 (OCPFDLNE_02582) (Fig. 2) and we observed expression of the anti-CRISPR genes in strain RO15.
As antibiotic resistance and pressure resistance were linked in a previous study [14], we also searched antibiotic resistance genes in all strains to identify relation of antibiotic resistance genes and pressure resistance. Four same antibiotic resistance gene families were detected in all the 10 strains (Table S4). Multiple sequence alignment of detected antibiotic resistance genes and following average distance tree generation showed that amino acid sequence of norB gene (encoding for a quinolone resistance protein) (OCPFDLNE_03068) was slightly different in barotolerant strain RO15 compared to other strains (Figure S3). Similarly, for lin gene (encoding for a lincomycin resistance protein) (OCPFDLNE_00980) amino acid difference was seen for strain 2HF33 and RO15 (Figure S3).
Roary pan-genome pipeline suggested that the 10 L. monocytogenes strains contain a total of 4825 orthologous gene clusters. Of these genes, 2250 were core genes found in all strains and 2575 genes were accessory genes found in at least one strain. The core genome was used for phylogenetic tree construction. This indicates that the two serovar 4b strains (ScottA and F2365) are closely related to each other and located closer to the serovar 1/2a strains in the phylogenetic tree. In addition, there was also a clear difference in accessory genome for serovar 4b strains compared to serovar 1/2a strains (Fig. 5).
To test any significant association between pressure resistance and genes (clusters) in the accessory genome, we performed pan-genome wide association analysis using Scoary. Based on pairwise comparison of reduction in log10 (cfu/ml) against strain RO15 using Student’s t-test (Table S5), statistically difference compared to RO15 was seen in all strains except MB5 and H2F33. Therefore, strain RO15, 2HF33, MB5 were used as barolerant strains for pan-genome association analysis. Scoary results showed that 13 gene clusters (Table S6) had p-values < 0.01 for the association with barotelerant strains (genes were only seen in barotolerant strains). Of these, 10 gene clusters were located in prophage regions but most of the genes were annotated as hypothetical proteins (Table S6).