Interpretable Machine Learning Accurately Classifies Endangered Austrian Cattle Breeds Using Genomic Markers

Researchers have developed an interpretable machine learning framework combining Light Gradient Boosting Machine and Random Forest classifiers with SHapley Additive exPlanations to accurately identify endangered Austrian cattle breeds and prioritize informative single nucleotide polymorphisms. Accurately classifying livestock breeds is an essential prerequisite for managing endangered populations, particularly when pedigree records are missing, because reliable identification enables targeted breeding programs aimed at population recovery. In traditional livestock conservation, distinguishing closely related populations using standard statistical pipelines often involves multi-step dimensionality reduction techniques such as principal component analysis or fixation index calculations that can miss complex non-linear interactions across genetic markers. By bypassing multi-stage pre-filtering and applying tree-based algorithms directly to comprehensive genomic profiles, the study demonstrated that machine learning can deliver high predictive accuracy while uncovering the specific genetic markers that underpin population differentiation.

The study analyzed genomic data from 6,850 cattle representing 11 endangered Austrian breeds, drawn from an initial cohort of 6,984 animals across commercial and locally conserved populations. Tissue and whole-blood samples were collected by certified veterinarians and local breeding organizations following standard animal ethics protocols, after which genomic DNA was extracted using standard phenol-chloroform protocols or commercial isolation kits. The animals were genotyped using the Illumina BovineSNP50 BeadChip, which originally yielded 53,218 autosomal markers. A quality control pipeline executed in PLINK version 1.9 removed individual samples with call rates below 90 percent as well as single nucleotide polymorphisms with call rates under 90 percent, minor allele frequencies below 0.05, or significant deviations from Hardy-Weinberg equilibrium at a threshold of p less than 1 multiplied by 10 to the power of negative 6. Following quality filtering, remaining missing genotypes were phased and imputed within each breed using Beagle version 4.6, yielding a final dataset of 47,418 high-quality autosomal single nucleotide polymorphisms across 6,850 individuals for downstream analysis.

To rigorously test model stability and mitigate data-partitioning bias, the researchers divided the complete genomic dataset into an 80 percent training and validation partition comprising 5,482 individuals and a 20 percent independent test partition containing 1,368 animals. This data-splitting protocol was repeated across six independent random seeds, numbered one through six, ensuring that model training, hyperparameter optimization, and evaluation were independently replicated. As a baseline benchmark, the team implemented a generalized linear model via the GLMNET framework in R version 4.5.1, utilizing elastic net regularization where the penalty parameter lambda was chosen through five-fold cross-validation based on minimum cross-validated multinomial deviance. Tree-based architectures, specifically Random Forest and Light Gradient Boosting Machine, were subsequently trained using Bayesian optimization combined with repeated five-fold cross-validation, using validation accuracy as the objective for Random Forest and multi-class log loss minimization for LightGBM.

The model training configurations addressed the severe dimensional imbalance characteristic of high-density genomic data, where feature counts vastly exceed sample sizes. For the Random Forest models, Bayesian optimization searched hyperparameter spaces spanning mtry values between 8,000 and 14,000, minimum node sizes from 1 to 10, and sample fractions ranging from 0.7 to 1.0, culminating in final ensemble models comprising 1,000 trees evaluated on the independent test sets. For LightGBM, the hyperparameter search covered leaf counts between 2 and 25, learning rates from 0.001 to 0.1, and L1 regularization strengths from 0.01 to 0.1. LightGBM incorporated algorithmic strategies including Gradient-based One-Side Sampling, which prioritizes instances with larger gradients that provide greater information gain, and Exclusive Feature Bundling, which groups mutually exclusive features. For LightGBM evaluation, five fold-specific models generated class probability estimates across the independent test partitions, which were averaged across folds before test accuracy was computed.

Evaluations across the six independent test sets revealed that tree-based algorithms consistently surpassed the linear baseline. The penalized linear model GLMNET attained a mean classification accuracy of 0.817, while the Random Forest classifier achieved a mean accuracy of 0.837 across the seeds. The tuned LightGBM model delivered the highest overall baseline performance, achieving an average classification accuracy of 0.842, with performance values ranging from 0.835 in seed 3 to 0.849 in seed 4. These findings confirmed that gradient-boosted decision trees maintain robust predictive capability across independent random data initializations, effectively handling the high dimensionality of tens of thousands of genomic markers without succumbing to fatal overfitting.

To interpret the classification mechanisms and evaluate the influence of individual genetic variants, the investigators implemented the TreeSHAP algorithm, an exact model-specific method grounded in cooperative game theory that computes Shapley values in polynomial time. By evaluating the marginal contribution of each marker across feature subsets, TreeSHAP accounts for non-linear interactions among genetic loci. Global marker importance was quantified by taking the mean absolute SHAP value across all samples to measure the overall magnitude of each marker’s contribution regardless of direction. SHAP values were calculated on the training-validation sets across all six repetitions, and markers were ranked by descending mean absolute SHAP value to isolate the top 500 informative variants, whose breed-specific relevance was statistically verified through Kruskal-Wallis testing at a significance threshold of p less than 0.001.

The researchers subsequently tested whether targeted panels of informative markers could match or exceed the predictive performance of the complete genomic array. When the top 500 single nucleotide polymorphisms selected by LightGBM-derived SHAP values were used as the input features, the LightGBM classifier reached an average test accuracy of 0.880, representing an improvement over the 0.842 accuracy obtained with the complete set of 47,418 markers. In parallel, evaluating the subset of markers with non-zero coefficients identified by the GLMNET model yielded an average classification accuracy of 0.857. These comparative tests demonstrated that pruning non-informative markers based on feature-importance metrics captures the primary predictive signals within the genome, reducing data dimensionality while simultaneously enhancing assignment accuracy on unseen test samples.

To assess the validity of the game-theoretic importance metrics, the team compared the SHAP-derived rankings against traditional population-genetic statistics calculated using Weir and Cockerham’s fixation index. Spearman’s rank correlation analysis between locus-specific fixation index values and SHAP importance scores demonstrated substantial concordance. When comparing the top 50 markers identified by each respective method, the researchers found an overlap of 19 single nucleotide polymorphisms, yielding a correlation of 0.829. Broader rank correlation analyses across marker subsets revealed a Spearman correlation of 0.571 for the top 50 SHAP-identified variants and 0.598 for the top 50 fixation index variants, confirming that SHAP values reliably identify biologically and genetically informative loci while accounting for complex multi-marker effects.

Detailed examination of the top-ranking single nucleotide polymorphisms revealed both global genetic discriminators and highly specialized breed-specific markers. A key global marker, ARS-BFGL-NGS-97995, showed high discriminatory power across several populations, stratifying major breed clusters and distinguishing groups such as Waldviertler Blondvieh, Tuxer, and Murbodner. In contrast, other loci displayed strong positive SHAP effects restricted to individual breeds, including DIAS-308 for Fleckvieh, ARS-BFGL-NGS-843 for Brown Swiss, ARS-BFGL-NGS-3513 for Original Braunvieh, and Hapmap40502-BTA-17,475 for Kärntner Blondvieh, all exhibiting extreme statistical significance at p less than 0.001. Pairwise overlap comparisons among the top 100 positive-impact markers per breed showed substantial shared genomic signals between specific pairs, such as Brown Swiss and Fleckvieh sharing 68 markers, Murbodner and Ennstaler Bergschecken sharing 67 markers, and Tuxer and Grauvieh sharing 49 markers.

Biological annotation using the biomaRt package in R version 4.5.1 connected top-ranking markers to annotated genes within 50-kilobase genomic windows in Ensembl release 112, based on the ARS-UCD1.2 bovine genome assembly. The globally influential marker ARS-BFGL-NGS-97995 on chromosome 3 mapped near OR6K4, a member of the olfactory receptor family linked to mammalian foraging, predator avoidance, appetite regulation, and livestock production traits such as feed intake and weight gain. In Grauvieh cattle, the discriminative variant ARS-BFGL-NGS-20141 was positioned adjacent to VDAC1, a gene involved in mitochondrial function and metabolic regulation. In Original Braunvieh, ARS-BFGL-NGS-3513 was located near the dopachrome tautomerase gene DCT, which influences coat color variation. In Fleckvieh, the discriminative marker DIAS-308 mapped near SMIM18, while in Kärntner Blondvieh, the marker Hapmap40502-BTA-17,475 was located close to the olfactory receptor gene OR10C1.

The authors noted important methodological scope conditions and operational limitations regarding the interpretation of their findings. The empirical results reflect the specific dataset of 6,850 animals with confirmed breed records, which does not encompass admixed or crossbred individuals that frequently appear in practical agricultural settings. Additionally, because genotype imputation using Beagle version 4.6 was conducted within breeds prior to partitioning the data into training, validation, and test subsets, the authors cautioned that this imputation step may have introduced slightly optimistic classification accuracies. Future research will need to incorporate broader genomic datasets containing admixed animals and explore the precise functional mechanisms of the prioritized variants. Nevertheless, the combination of LightGBM and SHAP provides a transparent decision-support framework capable of assisting conservationists in managing genetic diversity and verifying breed purity in conservation herds.

Subject of Research: Machine learning and genomic marker prioritization for endangered cattle breed classification

Article Title: Interpretable machine learning for cattle breed classification and SNP prioritization

Article References: Atrian-Afiani, F., Mészáros, G., Sölkner, J., & Waldmann, P. (2026). Interpretable machine learning for cattle breed classification and SNP prioritization. Genetics Selection Evolution, 58(1), Article 37. https://doi.org/10.1186/s12711-026-01056-7

Image Credits: AI Generated

DOI: 10.1186/s12711-026-01056-7

Keywords: machine learning, LightGBM, Random Forest, SHAP, single nucleotide polymorphisms, genomic classification, livestock conservation, Austrian cattle breeds, Fleckvieh, Brown Swiss, Original Braunvieh, Grauvieh, TreeSHAP, fixation index, Beagle imputation, Illumina BovineSNP50, OR6K4, DCT gene