All completely sequenced Lactobacillus genomes (132) were obtained from the NCBI RefSeq database (O'Leary et al., 2016) on December 31, 2016. Genomic repeats were identified as described before (Koren et al., 2013) using Nucmer 3.1 (Kurtz et al., 2004). The number of repeats (longer than 500 bp, sequence identity of 95% or greater) was plotted vs. the maximal repeat length using Seaborn 0.7.1 (https://github.com/mwaskom/seaborn, Figure 1A); two strains were excluded as they were chimeric or could not be processed.

L. helveticus strains FAM8105, FAM22155, and FAM8627 were obtained from the Agroscope culture collection (Agroscope, Liebefeld, Switzerland) and grown under aerobic conditions in 10 mL MRS broth (De Man et al., 1960) at 37°C overnight (see Supplementary Figure 1 for light microscopy images). Bacterial cells were treated with lysozyme (50 mg/mL) for 1 h at 37°C. Genomic DNA (gDNA) was isolated as described elsewhere (Moser et al., 2017), and the concentration determined (Qubit dsDNA BR Assay kit). The gDNA purity was assessed by controlling for RNA contamination using the Qubit RNA HS Assay kit (Thermo Scientific, Massachusetts, USA).

The gDNA was sequenced on the PacBio RS II platform (three SMRT cells per strain, P6-C4 chemistry, size-selection step (10 kb inserts) with BluePippin; for details, see Supplementary Table 1). Subsequent de novo genome assembly using HGAP (Chin et al., 2013) and resequencing steps with Quiver were performed as described before (Remus-Emsermann et al., 2016). Terminal repeats were removed and the genome was circularized using Circlator 1.1.2 (Hunt et al., 2015). Additional rounds of sequence polishing resulted in one complete chromosome and one complete plasmid sequence per strain. FAM8105 and FAM22155 were also sequenced with Illumina MiSeq (paired-end, 2 × 300 bp), and the reads mapped to the polished PacBio assemblies using BWA-MEM (version 0.7.10-r789, Li, 2013). MiSeq reads which did not map to the respective chromosome and plasmid assemblies were assembled with SPAdes (Bankevich et al., 2012) to check for the existence of additional small plasmids.

We compared the genomes assembled only from Illumina MiSeq data vs. the complete PacBio assemblies (Table 2). The short read-based genome assemblies were generated with SPAdes (v3.11.0) (Bankevich et al., 2012), requiring a minimum read coverage cutoff of 4 and a minimum contig size of 400 bp (i.e., following recommendations of the SPAdes tutorial). The assembled contigs were mapped to the PacBio assemblies for FAM8105 and FAM22155 with BWA (option -x intractg), files were parsed and the mapping quality was assessed using Samtools (v1.3.1), Bedtools (v2.26.0), and Qualimap (v2.2.1). Subsequently, all coding sequences (CDSs) and pseudogenes in areas with zero coverage were counted and classified as core, accessory or unique gene (see section Comparative Genomics). The circular plot was created with Circos (v0.69-6) (Krzywinski et al., 2009).

The complete genome sequences were deposited at NCBI GenBank (Table 2) and annotated by their Prokaryotic Genome Annotation Pipeline version 3.3 (Tatusova et al., 2016). Putative CRISPR repeats (and their total number) were detected with CRISPRs finder (Grissa et al., 2007) and PILER-CR (Edgar, 2007). Putative Cas proteins were first searched in the NCBI annotation and secondly identified using a Hidden Markov Model (HMM) based search with HMMCAS (Chai et al., 2017). Putative prophages were identified using PHASTER (Arndt et al., 2016), potential genomic islands with Islandviewer 3 (Dhillon et al., 2015); see Supplementary Table 2 for their respective predicted genome positions. All CDSs (without pseudogenes) were compared against the EggNOG 4.5 databases “bactNOG” (bacteria), “bacNOG” (bacilli), “firmNOG” (firmicutes) (Huerta-Cepas et al., 2016), selecting only the hit with the smallest e-value (hits with e-values above 0.001 were not considered) and extracting the respective Clusters of Orthologous Groups (COG) category.

The GenBank records of 24 selected, completely sequenced LAB genomes (see Supplementary Table 3) were downloaded from NCBI RefSeq (March 3rd, 2017). This included 9 L. helveticus genomes, 14 other LAB genomes (mostly reference or representative complete genomes; 9 from genus Lactobacillus and 5 other LAB) plus Bacillus subtilis subsp. subtilis str. 168 as outgroup. The predicted CDSs of all strains, including the three L. helveticus strains of this study, were used to calculate a maximum likelihood phylogenetic tree using bcgTree (Ankenbrand and Keller, 2016), which uses HMM models of 107 known housekeeping genes (Dupont et al., 2012). bcgTree was parameterized to perform 100 bootstrap runs while executing RAxML (Stamatakis, 2014), the output (Figure 4) was generated using FigTree 1.4.2 (http://tree.bio.ed.ac.uk/software/figtree/).

DNA from a NWC was extracted as described elsewhere (Moser et al., 2017), DNA libraries were prepared (TruSeq DNA PCR-Free LT Library Prep Kit; insert size: 350 bp) and sequenced on an Illumina HiSeq 3000 (paired-end, 2 × 151 bp). Reads were assembled using SPAdes 3.9.0 (Bankevich et al., 2012). Performing a blastn search of the resulting contigs against NCBI RefSeq, we first determined the genomes (and plasmids, where available) with most hits, requiring that at least 95% of the raw metagenome reads were assigned. This implied 8 L. helveticus, 7 S. thermophiles, and one L. delbrueckii RefSeq strain. The reads were then mapped to these genome sequences plus our three FAM strains using BWA-MEM (version 0.7.15-r1140 Li, 2013; using option -a); the resulting SAM file was filtered to remove non mapping reads and supplementary alignments. To determine a qualitative species level distribution, we counted how many reads mapped to genomes of respective species; reads mapping to several targets were attributed in a proportional fashion. As few reads mapped to more than one species, we did not correct for different numbers of target genomes. As all genomes had roughly the same size we did not correct for genome size either. Next, we counted reads that mapped exclusively either to the three FAM or to NCBI L. helveticus genomes. To do this, the resulting SAM file from the previous steps was first parsed and a correction for the different number of target genomes (three FAM, eight NCBI strains) was performed. In a last analysis step, we determined how many of the uniquely mapping reads (i.e. reads just mapping to one target sequence) mapped to the FAM and NCBI L. helveticus genomes, respectively.

Aiming to increase the number of completely sequenced L. helveticus strains, we first explored the repeat structure of all Lactobacillus strains for which complete genome sequences have been deposited in NCBI's RefSeq. The total number of repeats and the length of the longest repeat of bacterial genomes are two parameters that have a profound impact on their genome assembly complexity (Koren et al., 2013).

Using an in-house software prototype, we calculated the overall number of repeats above 500 bp and >95% identity identified by Nucmer (Kurtz et al., 2004) vs. the length of the overall longest repeat (see Methods) for 130 completely sequenced Lactobacillus strains. This analysis resulted in the classification of LAB genomes into three classes with increasing demands for assembly due to the number and size of the repeats. Most genomes are classified as class I genomes (N = 57, 43.8%) that are straight-forward to assemble with PacBio long reads as they harbor few repeats, with the longest repeat representing multi-copy rDNA operons typically around 6–7 kb in length (Koren et al., 2013; Figure 1A). Furthermore, 32 (24.6%) class II genomes were observed that are characterized by the presence of a large number (more than 100) of repeats (500 bp to several kb), but none larger than the rDNA operon. Using data from the PacBio RSII platform, such genomes should also be straightforward to be de novo assembled into complete genome sequences; in contrast, relying on only Illumina short reads would produce tens to hundreds of contigs. Finally, we also noted a sizable fraction of class III genomes (41, 31.5%; Figure 1A). Due to their long, almost identical repeats well above 6–7 kb, these genomes can be extremely difficult to assemble. Such repeats can be resolved only with very long reads, which can be obtained by including a size selection step in the library preparation.

We next specifically assessed L. helveticus strains, for which nine complete genomes were available at NCBI's RefSeq (blue circles, Figure 1A). Eight of these strains were isolates from the dairy environment, while strain D76 was described as acting as an ingredient in nutritional supplements (Table 1). Eight of the genomes are class II genomes (Figure 1A), while strain H10 (Zhao et al., 2011) is classified as a class III genome and harbors long, nearly identical repeats >30 kb.

Finally, an analysis of the predominant assembly complexity classification for different Lactobacillus species indicated that L. helveticus strains had the second highest mean repeat number among all Lactobacillus species, outnumbered only by L. backii (Geissler et al., 2016) (Figure 1B). More importantly, due to the overall higher percentage of class III genomes, L. delbrueckii (73%), L. backii (60%), L. casei (50%), L. plantarum (43%), and L. fermentum (40%) represent Lactobacillus species that may pose considerable challenges for researchers aiming to carry out complete genome assembly projects similar to that we describe here for L. helveticus.

We selected three L. helveticus strains (FAM8105, FAM22155, and FAM8627) from the Agroscope culture collection, which originated from different dairy products (including raw milk and natural whey cultures, Table 1), and sequenced them on the PacBio RSII platform. To obtain a high sequence coverage, we used three single-molecule, real-time (SMRT) cells per strain, and to possibly even completely assemble class III genomes we used the BluePippin size selection protocol (see Methods). A high coverage would allow us to rely on algorithms that remove the random errors of PacBio reads (Chin et al., 2013) and to obtain completely sequenced genomes of high quality.

This de novo genome assembly approach resulted in one completely assembled chromosome and one complete plasmid for each strain, with a PacBio sequence coverage above 150-fold except for the plasmid of strain FAM22155 (Table 2). This is most likely because this plasmid is relatively small (7.5 kbp) and thus selected against in the BluePippin size selection step (see Methods). MiSeq short read data (available for FAM8105 and FAM22155) were used to check the quality of the assembly and to eliminate potential single nucleotide mis-assemblies reported for PacBio data (Laehnemann et al., 2015). Although no evidence for any mis-assembly was found for strain FAM22155, FAM8105 harbored a single nucleotide deletion in a homopolymer stretch in the chromosome and the plasmid (data not shown), which were corrected. This confirms that our assemblies are of a high quality. In addition, we used the MiSeq data to search for small plasmids potentially missed due to the size selection step during library preparation. However, no evidence for additional small plasmids could be found. Finally, a repeat analysis of the complete genomes of FAM8105, FAM22155, and FAM8627 classified all three as class II genomes (red circles, Figure 1A).

The genomes were annotated at NCBI (see Methods). A detailed summary of their genome statistics is shown in Table 2. An overview of their genome features, including CDSs, rRNAs, and tRNAs, is shown for chromosomes (Figure 2) and plasmids (Supplementary Figure 2). The distribution of functional COG categories (see Methods; Table 3) was similar for all three L. helveticus strains. However, compared to annotation projects we carried out in the past (data not shown), we observed a high number of CDS classified in category “L” (“Replication, recombination and repair”). This has been reported previously for the genus Lactobacillus (Lukjancenko et al., 2012) and can, at least in part, be attributed to the many transposase sequences in the L. helveticus genomes which get classified as class “L.” For FAM8105, 56.9% of the class “L” hits were transposases (194/341), for FAM22155, 56.3% (192/341) and for FAM8627, 47.6% (130/273) (Table 2). Accordingly, among the COG categories for which a function could be assigned, category L contained the largest number of genes (Table 3).

Despite the value of complete genome sequences, they are still highly underrepresented in public sequence repositories. This is also true for LAB: Among the 213 strains recently analyzed, almost all the genome assemblies (> 96%) were fragmented (Sun et al., 2015). The number of contigs reported for the LAB strains ranged from 8 to 964 (a median of about 90). A comparison of the FAM8105 and FAM22155 genomes assembled only from Illumina MiSeq data using SPAdes vs. the respective final PacBio-based reference genomes provided an insight into the extent of such differences (Figure 3): About 230 regions with an average size of roughly 1,000 bp were not covered by the MiSeq assemblies, which amounted to 10–12% of the actual genome of these two strains (genome size around 2.2 Mbp; Table 2). Notably, for the repeat-rich L. helveticus genomes, around 10% of the annotated CDSs would be missed in addition to the rRNA operons that represent the largest repeats in the genome, i.e., roughly 200 genes (Table 2). For pseudogenes, the percentage of missed cases was even higher, surpassing 20% of all pseudogenes.

To ensure reliable phylogenetic placement of the newly sequenced strains, we selected several taxonomically diverse LAB genomes on top of all completely sequenced L. helveticus genomes from NCBI (Supplementary Table 3). The maximum likelihood phylogenetic tree generated from these genomes is based on 107 known housekeeping genes and has very good bootstrap support (Figure 4). The three newly sequenced strains form a monophyletic clade with the other nine L. helveticus strains (Figure 4A), which is clearly separated from the other groups and is characterized by a relatively low degree of intra-clade variation. As reported previously, the genus Lactobacillus is not monophyletic, but paraphyletic (Mayr and Bock, 2002; Sun et al., 2015); i.e., the Lactobacillus clade also includes Leuconostoc, Oenococcus, and Pediococcus species.

The three FAM strains fall into two different subgroups within the L. helveticus clade (Figure 4B). FAM8627 is close to CNRZ 32 (Christiansen et al., 2008; Broadbent et al., 2013), a strain used as a starter culture and for the production of bioactive peptides in milk. FAM8105 and FAM22155 are more closely related to DPC 4571 (Hannon et al., 2003; Callanan et al., 2008), another strain known to be beneficial in cheese production. Of note, strains FAM22155 and CNRZ 32 seem to lack a CRISPR/Cas system as neither CRISPR repeats nor Cas proteins were detected (Figure 4B).

To the best of our knowledge and as noted in a recent review (Stefanovic et al., 2017), no pan-core genome study has been reported for L. helveticus. Thus, we carried out such an analysis on the 12 complete genomes using Roary (Page et al., 2015). The pan genome is generally defined as the sum of all genes in a species, whereas the core genome is defined as the orthologous genes that are present in all strains of a species (Medini et al., 2005). Furthermore, the accessory genome comprises orthologous gene clusters (orthologous genes from here on are referred to as “gene clusters” or “genes,” depending on the context) that, in our example, are found in at least two and up to 11 strains. Finally, we also determined genes that occur in only one of the strains; these genes represent the unique or “strain-specific” genome.

As more genomes were added, the size of the core genome diminished and reached 988 gene clusters for all 12 genomes (Figure 5A). The curve of the pan genome hints at a still “open” pan genome, comparable to results from a similar-sized pan-core genome study for 17 L. casei genomes (Broadbent et al., 2012). Overall, we identified a pan genome of 3,335 gene clusters (Figure 5B). They could be further divided into a core genome of 988 gene clusters (29.6%) present in all 12 strains, an accessory genome of 1,278 gene clusters (38.3%) present in a subset of the 12 strains (Supplementary Figure 3), and a total of 1,069 strain-specific genes (32.1%). The number of strain-specific gene clusters ranged from 29 genes only found in strain H9 up to 225 gene clusters unique to strain R0052 (Figure 5B). Notably, among the CDSs missed in the two MiSeq assemblies (Figure 3), not only accessory genes were missed, but also around 30 core genes, i.e., roughly 3% of the close to 1,000 core genes identified.

To identify whether gene clusters of certain functional groups are preferentially found in all 12 strains (core genome) or if they diverge more between the strains, COG functional categories for core, accessory and unique genome of all 12 L. helveticus strains were explored (Table 3; see Methods). As expected, the COG functional categories of some highly conserved biological processes, such as “Translation, ribosomal structure, and biogenesis” (Class J), and “Intracellular trafficking and secretion” (Class U), were enriched among the core genes and depleted in the unique genome. In contrast, the category K “Transcription” was equally present in the core, accessory and unique genomes, indicative of specific transcriptional regulatory mechanisms in some of the strains. For the class “Replication, recombination, and repair” (L), enrichment in the accessory genome was observed. As mentioned, this is mainly caused by transposases which are often found in more than one of the 12 strains but not in all. Enrichment among accessory and unique genes was observed for “Defense mechanisms” (Class V). Importantly, for the unique genes, we observed enrichment of genes in the COG functional category “Cell wall/membrane biogenesis” (Class M). Finally, genes of unknown function (Class S) were enriched in the unique genome, as well as genes not assigned to any COG class. For the latter, the largest difference between their presence in the core genome (very low) vs. the unique genome could be detected (Table 3).

For the production of Gruyère cheese, the use of NWCs, where L. helveticus is a predominant species, is regulated. For quality management, therefore, determining the composition of NWCs at the species, and preferably at the strain level, is very relevant. To test whether the complete genomes of the three newly sequenced strains could contribute to a better assignment of the species composition of NWCs, and in particular, whether the genomes could help to distinguish different strains of a species (Smid et al., 2014), a metagenomic analysis of a Gruyère whey starter culture from a cheese factory in Switzerland was performed using the Illumina HiSeq platform.

First, the species composition of the NWC was determined. Using a reference-based approach, the Illumina HiSeq reads were mapped to the genomes of the strains that were able to explain 95% of the reads of the metagenome sample (see Methods). About 57% of the reads originated from L. helveticus, 34% from S. thermophilus and 5% from L. delbrueckii, respectively (Figure 7, pie chart). Another 5% of the reads could not be mapped. However, this last percentage might possibly decrease if the—based on the mapped reads—“rare” species in the sample were also considered, which was beyond the scope of this analysis. The percentages from this whole genome sequencing-based metagenomics approach are comparable to those reported in a reverse transcriptase length heterogeneity PCR-based analysis of the composition of NWCs in Grana Padano (Rossetti et al., 2008), where the domain A of the variable 16S rRNA gene was assessed. L. helveticus strains were always dominant, while the percentages of S. thermophilus and L. delbrueckii seemed to vary in these cheese whey starters.

Next, the relevance of the three FAM strains in relation to existing NCBI RefSeq L. helveticus strains was assessed. To achieve this, the mapping information above was filtered for two groups: (1) reads mapping only to FAM strains and (2) reads mapping only to NCBI RefSeq strains (Figure 7, right upper panel). Although this information provides only a “semi-quantitative” picture, nevertheless, it emphasizes the relevance of the FAM strains for this NWC from a Gruyère cheese.

Finally, we determined the numbers of all uniquely mapped HiSeq reads, i.e., reads that mapped exclusively to one target sequence among the L. helveticus strains. The result of this analysis showed a high relevance of FAM8105 and FAM22155 compared to the RefSeq L. helveticus reference strains. In contrast, strain FAM8627 had virtually no unique mappings (Figure 7, right lower panel).

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.