Genomic DNA (gDNA) was extracted from bacterial cultures, NWCs, and cheese samples as described by Moser et al. (2017). DNA quantity was determined fluorometrically using the Qubit dsDNA BR Assay Kit (Thermo Fisher Scientific, Baar, Switzerland).

The genomes of 57 L. helveticus strains from the Agroscope culture collection and that of one L. helveticus strain from the Direct Vat Set culture LH-32 from Christian Hansen (Copenhagen, Denmark) were sequenced (Supplemental Table S1). The “TruSeq DNA PCR-Free LT Library Prep” (FC-121-3003, Insert size option: 350 bp) was used to prepare the DNA libraries. The libraries were indexed and pooled, before being sequenced in one lane on an Illumina HiSeq 3000 instrument to produce paired-end reads (151 × 151). Trimmomatic (version 0.33, options: SLIDINGWINDOW:4:8 MINLEN:127; Bolger et al., 2014) was used to trim the raw reads. The remaining reads were assembled using SPAdes (version 3.6.1, options: —careful—mismatch-correction—k 21, 33, 55, 77, 99, 127) and the assembly was scaffolded using SSPACE (version 3.0, default options; Boetzer et al., 2011). The scaffolds were selected for length and coverage as follows: scaffolds that were shorter than 200 bp, and those with a lower median coverage than 20% of the median read-depth of all scaffolds larger than 5,000 bp, were excluded. The resulting assemblies were annotated using Prokka (version 1.11; Seemann, 2014). To determine the orthologous gene clusters (OGCs) between the genomes, all predicted protein-coding sequences were translated and compared using BLASTP (version 2.2.29+, default parameters; Altschul et al., 1990) and clustered using Ortho-MCL (version 2.0.9, default parameters; Li et al., 2003). Single copy OGCs that were present in all genomes were ranked according to the number of unique sequences per OGC. Afterwards the corresponding nucleotide sequences were extracted and aligned using Clustal Omega (version 1.2.1). The alignments were then inspected and not aligned sequences were removed manually using the CLC Main Workbench software (Version 7.5.1, Qiagen, Switzerland). Finally, the remaining aligned sequences were checked by eye for nucleotide polymorphisms. Thereby, a CDS encoding a surface layer protein, named slpH in this report, was found to carry a high amount of polymorphisms and could be used as a target for DNA-based strain typing.

For the assignment of sequence types (STs), additional orthologous genes were retrieved from the NCBI GenBank database. Based on the nucleotide polymorphism present in these genes, the nucleotide sequences were assigned to STs (Table 1). SlpH genes with newly identified polymorphisms were deposited in the NCBI database and their accession numbers are listed in Table 1. The sequences of all of the different slpH genes are presented in the nucleotide sequence alignment in Figure S1. The nucleotide sequences of the ORFs from all STs were translated using the CLC Main Workbench version 7.5.1 (CLC bio, Aarhus, Denmark) and the resulting amino acid sequences were analyzed using InterProScan (Jones et al., 2014). Furthermore, the genomic context of the slpH gene was analyzed in 10 complete L. helveticus genomes from the NCBI database (CP002081, CP000517, CP003799, CP002429, CP002427, CP009907, CP011386, CP012381, CP016827, and CP020029). Therefore, the genomes were reannotated using Prokka (version 1.11; Seemann, 2014) and the orthologous genes were determined using Ortho-MCL (Li et al., 2003). The slpH genes and the surrounding genes of all of the 10 strains were visualized using the package ggplot2 in R (Wickham, 2009).

The slpH sequence from L. helveticus CNRZ32 (Table 1) served as reference sequence for primer design. Based on the reference sequence, the primer pair LHslpF (5′-CAAGGAGGAAAGACCACATGA-3′) and LHslpR (5′-TGTACTTGCCAGTTGCCTTG-3′) that amplifies a 1,116-bp region was designed. Primers were designed using Primer3 (version 0.4.0; Untergasser et al., 2012). PCR was carried out on a Veriti® Thermal Cycler (Thermo Fisher Scientific, Baar, Switzerland) with the following conditions: 95°C for 2 min followed by 30 cycles of 95°C for 20 s, 60°C for 10 s, and 70°C for 30 s, and a final extension of 70°C for 7 min. If not otherwise specified, each PCR (25 μL) contained 50 ng of gDNA for cheese or NWC samples or 1 ng of gDNA for pure cultures, 5 pmol of each primer, 5 nmol of each dNTP, 31.25 nmol Mg₂SO₄, 0.02 U of KOD Hot Start Polymerase (Merck), and 1X buffer for KOD Hot Start Polymerase. PCR products were analyzed using the Agilent DNA 7500 kit on an Agilent 2100 bioanalyzer (Agilent Technologies, Waldbronn, Germany).

An NWC sample (100 μL, sample p in Table 2) was plated on modified MRS plates containing 20 g L⁻¹ of lactose instead of glucose (De Man et al., 1960). After incubation at 42°C for 48 h, 96 colonies were picked from the agar plates and suspended individually in 100 μL of TE buffer (pH 8.0) containing 10 mM Tris-HCl and 10 mM EDTA. A heat treatment (100°C, 8 min) was used to extract gDNA. After the heat treatment, the suspensions with the lysed cells were centrifuged at 5,000 g for 10 min. The supernatant containing the DNA was diluted tenfold in 10 mM Tris-HCl (pH 8.0) and examined by quantitative real-time PCR (qPCR) to identify colonies of L. helveticus. DNA from L. helveticus positive colonies was further examined by slpH PCR as described above. The PCR products were sent to Microsynth (Balgach, Switzerland) for amplicon purification and Sanger sequencing.

Identification of L. helveticus isolates and quantitation of the population density of L. helveticus representatives in NWCs and cheese was performed by a qPCR method targeting the single copy pheS gene that encodes the alpha-subunit of the phenylalanyl-tRNA synthetase as described elsewhere (Moser et al., 2017). Briefly, a plasmid containing the target sequence is used for absolute quantification. The equation of the standard curve was used to estimate the copy number for cheese (gene equivalent per g) and NWCs (gene equivalent per mL). Similarly, the pheS gene of Lactobacillus gallinarum was used as a target to develop a PCR-based detection method for this species. Based on the pheS nucleotide sequence from L. gallinarum DSM 10532 (NCBI GenBank Accession Number AZEL01000044:72250…73299), the primers 5′-TCAGGACCTTGTACTACCTTGTAA-3′ and 5′-TGCTACTAAGGCTGAAATCGT-3′ were designed, which enable amplification of a 180-bp fragment. The assay was used to analyze NWC samples for the absence or presence of L. gallinarum. PCR assays were conducted in a final volume of 25 μL containing 300 nM of each primer, 200 μM of each dNTP, 1.5 mM Mg₂SO₄, 0.02 U of KOD Hot Start Polymerase (Merck), and 1X buffer for KOD Hot Start Polymerase (Merck). Each PCR contained 50 ng of genomic DNA for NWC samples and 1 ng of DNA for pure cultures. The amplicons were amplified under the following conditions: 95°C for 2 min, followed by 30 cycles of 95°C for 20 s, 59°C for 10 s, and 70°C for 4 s, and a final extension of 70°C for 5 min. The PCR products were examined with an Agilent DNA 1000 kit on an Agilent 2,100 Bioanalyzer (Agilent Technologies, Waldbronn, Germany). To evaluate the lower limit of detection, DNA was extracted from a ten-fold dilution series of a pure culture from L. gallinarum DSM 10532 that had been grown overnight in MRS broth at 37°C and for which the population density had been determined by plate-counting.

After slpH-specific PCR, amplicons were purified using the Qiaquick PCR Purification kit (Qiagen, Hombrechtikon, Switzerland). Purified amplicons were fragmented by sonication of 1 μg of DNA in 130 μL nuclease-free water at 50 W, 200 cycles/burst, and 20°C for 90 s on a Covaris M220 instrument (Covaris, Brighton, U.K.). The fragmented DNA was end-repaired and adapter-ligated using the Ion Xpress Plus Fragment Library Kit (Thermo Fisher Scientific) according to the manufacturer's instruction. Samples were barcoded using the Ion Xpress Barcode Adapter 1-16 kit (Thermo Fisher Scientific). After barcoding, the DNA was size-selected for 400 bp reads using an E-Gel® SizeSelect™ Agarose Gel (Thermo Fisher Scientific). The libraries were quantified by qPCR using the Ion Library Quantitation Kit (Thermo Fisher Scientific). Before preparing the template-positive ion sphere particles with the Ion PGM Hi-Q OT2 kit (Thermo Fisher Scientific), each library was diluted to 100 pM in 10 mM Tris-HCl (pH 8.0) and 0.1 mM EDTA, before being pooled. The sequencing was performed using either an Ion 314™ or 316™ chip and the Ion PGM Hi-Q Sequencing kit on an Ion Torrent sequencer (Thermo Fisher Scientific).

First, fastq files were generated via the Torrent Suite Software (Thermo Fisher Scientific) using the default parameters; then, the reads were quality trimmed with Trimmomatic (version 0.36, options: SLIDINGWINDOW:4:20 MINLEN:101; Bolger et al., 2014). Finally, reads containing the identifying subsequences for slpH group1 (GGCTACACT, GATCAATTAA, AGTGTAGCC, and TTAATTGATC), for slpH group2 (CCTTAATGTA, CTGACGATGT, TACATTAAGG, and ACATCGTCAG), and for slpH group 3 (ATTGGTTCAG, GGTGTTGCTA, CTGAACCAAT, and TAGCAACACC) were extracted, trimmed, grouped based on 100% sequence identity, and mapped against a reference database containing known STs. Sequences that were not identical to one of the reference sequences were assigned to a new ST. For this procedure, a Python script was developed. The script, its usage and the database for the slpH STs are available at: https://github.com/danielwuethrich87/helveticus_strain_typing.

The gDNA from 10 different L. helveticus STs (Table 3) were pooled in equal amounts. The amplicons obtained using the slpH-specific PCR were subsequently sequenced on the Ion Torrent sequencer as described above. The experiment was repeated three times.

We manually inspected the single copy orthologs of 58 own-sequenced genomes of L. helveticus. One of the gene clusters encoding a putative surface layer protein, named slpH in this study, stood out from the others in terms of nucleotide sequence diversity and differentiated 15 STs. Other OGCs differentiated between 8 and 11 STs and were excluded from further analysis. The primer pair LHslpF/R was designed to amplify the slpH locus in L. helveticus. The specificity of the primer pair was tested using the gDNA of 18 L. helveticus strains and 20 other species of LAB (Table 4). By using the primer pair LHslpF/R, an amplicon with a size of approximately 1,252 bp (±60 bp) was obtained for all the L. helveticus and L. gallinarum tested. No amplicon was observed using gDNA from Streptococcus thermophilus and other LAB, including the more closely related Lactobacillus kefiranofaciens kefiranofaciens DSM 5016, L. kefiranofaciens kefirgranum DSM 10550, Lactobacillus crispatus DSM 10532, Lactobacillus acidophilus DSM 20079, and Lactobacillus amylovorus DSM 20531 (Table 4).

The nucleotide sequences of the L. helveticus amplicons were determined by Sanger sequencing and subsequently compared to slpH sequences that had been extracted from 79 bacterial genomes consisting of 58 genomes of own-sequenced strains and 21 from the NCBI GenBank. We found that the sequences were identical to the one derived from the illumina read assembly (data not shown) and that 30 STs could be unambiguously discriminated (Supplemental Figure S1). Additionally, the amplicon obtained from the L. gallinarum DSM 10532 was determined. A BLAST search revealed that the amplified region showed 99% identity to the lgsB gene (GenBank accession number AY597262) of L. gallinarum, which encodes a surface layer protein.

A pairwise sequence comparison of 30 STs showed that the slpH sequences from L. helveticus clustered into three groups (Figure 1). Strains that clustered within a group shared, on average, 91.4 (±4.6)% sequence identity, whereas the identity between groups was, on average, 57.1 (±1.7)%.

When the deduced amino acid sequences were compared, the primary sequences of ST1 and ST2, ST18 and ST27, and ST19 and ST28, respectively, were identical. All other STs could still be discriminated based on amino acid sequence and especially a high level of variation was observed between the three slpH groups (Figure 2). When the primary sequences were analyzed for domains, InterProScan analysis assigned all STs but ST8, ST9, ST10, and ST11, which belong to the slpH2 group, to the Lactobacillus surface layer protein family (InterPro accession number PIRSF037863) using the PIRSF family classification system (not shown). Furthermore, InterProScan identified in all STs the eight motif fingerprints of Lactobacillus surface layer protein family (InterPro accession number PR01729) based on the PRINTS database for protein fingerprints (blue bars, Figure 2). Additionally, the analysis predicted the presence of a tandem SLAP domain (InterPro accession number IPR024968) in all STs (red bars, Figure 2). Finally, a signal peptide with a length of 30 amino acids was predicted for all the STs (green bars, Figure 2).

When we compared the surrounding context of the slpH locus in 10 strains possessing different STs and slpH groups, we found that up- and downstream various OGCs were co-localized with the slpH locus between the strains (Figure 3). Nine of the genes surrounding slpH were present in all of the analyzed strains (Figure 3).

Since the primer pair LHslpF/R can also amplify a region from the gDNA of L. gallinarum strains, we established a PCR to assay dairy products for the presence of this species. A primer pair targeting the pheS gene of L. gallinarum showed species-specificity, as only a 180-bp amplicon was amplified from the gDNA of L. gallinarum, but not from closely related LAB (Table 4). The assay's limit of detection was determined by analyzing gDNA extracted from a dilution series of L. gallinarum DSM 10532. The lower limit of detection was determined to be at 2.45 × 10² colony forming units per mL broth. The method was generally applied to all samples used for semiconductor sequencing to confirm the absence of L. gallinarum.

A mixture of 10 L. helveticus gDNAs, each of which represented a different ST (Table 3), was prepared. The amplicons obtained after slpH amplification were sequenced on the Ion Torrent PGM sequencer. After extraction of the reads containing the slpH group specific identifying subsequences, the proportion of reads assigned to the 10 known STs ranged from 5.85 (±0.12)% to 15.73 (±1.00)% (Table 3). Each slpH group also comprised reads that were not assigned to the known STs and indicate the presence of new STs. The abundance of those reads per new ST did not exceed 1.7 (±0.3)% of all reads per slpH group. Based on this observation, we decided that at least 3% of all reads within a slpH group must be assigned to an ST to be considered a true positive result.

We also assessed the applicability of the slpH locus for typing, by comparing the slpH loci obtained after strain isolation using the culture-dependent approach with those obtained with the culture-independent approach. To effect this, the NWC from La Praz (sample p in Table 2) was plated on MRS agar plates containing lactose as a carbohydrate source. Of the 96 colonies picked, 81 were identified as L. helveticus using a species-specific qPCR. The slpH locus was amplified and sequenced from each L. helveticus isolate. In addition to two previously known STs, we identified two new STs, 31 and 32. The four STs were distributed as follows: 54.3% of the colonies were assigned to ST13 (slpH group 1), 8.6% were assigned to ST31 (slpH group 2), and 34.6 and 2.5% were assigned to ST17 and ST32, respectively (both slpH group 3; Figure 4).

This result was then compared to the sequence analysis obtained from the culture-independent approach. We identified the same STs as found by the culture-dependent approach. All assigned reads were distributed as follows: 21.66% were assigned to ST13, 49.04% were assigned to ST31, and 28.23 and 1.07% were assigned to ST17 and ST32, respectively (Figure 4).

The diversity of L. helveticus strains present in 15 NWCs, collected at variously located Gruyère cheesemaking factories (Table 2), was analyzed using our established culture-independent method. The quantification of L. helveticus using qPCR revealed an average population density of 4.42 (±3.79) × 10⁷ cells per mL of NWCs according to pheS copy numbers. The presence of L. gallinarum was not detected by the specific PCR assay. When we analyzed the slpH loci using the culture-independent approach, we detected a total of 10 STs in the 15 NWC samples (Figure 5). The presence of STs ranged from one to five per NWC. Therefore, for example, the NWC from La Brévine contained only one strain, whereas the specimen from Orsonnens contained five STs. The majority (10 NWCs) contained three different STs of L. helveticus. Remarkably, ST13 was present in all the NWC samples (Figure 5, blue bar).

We also monitored the development of L. helveticus strains during cheese ripening by analyzing samples collected from a cheese that had been manufactured with three L. helveticus STs, namely ST10 (FAM13019), ST 23 (FAM23236), and ST13 (FAM23237). QPCR estimated the population density of L. helveticus to be 10⁸ copies g⁻¹ cheese after 24 h of ripening and 10⁶ copies g⁻¹ cheese after 80–120 days (Table 2). By using the typing method, we detected only two of the three added L. helveticus STs in the cheese ripened for 24 h (Figure 6). After 80 and 120 days, all three of the STs used were detected. The relative abundance of ST13 increased from 27.36% after 24 h of ripening to 31.42% after 80 days and 67.47% after 120 days. For ST23, the relative abundance decreased from 72.64% after 24 h of ripening to 56.34% after 80 days and 25.92% after 120 days. ST10 was not detected in the cheese after 24 h of ripening and had a relative abundance of 12.24% after 80 days, which decreased to 6.6% after 120 days (Figure 6).

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.