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 tolerance differs between strains. The variance in colony forming units (cfu) of the treated samples was significantly larger (p<0.02, see ANOVA results below) compared to control, indicating a high variance in the level of tolerance 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 results were obtained 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).
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) (Figure 1a). The same statistical analysis for HPP at 600 MPa indicated that all strains including strain RO15 had significantly (p<0.01) lower cfu/ml in treated samples compared to controls (Figure 1b).
Genome sequencing, general features and RNA-Seq
We sequenced the genomes of seven L. monocytogenes strains for which no genome sequences were available in public databases. Genomes of strains RO15, 2HF33, MB5, AB199, AB120, C7, and RO4 were sequenced using Illumina MiSeq instrument and strain RO15 was sequenced using Pacbio RSII. Sequences of these seven strains were assembled and three additional strains (ScottA, F2365 and EGD-e, genome sequences of which are available in public databases) were used for comparative genome analysis. Assembly of PacBio long reads of strain RO15 resulted in one continuous 3,042,507 bp sized chromosome (Figure 2) at an average 308-fold 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). All genomes had a similar size, GC content and number of coding sequences (CDS) (Table 2).
The serotype of the strains was reported in previous studies [21, 22, 23, 24, 25] (Table 2). Although RO15 was reported as a serovar 4b strain using multiplex PCR in a previous study [21], 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 [24] were not found in the genome of strain RO15 based on our sequence match analysis, whereas sequences for primers targeting lmo0737, which are indicators of serovar 1/2a strains [24] were present. Multilocus sequence typing (MLST) based on 7 loci, core genome MLST (cgMLST; based on Moura scheme [26]), and clonal complex (CC) of the strains were also assigned based on the genome sequences (Table 2). Lineage assignment showed that only ScottA and F2365 belonged to lineage I and all other strains were members of lineage II (Table 2).
In addition to genome sequencing, transcriptome analysis using RNA-seq for strains RO15 and ScottA provided a basic view of transcriptional activity of the genomes. Differential expression analysis indicated that virulence genes and heat shock genes were upregulated after high pressure treatment in strain ScottA (Table S2).
Methylated DNA Motifs
High PacBio sequencing coverage allowed us to analyze DNA methylation modification motifs in strain RO15 in addition to assembling the genome. Using the long read modification data, 11 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 [29] 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.
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 the 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 significant (e-value <1E-50) BLASTP hits in the REBASE protein sequences database [29]. 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 (Figure 2). Prophage (region) 5 had the same sequence as the separate circular phage. RNA-Seq count data suggested that most of the genes of prophage 1,3 and 5 in RO15 were transcribed. Subsequently, the prophage prediction was performed for all strains, which showed that all strains had at least one or more prophage regions 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.
Genome alignment of the selected nine strains against the strain RO15 (Figure 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 strains ScottA and F2365 with other strains were lower than 95 (Table S3). To study how the large sequence gaps are localized in the other publicly available genomic sequences of L. monocytogenes, all 200 complete L. monocytogenes genomes in RefSeq database [30] were downloaded and aligned to strain RO15 (Figure S2). The visualization of the alignment showed that prophages 2, 3 and 4 were seen only in the minority of the L. monocytogenes strains, while prophages 1 and 5 were seen in the majority of the L. monocytogenes strains. The comparison indicated that these prophage regions were common locations of variance between L. monocytogenes strains.
Annotation of CRISPR systems revealed that the RliB-CRISPR system was seen in all the strains, which is in line with previous RliB-CRISPR system studies [31, 32]. Interestingly, CRISPR I or CRISPR II system genes were present only in strains RO15, 2HF33, and MB5 (Figure 4), which exhibited the lowest reduction in log10 cfu/ml with the 400 MPa pressure treatment (Table 1). Number of spacers in these strains were also significantly (p < 0.001) higher compared to barosensitive strains. The alignment of the spacer sequences of the CRISPR systems back to the genome itself revealed that only strains 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, which aligned to the addB gene (OCPFDLNE_02427) located in the chromosome encoding an ATP-dependent helicase.
Using prediction tools we detected anti-CRISPR genes in the prophage regions (Table 4) of strains RO15, 2HF33, MB5, RO4, and EGD-e. Homologues of all four previously annotated Listeria anti-CRISPR genes, i.e. acrIIA1, acrIIA2, acrIIA3, and acrIIA4 [33], 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) (Figure 2) and we observed expression of these anti-CRISPR genes in strain RO15.
As antibiotic resistance and pressure tolerance were linked in a previous study [19], we also searched antibiotic resistance genes in all strains to identify relation of antibiotic resistance genes and pressure tolerance. The same antibiotic resistance gene families were detected in all the 10 strains (Table S4). Multiple sequence alignment of the detected antibiotic resistance genes and the following average distance tree generation showed that amino acid sequence of the norB gene (encoding for a quinolone resistance protein) (OCPFDLNE_03068) was slightly different in barotolerant strain RO15 compared to the other strains (Figure S3). Similarly, for the lin gene (encoding for a lincomycin resistance protein) (OCPFDLNE_00980) amino acid difference was seen for strains 2HF33 and RO15 (Figure S3).
Roary pan-genome pipeline suggested that the 10 L. monocytogenes strains contain a total of 4820 orthologous gene clusters. Of these genes, 2247 were core genes found in all strains and 2573 genes were accessory genes found in at least one strain. The core genome was used for constructing a phylogenetic tree, which indicated that the two serovar 4b strains (ScottA and F2365) are closely related to each other and cluster separately from 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 (Figure 5).
To test any significant association between pressure tolerance and (clusters of) genes in the accessory genome, we performed a pan-genome wide association analysis using Scoary. Based on the pairwise comparison of reduction in log10 (cfu/ml) against strain RO15 using Student’s t-test (Table S5), a statistically significant difference compared to RO15 was seen in all strains except MB5 and 2HF33. Therefore, strains RO15, 2HF33, MB5 were used as barotolerant strains for the pan-genome association analysis. Scoary results showed that 13 gene clusters (Table S6) had p-values < 0.01 for the association with barotolerant strains, i.e., the genes were only seen in barotolerant strains. Of these, 10 gene clusters were located in the prophage regions, but most of the genes were annotated as hypothetical proteins (Table S6). These prophage genes were also searched for in 200 complete L. monocytogenes strains and found in some strains (Table S6). Interestingly, genes OCPFDLNE_02579 and OCPFDLNE_02580 were only seen in strains that harbored anti-CRISPR genes, Cas type I or II CRISPR system and self-targeting spacers (Table S6).