Idli batter was prepared from parboiled rice (Oryza sativa) and dehulled black gram (Phaseolus mungo) as described as before (Chandrasekar Rajendran et al., 2017). W. paramesenteroides strains W31, B1, and C8 were isolated from the overnight fermented batter on De Man, Rogosa and Sharpe (MRS) agar plates supplemented with 1% sucrose using dilutions and plating methods (Sanders, 2012). Strains SJ21-9, SJ27-4, and SJ27-5 were isolated from air dried onion powders (Säde et al., 2016). All isolates were grown in MRS liquid media supplemented with 1% sucrose.

For genome sequencing, total bacterial DNAs from 2 ml cultures of W. paramesenteroides isolates grown overnight in MRS at 30°C were extracted using MagAttract HMW DNA kit (Qiagen, Hilden, Germany) according to the supplier’s protocol.

Conventional spot-on-lawn method was used to test the antimicrobial activity of W. paramesenteroides strains. The indicator strain Micrococcus luteus NCIMB 86166 was grown in LB broth at 37°C overnight. Four hundred microliters of the indicator culture was mixed with 1 ml of sterile water, poured onto a LB agar plate and let dried. Ten μl of Weissella overnight cultures were spotted onto the surface of the M. luteus indicator plate, and the plate was incubated overnight at 28°C.

Genomes of the six Weissella paramesenteroides isolates were sequenced at the Institute of Biotechnology (University of Helsinki, Finland) using next-generation sequencing platforms. Genomic DNA (0.5 mg) was sheared using a Bioruptor NGS Sonicator (Diagenode) to approximately 600 bp fragments. The fragments were polished, A-tailed, and ligated to a TruSeq truncated adapter. Purification of the ligation reaction was done using AMPure XP beads (Agencourt, Beckman Coulter). PCR of the libraries were done using Phusion Hot Start II DNA Polymerase (Thermo Fisher Scientific, Waltham, MA, USA) and index P7 primers and full-length P5 adapter primers. The reactions were pooled and purified with AMPure XP beads. Size selection of the pool was done according to Lundin et al. (2010). The obtained library pool was paired-end sequenced on a MiSeq Sequencer using the v3 600 cycle kit (Illumina, San Diego, CA, USA).

The genomes were assembled using the SPAdes assembly pipeline v3.15.5 with default setting (Prjibelski et al., 2020). Plasmids were identified using the plasmid option of SPAdes (Antipov et al., 2016). The quality of the assemblies was determined with Quast (v5.2.0, Mikheenko et al., 2018) using default settings and the genome of W. paramesenteroides STCH-BD1 as reference. Contigs with length shorter than 500 bp were discarded.

Genomes of the six newly sequenced W. paramesenteroides strains were deposited in GenBank under the accession numbers listed in Table 1. The annotation was performed using the assembled DNA sequences of the six new draft genomes from these isolates. The genomes were run through an automatic annotation pipeline PGAP (NCBI Prokaryotic Genome Annotation Pipeline) (Tatusova et al., 2016), followed by manual curation in few cases.

To find out the phylogenetic relation between our isolates and other W. paramesenteroides strains, we reconstructed the phylogeny using all available 34 W. paramesenteroides genomes retrieved from NCBI genome (dataset downloaded by latest 6th July 2022, Table 1) together with the six newly constructed genomes. The pangenomes of all 42 W. paramesenteroides strains were analyzed using Roary pipeline version 3.13 (Page et al., 2015). First, the nucleotide fasta files and GFF files from the published 34 strains were downloaded from NCBI and converted to GFF3 using a script (Supplementary Datasheet 1). For strains which GFF files were missing in NCBI or could not pass the preliminary analysis in Roary pipeline, their nucleotide fasta files were annotated using Prokka (v1.14.6) to obtain a compliable GFF3 file (Seemann, 2014). The 42 GFF3 annotations were then fed to Roary to calculate the pangenome of the dataset and produce a multiple sequence alignment (MSA) of the core genes (present in > 99% strains) using MAFFT. The obtained MSA of the core genome was used to generate a best-fit maximum likelihood phylogeny using IQTREE version 1.6.12 using ModelFinder optimization (Kalyaanamoorthy et al., 2017). Finally, the trees were visualized in an R package ggtree (Yu, 2020). A python script roary_plots.py¹ was used to show the presence and absence of core and accessory genes. All genomic fasta files were quality controlled to eliminate shorter sequences under the length of 300 nucleotide prior to Average Nucleotide Identity (ANI) analysis. The ANI of all genomes was calculated using a python module pyani v0.2.7 (Pritchard et al., 2016).

To predict the presence of bacteriocin encoding genes, 42 W. paramesenteroides genomes were analyzed by online prediction tool BAGEL4.² Any resultant genomic area of interest from BAGEL4 was analyzed further using BLASTP (Altschul et al., 1997). Virulence and the presence of any antibiotic resistance genes were predicted by VirulenceFinder 2.0³ and in RGI tool,⁴ respectively. Prophage sequences were searched using PHAST (Zhou et al., 2011), while the CRISPRs and Cas genes developed by bacteria for phage defensing were searched on CRISPRCasFinder (Couvin et al., 2018). To further understand the function of genomes even with novel sequences, the core genome of the W. paramesenteroides was annotated through eggNOG-mapper v2.⁵ The metabolic pathways were reconstructed using KEGG mapper.⁶

In this study, we sequenced the genomes of six W. paramesenteroides isolates. The genome of strain W31 was complete, while the other five draft genomes contain 2 to 36 contigs. For all the assemblies, 1−3 contigs together constituted more than 50% of the genome size. The six newly assembled genomes had a mean size of 2.035 Mbp, from 1.92 Mbp of strain C8 to 2.1 Mbp of the three SJ strains. The G + C content of these strains was 37.9%, except W. paramesenteroides W31, which contains higher GC percentage of 38.2% (Table 1). An over 75-fold coverage was obtained for the six newly sequenced genomes. One plasmid was identified from W. paramesenteroides strains B1 and W31 each, and strain C8 seemed to carry three plasmids. Since the genome assemblies of the three SJ strains contained more fragmented contigs, plasmidSPAdes could only detect short contigs from their genomes as putative plasmids, while the plasmid status of these detected contigs were still questionable (Supplementary Table 1). To confirm the existence of extrachromosomal DNA of these strains, plasmid isolation, and sequencing would be needed. Closing gaps between contigs e.g., using hybrid assembly would also provide a better view on the plasmids of these newly assembled genomes. Having plasmid DNA is an uncommon feature of W. paramesenteroides, as among the previously published 34 W. paramesenteroides genomes only W. paramesenteroides STCH-BD1 and WpK4 contain plasmids with the latter having extraordinary seven plasmids. The genomes of two strains DmW_118119 and DmW_118 were excluded from several further analyses, due to their failure in taxonomy check by NCBI and share less than 50% symmetric identities with the other W. paramesenteroides strains.

The six new genomes contain coding sequences (CDSs) ranging between 3389 and 3469, with an average of 3,407 genes. To identify the orthologous genes, all the CDSs were compared by blastp by all-against-all strategy. The pangenome of W. paramesenteroides without the outliers DmW_118 and DmW_118119 included 4,526 orthologous genes and the core genome contained 1,245 genes in addition to the softcore genome (138 genes) which together take the share of 30.56% of the total genes. The accessory genome includes 1,153 shell genes and 1,990 cloud genes (Figure 1D). The cumulative gene curve of 40 W. paramesenteroides pangenome fit the Heap’s law (n = κN^γ). The value of γ was calculated to be 0.85 (< 1.0) for the 40 genomes, using the result of the present and absent genes from roary pipeline, with the heaps () function in R package micropan (Snipen and Liland, 2015). Therefore, W. paramesenteroides pangenome was considered open (Figure 1A; Tettelin et al., 2008). The numbers of accessory genes in individual strains and numbers of unique genes varied between genomes (Figure 1C), while the number of unique genes stabilized below 900 (Figure 1B).

The roary matrix based on the presence and absence of the genes among 40 W. paramesenteroides genomes was generated and is shown in Figure 2A. Two genomes of strains DmW_118119 and DmW_118 were excluded from the matrix here. Nevertheless, we still performed the phylogenetic analysis on all 42 strains (Supplementary Figure 1), in which W. paramesenteroides DmW_118 and DmW_118119 apparently stood out from the rest of the Weissella strains, indicating their phylogenetic differences from the other strains.

The best assembled genomes of strain W31 in our study is closely related to and share the immediate evolutionary origin with strain WpK4. Two other isolates W. paramesenteroides SJ27-4 and SJ27-5 are close to each other, as expected from their origin of isolation, which also is reflected the relatedness of SJ21-9 strain to them.

A comparative analysis of genomic sequences was carried out for the 40 W. paramesenteroides genomes excluding DmW_118 and DmW_118119. Of the 1,245 genes of the core genome, 1,164 were successfully scanned by the eggNOG-mapper and ascribed to 20 COG families (Figure 3). Functional annotation of the COGs of the core genome showed that the majority encoded proteins responsible for the information processing systems associated to translation, ribosomal structure, and biogenesis (11.9%), and for amino acid and nucleotide metabolism and transport (7.8% each). The core genome also contained a large number of genes of unknown function (230, 19.8%) and 28 genes which could not be assigned to any functional category in the tested database (2.4%).

The metabolic pathways of the core genome were reconstructed using KEGG gene annotation predicted by eggNOG-mapper. The core genome of W. paramesenteroides is heavily involved in biosynthesis of secondary metabolites, cofactors, and amino acids. As shown in Supplementary Table 2, 96 orthologous enzymes are responsible for biosynthesis of secondary metabolites. W. paramesenteroides strains utilize 39 core proteins for amino acid biosynthesis, while dedicating an additional 87 orthologous enzymes to amino acid metabolic pathways. Among the 44 essential proteins needed for the synthesis of cofactors, these Weissella strains seem to be particularly good at biosynthesis of folate.

As shown in Figure 4, the pairwise comparison between all the 42 publicly available W. paramesenteroides genomes revealed that the 40 true W. paramesenteroides genomes had ANI scores 97.7% higher than the threshold for species demarcation (95%, Kim et al., 2014). Two strains DmW_118119 and DmW_118 significantly differ from the other strains, with only 84.4% ANI, as already noticed in the phylogenetic analysis. Interestingly, ANI analysis between strains DmW_118, DmW_118119 and all ten publicly available Weissella hellenica genomes could possibly demarcate the two outlier strains to W. hellenica, as their ANI were over 99% when compared with six of the W. hellenica strains (0916-4-2, 1.2.50, CBA3632, MBEL1842, PC1A, and Wikim14, Supplementary Table 3). Species identification using 16S rRNA sequences also agreed with this observation (data not shown).

The presence of prophages in all W. paramesenteroides genomes was investigated using search tools PHAST. All strains contained 2−7 phage sequences (Supplementary Table 4). Strains FDAARGOS_414, DRD-142, DmW_118, and DmW_118119 did not carry any intact prophages, but the others contained at least one intact and a few incomplete or questionable prophages per genome. Intact prophage regions spanned from 18.8 kb in strain DRD-187 genome to 57.9 kb in DRD-137 genome. The strain DmW_107 carried extraordinarily seven prophage regions, of which two were intact harboring necessary genes for phage assembly. Strains C8, 3.2.24, A37, N43, and STCH-BD1 carried only intact prophages without any incomplete phage fragments.

Since the prophages were shown to be prevalent among all W. paramesenteroides genomes, we wondered whether the phage defensing system CRISPR/Cas was also widespread among those bacteria. Nevertheless, according to the online CRISPRCasFinder tool, none of the 42 W. paramesenteroides strains showed strong evidence (evidence level three or four) of the presence of CRISPR sequences and cas genes in their genomes (Supplementary Table 5).

Weissella share similarities with Leuconostoc species, which have different roles in foods, for instance as dairy starters, meat spoilage bacteria, and meat protective cultures. The last-mentioned role is based on the production of antimicrobial peptides, bacteriocins. Previously, we had isolated and identified 50 Weissella strains from spontaneously fermented idli batter. When testing their antimicrobial potential, we found that W. paramesenteroides strain W31 kills other Gram-positive bacteria. As exemplified on M. luteus indicator plate with ten Weissella idli isolates (Figure 5A), unlike other slightly bactericidal Weissella strains 33, 38, and 39, W31 produced larger inhibition zone indicating the production of a bacteriocin (Figure 5A). For this reason, in addition to the strain W31, we screened all 42 W. paramesenteroides genomes for the presence of bacteriocin related genes using BAGEL4. The distance of DmW_118 and DmW_118119 from the rest of Weissella strains was also shown in the presence of bacteriocin related genes, as based on BAGEL4, only these two strains carried bacteriocin genes. Both strains harbored identical cassette for the production of a leucocin A/sakacin P family class II bacteriocin, including the genes encoding bacteriocin precursors, a putative immunity protein and an ABC type transporter with a bacteriocin secretion accessory protein (Figure 5C and Supplementary Table 6). In the cassette, there are also a Blp family class II bacteriocin, an MFS transporter and a PTS sugar transporter subunit IIA.

BAGEL4 did not find bacteriocin cassette in antimicrobial W. paramesenteroides W31. However, when comparing the bacteriocin sequences from DmW_118 and DmW_118119 to the proteome of W31, a similar 45-aa leucocin A/sakacin P class II bacteriocin with 23-aa double glycine signal peptide was found but annotated as hypothetical protein (80.95% identical, excluding the signal peptide, Figure 5B). When this putative W31-bacteriocin was compared with other W. paramesenteroides proteomes, identical or nearly identical proteins were found from all of them, including the idli strains B1 and C8, which have not shown any antimicrobial activity. All these putative bacteriocins shared a YWGNGV consensus motif, resembling the YGNGV “pediocin box” motif of class IIa bacteriocins. The 297-bp gene immediately downstream of the leucocin A/sakacin P family bacteriocin encoding gene (Figure 5C) was also found from other Weissella genomes and shared 59% identity with that from DmW_118, DmW_118119, and W31. The location and the size of this gene suggested its role as a bacteriocin immunity protein. Other bacteriocin related genes, e.g., those encoding bacteriocin transporters, were not found from the rest 39 strains.

The genomes of all 42 strains were screened for antibiotic resistances in RGI online tool. All showed “Strict” hits to vanT gene from glycopeptide resistance gene cluster vanG. Three SJ strains also returned “Perfect” hits on mefF major facilitator superfamily (MFS) antibiotic efflux pump, which causes macrolide antibiotic resistance, and “Strict” hits on msrG msr-type ABC-F protein, macrolide and streptogramin antibiotic resistance. DmW_118 and DmW_118119 had the Strict hits to qacJ, encoding small multidrug efflux pump (Supplementary Table 7). Since some Weissella spp. like W. cibaria and W. confusa have caused infections to human, we also screened the pathogenicity of all the 40 W. paramesenteroides strains with VirulenceFinder. But none of the genomes was found to carry genes associated with virulence.