The L. parabuchneri strains (Table 1) used in this study were cultivated in MRS broth (de Man et al., 1960) at 30°C. For longterm storage, strains were stored at −80°C in sterile resonstituted skim milk powder (10%, w/v). To determine the formation of histamine, the strains were cultivated in MRS broth supplemented with 0.3% L-histidine for 3 days at 30°C. Afterwards, the culture supernatants were diluted 1:100 in methanol and 10 μL was applied to high-performance thin-layer chromatography (HPTLC) cellulose plates. After the plates had been developed with 2-propanol, 25% ammonia (3:1, v/v), the imidazole ring of histamine and histidine was visualized by dipping the plates into Pauly's reagent.

DNA was extracted from L. parabuchneri with ZR Fungal/Bacterial DNA MiniPrep Kit (Lucerna Chem AG, Lucerne; Switzerland) and concentrated with Genomic DNA Clean & Concentrator Kit (Lucerna Chem AG) according to the manufacturer's instructions.

For the PacBio sequencing, 5 μg of high molecular weight DNA from L. parabuchneri FAM21731 were sheared in a Covaris g-TUBE (Covaris, Woburn, MA, USA) to obtain 20 kb fragments. After shearing, the DNA size distribution was analyzed with a FragmentAnalyzer (Advanced Analytical Technologies, Ames, IA, USA). A SMRTbell library was prepared using the PacBio DNA Template Prep Kit 2.0 (Pacific Biosciences, Menlo Park, CA, USA) according to the manufacturer's recommendations. The library was sequenced using two SMRT cells with P4/C2 chemistry on a PacBio RSII system with a movie length of 120 min. The sequencing yielded 118,000 post filter reads corresponding to 629 Mb with a mean read length of 5,326 bases.

FAM21731 was also sequenced using Illumina technology. The library preparation of the DNA from Lactobacillus parabuchneri FAM21731 using “TruSeq DNA Sample Preparation Kit” (15025064) was followed by gel size selection for 400–500 bp fragments. The library was paired-end sequenced (2 × 101 bp) in a fraction of a lane on an Illumina HiSeq 2000 instrument.

Library preparation of the DNA from the remaining 17 L. parabuchneri strains was performed using “TruSeq DNA PCR-Free LT Library Prep” (FC-121-3003, Insert size option: 350 bp). The libraries were paired-end sequenced (2 × 151 bp) in a fraction of a lane on an Illumina HiSeq 3000 instrument.

The strains FAM21809, FAM23164, FAM23165, FAM23166, and FAM23167 were additionally sequenced using an Ion Torrent device. Library preparation, amplification and sequencing was performed using Ion Xpress Plus Fragment Library Kit for AB Library Builder System, Ion PGM Template OT2 200 Kit and Ion PGM Sequencing 200 Kit v2 (Life Technologies Europe BV, Zug, Switzerland) according to the manufacturer's instructions. Libraries were barcoded using Ion Xpress Barcode Adapters 1-16 Kit (Life Technologies Europe BV). Five Libraries were pooled and sequenced on a single Ion 318 chip (Life Technologies Europe BV).

The Pacbio reads were assembled using the HGAP 3 (smrtanalysis-2.2.0) (Chin et al., 2013) standard procedure. To close the cyclic DNA of the chromosome and the plasmids, the assembly was performed twice. The resulting scaffolds of both assemblies were aligned using mauve (Darling et al., 2010) and the ends of the scaffold were replaced by the contiguous sequence of the other assembly. The assembly resulted in one chromosome (2,600,578 bp) and two plasmids (58,093 and 67,905 bp). Remapping and variant calling of the PacBio reads was performed using quiver (smrtanalysis-2.2.0) in which GATK (McKenna et al., 2010) is embedded. Finally, the kineticsTools (smrtanalysis-2.2.0) was used for the detection of modified bases in the genome. To find modified motifs MotifMaker (smrtanalysis-2.2.0) was used.

The Illumina reads were quality trimmed using Trimmomatic (version 0.33, options: SLIDINGWINDOW:4:8 MINLEN:127) (Bolger et al., 2014). The trimmed reads were assembled using SPAdes (version 3.6.1, options: –careful –mismatch-correction -k 21,33,55,77,99,127). The resulting sequences were scaffolded using SSPACE (version 3.0, default options) (Boetzer et al., 2011). Scaffolds with a lower median coverage than 20% of the median read-depth of the whole genome and scaffolds shorter than 200 bp were excluded.

The reads were mapped to the assembly of PacBio reads using BWA (version 0.7.10) (Li and Durbin, 2009). The variants were determined using haplotype-caller integrated in GATK (version 3.3.0) (McKenna et al., 2010). Manual curation was performed using IGV (Thorvaldsdóttir et al., 2013).

To locate the position of Lactobacillus parabucheri the protein sequence of all the proteins that are present as single ortholog in all strains were used for the construction of the phylogenetic tree. The 18 L. parabuchneri (Table 1) strains and also the strains downloaded from Genbank (Table S1) were annotated with Prokka (version 1.11) (Seemann, 2014) to avoid differences derived from different ab initio annotation tools. The orthologs were determined using OrthoMCL (version 2.0.9, default parameters) (Li et al., 2003). The amino acid sequences of the OGCs were aligned separately using clustal-omega (version 1.2.1) (Sievers et al., 2011). The tree was calculated using RAxML [version 8.1.2, (-m PROTGAMMAWAG), Stamatakis, 2006] based on the concatenated alignments.

The phylogenetic tree of the L. parabuchneri and L. buchneri strains was constructed with the same pipeline, with the difference that the nucleotide sequences of the OGCs were used and the model of RAxML was adjusted (-m GAMMAWAG).

The phylogenetic trees were visualized using iTol (Letunic and Bork, 2007).

To estimate the trend of the pan and core genome size the strains were randomly subsampled without replacement, from one to the total number of strains. For each of the 1,000 subsampling, the core and the pan genome was calculated, using an in-house python script.

The T-boxes of the HDC island that were identified using Prokka (Seemann, 2014) were aligned against the T-boxes from Lactococcus lactis (strains SK11 and Il1403, downloaded from http://regprecise.lbl.gov) (Novichkov et al., 2013) using the MUSCLE (Edgar, 2004) web interface. The conserved regions were adapted from a previous study (Lebeer et al., 2007).

GO terms were assigned to the CDSs of FAM21731 using blast2go-pipe [version 2.5, BLASTP (version 2.2.29+, NCBI non-redundant database from August 2015), Conesa et al., 2005] and InterProScan (version 5.14-53.0) (Quevillon et al., 2005).

The GO enrichment analysis was performed using the elim algorithm of the topGO package (Alexa and Rahnenführer, 2010). We took the GO terms of the genes from the HDC island as test set, and the GO terms of the remaining genes of the genome as background.

The SynOrF algorithm was implemented in python. The NetworkX (version 1.9.1, https://networkx.github.io/) library was used for graph analysis.

Two separate whole genome graphs were constructed based on the coding sequences (CDSs) and their positions in the genomes of the two compared strains. The CDSs were set as nodes in the graphs, and the edges were set between each CDS and its neighboring CDSs from the first to third degree in relation to genomic location. To determine homologs, all CDSs of one strain were aligned with all CDSs of the compared strain using BLASTP (version 2.2.28+) (Altschul et al., 1990). Two CDSs with an alignment similarity of more than 50% and an alignment that covered at least 50% of the longer CDSs were considered as homologs.

The selected homologs were used as edges to connect two whole genome graphs of the compared strains. Nodes of the whole genome graphs were removed if they had no edges to the other genome graph (i.e., if the CDS had no homologs in the compared genome). This removal resulted in a genome graph that consisted of unconnected subgraphs.

With the chosen homology cutoffs, it is possible for a CDS to have several homologs in the compared genome. To select the best synteny ortholog, the neighboring CDSs of the homologs in the genome graph are also taken into account (i.e., checked whether the CDS of the two genomes are also homologs). The homolog with the highest number of homolog neighbors (at least two are required) is selected as the best synteny ortholog. The whole workflow is depicted in Figure S1.

The script is available for download at https://github.com/danielwuethrich87/Synteny_ortholog_finder.git.

Between the genome of FAM21731 and the other 17 L. parabucheri strains and the L. buchneri stains (CD034, NRRL B-30929) the synteny orthologs were identified using the SynOrF algorithm. Not conserved regions became visible after visualization of the conserved synteny othologys in a Gantt plot using ggplot2 (Wickham, 2010).

The phylogenetic trees in Figures 1A, 2A, 3C, 4B were visualized using iTol (Letunic and Bork, 2007). The plots in Figures 1B, 2B, 3C, 4A were created with ggplot2 (Wickham, 2010). Figure S1 was drawn using Microsoft office PowerPoint 2013. Figure 3B, Figure S3 are copied text files. Figure 3A was created using genoPlotR (Guy et al., 2010).

Histamine in cheese is mainly formed by bacterial decarboxylation of histidine. When we analyzed cheese with elevated histamine levels (>25 mg per kg cheese), we repeatedly isolated histaminogenic strains of L. parabuchneri (Berthoud et al., 2016). Because a complete genome of this species was not available at the time we started the study, we selected the single molecule long read technology of PacBio to sequence the L. parabuchneri strain FAM21731 and to assemble a complete genome. The de novo assembly of the PacBio reads resulted in three contigs that we could close to circular molecules using manual correction. The largest contig is the chromosome (2,600,578 bp) and the two short contigs plasmids (58,093 bp; 67,905 bp, Table 1). Additionally, we sequenced FAM21731 using the short read technology of Illumina for quality control. Re-mapping of the PacBio and Illumina reads showed that the assembly does not have indels or base mismatches. We annotated 2,566 protein coding sequences (CDS) in the genome using Prokka (Seemann, 2014).

Besides the genome assembly, we also used the PacBio reads to identify methylated nucleotides in the genome of the strain FAM21731. We detected that adenine is methylated in the four motifs ACANNNNNNNTTYG, CACCAG, GATC and CGCCAT in 99.9, 99.9, 99.6, and 28.8% of the motifs in the genome, respectively. While the motif GATC was described in E. coli (Barras, 1989) and other species, so far no report about the three other motifs is available. Furthermore, all adenine in the motive ACANNNNNNNTTYG and its reverse complement CRAANNNNNNNTGT (97.0%) are methylated. The distribution of methylated motifs within the genome did not show any accumulation at specific locations.

Additionally, to the strain FAM21731, we sequenced 17 L. parabuchneri strains isolated from cheese and milk samples (Table 1) using Illumina technology. The assembly sizes varied between 2,562,122 and 2,799,109 bp. We annotated the genomes with Prokka and clustered the genes using OrthoMCL. Based on these gene cluster, we constructed a phylogenetic tree including—in addition to our L. parabuchneri strains - all complete Lactobacillus genomes available (57 strains, Table S1) in the Refseq database based on the OrthoMCL clusters. In total, we found 388 orthologous gene clusters (OGCs) that are present in all L. parabuchneri and in other Lactobacillus strains as single copy genes. The tree shows that the L. parabuchneri strains sequenced in this study cluster together with the L. buchneri strains (Figure 1A).

To get an insight into the gene repertoire present in L. parabuchneri, we constructed the pan and the core genome of the 18 sequenced strains. The core genome has a size of 1,947 OGCs and the pan genome a size of 3,461 OGCs. We applied random subsampling to the strain order to show the trend curves of the pan and core genome (Figure 1B). The size of the core genome seems to reach a plateau at 18 strains whereas the pan genome does not reach a plateau at 18 strains, i.e., additional strains would significantly increase the number of OGCs.

The phylogenetic tree shows that L. buchneri and its subspecies L. parabuchneri are very closely related (Figure 1A). We determined the size of the pan genome of the 18 L. parabuchneri to have 3,461 OGCs and with additional two L. buchneri (CD034 and NRRL B-30929) strains to have 3,910 OGCs. The small difference of 396 OGCs between the pan genomes shows that the gene content is very similar in these two species.

We screened the 18 strains for their ability to synthesize histamine. Therefore, the culture supernatants of the strains were analyzed for the presence of histamine after 3 d incubation using thin-layer chromatography. We found that 12 strains produced histamine (Table 1) and that the band intensities of histamine were very similar. In these samples histidine could not be detected anymore (an example is shown in Figure S2). The other six strains, namely FAM21835, FAM21838, FAM23279, FAM23280, FAM23281, and FAM23282, did not form histamine whereas histidine was still present. This indicated that the histamine-forming strains converted histidine completely to histamine.

To link the phenotype of histamine synthesis with a genotype, we assessed the phylogenetic relationship between the histamine producing and the non-producing L. parabuchneri strains. These 18 strains and the two L. buchneri strains share 1,685 OGCs that are present in every strain as single copy genes. Based on this, we constructed a phylogenetic tree using the nucleotide sequence of these OGCs and found that the L. parabuchneri strains are separated into three branches (Figure 2A). We show that the phylogenetic relation and the capability to synthesize histamine correlate perfectly. Furthermore, the 12 histamine positive (HDC+) strains were separated into two branches (branch 1, branch 2) in the phylogenetic tree (Figure 2A) whereas the six histamine negative (HDC−) strains were all located on the same branch (branch 3).

The HDC gene cluster found in the chromosome of FAM21731 encodes the following proteins: A putative histidine/histamine antiporter (FAM21731_00901), a histidine decarboxylase (FAM21731_00900), a HdcB protein (FAM21731_00899) and a histidyl-tRNA synthetase (FAM21731_00897, synonym prokka annotation: Histidine–tRNA ligase). We found orthologs of these four genes in all 12 strains with the capability to synthesize histamine. The other six HDC strains lack this gene cluster completely. Sequence alignments (Figure 3C) revealed that the HDC gene cluster nucleotide sequence from L. parabuchneri is identical to the HDC operon of L. buchneri B301 (GenBank AJ749838) that was characterized by Martín et al. (2005). This indicates that strain B301 probably belongs to the L. parabuchneri species. Meanwhile, the nucleotide sequences of pheS (GenBank LT547856), rpoA (GenBank LT547858), and the 16S rRNA gene (GenBank LN877766) from strain B301 were deposited confirming that this strain indeed is a L. parabuchneri strain. Remarkably, also the histamine-forming strains LMG11773 (= St2A) and DMS 5987 (= Lb14) that originally were deposited as L. buchneri were re-identified as L. parabuchneri.

We searched the genome of FAM21731 for regions that are not conserved within the L. parabuchneri and L. buchneri strains. For this purpose, we developed an algorithm termed “Synteny Ortholog Finder” (SynOrF) (see Figure S1 and the methods section). It determines orthologs between two genomes not only based on sequence homology but also in consideration of the syntenic relationship of the genes. This pairwise comparison to identify synteny orthologs is carried out between the genome of interest and each of the remaining genomes separately. We found one 62.4 kbp long genomic region (FAM21731_00889–FAM21731_00946, Table 2) in FAM21731 that is not present in the HDC- L. parabuchneri and L. buchneri strains. This region shows a significant difference in GC-content compared to other regions of the genome (39.5 vs. 43.7%), indicating that it is a genomic island that was gained by horizontal gene transfer. We refer to this region as HDC island (HDC island, Figure 2A).

We applied gene ontology (GO) enrichment with the genes coded in the HDC island and found that the GO terms gluconate transmembrane transport (GO:0035429, p = 1.7e-06), D-ribose metabolic process (GO:0006014, p = 2.8e-04) and histidine metabolic process (GO:0006547, p = 8.0e-04) are enriched with a high significance (p < 0.001). The genes related to histidine metabolic process are the four CDSs of the HDC gene cluster. The analysis of the gluconate transmembrane transport revealed that a cluster of three neighboring gluconate permeases on the HDC island is directly adjacent to the HDC gene cluster and a cluster of five genes involved in gluconate and galactonate (a stereoisomer of gluconate) metabolism (Figure 3A). The functional annotation shows that the genes involved in the ribose metabolism are arranged in a cluster of seven neighboring genes (Figure 3A). Finally, we identified two T-boxes on this genomic island that showed high similarity to two histidine-specific T-boxes of Lactococcus lactis (Figure 3B). One of these T-boxes (FAM21731_00898) is upstream of the histidyl-tRNA synthetase, the fourth gene of the HDC gene cluster. The second T-box (FAM21731_00905) is upstream of a gene cluster that consists of a peptidase (FAM21731_00903) and a peptide transporter (FAM21731_00904). We also found a second gene cluster consisting of a peptide transporter (FAM21731_00939) and a peptidase (FAM21731_00940) within the HDC island (Figure 3A).

We determined single nucleotide polymorphisms (SNPs) between the different strains and FAM21731. We found a 400 kbp region of the FAM21731 genome (FAM21731_00682–FAM21731_01148) that was conserved within the strains of branch 1 and branch 2 but not of branch 3 (Figure 4A). Interestingly, the HDC island is located within this 400 kbp region. Phylogenetic analysis of the L. parabuchneri strains only including the CDSs from 400 kbp region showed a clustering of the strains from branch 1 and branch 2 into a single branch (Figure 4B). Contrary, the phylogenetic analysis based on the CDS of the core genome, excluding the genes of the 400 kbp region, showed a clustering of the strains from branch 2 and branch 3 into a single one (Figure 4B).

L. parabuchneri is not the only species that can synthesize histamine. The HDC gene cluster found in L. parabuchneri shows the same gene order as the ones in Tetragenococcus halophilus, Tetragenococcus muriaticus, Lactobacillus hilgardii, Lactobacillus sakei, Lactobacillus saerimneri, Lactobacillus vaginalis, Lactobacillus fermentum, Oenococcus oeni and Lactobacillus reuteri. Gene clusters from species with different gene order (e.g., S. thermophilus) were not included. We observed that the nucleotide sequence of the gene encoding the histidine decarboxylase is closely related to the ones of these species (Figure 3C) (74.7–89.2% sequence identity). Comparing these results to the phylogeny of the core genome of the lactobacilli, we observed that the histidine decarboxylases are more closely related (e.g., L. sakei) to one another than the core genomes (Figure 1A).

Our analyses showed that the HDC gene cluster is located in the chromosome of L. parabuchneri as well as L. reuteri (GenBank NC_009513). In L. vaginalis, the HDC locus is in a 693-kb scaffold (GenBank NZ_GG693412). The scaffold size indicates that the HDC is probably also chromosomal. This hypothesis is supported by the observation that no plasmids were found in histamine-producing L. vaginalis strains (Diaz et al., 2015). On the contrary, the HDC gene cluster was found in plasmids in strains of O. oenii (Lucas et al., 2005), L. hilgardii (Lucas et al., 2008), L. fermentum (GenBank NZ_AVAB01000110), T. halophilus (Satomi et al., 2008), and T. muriaticus (GenBank NC_918355). The sizes of these plasmids vary between 100 kb (O. oenii) and 23 kb (T. muriaticus). In the draft genome sequence of L. saerimneri 30a, one of the first lactobacillus species described to produce histamine (Rosenthaler et al., 1965), the HDC gene cluster is present in a 92-kb contig (GenBank ANAG01000006). As we did not find the hdc genes in the draft genome of L. saerimneri DSM 16094, we assume that the HDC locus is also encoded on plasmids in this species.

The alignment of these HDC gene cluster (Figure 3C) revealed that the sequences are highly conserved in the species L. parabuchneri and the other species that can also be present in food (T. halophilus, T. muriaticus, L. hilgardii, L. sakei and O. oeni). Interestingly, we found that the nucleotide sequence of the hdcB gene, which is predicted to catalyze the maturation of pyruvoyl-dependent histidine decarboxylase HDC (Trip et al., 2011) of L. parabuchneri, is highly different in L. saerimneri, L. vaginalis, L. fermentum, and L. reuteri. Especially, the first two-thirds of the hdcB gene show considerable less sequence identity than the last third.

Finally, we compared the conservation of the HDC gene cluster from different species to the general sequence conservation between L. parabuchneri and L. buchneri. We calculated the mean similarity of the nucleotide sequence between the four genes of the HDC gene cluster of L. parabuchneri and the HDC gene cluster from T. halophilus, T. muriaticus, and L. hilgardii. We found that the average DNA sequence identity is 82.5% of the genes of the HDC gene cluster. Compared to this the average sequence identity of the core genome of L. parabuchneri FAM21731 and L. buchneri NRRL B-30929 is 77.4%. We assume that the events of the separation of the species L. parabuchneri and L. buchneri and the gain of the HDC gene cluster lie chronologically rather close to each other.

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.