Does habitat shape the gene content of dairy Lactobacillus helveticus?

A comparative pangenomics analysis of nine Lactobacillus helveticus genomes isolated from Swiss-type cheese, other cheeses and fermented milk, asking whether strains that share a fermentation substrate also share accessory gene content — and whether that signal can be separated from shared ancestry.

Short answer: gene content groups the strains almost perfectly by isolation source, but the grouping is explained by phylogeny rather than by habitat, and with nine genomes no gene reaches genome-wide significance for niche association. The interesting result here is a negative one, and the analysis is built to show why.

---

Contents

• Background
• Data
• Methods
• Results
• Discussion
• Limitations
• Reproducing the analysis
• Repository layout
• Provenance and corrections
• References

---

Background

Lactobacillus helveticus — one of the species that kept the genus name Lactobacillus when the Lactobacillaceae were reorganised into 25 genera (Zheng et al. 2020) — is a thermophilic lactic acid bacterium used as a starter culture in Swiss-type and Italian hard cheeses and in fermented milks. Its industrial value rests almost entirely on its proteolytic system: a cell-envelope proteinase cleaves casein into oligopeptides, an oligopeptide transporter imports them, and a large battery of intracellular peptidases degrades them to free amino acids (Savijoki et al. 2006). Those amino acids both drive acidification and feed the catabolic reactions — transamination, methionine demethiolation, decarboxylation — that generate cheese aroma compounds and reduce bitterness.

Dairy lactobacilli are a textbook case of reductive genome evolution under domestication: in a nutrient-rich, low-competition environment, genes for functions the milk provides are inactivated by IS-element insertion and eventually lost, while proteolytic and transport capacity is retained or expanded (Callanan et al. 2008). L. helveticus DPC 4571 carries over 200 IS elements and a correspondingly decayed genome.

This raises a straightforward comparative question, and the one this project asks:

If domestication tunes gene content to substrate, do strains isolated from the same kind of dairy product carry the same accessory genes — over and above what their shared ancestry already predicts?

Data

Nine complete or draft L. helveticus assemblies from NCBI GenBank, spanning three isolation sources:

Niche	n	Strains	Source description
Swiss-type cheese	5	Lh11051, Lh11961, Lh21456, Lh21462, Lh21463	swiss cheese, cheese
Other cheese	2	CNRZ32, DPC 4571	milk, cheese / dairy, cheese
Fermented milk	2	H10, TJA10	milk, fermented milk / milk, dairy, fermented milk

The five Lh-prefixed strains come from a single 2024 study of starter-culture domestication (PMID 39366947); H10 was isolated from Tibetan fermented milk, TJA10 from a Tajik fermented dairy product, and CNRZ32 and DPC 4571 are long-established cheese starters. Full metadata, including assembly accessions, is in results/tables/strains.csv.

The upstream pangenome was built with Prokka 1.14 for annotation and Roary 3.13 for clustering:

prokka --outdir "$out_dir" --prefix "$name" --locustag "${name:0:20}" --cpus 4 "$fna"
roary -e -n -v prokka_output/*/*.gff -f roary_output

Functional annotation of the pangenome reference sequences was done with eggNOG-mapper. scripts/fetch_data.sh downloads these precomputed intermediates; this project consumes them and does not re-run Prokka or Roary.

Methods

Pangenome partition. Gene clusters are classified core / soft-core / shell / cloud using Roary's own occupancy thresholds, reimplemented so the partition is explicit and testable. As a check on the whole loading path, the reimplementation reproduces the counts in the upstream summary_statistics.txt exactly (1261 / 0 / 1432 / 1608).

Distances. Strains are compared on the accessory genome only. Including the 1261 core clusters would add a constant to every genome pair and compress all distances toward zero. Jaccard rather than Hamming distance is used, because Hamming counts shared absences as similarity: two genomes that both happen to lack the same large block of accessory genes would look related for no biological reason.

Phylogenetic control. Gene-content distances are compared against the core-genome maximum-likelihood tree shipped with the data, in two ways: an unrooted Robinson–Foulds comparison of the core-gene and accessory-content trees, and a Spearman correlation between the two distance matrices over all 36 strain pairs. This is the step that separates "these strains share genes because they share a habitat" from "…because they are close relatives".

Niche association. Every informative cluster is tested for association with the Swiss-type cheese niche by a two-sided Fisher exact test, with Benjamini–Hochberg correction across all 3040 tests. Alongside it, the descriptive prevalence filter used in the original coursework (present in ≥80% of in-group genomes, ≤20% of the rest) is retained for hypothesis generation and reported as exploratory throughout.

Functional sets. Peptidase and flavour-forming gene sets are assembled by keyword matching against Prokka annotations. This is a recall-oriented screen, not a classification: it deliberately pulls in housekeeping proteins that share vocabulary (signal peptidases, Clp proteases) and those are visible in the output tables rather than silently filtered.

Results

1. The pangenome is open, and mostly not shared

4301 gene clusters across nine genomes: 1261 core (29%), 1432 shell, 1608 cloud. The U-shaped occupancy spectrum — a large core, a larger set of near-private genes, little in between — is the standard signature of an open pangenome with active gene gain and loss. More than a third of all clusters are found in a single genome.

2. Gene content groups strains by isolation source, with one exception

The dendrogram at the top of this README shows four of the five Swiss-cheese strains (Lh11051, Lh21456, Lh21462, Lh21463) forming a tight group at Jaccard distance < 0.3. Lh11961, despite being a Swiss-cheese isolate, sits with the cheese starter CNRZ32 instead. The two fermented-milk strains, H10 and TJA10, do not group together at all.

Taken alone, this looks like the classic story: gene content tracks habitat, with one anomaly to be explained by convergent adaptation. The next result shows that reading is wrong.

3. Gene content tracks ancestry, not habitat

Accessory gene-content distance correlates strongly with core-genome patristic distance (Spearman ρ = 0.84 over 36 strain pairs). The core-gene and accessory-content trees share 4 of their 6 non-trivial splits (Robinson–Foulds distance 4).

Critically, the core-genome tree already places Lh11961 as the sister taxon of CNRZ32, nested inside the non-Swiss clade rather than with the other Swiss isolates. Lh11961 does not resemble CNRZ32 because the two converged on a shared substrate; it resembles CNRZ32 because they are close relatives. The apparent anomaly in gene content is simply vertical inheritance showing through.

The same logic undercuts the positive signal: the four tightly clustered Swiss strains are also nearly identical in the core tree (patristic distance ≈ 0). In this panel, niche and ancestry are confounded by design — the Swiss strains are all from one study and one clonal lineage. No amount of analysis of these nine genomes can attribute their shared gene content to their shared habitat.

4. No gene is significantly associated with niche

Quantity	Value
Informative clusters tested	3040
Smallest achievable p for this design	0.0079
Smallest observed p	0.0079 (15 clusters, perfect 5/5 vs 0/4 split)
Smallest BH-adjusted p	0.674
Clusters significant at FDR 0.05	0

With 5 in-group and 4 out-group genomes, the most extreme table available — a gene present in all five Swiss strains and none of the others — yields p = 1/126 ≈ 0.0079. Against 3040 tests, that cannot survive multiple-testing correction: the study is underpowered by construction, not merely negative in outcome. Reporting the achievable floor alongside the result makes this checkable rather than a matter of interpretation.

The exploratory prevalence filter nevertheless returns 105 clusters that are common in Swiss strains and absent elsewhere, 58 of them with informative annotation (figure, table). Their dominant functional category is L — replication, recombination and repair (18 clusters), which in this organism largely means transposases and IS elements rather than DNA metabolism proper. That is consistent with mobile-element expansion in the domesticated lineage, and it is equally consistent with those five near-clonal strains simply inheriting the same insertions. These are hypotheses for a larger panel, not findings.

5. The proteolytic system is near-universal; the cell-envelope proteinase is not

Of 50 peptidase/protease clusters, 41 are core — every strain carries essentially the complete intracellular peptidase battery (PepC, PepN, PepX, PepT, PepO, PepQ, PepV, the Pep D/E/DA families). This is expected: casein degradation is the shared requirement of every dairy niche, so the proteolytic core is not where niche differences would appear.

The nine variable clusters are more informative, and the striking pattern is in the cell-envelope proteinases, the enzymes that perform the first cleavage of casein:

• prtP (PI/PII-type proteinase) — CNRZ32 only
• group_2622 (PIII-type proteinase) — H10 only
• group_4001 (Dipeptidase A) and yfeW — H10 only
• group_1343 (Aminopeptidase C) — CNRZ32 only

None of the five Swiss-cheese strains carries an annotated cell-envelope proteinase in this pangenome. In Swiss-type cheese manufacture L. helveticus is used together with Streptococcus thermophilus and propionibacteria, and proteinase-negative strains can subsist on peptides liberated by their partners — but with five clonal genomes this is a plausible reading, not a demonstrated one.

6. H10 carries a large, mobile-element-rich private gene set

H10, the Tibetan fermented-milk isolate and the phylogenetically most isolated strain in the panel, carries 475 strain-specific clusters. Excluding the uninformative "function unknown" category (102), the largest category is replication, recombination and repair (96) — again dominated by transposases and recombinases — followed by transcription (34), amino acid transport and metabolism (26) and carbohydrate transport and metabolism (26). The amino acid and carbohydrate transport genes are the biologically interesting ones for a strain growing on a casein-rich substrate, and H10's private PIII-type proteinase and dipeptidase fall in this set.

Of 49 flavour-related clusters (amino acid transaminases, cystathionine lyases, decarboxylases, methionine-sulfoxide reductases), 29 are core and 20 variable (figure) — the flavour-forming potential is, like proteolysis, largely a species-level property here.

Discussion

The headline claim this analysis was set up to test — that dairy L. helveticus strains cluster by fermentation substrate — is observed but not supported. Gene content does partition the panel by isolation source, and it would have been easy to stop there and write up a story about convergent adaptation to Swiss cheese. The core-genome comparison shows that the partition is a phylogenetic one: the Swiss strains are a clonal group from a single study, and the one strain that breaks the pattern breaks it in exactly the way its position in the core tree predicts.

This is the general hazard in comparative genomics of domesticated microbes. Isolates from the same product are frequently also close relatives, because they descend from the same commercial starter lineage. Any niche-association analysis that does not condition on phylogeny will recover the niche signal and misattribute it.

What the data do support is narrower and still worth stating: the proteolytic and flavour-forming machinery of L. helveticus is largely core, with strain-level variation concentrated in the cell-envelope proteinase and a handful of accessory peptidases; and accessory genome content in this species is dominated by mobile elements, consistent with the IS-driven genome decay documented for DPC 4571.

A design that could actually answer the original question would need isolates sampled across the phylogeny within each niche — several independent lineages per substrate — and would test niche association with the tree as a covariate (a phylogenetic logistic regression or a gene-wise phylogenetic ANOVA, e.g. the approach implemented in Scoary or treeWAS). That is the natural next step and is out of reach with nine genomes.

Limitations

• n = 9, and not independent. Five of nine genomes are near-clonal, so the effective sample size is closer to five lineages. Niche and phylogeny are confounded; the FDR result quantifies this rather than working around it.
• Annotation-driven gene sets. Peptidase and flavour sets come from keyword matching on Prokka annotations. Recall is incomplete (unannotated hypothetical proteins are invisible) and precision is imperfect (housekeeping proteases are included). Absence of a cell-envelope proteinase in the Swiss strains is absence of an annotation, and would need HMM-based confirmation.
• Draft assemblies. A gene called absent may be missing from an assembly rather than from the genome. Cloud genes in particular are sensitive to assembly fragmentation, and no assembly-quality filtering was applied upstream.
• Roary clustering thresholds were left at their defaults; paralogous families can be split or merged, which affects private-gene counts such as H10's 475.
• No phylogeny-aware association test was run, because with this panel it would be uninformative. The correlation and Robinson–Foulds comparison are descriptive.

Reproducing the analysis

Requires Python ≥ 3.10 and git. The whole pipeline runs in well under a minute.

make setup     # create .venv, install the package and dev dependencies
make data      # download the Roary / eggNOG inputs into data/raw/
make results   # regenerate every figure, table and summary.json
make test      # 33 unit tests

Or without make:

python -m venv .venv && .venv/bin/pip install -e ".[dev]"
./scripts/fetch_data.sh
.venv/bin/lhpan --verbose
.venv/bin/pytest -q

Every number quoted in this README is written by the pipeline to results/summary.json; every figure and table under results/ is regenerated from scratch on each run. Raw inputs are not committed — make data fetches them.

Repository layout

src/lhpan/
  config.py      paths, niche mapping, COG table, keyword sets, plotting theme
  data.py        loading and the three-way identifier reconciliation
  pangenome.py   binary matrix, occupancy classes, distances, gene tables
  enrichment.py  Fisher/BH niche association, prevalence filter, COG summaries
  phylogeny.py   Newick parsing, bipartitions, Robinson-Foulds, patristic distances
  figures.py     all figures, one shared theme
  pipeline.py    end-to-end run; writes results/ and summary.json
tests/           33 unit tests on hand-built matrices and toy trees
scripts/         data download
results/         figures, tables, summary.json (committed)

Provenance and corrections

This project is a rewrite of a first-year MSc bioinformatics assignment (HSE University, 2025) into a reproducible analysis. The original was a single 82-cell notebook with Russian-language commentary; the biological questions and the choice of dataset are from that assignment. The following changes were made, and are listed because they affect the results:

1. Identifier join fixed. The original mapped eggNOG query names to Roary clusters with query_name.rstrip('_1'). str.rstrip strips a character set, not a suffix, so every identifier ending in 1 was truncated (..._00001 → ..._0000) and lost its annotation — 373 of 3934 rows, 9.5% of the functional table. eggNOG query names match the FASTA record ids exactly, so the join is now direct. Covered by a regression test.
2. "Swiss-specific" genes were one strain's genes. The original selected the Swiss group with meta['short_name'].str.contains('Lh')...values[0], which returns the first matching strain. What was reported as Swiss-cheese-unique gene content was Lh21462-unique gene content. Replaced with an explicit group contrast.
3. A heatmap labelled "peptidase genes" plotted all 4301 clusters — the filtered variable was computed but not used. Removed.
4. Hamming → Jaccard for gene-content distance, so that shared absences no longer count as similarity.
5. A significance test was added. The original reported the ≥80%/≤20% prevalence filter as a result. That filter is kept as an exploratory screen, and the Fisher/BH test plus the achievable-p-floor calculation were added; they show the design cannot support the claim.
6. The phylogenetic comparison is new. The original attributed Lh11961's placement to convergent evolution and transposon activity. The core-genome tree, which was present in the input data but unused, shows Lh11961 and CNRZ32 are sister taxa, which explains the placement without invoking convergence.

References

1. Callanan M, Kaleta P, O'Callaghan J, et al. Genome sequence of Lactobacillus helveticus, an organism distinguished by selective gene loss and insertion sequence element expansion. J Bacteriol 190:727–735 (2008). doi:10.1128/JB.01295-07
2. Savijoki K, Ingmer H, Varmanen P. Proteolytic systems of lactic acid bacteria. Appl Microbiol Biotechnol 71:394–406 (2006). doi:10.1007/s00253-006-0427-1
3. Page AJ, Cummins CA, Hunt M, et al. Roary: rapid large-scale prokaryote pan genome analysis. Bioinformatics 31:3691–3693 (2015). doi:10.1093/bioinformatics/btv421
4. Seemann T. Prokka: rapid prokaryotic genome annotation. Bioinformatics 30:2068–2069 (2014). doi:10.1093/bioinformatics/btu153
5. Cantalapiedra CP, Hernández-Plaza A, Letunic I, et al. eggNOG-mapper v2: functional annotation at the metagenomic scale. Mol Biol Evol 38:5825–5829 (2021). doi:10.1093/molbev/msab293
6. Zheng J, Wittouck S, Salvetti E, et al. A taxonomic note on the genus Lactobacillus. Int J Syst Evol Microbiol 70:2782–2858 (2020). doi:10.1099/ijsem.0.004107
7. Brynildsrud O, Bohlin J, Scheffer L, Eldholm V. Rapid scoring of genes in microbial pan-genome-wide association studies with Scoary. Genome Biol 17:238 (2016). doi:10.1186/s13059-016-1108-8
8. Benjamini Y, Hochberg Y. Controlling the false discovery rate. J R Stat Soc B 57:289–300 (1995). doi:10.1111/j.2517-6161.1995.tb02031.x