From Genome to Habitat: Hierarchical Machine Learning for Ecological Classification of Salmonella enterica Isolates

1. Introduction

Salmonella enterica is the causative agent of salmonellosis, a zoonotic infection transmitted through contaminated food or water [1]. Disease severity ranges from mild gastroenteritis to invasive systemic infection that can cause multiple-organ failure. The global emergence of antibiotic-resistant strains complicates treatment and poses a serious public health challenge, particularly in developing countries.

This species comprises numerous serovars, many of which were previously classified as separate species [2]. These diverse pathogenic lineages differ significantly in disease manifestation, host specificity, transmission routes, and environmental persistence. Thus, S. enterica serovars Typhi and Paratyphi spread primarily between humans through food or water contaminated with human feces. They commonly cause systemic enteric fever rather than gastroenteritis [3,4,5]. Sources of infection are often asymptomatic individuals who chronically carry the pathogen. By contrast, S. enterica serovar Typhimurium infects numerous animal species and is transmitted to humans through meat, particularly pork and poultry. It can also spread through contaminated water, fresh produce, contact with animals, and the broader farm environment [6,7].

Serovars Enteritidis, Heidelberg, Infantis, Kentucky, and Hadar are also strongly associated with eggs and poultry meat [8,9], whereas serovars Dublin, Derby, and Choleraesuis are more frequently associated with cattle, pork, and industrial meat-production chains [10,11]. Strains circulating in livestock and poultry production systems frequently develop antibiotic resistance because of the inappropriate use of antibiotics [12,13,14,15]. Serovars Montevideo, Senftenberg, Agona, Oranienburg and Typhimurium are notable for their persistence in low-moisture foods and food-processing environments, including nuts, seeds, spices, cereals, chocolate, and animal feed [16,17]. Their principal survival strategy is a strong tolerance of desiccation, which enables prolonged persistence without active growth [18,19].

Serovars Newport, Saintpaul, Javiana, and Braenderup have caused outbreaks associated with tomatoes, leafy vegetables, melons, sprouts, and other fresh produce [20,21]. Contamination may occur through irrigation water, manure, soil, wildlife, or handling [22]. Some of these strains can attach to plant surfaces, colonize the rhizosphere, or enter internal plant tissues [23].

Many apparently sporadic infections cannot be attributed to a specific source, commonly because the contaminated food or relevant environmental exposure is no longer available when the illness is investigated. Therefore, serovar identity may indicate probable reservoirs and survival strategies but does not, by itself, establish the source of an individual infection. For example, salmonellosis outbreak in Brazil caused by S. enterica serovar Newport was transmitted to humans from parrots [24].

Advances in next-generation sequencing (NGS) technologies have provided more detailed insights into the genetic diversity and plasticity of different S. enterica lineages [25]. Multilocus sequence typing (MLST) has revealed numerous sequence types (STs) with distinct geographical distributions [26,27]. For example, S. enterica serovar Typhimurium ST313 has been identified as a major cause of invasive non-typhoidal salmonellosis, particularly among young children and immunocompromised individuals in sub-Saharan Africa [28]. Another study reported that S. enterica serovar Enteritidis ST11 was the most prevalent sequence type among hospital isolates in Kazakhstan [29].

Although serotyping and genotyping may provide valuable information about possible sources of infection and support epidemiological surveillance, ecological classification has not been a primary objective of most such studies. We hypothesized that adaptation to particular survival and dissemination strategies creates characteristic patterns of genetic polymorphisms in bacterial genomes, producing distinguishable allelic profiles even among strains that evolved from different lineages. The availability of numerous whole-genome sequences in public databases, together with metadata describing isolation sources and phenotypic characteristics, enables the application of machine-learning methods and neural-network models to the ecological classification of isolates [30,31,32].

However, the application of these computational approaches is hindered by inconsistent standards for recording metadata in public databases, erroneous records resulting from contamination or inaccurate reporting at the time of sampling, and the subjective assignment of bacterial strains to ecological categories. These limitations impede the construction of reliable ground-truth datasets required for supervised neural-network training. Nevertheless, previous studies have demonstrated the feasibility of attributing Salmonella isolates to probable animal hosts or zoonotic sources using machine-learning methods, including neural networks [30,31]. Despite their promising results, these studies did not provide maintained, ready-to-use software that would enable non-specialist users to classify newly sequenced Salmonella isolates. In the present study, we implemented a semi-supervised, genotype-guided hierarchical learning network to predict the probable habitats and ecovariants of S. enterica isolates and to identify the major genes involved in the ecological adaptation of these pathogens.

2. Materials and Methods

In total, 5,893 complete genomes of S. enterica strains, together with all associated metadata records, were retrieved from the NCBI database.

Core genome alignment and a matrix of single-nucleotide polymorphisms were generated using Panaroo v.1.6.0 [33].

The Hierarchical Density-Based Spatial Clustering of Applications with Noise (HDBSCAN) algorithm, implemented in the Python module hdbscan v.0.8.44, was used for unsupervised clustering.

The chi-square test implemented in the scipy.stats module v.1.17.1 was used for feature scoring. To control for multiple testing, p-values were corrected using the Benjamini–Hochberg false discovery rate (FDR) procedure implemented in statsmodels v.0.14.6. Custom scripts were developed in Python 3.14.4. An automatic model-selection script was used to evaluate eight supervised machine-learning algorithms: logistic regression (LR), support vector classification (SVC), Random Forest (RF), decision tree (DT), multilayer perceptron (MLP), k-nearest neighbours (KNN), categorical Naive Bayes (NBayes), and XGBoost. Candidate models were implemented using the scikit-learn Python package v.1.8.0, with XGBoost provided by the xgboost package v.3.4.1. Categorical SNP states were one-hot encoded for all models except NBayes, which used ordinal encoding; MLP inputs were additionally standardized using sklearn.preprocessing.StandardScaler. Candidate models were ranked by mean cross-validation balanced accuracy, and the highest-scoring algorithm was selected. Model parameters were heuristically adjusted according to sample size, feature dimensionality, number of classes, data type, and class imbalance; class weighting was applied when the largest-to-smallest class ratio was ≥2. When stratified cross-validation was not possible, NBayes was used as a conservative fallback for nodes with < 50 samples and RF for larger nodes.

The OpenAI ChatGPT service was used to assist with custom Python scripting and the design of the architecture of the model PKL files.

3. Results

This section may be divided by subheadings. It should provide a concise and precise description of the experimental results, their interpretation, and the experimental conclusions that can be drawn.

3.1. Machine-Learning Algorithm

3.1.1. Genomic Data Collection and Ontology-Based Grouping

The ontology-based classification scheme shown in Table 1 was developed to assign S. enterica strains to two-level ecological categories.

The resulting two-level hierarchy was used as input for the ontology-guided supervision of the ML workflow by assigning composite labels to genome sequences. These labels do not represent genomic ground truth but rather prior biological hypotheses inferred from metadata according to predefined rules. Of the 5,893 S. enterica genomes retrieved from NCBI, 2,406 were assigned to two-level ecological categories through keyword searches. A genome could be assigned to multiple categories when its metadata contained keywords associated with different categories. The resulting dataset comprised 2,408 records, including only a few genomes assigned to more than one category.

3.1.2. Genomic Feature Selection for Machine Learning

The aim of this step was to construct a matrix of single-nucleotide polymorphisms (SNPs) suitable for distinguishing among S. enterica habitats. The selected features needed to be sufficiently diverse to ensure their broad applicability. A feature-discovery cohort was created by randomly selecting five reference strains from each composite category. Because Table 1 defined 20 unique two-level composite categories, 100 strains were selected for the feature-discovery cohort.

The selected genomes were analyzed using Panaroo, a pangenome analysis tool that generates a core-genome SNP matrix for a given set of genomes. SNPs with a minor allele frequency below 5% were removed to exclude variants specific to individual genomes and potential sequencing errors. SNP positions containing non-canonical nucleotide symbols were also removed. In total, 172,397 SNPs were retained.

In parallel, an SNP annotation file was created containing the following fields: a unique SNP ID; the SNP position in the core-genome alignment generated by Panaroo; its position and DNA strand in the reference genome; the observed allelic states (e.g., A/T); the locus tag of the corresponding reference-genome gene, if the SNP was located within a coding sequence; the protein-coding gene annotation, where applicable; the resulting amino acid substitution for nonsynonymous mutations (e.g., N131K); and an 81-bp context sequence containing the variable nucleotide at the central position, flanked by 40-bp upstream and downstream DNA sequences.

Genome assembly GCA_053410235.1 was selected as the reference genome. There was no specific biological rationale for selecting this particular assembly beyond its high level of completeness. Because the reference genome was used exclusively for SNP annotation, its selection did not influence any subsequent ML steps.

An additional file was generated containing the most common allelic state at each selected SNP position, as determined from the core-genome alignment of all genomes. Using this information, the SNP dataset was converted from a nucleotide matrix (A/T/G/C) into a binary matrix (0/1), in which 0 represented the most common allelic state and 1 represented the minor allelic state or states. This conversion resulted in little information loss because only 5,111 positions contained more than two allelic states (approximately 3% of all selected positions). At the same time, representing allelic states numerically facilitated subsequent pattern comparisons.

The number of potentially diagnostic SNPs was substantially reduced during subsequent processing steps. At each reduction step, the SNP annotation and allelic-state datasets were updated to keep their records synchronized. During the first reduction step, the annotation data were used to remove SNPs located in non-coding regions and synonymous substitutions within coding sequences. This step reduced the SNP matrix to 35,838 positions.

3.1.3. Unsupervised Clustering and Further Dimensionality Reduction

It can be assumed that S. enterica strains sharing the same habitat do not necessarily belong to phylogenetically related lineages. Nevertheless, adaptation to similar environmental conditions may result in similar patterns of adaptive genetic changes, while lineage-specific SNP patterns are also retained. These lineage- and environment-specific signals must be distinguished because they may produce conflicting patterns. To address this issue, unsupervised clustering was applied to the matrix of allelic states at 35,838 SNP positions in the feature-discovery cohort of genomes.

For unsupervised clustering, the Hierarchical Density-Based Spatial Clustering of Applications with Noise (HDBSCAN) algorithm was selected because it is well suited to high-dimensional binary datasets for which the number of clusters is not predefined. Its ability to identify noise is particularly relevant because metadata-derived ecological categories are not necessarily expected to correspond closely to genomic distances. The HDBSCAN parameters were configured as follows: Jaccard distance was used because it is appropriate for binary SNP data; the minimum cluster size was set to 10 genomes; and the min_samples parameter was set to 5. The cluster identifiers generated by HDBSCAN were added as the first level of the composite labels, as illustrated below. Vertical bars separate the cluster identifier, hierarchical ecological categories, and genome accession:

1|Environment|Environmental swabs and sponges|GCA_053410235.1

2|Environment|Environmental swabs and sponges|GCA_053411835.1

5|Environment|Environmental swabs and sponges|GCA_053410975.1

In this example, three genomes identified from their metadata as environmental isolates belonging to the same second-level category, “Environmental swabs and sponges,” were assigned to three different genomic clusters (1, 2, and 5), because they represented distinct genetic lineages.

The next step assessed whether SNP allele distributions were associated with the proposed three-level hierarchical classification of S. enterica genomes. For each SNP, a contingency table was constructed from its observed allelic states and the genome groups defined at each hierarchical level. A chi-square test of independence was then performed separately for every SNP at every hierarchy level. To account for multiple testing, the resulting p-values were adjusted using the Benjamini-Hochberg procedure. The script calculated two types of adjusted p-values: a level-specific FDR-adjusted p-value across all SNPs tested at a particular hierarchy level and a global FDR-adjusted p-value across all valid SNP-by-level tests. The global correction was used for feature selection because a SNP was retained if it was significantly associated with the grouping at least at one hierarchical level. SNPs with a global FDR-adjusted p-value below 0.05 were classified as potentially diagnostic features. This analysis retained 21,290 SNPs associated with at least one level of the proposed hierarchical classification.

A small number of significant SNPs would indicate limited genomic support for the proposed classification. Conversely, a large number of retained SNPs indicates a strong statistical association between allelic distributions and at least one hierarchy level.

Further dimensionality reduction was performed by controlling for redundancy among allelic-state patterns and for the genomic proximity of selected SNPs. Although neural-network models do not strictly require independent input features, highly correlated SNPs provide little additional information, may disproportionately weight particular genomic signals, and increase model complexity.

First, SNPs were grouped into genomic neighborhoods according to their coordinates in the reference genome. SNPs were assigned to the same neighborhood when consecutive SNP positions were separated by no more than 80 bp. Because each SNP is subsequently identified using an 81-bp BLASTN query sequence containing 40-bp upstream and downstream flanking regions, closely located SNPs may have overlapping query sequences and may interfere with reliable allele calling. Therefore, only the statistically highest-ranked SNP was retained from each genomic neighborhood.

Redundancy among the remaining SNPs was then assessed by comparing their allelic-state patterns across all genomes in the feature-discovery cohort. Pairwise dissimilarity was calculated as the proportion of genomes with mismatching allelic states, equivalent to the normalized Hamming distance. The SNPs were grouped using complete-linkage hierarchical clustering with a maximum mismatch threshold of 10%. Consequently, every pair of SNP patterns within a redundancy cluster was at least 90% identical. Up to five SNPs were retained from each cluster, prioritizing those with the lowest global FDR-adjusted p-values; level-specific FDR-adjusted and raw p-values were used as successive tie-breakers. Retaining several representatives from each redundancy cluster provides robustness when genomic regions containing otherwise informative SNPs are absent or cannot be reliably analyzed in newly sequenced genomes.

These filtering steps reduced the feature set to 183 SNPs, providing a manageable input matrix for subsequent ML analyses. However, the final classifier determines the allelic state of each SNP through a separate BLASTN search, making the analysis of all 183 loci time demanding. The remaining SNPs were therefore ranked according to their statistical significance, and the 100 highest-ranked SNPs were retained for model development.

3.1.4. SNP Calling and Machine Learning

Allelic states at the 100 selected diagnostic SNP positions were determined for all 2,406 S. enterica genomes assigned to ecological categories by aligning SNP-context sequences against the genome assemblies using BLASTN. An additional hierarchy level representing genomic lineage was generated by unsupervised clustering with HDBSCAN, as described above. The resulting matrix was divided into non-overlapping training and test datasets. The test dataset comprised 554 genomes randomly selected to represent the different ecological categories, whereas the remaining genomes were assigned to the training dataset.

A separate supervised classification model was trained at each internal node of the hierarchical classification tree to distinguish among its immediate child nodes. Candidate models were compared using mean balanced accuracy obtained through stratified cross-validation, with a maximum of five folds. The number of folds was reduced when required by the size of the smallest child class. Balanced accuracy was used to reduce bias caused by unequal representation of ecological categories. The highest-scoring algorithm was selected independently at each eligible node and was subsequently trained using all training genomes assigned to that node. Supervised models were trained only at splitting nodes containing at least 10 genomes and two or more child classes. Nodes containing fewer than 10 genomes were classified using allele-profile similarity rather than a fitted supervised model. For nodes with only one possible child, that child was assigned directly without model fitting. Allele-frequency profiles were calculated and stored for every node and each of its immediate children. These profiles provide a fallback when supervised-model training fails or when a prediction is considered unreliable. In particular, a model prediction is rejected when its normalized Shannon entropy exceeds 0.80, indicating an uncertain probability distribution across the candidate child classes; classification then proceeds according to similarity between the query SNP profile and the stored child profiles. The same profile-based procedure may be used when missing SNP calls prevent reliable model-based prediction.

The resulting hierarchical ensemble of supervised machine-learning models was packaged in a single model.pkl file. In addition to the fitted node-specific models, the file contained the hierarchical network structure; SNP-context sequences used for BLASTN searches; common allelic states; SNP annotations; allele-frequency profiles for all nodes and their immediate children; model-selection parameters and scores; training-set accuracy values; the reference-genome identifier; the original training matrix; and descriptive metadata. This self-contained model package was used to classify genomes in the test dataset and can also be applied to newly sequenced S. enterica genomes.

Evaluation using the independent test dataset resulted in correct habitat classification for 52% of the 554 genomes. The initial model was subsequently applied to all 2,406 ecologically categorized genomes to identify recurrent patterns of cross-identification. For each genome, all successful terminal predictions were collected from the classification reports. Within each original HDBSCAN cluster, genomes sharing the same complete set of terminal predictions were assigned to the same cross-group. Genomes without a successful terminal prediction were assigned to a separate cross-group with the signature “Not identified.” Cross-group identifiers were numbered independently within each HDBSCAN cluster and inserted between the top-level HDBSCAN clusters and the two ontology-derived ecological levels creating a four-level hierarchy:

HDBSCAN cluster | cross_group | first-level category | second-level category

The hierarchical ensemble was then retrained using the training dataset and the revised four-level labels. Incorporating the cross-group level increased the proportion of genomes receiving correct terminal habitat predictions from 52% to 62%. The architecture of the final four-level hierarchical classification network stored in the model.pkl file is shown in Figure 1. The figure indicates the supervised machine-learning algorithm selected at each internal node and its training accuracy. This value represents the proportion of genomes assigned to the correct immediate child node when the fitted model was applied to the same data used for its training. An accuracy of 1.0 indicates complete separation of the child classes in the training dataset, whereas lower values indicate that their SNP profiles overlap and are more difficult to distinguish. Nodes without a fitted model were classified using allele-profile similarity and therefore have no model-training accuracy value. Because these values were calculated from the training data, they describe model fit rather than predictive performance on independent genomes and may overestimate generalization accuracy.

3.2. Salmonella Enterica Habitat Classifier

3.2.1. Program Overview

The trained model was incorporated into a software framework developed in Python 3. The resulting program, Salmonella_classifier, is available from the SeqWord Classifiers repository (https://github.com/SeqWord/Classifiers).

The program does not require a formal installation and can be run after been downloaded and the required Python dependencies have been installed. The program directory contains the main script, classifier.py, together with the ‘input’, ‘output’, ‘model’, and supporting program directories. Genome assemblies in FASTA or GenBank format may contain either a single chromosome sequence or multiple contigs. GZ-archives of these files can also be used as an input. Input files should be placed in a project subfolder newly created by the user in folder ‘input’. Classification can then be initiated using the following command:

python3 classifier.py new_project_folder

The program saves the classification results as text files in folder ./output/new_project_folder created by the program.

Each query genome is listed in the output file, followed by an indented list of predicted ecological categories ranked according to their classification scores. A query genome may reach multiple terminal nodes, producing a composite profile of predicted ecological associations. If none of the possible terminal categories satisfies the classification criteria, the query is reported as “Not identified”.

Users can adjust several additional command-line parameters that control the sensitivity and specificity of predictions, as described on the project website.

3.2.2. Validation of the Program Using Newly Sequenced Genomes

We applied the program to several recently sequenced S. enterica genomes described in our previous publication [29]. These strains were isolated from clinical samples collected in Kazakhstan, and most belonged to MLST ST11. The ST11 isolates were genetically homogeneous; however, one strain, 19S, exhibited a distinct genome-wide adenine-methylation pattern. We hypothesized that this alternative methylation pattern might be associated with adaptation to a different habitat or life-cycle strategy. The classifier results were consistent with this hypothesis. 19S was assigned a combination of probable ecological categories that differed from those predicted for the other ST11 isolates, which were grouped together by the classifier as likely waste water and clinical isolates. Conversely, the program suggested an association of the strain 19S livestock and food chain environment. The predictions for 19S and a representative ST11 strain, 1S, are compared in Table 2.

These results indicate that the classifier distinguished 19S from the otherwise genetically homogeneous ST11 group. Although this observation does not establish a causal relationship between methylation and ecological adaptation, it provides a basis for further investigation of the distinctive biological characteristics of the isolate 19S

As a specificity control, sequences from the closely related species Escherichia coli, as well as datasets containing mixtures of Salmonella and Escherichia contigs, were submitted to the classifier. In both cases, the program returned “Not identified,” indicating that these sequences did not produce a supported terminal classification within the trained S. enterica hierarchy.

3.3. Functional Analysis of Genes Bearing Marker Mutations

Assuming that the patterns of habitat-specific mutations are formed due to adaptive changes to specific environments and pathogen survival strategies, an analysis of functions of genes bearing these mutations can tell great deal about molecular mechanisms of adaptation. The analysis of the distribution of 183 top-selected marker SNPs representing non-synonymous mutations in protein-coding genes showed non-random scattering of these polymorphisms across housekeeping proteins summarized in Table 3.

4. Discussion

4.1. Classification Model Design

This study presents a semi-supervised algorithm for constructing a hierarchical framework of trained models to predict the probable ecological habitats of S. enterica isolates from their SNP genotypes. The workflow begins with human-interpretable ecological concepts, which are translated into provisional two-level labels. The genomic data are then used to determine whether the proposed categories are genetically distinguishable. The initial two-level hierarchy therefore provides ontology-guided weak supervision for genome classification. An additional level generated through unsupervised HDBSCAN clustering accounts for genomic relatedness and helps stratify lineage-associated variation before ecological associations are evaluated.

The resulting three-level hierarchy should not be regarded as ground truth. Rather, it represents a prior biological hypothesis inferred from metadata according to expert-defined rules and supplemented with genomic clustering [34] The proposed workflow circumvents the requirement of many supervised machine-learning methods for a fully validated ground-truth dataset at the beginning of the analysis. Instead, it evaluates the initial hypothesis by testing whether SNP allele distributions are significantly associated with the proposed hierarchical categories. The number of SNPs passing the global FDR threshold, as well as their level-specific FDR-adjusted and unadjusted p-values, depends on the extent to which the proposed hypothesis corresponds to reproducible genomic patterns. An arbitrary or biologically uninformative grouping would be expected to yield none or few significant features at the ecological hierarchy levels, whereas enrichment of significant associations provides evidence that the proposed categories capture reproducible genomic structure.

After establishing that the proposed hierarchy was associated with a substantial set of genomic features, we performed further dimensionality reduction. Principal component analysis (PCA) is commonly used to reduce a large set of correlated features to a small number of composite variables [31]. However, the principal component eigenvalues cannot usually be attributed directly to individual polymorphisms, genes, or other readily interpretable biological characteristics. PCA was therefore unsuitable for the aims of this study, which included identifying candidate molecular mechanisms underlying the adaptation and dissemination of Salmonella pathogens, in addition to assigning isolates to probable ecological niches.

Instead, a series of feature-filtering procedures was used to select informative SNPs while controlling for redundant allele patterns and the genomic proximity of polymorphic sites. Allelic states at the resulting subset of informative SNP positions were then determined in the complete collection of categorized genomes. This targeted approach substantially reduced the computational cost of feature extraction while preserving the biological interpretability of individual SNPs and their associated genes.

The resulting matrix of allelic states for the training genomes was used to construct a hierarchical ensemble of supervised classifiers, with a separate model trained at each eligible internal node (Figure 1). Different machine-learning algorithms were evaluated independently at each node, and the best-performing algorithm was selected to distinguish among its immediate child nodes. Nodes containing insufficient training data, nodes for which model fitting was unsuccessful, and uncertain predictions were handled using stored allele-frequency profiles. Efficacy of this approach was demonstrated previously in the model trained for identification of patterns of antibiotic resistance of Mycobacterium tuberculosis isolates [35].

Cross-validation and subsequent application of the initial model revealed recurrent cross-identification among some ecological categories. These patterns were used to introduce an additional cross-group level between the top-level genomic clusters and the two ontology-derived ecological levels. Similar approaches of model improving were discussed previously [36]. This post-training hierarchy refinement aimed at latent class discovery resulted in four-level hierarchy with an improved correspondence between the predicted and metadata-derived terminal categories.

The resulting program, which is freely available from a public GitHub repository, returns a ranked set of ecological categories assigned to each query genome, as illustrated in Table 2. These categories should not be interpreted merely as mutually exclusive alternatives. Collectively, they constitute a predicted ecological profile that may reflect multiple habitats, reservoirs, transmission routes, or clinical sources associated with genetically similar isolates. For example, the output for strain 1S in Table 2 associates its SNP profile with categories representing sediment and waste, livestock, food-chain products, and several types of clinical samples. These predictions suggest possible ecological and epidemiological associations; however, do not, by themselves, demonstrate a specific transmission pathway or clinical outcome. In Figure 1, the terminal ecological categories are represented by circles of different colors. The combination of terminal categories reached by a query genome describes a predicted survival and dissemination strategy associated with specific genomic profiles.

4.2. Molecular Mechanisms of Adaptation of S. enterica to Versatile Ecological Niches

This study presents a semi-supervised algorithm for constructing a hierarchical framework of trained models to predict the probable ecological habitats of S. enterica isolates from their SNP genotypes. The workflow begins with human-interpretable ecological concepts, which are translated into provisional two-level labels. The genomic data are then used to determine whether the proposed categories are genetically distinguishable.

The analysis of functional genes bearing mutations associated with ecological adaptation provides important insights into the strategies and molecular mechanisms exploited by these pathogens (Table 3). Of particular interest are the genes:

iclR transcriptional regulator: glyoxylate shunt ← aceA ⇄ icd → aceK → TCA

which form a tightly connected regulatory and metabolic module controlling the balance between the TCA cycle and the glyoxylate shunt, particularly when bacteria utilize acetate or fatty acids as major carbon sources. The glyoxylate shunt bypasses the two CO₂-generating steps of the TCA cycle, thereby conserving carbon skeletons that can be channeled into anabolic pathways through gluconeogenesis [37]. Microorganisms with an inactive or absent glyoxylate shunt may utilize lipids for energy production but have a limited capacity to convert lipid-derived acetyl-CoA into biomass. This metabolic capability may therefore contribute to the preferential colonization and persistence of pathogens in particular host organs and tissues [38]. For example, the ability of Mycobacterium tuberculosis and Pseudomonas aeruginosa to persist in lipid-rich host environments, including the lungs, has been associated with the activity of the glyoxylate shunt and utilization of host-derived fatty acids [39,40]. Other host compartments, such as lymphoid tissues and the brain, can also provide lipid-rich nutritional niches, whereas the intestinal tract and bloodstream generally offer greater access to carbohydrates and other readily metabolizable carbon sources. Environmental habitats such as wastewater and sediments may likewise be relatively depleted in readily available carbohydrates while containing diverse organic acids and lipid-derived compounds. Adaptation to these contrasting nutritional environments requires precise regulation of catabolic and anabolic pathways, particularly the balance between glycolysis, the TCA cycle, the glyoxylate shunt, and gluconeogenesis. Strain-specific mutations in proteins associated with anaerobic respiration, including the NarX nitrate/nitrite sensor, the NarK nitrate/nitrite transporter, and nitrate reductases, which catalyses a key reaction in nitrate respiration, contribute to the versatility of Salmonella adaptation to diverse environments [41,42].

There is a broad collection of carbohydrate-related proteins bearing strain-specific mutations in different S. enterica isolates. These proteins are involved in the metabolism of mannose, sorbitol/glucitol, chitobiose/N-acetylglucosamine, fucose, and β-glucosides, as well as in central carbohydrate metabolism. This finding suggests that a substantial fraction of strain-specific diversification is associated with nutritional niche specialization rather than merely with modifications of core glycolysis. Salmonella pathogens disseminated via contaminated plants can actively degrade plant-derived polysaccharides including cellulose, providing substrates for energy production and potentially facilitating access to internal plant tissues [23]. On the other hand, mutations in fucK and nagK/chbC are particularly interesting in Salmonella, because fucose and amino sugars are abundant in warm-blooded host-associated environments, where they occur in host glycans and mucus [43,44].

The type III secretion system is one of the major virulence factors of Salmonella, used by the pathogen to inject effector molecules directly into the cytoplasm of host cells, including those of warm-blooded animals and plants. The large number of type III secretion system components bearing environment-associated mutations (Table 3) highlights the importance of this secretion system in the adaptation of the pathogen to different ecological niches [45].

Several other adaptation mechanisms can be inferred from the list of genes bearing strain-specific mutations in Table 3, including those encoding enzymes involved in amino acid utilization and nitrogen metabolism, aromatic compound utilization, the synthesis and salvage of important vitamins and cofactors, as well as global regulators of gene transcription.

4.3. Prospects for Improving the Accuracy of the Classifier

The performance of the current version of the classifier can be improved through reciprocal analysis and refinement of the initial grouping of genomes into two-level ecological categories. Although the accuracy values assigned to individual nodes of the classification framework (Figure 1), which describe model fit, should not be interpreted as general estimates of overall system performance, they can identify parts of the framework that require improvement. Such improvements may be achieved by adding more genomes from the relevant ecological categories to the training dataset, restructuring the ecological grouping of these genomes, or searching for additional latent layers within the classification framework.

An additional source of information for model training is the matrix of accessory genes present in subsets of genomes. This matrix is automatically generated by Panaroo, a program designed for pangenomic analysis, and is therefore readily available. It uses a binary 0/1 encoding similar to that employed in the current version of the classifier for alternative SNP alleles. Thus, the SNP and accessory-gene matrices can, in principle, be combined. The applicability of accessory-gene presence/absence matrices to bacterial classification has been demonstrated in previous studies (Lupolova et al., 2019). However, we decided against using this matrix in the current study and reserved its incorporation for future updates for several reasons.

The major concern is the difficulty in distinguishing the true absence of an accessory gene from a false-negative result caused by incompleteness of the query genome or by a high level of sequence divergence within the sequence context used for BLASTN alignment. The current classifier, which is based on patterns of SNP alleles, treats BLASTN alignment failures as missing data rather than as a specific allelic state. In contrast, when accessory-gene presence/absence is encoded as a binary variable, an alignment failure could be incorrectly interpreted as gene absence. Nevertheless, we do not exclude the possibility that careful selection of conserved, reliably detectable, and ecologically informative accessory genes could improve the accuracy of future versions of the classifier.

5. Conclusions

This study demonstrates that ecological information encoded in heterogeneous and potentially incomplete genome metadata can be combined with genomic variation to construct an interpretable hierarchical framework for predicting the probable habitats and survival strategies of S. enterica. Rather than treating metadata-derived ecological categories as ground truth, the proposed semi-supervised approach uses them as an initial biological hypothesis and subsequently evaluates and refines this hypothesis according to patterns present in the genomic data. Unsupervised clustering, statistical selection of habitat-associated SNPs, redundancy reduction, and node-specific supervised classification together provide a framework capable of accommodating both phylogenetic structure and ecological diversification.

An important outcome of the study is the biological interpretability of the selected genomic markers. Genes bearing habitat-associated nonsynonymous mutations participate in interconnected processes including carbon and nitrogen utilization, the balance between the TCA cycle and glyoxylate shunt, anaerobic respiration, redox metabolism, vitamin and cofactor metabolism, environmental sensing, motility, stress responses, transcriptional regulation, and type III secretion. Their functional diversity is consistent with adaptation of S. enterica to markedly different nutritional, physicochemical, environmental, and host-associated niches.

The developed Salmonella_classifier translates this approach into a freely available tool applicable to newly sequenced genomes. Its architecture also provides opportunities for iterative improvement as additional genomes and better ecological metadata become available. More broadly, the approach illustrates how weak ecological supervision combined with unsupervised genomic structure and interpretable molecular markers can be used to investigate microbial ecological adaptation when a reliable ground-truth classification cannot be established in advance.