A total of 121 strains from lineage I and 104 from lineage II were included in the study. These strains come from our internal collection and have been phenotypically characterized and used in previous studies (Balandyte et al., 2011; Rupp et al., 2015; Dreyer et al., 2016). The strains were isolated by enrichment in Oxoid Novel Enrichment Broth at 30°C and subsequently growth on Brilliance Listeria agar (Oxoid, Ltd., Basingstoke, United Kingdom) at 37°C for 24 h. Single colonies suggestive for L. monocytogenes, were then transferred to Tryptic Soy Agar (TSA) containing 5% (v/v) sheep blood (BD, Becton Dickinson and Company, Sparks, U.S.A.) and incubated at 37°C for another 24 h. For a few strains, colonies presenting haemolysis on the TSA were then applied to the VITEK Compact 2 phenotypic analysis identification system, using Gram-positive identification cards (Biomerieux, Geneva, Switzerland) for the phenotypic identification of the species Listeria monocytogenes. For the rest of the strains, the species and lineage were defined according to Matrix-Assisted Laser Desorption Ionization-Time Of Flight Mass Spectrometry (MALDI-TOF MS) (Dreyer et al., 2016), after confirming the equivalence between the two methods. Two strains from lineage I (JF5203 and JF5861) and two strains from lineage II (JF4839 and LMNC088) were selected as internal reference strains (Table 1). The JF5203 strain belongs to sequence type (ST) 1 and was isolated from a rhombencephalitis case in cattle. The JF5861 strain belongs to ST4 and originated from a human CNS infection. The LMNC088 strain belongs to ST412 and came from the farm environment. The JF4839 strain belongs to ST9 and originated from food not related to any listeriosis outbreak.

Listeria monocytogenes strains that were re-sequenced in this study were grown overnight at 37°C on TSA supplemented with 5% (v/v) sheep blood (Becton Dickinson GmbH, BD^TMTrypticase^TM, PA-254053.07). Colonies were picked and directly treated with lysozyme (at a final concentration of 0.4 μg/μL). Thereafter, the bacterial cells were lysed in guanidium buffer (60% w/v) (Pitcher et al., 1989) and genomic DNA (gDNA) was extracted according to a previously published phenol-chlorofom-isoamyl alcohol protocol (Wilson, 1987).

Ninety-one of the total number of strains used in this study were sequenced in a previous study (Accession numbers PRJEB15123 and PRJEB15195; Dreyer et al., 2016). Some of them were re-sequenced to improve the data quality and coverage. The remaining 134 strains were sequenced specifically for this study (Table S1).

All strains were sequenced using the Illumina® technology (https://www.illumina.com/), either on MiSeq (300 bp paired-end reads) or HiSeq 2500/3000/4000 (95–150 bp paired-end reads) platforms, according to the manufacturer's protocols. Genome coverage varied from 19 x to more than 1000 x (Table S1). The four internal reference strains (JF5203, JF5861, JF4839, and LMNC088) were also sequenced using the Pacific Biosciences® (PacBio) technology (http://www.pacb.com/) in order to combine the high level of accuracy from the short-reads sequencing generated by Illumina technology with the long fragments from PacBio sequencing technology.

De novo assembly of PacBio data was done using HGAP v3.0 (Chin et al., 2013) from the SMRT® Analysis package v2.3.0. Quality control of the assembly was performed by mapping the Illumina reads to the obtained contigs and then performing an analysis with Qualimap v.2.2.1 (Okonechnikov et al., 2016). Circularization of the single contigs was carried out using the AMOS package v.3.1.0 (Treangen et al., 2011). Genomes were compared using the BRIG application v.0.95 (Alikhan et al., 2011) and ANI server (Rodriguez-R and Konstantinidis, 2016). Annotation of the whole genomes was made using Prokka v1.12 (Seemann, 2014) and MicroScope (Vallenet et al., 2013).

To detect variations between the two lineages, Illumina reads of all the sequenced strains were mapped using BWA v0.7.13 (Li and Durbin, 2009) to the whole genome sequence of JF4839 (internal reference strain of lineage II). This was done to obtain the same position in all genomes relative to the same position in the reference genome. Reads with quality values below 20 in Sanger scale (Phred+33) were excluded from the analysis using sickle (https://github.com/najoshi/sickle). The total number of genomic variants including short insertions/deletions (INDELs) and single nucleotides variants (SNVs) were identified per strain using SAMtools v0.1.19 (Li et al., 2009). For variant calling files (vcf) filtering and manipulation SAMtools and vcflib (https://github.com/vcflib/vcflib) were used. Variants with mapping and assertion quality values lower than 30 in Phred-scale and with less than 20 reads supporting the alternate allele were filtered out from the vcf files in individual genomes. A Mann-Whitney-Wilcoxon test (“stats” R-package) was performed to identify differences between the numbers of SNVs per lineage (Neuhäuser, 2011). The vcf files per strain were then combined into a single merged vcf file using VCFtools v0.1.14 (Danecek et al., 2011). Manipulation and annotation of the merged vcf file was done with VCFtools, SnpEff v4.3i, (Cingolani et al., 2012) and in house bash scripts.

A phylogenetic tree based on the multiple alignment of the SNVs found (ignoring heterozygous sites) was built using RAxML v8.2.9 (Stamatakis, 2014), with a generalized time reversible (GTR) substitution model and using a strain from lineage III as an outgroup. Bootstrap scores (350 replicates) were calculated. The tree was re-rooted to the outgroup genome using FigTree v.1.4.3 (http://tree.bio.ed.ac.uk/software/figtree/), edited and displayed adding metadata information with CLC Genomics Workbench v.9.5.2 (https://www.qiagenbioinformatics.com/products/clc-genomics-workbench).

The list of variants present in lineage I with respect to lineage II was filtered by excluding the SNVs with low impact according to SnpEff (http://snpeff.sourceforge.net/SnpEff_manual.html). The same approach as detailed above for read mapping and SNVs filtering was repeated with JF5203 (lineage I, CC1 and ST1) as the reference to identify SNVs only present in lineage II strains with respect to lineage I strains.

Using an in house python script the variants per gene were counted in the CNS-related strains, in order to look for genes with multiple SNVs. The BED file of the annotated JF4839 genome and a filtered merged vcf file containing the variants private to CNS cases (according to SnpSift) were used as input files.

Thirty-six published genomes with sufficient information about their lineages (18 belonging to lineage I and 18 to lineage II), along with the genomes of the four internal reference strains underwent pan-genome analyses (Table S2). Using the MicroScope platform (Vallenet et al., 2013), the core-genome of lineage I excluding the pan-genome of lineage II was calculated (at the protein level) for 80% sequence identity and 80% length coverage. As a result, the protein-coding genes shared by the 20 strains of lineage I but absent in any of the 20 strains of lineage II were predicted. A putative function was assigned to proteins with no described function according to InterPro (Apweiler et al., 2000) and BLASTp search (Altschul et al., 1990) against UniProtKB (Bairoch et al., 2005). The list was further filtered by taking into account the presence of certain amino acid motifs and domains potentially related to surface proteins and virulence factors (Bierne and Cossart, 2007).

To check for the presence of the previously (section Core- pan- genome analyses based on a reduced set of strains) selected genes present in all lineage I strains of our set but absent in all lineage II strains, we used the information recorded in the bam files (all Illumina reads of each strain mapped to the JF5203 genome) to calculate the Reads Per Kilobase per Million mapped (RPKM) values. RPKM analysis is an established method for RNA sequencing (RNA-seq) data examination (Mortazavi et al., 2008; Deng et al., 2012; Tonner et al., 2012). It allows the RNA-seq gene expression quantification by normalizing for total read length and the number of sequencing reads. RPKM values were obtained according to the following equation: RPKM= numReadsgeneLength1000*totalNumReads1000000

numReads is the number of reads mapped to a gene sequence, geneLength is the length of the gene sequence (in bases) and totalNumReads is the total number of reads mapped to the genome.

Two housekeeping genes, dnaA and gyrB genes were used as controls. The Mann-Whitney-Wilcoxon test (“stats” R-package) was used to check for differences between the RPKM values for each of the previously identified genes in the two lineages (Neuhäuser, 2011).

To corroborate that there were differences between the RPKM values only in the selected genes with respect to the controls, a pairwise comparison between elements of different lineages by calculating the RPKM-difference values was performed (subtracting each RPKM value of lineage II to each RPKM value of lineage I). The post hoc Dunn's test for the Kruskal-Wallis multiple comparison test was performed with “dunn.test” R-package v.1.3.4 to assess the significance among the groups (Dinno, 2017).

RPKM values were calculated in the 121 sequenced lineage I genomes and in the 104 sequenced lineage II genomes using the 2981 genes of the annotated JF5203 genome as references. Thereafter a similar procedure as described in section Reads Per Kilobase per Million analysis was performed to calculate the RPKM-difference values amongst the different lineages for all genes. The median of the RPKM-difference per gene was calculated and genes with a median greater than or equal to 233.77 (2^*standard deviation) were retained to obtain a list of genes predominantly absent in lineage II. A matrix was built based on the RPKM values of the selected genes. A heatmap was calculated and plotted in R (“gplots” R-package) (Warnes et al., 2016) using hierarchical clustering algorithms (average linkage clustering) based on Euclidean distance.

The same RPKM approach described above for the 2981 genes in the genome was used for comparing the CNS related strains with the non-CNS infection associated ones. Genes with a median greater than or equal to 103.38 were kept. These genes are predominantly absent in food, environmental and non-neurolisteriosis strains.

Finally, the method was applied to compare strains from the clinical group (AD) with the ones present in the non-clinical group BC (Section RPKM analyses at whole genome level and PCA analysis). Genes with a median greater than or equal to 154.01 were kept.

After every comparative analysis, the selected genes underwent a Gene Ontology (GO) enrichment step by Blast2GO (Conesa et al., 2005) and Interproscan v.5.2 (Jones et al., 2014).

A principal component analysis (PCA) of the RPKM values for the 2981 genes in the 225 strains was performed (“stats” R-package). A permutation multivariate analysis of variance (PERMANOVA) test (Anderson, 2001) was used to identify significant differences between the different clusters (“vegan” R-package) (Oksanen et al., 2017).

All the statistical analyses were done in R v3.3.2 (R Development Core Team, 2016) and all p-values < 0.0001 were considered as significant.

The four reference genomes with their annotations and the sequencing data for all the strains used in this study were submitted to the European Nucleotide Archive (ENA) under the Project number PRJEB22706 (See Tables S1, S3 for details).

Our first aim was to use WGS and comparative genomics to find characteristics distinguishing L. monocytogenes lineage I strains from lineage II strains. To this end, we selected two strains belonging to lineage I (JF5203 and JF5861) and two strains of lineage II (JF4839 and LMNC088) (Table 1) and obtained their entire genome sequences using PacBio and Illumina sequencing technologies. The one-contig assembly for each of the chosen reference strains was obtained with PacBio data and the quality control of the assemblies was performed using the short reads from Illumina. Only very few bases (less than 30) needed to be corrected. The sequences were further circularized to generate whole non-fragmented circular chromosomes. Each reference genome had a chromosome size of approximately 2.9 Mb and GC content of 38% which is in line with previously published Listeria genomes (Hain et al., 2006).

Additionally, episomal sequences were generated for strains JF5203 and JF4839. In strain JF5203, three low coverage contigs, named 1, 2, and 3, belonging to phages were sequenced concomitantly with the bacterial genome. In the samples sequenced by PacBio the coverage of these regions was 2.4–4.7 times below the genome coverage while, they were barely detectable in the samples processed by Illumina sequencing (Image S1). The phages identified by PHASTER (Arndt et al., 2016) are the following: in contig_1 and contig_2, two intact prophages LP-030-3 (GenBank accession number NC_024384.1) and vB_LmoS_293 (GenBank accession number NC_028929.1), respectively, while contig_3 contained an incomplete prophage. They are bacteriophages of the Siphoviridae family Orthocluster IV which have already been described as L. monocytogenes phages (Denes et al., 2014; Casey et al., 2015). Phages from this cluster are typically between 38 and 41 kb long and have GC contents of 35.5–36.6%. Indeed, the phages identified here are approximately 36, 33, and 3.5 kb in contig_1, contig_2, and contig_3, respectively and each have a GC content of between 35.3 and 37.1%. Most likely contig_2 and contig_3 are part of the same phage according to a Mauve comparison (Darling et al., 2004; data not shown). Therefore, the phages identified in JF5203 have approximately the same size as those previously reported for this cluster, as well as a similar GC content. Additional examinations to verify the contiguity of contigs_2 and 3 were not performed because it was not possible to re-isolate or re-identify the phages. Additionally, an incomplete prophage is also integrated into the chromosome of JF5203.

A 74 kbp plasmid was found in the lineage II strain JF4839, isolated from cheese unrelated to a listeriosis outbreak (Filiousis et al., 2009). The plasmid contains genes associated with metal transport and resistance to cadmium and camphor (Image S2). Cadmium is an important environmental pollutant and a potent toxicant to bacteria (Trevors et al., 1986). The metal transport and resistance genes are common in environmental strains to allow them to better adapt to the different environmental conditions. This plasmid shows 99% identity at the DNA level to the L. monocytogenes strain N1-011A plasmid (GenBank accession number NC_022045.1), representing approximately 79% of its length.

When using the ANI server (http://enve-omics.ce.gatech.edu/ani/) and Nucmer (Kurtz et al., 2004) to compare the reference strains from lineage I (JF5203 and JF5861) to the previously published strain F2365 from a listeriosis outbreak in California (Mascola et al., 1988), we found more than 99.6% identity at the DNA level (Image S3). On the other hand, a comparison between the reference strains from lineage II (JF4839 and LMNC088) to the EGD-e strain, showed more variation, resulting in 99% sequence identity (Image S4). The differences between the reference strains of lineage I and II were much larger representing 5.7% (Image S5).

The numeric summary of results of the annotation step using Prokka, as well as basic metrics of the genomes obtained are detailed in the Table 2. A similar number of internalin-like proteins were identified in all four reference sequences. Likewise, we examined the integrity and synteny of the LIPI-1 island in the four reference genomes and found a preserved co-localization and order of all the genes in this region (Image S6).

The Illumina sequencing data was used for the detection of variants that can distinguish between lineages I and II. Using the genome sequence of lineage II strain JF4839 as a reference to map all the strains in the study, the average number of SNVs in lineage II was 26,826 while in lineage I this value was 129,632 SNVs. The distribution of all SNVs in each lineage was represented in a kernel density plot (Figure 1, Table S4), showing more heterogeneity within lineage II (see the wider x-axis range). Significant differences with a p-value < 0.0001 were obtained between the two groups (Mann-Whitney-Wilcoxon test). The distribution of the total genomic variants was very similar to the distribution of the SNVs (results not shown).

A tree based on the number of SNVs at the whole genome level using JF4839 (lineage II) as reference is shown in Figure 2. A more distant strain, LMNC318 from lineage III, was used as an outgroup. In the resulting tree, three main branches are observed, corresponding to the different lineages. The topology of the tree confirms the clustering based on CC and ST classifications; interestingly a single branch of two strains corresponding to ST91 (CC14) are not grouped with other CC14 strains (see asterisks in Figure 2).

In lineage I, there are more clinical cases associated with CNS infections (75 compared to 30 cases in lineage II). One septicemia cases is also present in this group. In contrast, in lineage II, environmental strains are more common along with two strains of food origin. Other clinical manifestations such as mastitis, gastroenteritis and abortion are present in both groups. It is also relevant that lineage I has far less SNVs among them than lineage II, showing a closer distribution, while lineage II displays more diversity.

The total list of conserved variants in lineage I (present in ≥80% of the strains taking the JF4839 annotated genome of lineage II as a reference) allows the differentiation between the two lineages and all are documented in Supplementary Material DS1.

Addressing our next aim, the SNVs private to CNS infection cases were analyzed separately in order to find a pattern specific for strains of CNS origin. We looked for genes with high number of SNVs in this fraction of strains. However, we did not identify specific genes correlating to more variations in the neurolisteriosis strains. The maximum percentage of strains having common SNVs was approximately 20% (Supplementary Material DS2). Given this low percentage we think that the SNVs identified in our study are not related to neurovirulence in ruminants.

The SNVs analysis (calculation of number of SNVs per strain) was done in both directions because both lineages contain elements or genes not present in the other one. Thus, taking the JF5203 annotated genome (lineage I) as a reference, the SNVs were identified for lineage II strains. A list of variants of lineage II present in ≥80% of the strains with respect to lineage I was also created (Supplementary Material DS3).

Our whole genome variant analysis revealed variations in the prfA gene between lineage I and II. The PrfA protein is a master regulator essential for the activation of the transcription of many bacterial virulence factors within infected host cells. Specifically, we identified the presence of two substitutions of cytosine (C) in lineage II to thymine (T) in lineage I at positions 10 and 13 in the 5′UTR of the prfA gene in all of the 121 lineage I strains analyzed (Table S5). This specific untranslated region acts like a thermosensor in L. monocytogenes (Johansson et al., 2002).

Another pair of variants between the two lineages was found in the S-adenosylmethionine (SAM) riboswitch SreA. An Adenine (A) at position 83 in lineage II is substituted by a guanine (G) in lineage I, and a G at position 88 in lineage II is changed to an A in lineage I (Table S5). This SAM riboswitch participates in the negative regulation of PrfA translation, since it can bind and make and hybrid structure with the prfA transcript (Loh et al., 2009).

These two pairs of variants are given as interesting examples, however there are thousands of other potentially interesting variants to look at.

Genetic differences between the lineages I and II are not only due to single nucleotide differences, but also due a different gene composition. Therefore, in order to identify a group of genes shared by lineage I strains but not present in lineage II, L. monocytogenes genomes (20 of each lineage including our references strains) were analyzed in the MicroScope platform (http://www.genoscope.cns.fr/agc/microscope/home/). The lineage I differential pan-genome comprised a total of 5,730 genes (2,838 families according to MicroScope MICFAM parameters 80% amino acid identity and 80% alignment coverage) including the lineage I differential variable-genome of 5,290 genes (2,816 families) and the lineage I differential core-genome of 440 genes. This core-genome corresponds to 22 gene families specific to lineage I (Image S7). From these 22, a further filtering was done, taking into account the gene length (longer than 90 bp, according to Prodigal Hyatt et al., 2010) and the presence of certain domains/motifs (internalin-like domains, cell wall anchor protein domains, adhesion domains such as LPXTG). According to these criteria, we identified a reduced group of six genes specific for lineage I strains (LMOF2365_RS01905, LMOF2365_RS06250, LMOF2365_RS12245, LMOF2365_RS11140, LMOF2365_RS03470, LMOF2365_RS13380), which have the potential to be putative virulence attributes of L. monocytogenes (Figure 3).

In order to confirm the absence of these six selected genes in all the strains of lineage II and their exclusive presence in lineage I, we examined the remaining 221 sequenced strains. Specifically, RPKM values were calculated and compared between lineages (Figure 4A). Significant differences (p-value < 0.0001) were found between the lineages (Mann-Whitney-Wilcoxon test). All six genes were absent in lineage II and present in all lineage I strains analyzed, with the exception of the LMNC284 strain, where the gene LMOF2365_RS06250 is not present after verification with the Integrative Genomics Viewer (IGV) (Robinson et al., 2011; Thorvaldsdottir et al., 2013) (Image S8).

After RPKM pairwise analysis, box plots of the difference between lineage I and II gene by gene were generated (Figure 4B). For the control genes (dnaA and gyrB), the difference between the two lineages is close to zero whereas for the six selected genes, this difference increases to approximately 400, showing similar behavior for all the genes. A Dunn's test revealed significant differences (p-value < 0.0001) between the control group and the six selected genes.

The RPKM values-based analysis was extended to a whole genome level in order to explore the presence of the 2981 genes from the reference genome JF5203 in the 225 genomes of this study.

A heatmap showing the degree of presence/absence of the 167 genes predominantly absent in lineage II was created for the 225 strains (Figure 5, Table S6). For the gene filtering, a cut-off of two times the standard deviation of the RPKM-difference values was selected (Section RPKM analysis at whole genome level). In Figure 5, the heatmap shows that the two lineages are perfectly separated based on the genes selected and that the strains generally grouped according to the CC classification, except for the same two strains of ST91 (CC14) already mentioned in section Variant calling at whole genome level and phylogenetic relationship determination (Data not shown in the Figure 5; refer to Table S1 for CC/ST classification). Furthermore, it was possible to distinguish 28 genes that are only present in lineage I (Table 2). After Blast2GO analysis, only 5 proteins remain uncharacterized, while the remaining 23 are proteins of interest because of their classification (internalin-like proteins, cell wall anchor proteins, transcriptional regulators, ABC transporters). Notably, the majority of them are bacterial surface proteins.

Figure 6 shows the results of a PCA analysis with the whole RPKM matrix (RPKM values of all 2981 genes in all 225 strains). PCA 1 and 2 explains the 55% of the variance. Four groups are clearly defined. On one side the lineages are perfectly separated with significant differences (p-value < 0.0001; PERMANOVA test). On the other side, two other groups can be distinguished, one with the majority of the clinical strains (either from lineage I or II) and the other one with the majority of the environmental strains. Significant differences were also detected with the PERMANOVA test (p-value < 0.0001). Based on these results, the groups were defined as follows: A-lineage I clinical, B-lineage I non-clinical, C-lineage II non-clinical, and D-lineage II clinical. Only one strain (JF5593) was not classified as either clinical or non-clinical, since is located at the middle of the two groups.

In order to apply a clinically relevant filter to look for genes that could possibly be related to CNS infections, the RPKM method was applied, but this time to compute the differences of RPKM values between the CNS infection-related and non-related strains. We found that 77 genes are predominantly absent in the non-CNS group of the strains. According to our data, a single gene cannot perfectly separate the two groups and be assigned as CNS-infection causative, but the combination of various genes could be a signature of neurolisteriosis (Figure 7, Table S7).

There are 65 mutual genes from the 167 predominantly absent in the less virulent lineage II and the 77 genes predominantly absent in the non-CNS group. Of these 65 genes, 24 encode for membrane proteins, 5 for transcriptional regulators and the rest have other functions or are hypothetical proteins (asterisks in Table S6). These genes would be of further interest since are differentially present in strains that are pathogenic and have been associated to a CNS-infections.

We decided to do a comparison between the clinical AD group vs. the non-clinical BC group because of the PCA results and considering the fact that strains originating in the environment or food does not exclude the possibility of potential pathogenicity. In this case, the combination of the resulting 39 genes seem to be essential to visibly separate clinical from the non-clinical strains (Figure 8, Table S8).

Combining all the methods together and looking for common genes among the three comparisons, we have compiled a list of six genes that may represent interesting targets for future in vitro and in vivo studies to evaluate their neuroinvasion potential. These genes are: LMJF5203_02482 (transcriptional regulator), LMJF5203_01155 (ncRNA rli38) (Toledo-Arana et al., 2009), LMJF5203_00294 (ABC transporter ATPase), LMJF5203_00370 (uncharacterized conserved protein) and two short membrane proteins LMJF5203_00013 and LMJF5203_00470.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.