Repeatability and genetic architecture of milk urea nitrogen in a multi-breed Romanian dual-purpose cattle population: a genome-wide association study and the impact of phenotype definition

Milk urea nitrogen (MUN) indicates nitrogen utilization efficiency in dual-purpose cattle, and the definition of its per-cow phenotype represents a decisive choice in small two-step genome-wide association studies (GWAS). We characterized MUN repeatability and genetic architecture in 406 Romanian dual-purpose cows from multiple breeds (324 Simmental, 57 Brown Swiss, and 25 Holstein), genotyped at 42,440 single nucleotide polymorphisms (SNPs). MUN was moderately repeatable (intraclass correlation of 0.243 and 0.388 in the full and genotyped samples, respectively). Under a well-calibrated genomic mixed model (13,722 test-day records from 406 cows; λ = 1.034), no SNP reached genome-wide significance, and three were suggestive, none of which survived false discovery rate control (minimum q = 0.723). A chromosome 18 signal near ZNF423, significant under a conditional-residual phenotype, proved to be driven by phenotype definition rather than biology: switching to a marginal adjusted mean, with the covariate model held fixed, collapsed the signal from P = 4.42 × 10 -7 to 0.060 through standard-error inflation rather than a change in the effect size, with residual breed-farm structure serving as a secondary cause. Narrow-sense heritability was not separately identifiable in this mostly unrelated cohort, and established MUN loci from larger studies showed no detectable transferability, consistent with the limited statistical power of this cohort. We provide the first Romanian MUN data, a rigorous null scan, and a demonstrated caution regarding phenotype definition for two-step GWAS.

Dual-purpose and dairy cattle production systems face mounting pressure to improve nitrogen (N) utilization efficiency while minimizing environmental impact. The efficiency with which lactating cows convert dietary N into milk N is generally low, typically averaging approximately 25%, meaning that a substantial proportion of ingested N is excreted in urine and feces (Lavery and Ferris, 2021). Nitrogen in manure can be lost to the environment through ammonia volatilization, nitrate leaching into groundwater, and nitrous oxide emissions, all of which carry significant consequences for air quality, water quality, and climate change (Chadwick et al., 2018). Improving N utilization efficiency is therefore a central objective for sustainable cattle production, with implications for both the environment and economic performance through reduced feed protein costs.

Milk urea nitrogen (MUN) has emerged as a practical, non-invasive indicator of N utilization efficiency in cattle. MUN reflects the concentration of urea in milk, which equilibrates with blood urea and is closely related to the balance between dietary protein intake and fermentable energy supply in the rumen (Spek et al., 2013). When dietary crude protein exceeds the capacity of ruminal microorganisms to incorporate ammonia into microbial protein, the surplus ammonia is absorbed, converted to urea in the liver, and partitioned into blood, urine, and milk. Substantial between-cow variation in MUN persists even when cows are fed the same ration, indicating that individual physiological and genetic factors contribute to MUN concentration independently of diet (Aguilar et al., 2012). Between-cow differences in MUN appear to be governed substantially by renal urea handling rather than hepatic urea synthesis (Prahl et al., 2022, 2023), which bears on how a polygenic architecture for the trait should be interpreted. Because MUN is routinely measured in many national milk recording programs at no additional cost, it is an attractive target for both management interventions and genetic improvement.

In addition to its role as a nutritional indicator, MUN has been associated with reproductive performance; a meta-analysis demonstrated 43% lower odds of conception when blood urea exceeded 7.0 mM (or milk urea exceeded 420 mg/L) (Raboisson et al., 2017). MUN is a repeatable trait with low-to-moderate reported heritability that is strongly population-dependent. Published estimates span a wide range across populations and models, from approximately 0.06 in a single-step repeatability model in Chinese Holsteins (Ma et al., 2023) to approximately 0.20 in Italian Brown Swiss (Bobbo et al., 2020) and approximately 0.22 to 0.28 in grazing New Zealand cattle (Beatson et al., 2019), with a meta-analytic synthesis confirming low-to-moderate genetic parameters for the trait (Ghavi Hossein-Zadeh, 2024). Selection for lower MUN has been proposed as a route to reduce urinary N excretion in grazing systems (Beatson et al., 2019). These findings support the rationale for considering MUN in genetic improvement strategies, while indicating that the magnitude of the heritable component is both modest and population-dependent.

Efforts to elucidate the genomic architecture of MUN through quantitative trait locus (QTL) mapping and genome-wide association studies (GWAS) have identified several chromosomal regions of interest, though with limited overlap across studies. Bouwman et al. (2010) conducted a genome-wide scan in Dutch Holstein-Friesian cows and identified suggestive QTL on Bos taurus autosomes (BTAs) 1, 6, 21, and 23, each explaining only 2% to 3% of phenotypic variance. van den Berg et al. (2022), using a meta-GWAS with imputed whole-genome sequence data from up to 33,873 dairy cows in Australia and New Zealand, detected 14 QTL for MUN, with the most significant signals near the casein gene cluster (BTA 6), PAEP (BTA 11), and DGAT1 (BTA 14). Cecchinato et al. (2014) used a candidate-gene approach in Italian Brown Swiss and identified associations between milk urea and polymorphisms in genes related to milk composition. Honerlagen et al. (2021) identified genomic regions influencing N-metabolism traits in 371 lactating Holstein-Friesians, a sample comparable in size to that of the present study. Ariyarathne et al. (2021) performed a GWAS for milk urea in 490 mixed-breed New Zealand grazing dairy cows, a close multi-breed, 50K analogue to this work, and identified several candidate genes. Ma et al. (2023) performed a single-step GWAS in Chinese Holsteins and found significant single nucleotide polymorphisms (SNPs) on BTAs 11, 12, 14, 17, and 18. The limited concordance among these studies reflects the polygenic nature of MUN, the typically small effect sizes of MUN QTL, and the importance of investigating diverse populations and analytical approaches.

The majority of MUN GWAS to date have been conducted within single breeds, predominantly Holstein. Multi-breed analyses offer the potential to increase sample size and to detect shared QTL, but they also introduce methodological challenges related to population structure, differing linkage disequilibrium (LD) patterns between breeds, and the potential for breed-specific genetic effects (Hayes et al., 2009). Mixed-model approaches that incorporate a genomic relationship matrix (GRM) can account for population stratification and cryptic relatedness simultaneously in multi-breed settings, as the GRM captures both within-breed and between-breed genetic similarity (VanRaden, 2008; Yang et al., 2014). The leave-one-chromosome-out (LOCO) strategy further avoids proximal contamination, whereby the SNP being tested would otherwise be included in the relationship matrix and its effect absorbed into the polygenic term (Yang et al., 2014).

A further methodological choice in such studies is the definition of the phenotype itself. When repeated per-cow measurements are first condensed into a single value through a mixed model, the way the random animal effect is handled determines whether between-cow variation is preserved or inadvertently shrunk. That choice proved decisive here. A factorial experiment was therefore used to isolate its contribution, establish the mechanism through which it operates, and define the conditions under which it can arise in other two-step analyses.

Romania maintains a diverse cattle population that includes both dual-purpose breeds (Romanian Simmental, national name Bălţată Românească, and Romanian Brown Swiss, national name Bruna de Maramureş) and dairy breeds (Romanian Holstein, national name Bălţată cu Negru Românească), with breed composition varying across farms and regions. Hereafter, these are referred to as Simmental, Brown Swiss, and Holstein for brevity and international comparability. Our group has previously reported genome-wide association analyses in Romanian cattle for milk somatic cell score in Romanian Spotted and Romanian Brown cows (Ilie et al., 2021) and for milk protein content in Romanian dual-purpose cattle (Bratu et al., 2025). Despite the economic and environmental importance of N utilization efficiency, no MUN data or MUN GWAS have been reported for Romanian cattle to date.

Accordingly, this study characterizes the genetic architecture of MUN in a multi-breed Romanian cattle population comprising Romanian Simmental, Romanian Brown Swiss, and Romanian Holstein cows, with four objectives: (i) to estimate the repeatability of MUN and assess whether its narrow-sense heritability can be separately identified in this cohort; (ii) to test for genome-wide associations with MUN under a well-calibrated mixed model; (iii) to quantify how phenotype definition, covariate specification, and structural breed-farm confounding influence association and variance-component results in a small, multi-breed context, including through a factorial experiment that isolates the contribution of phenotype definition; and (iv) to test whether genomic regions previously associated with MUN in larger studies are detectable in this cohort. This study provides the first MUN repeatability estimate for a Romanian dual-purpose cattle population, a phenotype-derivation recipe transferable to any small longitudinal livestock GWAS, and a realistic power benchmark for the cohort sizes that future regional consortia would require.

Animals and phenotype data

The study population comprised lactating cows from two farms in Romania: SCDCB Arad (Research and Development Station for Bovine, Arad) and a private farm (Wagner, Timiş County). After intersecting genotyped and phenotyped animals, the final GWAS population comprised 406 cows: 324 Romanian Simmental (Bălţata Românească) (79.8%), 57 Romanian Brown Swiss (Bruna de Maramureş) (14.0%), and 25 Holstein (6.2%). Animals were distributed across farms as follows: SCDCB housed 347 cows (290 Simmental and 57 Brown Swiss), while Wagner housed 59 cows (34 Simmental and 25 Holstein). Only Simmental was present at both farms; Brown Swiss was exclusive to SCDCB Arad, and Holstein was exclusive to Wagner.

The cows from SCDCB Arad had unrestricted access to both water and outdoor paddocks. Feeding was carried out twice daily. Their diet consisted of a daily ration including 15 kg of alfalfa and 15 kg of fresh forage, supplemented with 12 kg of corn silage, 6 kg of alfalfa hay, and 6 kg of concentrates from spring through late autumn. During the winter period, the ration was adjusted to 15 kg of alfalfa and 25 kg of corn silage, along with 6 kg of alfalfa hay and 6 kg of concentrates. At the Wagner farm, cows were managed in an automatic milking system (DeLaval) and received a partial mixed ration based on maize silage, haylage, and brewer’s spent grains, with the concentrate fraction dispensed individually at the milking robot according to milk yield. The daily ration of a lactating cow comprised 25 kg of maize silage, 8 kg of brewer’s spent grains, 4 kg of rye haylage, 4 kg of alfalfa haylage, and 1.5 kg of alfalfa hay, together with a concentrate mixture of 2.2 kg of ground maize, 2 kg of rye, 1.8 kg of wheat, and 1 kg of sunflower meal, supplemented with 250 g of a vitamin-mineral-electrolyte premix, 150 g of sodium bicarbonate, and 50 g of salt.

The trait of interest was milk urea nitrogen (MUN), expressed in mg/dL. Urea is reported by the Romanian official milk recording program as urea in mg/dL, not as mmol/L or as urea nitrogen. Both farms are sampled under that scheme, and both sets of samples are analyzed using the same type of instrument, although sampling at the two farms is administered by different breeding associations, so the values from the two farms are on a common scale. Urea concentrations were converted to MUN using the molecular weight conversion factor MUN = Urea/2.14 (urea has a molecular weight of 60.06 g/mol and contains two nitrogen atoms [28.01 g/mol], leading to the ratio 60.06/28.01 = 2.14). The genotyped cows were extensively recorded: the 406 animals contributed 13,722 test-day MUN records, with a minimum of 2 records per cow, a median of 25, a mean of approximately 34, and a maximum of 147 records.

To account for systematic environmental effects on repeated MUN measurements, a linear mixed model (LMM) was fitted using the lme4 package in R (Bates et al., 2015):

where Lactation stage (three levels: Early < 100, Mid 100–199, Late ≥ 200 days in milk), Season (Winter, Spring, Summer, Fall), Year, and Parity (parities of 5 and above grouped into a single class, with parity capped at 5) were fitted as fixed effects, and Animal was fitted as a random effect assumed to follow Animal m ~ N(0, σ 2 a), with residuals e ijklm ~N(0, σ 2 e), independent of the random animal effects.

Breed and Farm were deliberately excluded from this model because they were handled in the subsequent GWAS through the GRM and a farm covariate. The model was fitted by restricted maximum likelihood (REML) with the bobyqa optimizer (maxfun = 1 × 10 5). A fixed random seed [set.seed(42)] was used throughout; outliers were removed within each breed using the 1.5 × interquartile range (IQR) method, and only animals with at least 2 records after fitting were retained. Records with a non-positive or non-finite urea value were excluded before this step. Because the interquartile rule is symmetric and the lower fence lies below zero in every breed, it removed only high outliers; therefore, a small number of very low values were retained, and no biological plausibility bound was imposed beyond the positivity requirement.

The way the per-cow phenotype was extracted from this LMM requires an explicit statement. The corrected and primary phenotype used throughout this manuscript is the fixed-effects-only adjusted mean, equivalent to the marginal residual. For each record, the fixed-effect prediction was obtained with predict(fit, re.form = NA), subtracted from the observed value, and added back to the grand mean; the resulting marginal residuals were then averaged within each cow. This marginal adjusted mean preserves each cow’s full environment-corrected deviation irrespective of how many records she contributed, and it is the appropriate two-step GWAS phenotype in the absence of a national evaluation from which de-regressed estimated breeding values could be constructed (Garrick et al., 2009; Twomey et al., 2019; Cai et al., 2023).

The conditional-residual phenotype, against which the marginal definition is compared in the factorial experiment (Section 2.6), instead uses the residuals returned by default from lme4. These implicitly subtract the estimated, shrunken animal effect (its best linear unbiased prediction, BLUP) in addition to the fixed effects. The LMM, fixed effects, filters, seed, and outlier rule are otherwise identical; only the residual type differs.

Model diagnostics confirmed the suitability of the adjustment procedure. The Shapiro-Wilk test was applied to the LMM residuals (a random subsample of at most 5,000 residuals, drawn using a fixed random seed) to assess the normality assumption. The intraclass correlation coefficient (ICC), interpreted here as the between-animal repeatability of MUN, was computed from the estimated variance components (Section 2.8).

Genotyping and quality control

Animals were genotyped using the Thermo Fisher Axiom Bovine Genotyping Array (approximately 65K SNPs). Genotype data from two batches were merged into binary PLINK format using PLINK v1.9 (Purcell et al., 2007; Chang et al., 2015), with automated strand orientation correction. SNP coordinates on this array correspond to the UMD3.1 reference assembly (bosTau6).

Quality control (QC) was performed on the 406-animal GWAS population (starting from 64,773 merged SNPs). The following sequential filters were applied using PLINK v1.9: (1) genotype call rate, SNPs with more than 10% missing genotypes were removed (879 SNPs); (2) minor allele frequency (MAF), SNPs with MAF < 0.05 were removed (19,461 SNPs); (3) Hardy-Weinberg equilibrium (HWE), SNPs deviating from HWE at P < 1 × 10 -6 were removed, with HWE tested only within the Simmental subpopulation (N = 324) (389 SNPs); and (4) chromosome filter, non-autosomal SNPs were removed, retaining only bovine autosomes 1 to 29 (1,604 SNPs). The final QC-passed dataset comprised 42,440 SNPs across 406 animals; 88.8% of the retained markers had a minor allele frequency of at least 0.10, and the remaining 11.2% lay between 0.05 and 0.10.

The decision to test HWE within Simmental only, rather than in the combined multi-breed population, was motivated by the Wahlund effect (Wahlund, 1928). When genetically distinct subpopulations with different allele frequencies are pooled, heterozygosity is systematically reduced relative to HWE expectations even if each subpopulation is individually in equilibrium, which can cause spurious HWE deviations and the inappropriate removal of informative loci (Pearman et al., 2022). Whether filtering thresholds should be applied across an entire dataset or separately within biologically meaningful groups is a recognized decision point in SNP filtering (O'Leary et al., 2018). Testing HWE on the largest breed subpopulation at a stringent threshold detects genuine genotyping artifacts while avoiding Wahlund-induced false rejections. The threshold of P < 1 × 10 -6 was chosen so that the filter removes gross genotyping failures without discarding informative variation. At N = 324, a one-degree-of-freedom Hardy-Weinberg test at this threshold has 80% power only against departures of |F| >= 0.32, that is, a heterozygote deficiency or excess of roughly one third relative to expectation, and 50% power at |F| >= 0.27. Departures of that size are characteristic of array cluster-calling failures rather than any demographic or selective process acting within a single breed, so the threshold is calibrated to flag technical artifacts while leaving informative loci in the panel.

Population structure analysis

Population structure was assessed using principal component analysis (PCA) computed on a linkage disequilibrium (LD)-pruned SNP set. LD pruning was performed using PLINK v1.9 (--indep-pairwise 50 5 0.5), and the first 10 principal components (PCs) were extracted from the pruned set; the first 4 were used as fixed covariates in the association model (Section 2.5), and the leading PCs were used to visualize population structure.

Genomic relationship matrix

A genomic relationship matrix (GRM) was computed using GCTA v1.95.0 (Yang et al., 2011) on all 42,440 post-QC autosomal SNPs, following the marker-based definition of VanRaden (2008). The GRM captures realized genetic relationships among individuals, including breed structure, familial relatedness, and cryptic kinship; the interpretation of such SNP-based relatedness and its use for variance-component estimation follow the framework of Speed and Balding (2015).

Genome-wide association analysis

Association analysis was performed using the mixed linear model association method with leave-one-chromosome-out (MLMA-LOCO) implemented in GCTA v1.95.0 (Yang et al., 2011, 2014). In MLMA-LOCO, for each SNP tested on a given chromosome, the GRM is recomputed excluding all SNPs on that chromosome, avoiding proximal contamination (Yang et al., 2014). The primary model, designated Model GS, was specified as:

where y is the vector of corrected (marginal adjusted mean) MUN phenotypes; μ is the overall mean; x j is the vector of genotypes at the SNP being tested, and b j is its allele-substitution effect; W is the matrix of fixed covariates, namely, a binary farm indicator (SCDCB = 0, Wagner = 1) and the first four PCs, with effects c; g is the random polygenic effect with g ~ N(0, G·σ 2 g), where G is the LOCO-GRM; and e is the residual, e ~ N(0, I·σ 2 e).

The method tests b j separately for each SNP in turn. Breed was not included as a fixed covariate because the GRM already captures breed structure through allele sharing patterns (VanRaden, 2008; Yang et al., 2011). The PCs were added to absorb residual breed axes not fully removed by the LOCO-GRM, following the established use of principal components to control residual population stratification in association studies (Price et al., 2006). The number of PCs was fixed at four a priori rather than tuned: with three breeds, the two leading PCs separate the breeds, and the next two capture residual within-breed substructure, so four axes span the structure present in the sample. This a priori choice was examined by a sensitivity analysis over 0, 2, 4, 6, and 10 PCs. This primary model differs deliberately from the farm-only specification: handling breed through the GRM (Misztal et al., 2020) and accounting for residual structure with PCs (Price et al., 2006) is the better-specified choice in this multi-breed sample.

Sensitivity analysis

The robustness of the association results was assessed through four complementary checks. First, genomic inflation was monitored for the primary model and for the Simmental-only model; the genomic inflation factor (λ) was computed as the median observed χ² test statistic divided by the expected median under the null (Devlin and Roeder, 1999). Second, and most importantly, a within-breed sensitivity GWAS was performed on the largest homogeneous breed subgroup (Romanian Simmental, N = 324, comprising 290 cows from SCDCB Arad and 34 from Wagner) using the same Model GS configuration, restricted to Simmental animals, with the Simmental-only GRM recomputed after subsetting the post-QC binary PLINK files. Because Simmental is the only breed present at both farms and this analysis deliberately retains the Simmental cows from both farms, it is the primary deconfounding check: it tests whether any full-cohort signal survives once between-breed allele-frequency differences are removed while a single breed spanning both farms is retained.

Third, to isolate the contribution of phenotype definition to the chromosome 18 signal, a phenotype-by-model factorial (2 × 2) experiment was run on the identical genotypes and animals. A factorial design varies two factors across all their combinations so that the effect of each can be separated from the other. The two factors here were the phenotype (conditional-residual versus corrected marginal adjusted mean) and the covariate model (the farm-only specification versus the corrected Model GS, farm plus four PCs, with breed carried by the LOCO-GRM), giving four combinations: (1) conditional phenotype with the farm-only model; (2) conditional phenotype with Model GS; (3) marginal phenotype with the farm-only model; and (4) marginal phenotype with Model GS, the corrected primary analysis. Each of the four combinations was run as an independent MLMA-LOCO analysis, and the association statistics for AX-124375815 (BTA 18, near ZNF423), along with the genomic inflation factor, were compared across them. Crossing the analytical choices factorially in this way, rather than changing several at once, is the recommended approach for revealing which choice a result actually depends on (Steegen et al., 2016), and it allows the collapse of the signal to be attributed to a single cause rather than confounded across changes. Fourth, the influence of the most closely related pairs was tested by refitting Model GS after removing one member of the pair with the highest relatedness (N = 405) and one member of each of the three pairs with a genomic relationship above 1.0 (N = 403).

Multiple testing correction

Multiple testing was addressed using the SimpleM method (Gao et al., 2008, 2010; Gao, 2011), which estimates the effective number of independent tests (M eff) by accounting for LD among markers. The composite LD correlation matrix was decomposed per chromosome in 1,000-SNP blocks using eigenvalue decomposition, retaining the leading eigenvalues that together explained a declared proportion C = 0.995 of the variance, following the SimpleM calibration of Gao et al. (2008), validated on a genome-wide scale by Gao et al. (2010). This gave M eff = 13,988, a genome-wide threshold of P < 0.05/13,988 = 3.574 × 10 -6, and a suggestive threshold of P < 1/13,988 = 7.149 × 10 -5. As a parameter-free reference, the Bonferroni bound on all 42,440 SNPs (P < 0.05/42,440 = 1.178 × 10 -6) is in fact stricter than the SimpleM threshold; both bounds are applied in Section 3.4.

False discovery rate (FDR) was additionally controlled using the Benjamini-Hochberg (BH) procedure (Benjamini and Hochberg, 1995) applied to the full vector of per-SNP MLMA-LOCO P-values using the p.adjust function in R; under this criterion, a SNP would be declared significant only if its BH-adjusted P-value (q-value) fell below 0.05. FDR controls the expected proportion of false positives among any declared discoveries and complements the genome-wide significance threshold.

Repeatability and heritability

Repeatability was estimated as the ICC from the phenotype LMM (Section 2.1) for both the full milk-recording cohort and the matched genotyped cows. Because the genotyped cohort is only weakly related, the additive genetic and permanent-environment components of variance cannot be separated; we therefore do not estimate narrow-sense heritability and instead treat repeatability as its upper bound (0 ≤ h 2 ≤ repeatability) (Wilson et al., 2010; Huţu et al., 2020). For completeness, supporting genomic variance-component analyses [genomic restricted maximum likelihood (GREML) on the per-animal mean and a record-level animal model in the sommer package] are reported in the Supplementary Materials (Supplementary Tables S2, S3); neither is interpreted as a narrow-sense heritability.

Gene annotation and coordinate conversion

GWAS results were obtained in UMD3.1 (bosTau6) coordinates, corresponding to the genotyping array design. To annotate genes nearest to the leading SNPs using the current bovine reference assembly, SNP positions were converted from UMD3.1 to ARS-UCD1.2 (bosTau9) using the UCSC liftOver tool with the bosTau6ToBosTau9 chain file (Kent et al., 2002; Hinrichs et al., 2006) through the rtracklayer package in R/Bioconductor (Lawrence et al., 2009). All 15 top SNPs were successfully lifted over. Gene queries were performed using the Ensembl REST API (Cunningham et al., 2022), searching for annotated features within ±500 kb of each leading SNP. This was used as a positional rather than functional annotation (Supplementary Table 4). Reference SNP (rs) identifiers were retrieved for the leading SNPs from their ARS-UCD1.2 positions through Ensembl Variation, accepting a match only when the reported alleles agreed (Supplementary Table 4).

Transferability of established MUN loci

To assess whether loci previously associated with MUN in larger populations show any signal in the present cohort, a pre-specified, hypothesis-driven transferability analysis was performed against a fixed reference set assembled before any lookup. The reference set comprised 55 loci or regions compiled from prior MUN studies (Bouwman et al., 2010; Cecchinato et al., 2014; Ariyarathne et al., 2021; Honerlagen et al., 2021; van den Berg et al., 2022; Ma et al., 2023) together with the six canonical urea-cycle genes (CPS1, OTC, ASL, ASS1, ARG1, and ARG2). Of these, 41 had exact source-table base-pair positions, whereas 14 were coarser: the Bouwman et al. (2010) QTL were reported only in linkage centimorgans with wide confidence intervals, the Ariyarathne et al. (2021) loci were 1-Mb Bayes C windows, and the suggestive Honerlagen et al. (2021) SNPs lacked fine positions.

Coordinates were harmonized to a single assembly before lookup. Positions reported in earlier assemblies were lifted to ARS-UCD1.2 (bosTau9) with the rtracklayer package and the UCSC chain file (Hinrichs et al., 2006; Lawrence et al., 2009); 99.7% of our SNPs mapped, and the liftover of the two anchor loci (the BTA 18 lead SNP and a DGAT1 -region SNP on BTA 14) reproduced their expected ARS-UCD1.2 positions to the base pair. For each harmonized locus, the best-proxy SNP was defined as the SNP with the lowest P-value within ±250 kb of the locus position in the Model GS results, thereby providing a per-locus association test. A set-based test was then applied using GCTA-fastBAT (a fast set-based association test that aggregates the signal from all SNPs in a region while accounting for their linkage disequilibrium) (Yang et al., 2011; Bakshi et al., 2016), defining each locus as the SNPs within ±250 kb of its position and using our own post-QC genotypes as the linkage-disequilibrium reference. Finally, the genome was partitioned into two genomic relationship matrices, one built from SNPs lying in the known-MUN regions and one from all remaining SNPs, and GREML was fitted on the corrected per-animal phenotype with the farm covariate to estimate the share of genomic variance attributable to the known-MUN regions, with a likelihood-ratio test against the single-GRM model.

A polygenic score and an allele-aligned directional sign test were additionally considered but were not feasible. We examined the full texts, Supplementary Materials, and data-availability statements of the prior MUN studies, together with PubMed, the GWAS Catalog, and Animal QTLdb: none provided genome-wide per-SNP effect sizes with effect-allele coding. The studies report only top-hit or QTL-window summaries, and the largest datasets are industry-owned and not publicly available. A weighted polygenic score and a sign-concordance test therefore could not be aligned to our genotypes.

Software

All analyses were performed using PLINK v1.9 (Purcell et al., 2007; Chang et al., 2015) for genotype QC, LD pruning, PCA, and supplementary association analyses; GCTA v1.95.0 for GRM computation and GREML (Yang et al., 2011) and for MLMA-LOCO association testing (Yang et al., 2014); and R v4.4.0 for phenotype preparation, the record-level repeatability model, data processing, and visualization. The R packages used were lme4 v1.1-37 (Bates et al., 2015) and lmerTest v3.1-3 (Kuznetsova et al., 2017) for linear mixed-model fitting, sommer v4.4.1 for the record-level genomic repeatability model, data.table v1.17.4 for data handling, ggplot2 v3.5.2 for visualization, and rtracklayer v1.66.0 (Bioconductor 3.20) (Lawrence et al., 2009) for coordinate conversion. Gene annotation was performed by querying the Ensembl REST API (Cunningham et al., 2022) over HTTP from R.

This section reports, in turn, the distribution of the corrected MUN phenotype, the population structure and GRM, the calibration of the primary association model, the repeatability and heritability bounds, and the association results under SimpleM, the Bonferroni bound, and FDR control. Particular attention is given to the collapse of the chromosome 18 signal under the corrected phenotype.

Corrected phenotype distribution and the contribution of farm

The corrected MUN phenotype (fixed-effects-only marginal adjusted mean, averaged per cow) for the 406-animal GWAS population had a mean of 9.33 mg/dL (standard deviation (SD) = 3.74; range 3.61–19.78 mg/dL; n = 406; Supplementary Figure 1). The grand mean of the underlying records was approximately 8.5 mg/dL. For comparison, the conditional phenotype had a mean of 8.92 mg/dL, with an SD of only 1.44 and a range of 6.36–15.21 mg/dL: the conditional-residual procedure had compressed the between-cow SD roughly 2.6-fold. The larger between-cow spread under the marginal definition, closer to the biological range, is the first indication that the conditional-residual procedure, not the marginal one, had distorted the between-cow variation.

Before considering genetic effects, we quantified how much of the variation in the corrected per-animal phenotype was attributable to farm and breed. A one-way model with farm alone explained 84.7% of the phenotypic variance. Mean corrected MUN was 7.92 ± 1.54 mg/dL at SCDCB Arad (n = 347) versus 17.67 ± 0.86 mg/dL at Wagner (n = 59). Breed alone explained 36.2%, but almost all of that reflected farm rather than breed: adding breed to a model already containing farm raised the explained variance by only 0.44 percentage points, and although this increment was statistically significant, it was small (Type II F(2,402) = 6.03, P = 0.0026), whereas the effect of adding farm to a model already containing breed was large (F(1,402) = 1328.4).

The comparison that separates the two is within Simmental, the only breed represented at both farms and therefore the contrast in which breed is held constant. There, farm still explained 80.8% of the variance (SCDCB Arad, 8.03 ± 1.48 mg/dL, n = 290; Wagner, 17.57 ± 0.84 mg/dL, n = 34; P = 1.7 × 10-117). Within the Wagner farm, by contrast, the two breeds present were almost indistinguishable (Holstein, 17.79 ± 0.88 mg/dL; Simmental, 17.57 ± 0.84 mg/dL). The bimodality visible in the Simmental panel of Supplementary Figure 1 is therefore a farm split rather than breed substructure.

Environmental and management factors, principally ration composition but also herd management, thus account for roughly an order of magnitude more of the variation in MUN than any genetic effect this scan was powered to detect. The null genome-wide scan is not only a consequence of limited power; it is also consistent with a trait whose expressed variation in this cohort is largely environmental.

Population structure and model calibration

Principal component analysis confirmed the expected population structure among the three breeds (Figure 1). PC1 separated Simmental from Brown Swiss and Holstein, while PC2 distinguished Holstein from Brown Swiss. Within Simmental, the only breed represented at both farms, animals from the two farms overlapped on the leading components but separated on the lower ones (standardized mean difference of 1.3 on PC2 and 1.8 on PC3), consistent with the Wagner Simmental cows forming a sub-line within the breed. Breed remains the primary axis of genetic structure, and this residual within-breed separation is what the third and fourth components absorb (Section 2.5).

The genomic relationship matrix computed from the 42,440 post-QC autosomal SNPs showed a mean diagonal (self-relatedness) of 1.021 (SD = 0.137; range 0.884–1.568) and a mean off-diagonal of –0.0025 (range –0.127–1.155). These values are consistent with what is expected for a SNP-based GRM: diagonal elements close to 1 reflect expected self-relatedness, a mean off-diagonal near 0 indicates that the cows are, on average, only weakly related, and the few larger off-diagonal values correspond to a small number of closely related, twin, or near-duplicate pairs; the single largest (relatedness approximately 1.15, between two animals with near-identical herd-book identifiers) is consistent with a monozygotic twin or duplicate-sample pair. The GRM is therefore close to an identity matrix away from those related pairs. Because separating the additive-genetic from the permanent-environment variance relies on relatedness contrasts among animals, this weak relatedness leaves the two components confounded, so they are not separately identifiable in this cohort (Section 3.3 and Section 4.4) (Speed and Balding, 2015). Three pairs had an off-diagonal coefficient above 1.0 (1.155, 1.120, and 1.007); seven exceeded 0.90, out of 82,215 pairs. All animals were retained in every analysis. The mixed model conditions on the genomic relationship matrix, so closely related or duplicated animals contribute to the polygenic term rather than inflating the association statistics, which is the situation MLMA-LOCO is designed to handle; excluding one member of such a pair would be appropriate under a fixed-effects or unrelated-samples framework but is not required here. A sensitivity analysis confirmed that the choice was immaterial: removing one member of the highest pair (N = 405) left the genomic inflation factor at 1.037 versus 1.034, and removing one member of all three pairs above 1.0 (N = 403) left it at 1.036, with no genome-wide-significant SNP under either scenario (Spearman correlations of P-values with the full-cohort scan were 0.999 and 0.989, respectively). The same three suggestive leads were recovered at N = 405; at N = 403, a fourth marker also crossed the suggestive threshold (AX-106727300, from P = 7.6 × 10 -5 to 6.2 × 10 -5), which we do not interpret.

Under the corrected primary model (Model GS: LOCO-GRM + Farm + 4 PCs), the genomic inflation factor was λ = 1.034, indicating good calibration. The within-Simmental model gave λ = 1.027, with its own quantile-quantile plot provided as Supplementary Figure 4 on the same axes as Figure 2. The quantile-quantile plot showed well-calibrated P-values across the distribution (Figure 2).

Repeatability and heritability

MUN was moderately repeatable. The ICC from the phenotype LMM was 0.243 for the full milk-recording cohort (954 cows) and 0.388 for the matched 406 genotyped cows (13,722 records); these are two estimates from two different samples, not a confidence interval. The higher value in the genotyped subset is expected because those cows skew toward the more densely recorded SCDCB Arad animals, which lowers within-cow variance. We take 0.243 as the more representative population repeatability and 0.388 as a recording-enriched upper estimate; both lie within the moderate range reported for MUN.

Narrow-sense heritability could not be separately identified. Because the genotyped cohort is mostly unrelated (near-identity GRM, Section 3.2), the additive-genetic and permanent-environment variances are confounded, and the supporting genomic models cannot resolve them: the GREML estimate for the per-animal mean exceeds the repeatability ceiling and is therefore not a heritability estimate, and the record-level animal model drives the permanent-environment variance to its boundary (Supplementary Tables S2, S3; Supplementary Figure S2). We therefore treat repeatability as the upper bound for heritability (h 2 ≤ repeatability) (Nakagawa and Schielzeth, 2010; Wilson et al., 2010), a bound that holds under the usual variance decomposition, although not in every case (Dohm, 2002), and assign h 2 no point or interval value; resolving it would require a more related, better-balanced sample, ideally with pedigree information.

Genome-wide and suggestive associations

Under the corrected primary model, no SNP reached the SimpleM genome-wide threshold (P < 3.574 × 10 -6) or the stricter parameter-free Bonferroni bound (P < 1.178 × 10 -6). Three SNPs reached the suggestive threshold (P < 7.149 × 10 -5). No SNP survived FDR control: the minimum BH-adjusted P-value across all 42,440 SNPs was 0.723. The scan is therefore null at genome-wide significance and underpowered for FDR discovery (Figure 3).

The three suggestive leads are listed in Table 1. None survived FDR control (minimum BH q = 0.723), and the three leads at the suggestive threshold match the number expected by chance almost exactly (42,440 SNPs x 7.149 × 10 -5 ≈ 3.0 expected; Section 4.2). The nearest genes, after LiftOver to ARS-UCD1.2 and Ensembl annotation within ±500 kb, were DGKI (BTA 4, within the gene), PLPP3 (BTA 3, approximately 247 kb), and OXR1 (BTA 14, within the gene). Additional non-significant signals (for example, on BTA 8, 13, 18, and 19) are tabulated in Supplementary Table 1.

The chromosome 18 signal is an artifact of phenotype definition

Under a conditional-residual phenotype, the SNP AX-124375815 (BTA 18, UMD3.1 18,238,978 bp), located near ZNF423, reached genome-wide significance (P = 4.42 × 10 -7 with the farm-only model). Under the corrected marginal phenotype, it did not survive: in the full cohort under the primary model (Model GS) it gave P = 0.065 (beta = 0.393, MAF = 0.064), and in a separate within-Simmental analysis, it gave P = 0.795 (beta = –0.071), with the effect direction reversing (Supplementary Figure 3). The two phenotype derivations, conditional-residual versus marginal adjusted mean, are the decisive difference; the covariate model (farm-only versus the GRM-plus-four-PCs Model GS) is secondary, as the factorial experiment below shows. The locus therefore provides no evidence of association once the phenotype artifact is removed; the within-Simmental result is consistent with the residual full-cohort signal being tied to the breed-farm structure rather than to a within-breed effect (the BTA 18 minor allele is concentrated in the Brown Swiss and Holstein animals, each confined to a single farm), although at N = 324 and this minor-allele frequency, this result is corroborative rather than definitive.

The phenotype-by-model factorial experiment (Section 2.6) shows that this collapse is driven by the phenotype definition, not by the model change (Table 2). Switching only the phenotype from the conditional residual to the marginal adjusted mean while holding the model at the original farm-only specification was sufficient to collapse the association at AX-124375815 (dbSNP rs110407992) from P = 4.42 × 10 -7 to P = 0.060. Adding the four PCs and moving breed into the GRM (the full Model GS) then shifted it only marginally further, to P = 0.065, whereas the same model change applied to the unchanged conditional phenotype left the SNP genome-wide significant (P = 2.68 × 10 -6). The mechanism lies in the standard error, not the effect size: the estimate was stable in sign and magnitude across all four runs (beta approximately 0.35–0.41; minor allele frequency 0.064 throughout), while the standard error increased approximately 2.6-fold, from approximately 0.082 under the conditional phenotype to approximately 0.214 under the marginal phenotype (Figure 4). This 2.6-fold increase is, to within rounding, the same factor by which the corrected phenotype’s between-cow standard deviation exceeds that of the conditional phenotype (3.74 versus 1.44 mg/dL): restoring the full between-cow variance restores a correspondingly larger standard error. The genomic inflation factor is consistent with this finding and clarifies the mechanism: under the conditional phenotype, the bulk of test statistics is deflated (λ = 0.954, falling to 0.840 with PCs), so the overshrinkage does not inflate all standard errors uniformly but instead redistributes apparent signal toward the few loci whose allele-frequency structure aligns with the record-count gradient, here the BTA 18 allele concentrated in the single-farm breeds, whereas under the marginal phenotype, the scan is well calibrated (λ = 1.055 and 1.034).

The record-count dependence that drives the artifact is absent from the corrected phenotype. We computed the Spearman correlation between the number of test-day records per cow and the absolute deviation of her per-cow phenotype from the mean. Because the two farms differ both in recording depth (median 28 records at SCDCB Arad versus 5 at Wagner) and in mean MUN (Section 3.1), deviations were taken from the farm mean so that the diagnostic would reflect record-count dependence rather than the between-farm difference. On this basis, the conditional-residual phenotype gave a Spearman correlation of –0.458 (P = 2.0 × 10 -22), whereas the corrected marginal adjusted mean gave +0.025 (P = 0.62). Computed without farm centering, the corresponding values were –0.691 and –0.554; the residual negative correlation of the marginal phenotype in that uncentered form is attributable to the farm structure, since the Wagner cows have both the fewest records and the most extreme MUN values, and it is absent within SCDCB Arad, the farm holding 347 of the 406 cows and essentially all of the variation in recording depth (conditional –0.302, P = 9.3 × 10 -9; marginal –0.078, P = 0.15).

Comparison with known QTL

None of the three suggestive SNPs overlapped previously reported QTL regions for MUN. Consistent with the null genome-wide result, targeted examination of regions implicated in prior MUN GWAS [for example, the casein cluster on BTA 6 and DGAT1 on BTA 14 (van den Berg et al., 2022)] produced only nominal P-values that did not reach the suggestive threshold.

Transferability of established MUN loci

The pre-specified transferability analysis (Section 2.10) asked, more strictly, whether loci established for MUN in larger populations are detectable in this cohort. Of the 55 reference loci, 45 were testable in our autosomal panel and 10 were not, for the coordinate and assembly reasons already described (Section 2.10). The full per-locus results are given in Supplementary Table 5.

No locus showed a transferable signal under either test. In the per-locus best-proxy test, none of the 45 loci passed the Bonferroni threshold (0.05 / 45 = 1.1 × 10 -3); the smallest P-value was 0.012 (the SLC24A4/LGMN region on BTA 21), and only three loci were nominally below 0.05 (LGMN on BTA 21, P = 0.012; the ADGRB1 region on BTA 14, P = 0.020; and a BTA 12 region, P = 0.025), close to the 2.25 hits expected by chance across 45 such tests. The set-based test agreed: no region passed Bonferroni, and the lowest set P-value was for the DGAT1 region on BTA 14 (GCTA-fastBAT set P = 0.026, 8 SNPs), the canonical milk-fat-linked MUN signal, which did not survive correction. The other canonical regions were also non-significant: the casein cluster (BTA 6, set P = 0.513), the PAEP region (BTA 11, set P = 0.173), and SCD (BTA 26, set P = 0.551). The six urea-cycle genes were likewise non-transferring (set P ranging from 0.40 for CPS1 to 0.999 for ASL); the absence of any hepatic urea-cycle signal is concordant with the renal-transport interpretation of MUN variation (Section 4.7). The headline regions are summarized in Table 3.

The two-GRM genomic partition was uninformative rather than indicative of absence. The known-MUN regions contained 0.97% of the SNPs and captured a non-significant share of genomic variance (point estimate 0.039, standard error 0.052; likelihood-ratio test P = 0.22; total genomic variance of the adjusted mean approximately 0.62, which is the genomic variance of a per-animal mean under a near-identity GRM, inflated by the averaging and by the confounding of additive-genetic with permanent-environment variance; it is not a narrow-sense heritability estimate, and it exceeds the repeatability ceiling of 0.388; see Section 3.3). The point estimate is approximately four times the SNP share of these regions, but because its standard error exceeds the estimate, the analysis cannot distinguish a fourfold enrichment from zero at N = 406. The planned polygenic score and allele-aligned sign test could not be computed because no prior MUN study provides genome-wide per-SNP effect sizes with effect alleles, only top-hit or QTL summaries (Section 2.10). Taken together, the established MUN loci show no detectable transfer to this cohort, whether tested per locus, by region, or as an aggregate genomic partition.

This study characterized the repeatability and genetic architecture of milk urea nitrogen in a multi-breed Romanian dual-purpose cattle population. The trait was moderately repeatable, and the association scan was null: no SNP reached genome-wide significance, no major-effect locus was detectable, and the genome-wide pattern was polygenic, consistent with prior MUN work. Narrow-sense heritability was bounded rather than point-identified because the cohort is mostly unrelated. The disappearance of the chromosome 18 signal under the corrected phenotype is explained by two threads discussed below: a phenotype-definition artifact and irreducible breed-farm confounding. The phenotype-definition issue is also an easily overlooked pitfall of general relevance. We discuss each in turn, then situate the findings within prior MUN work and the environmental rationale for the trait.

A phenotype-definition pitfall in two-step GWAS

A practical, easily overlooked pitfall in this study concerns how the per-animal phenotype is defined. In a two-step GWAS, repeated per-cow measurements are commonly condensed into a single adjusted phenotype through a mixed model. This two-step approach is standard in livestock GWAS: where a national genetic evaluation exists, de-regressed breeding values are used (Cai et al., 2023), and where one does not exist, records adjusted for fixed effects are used instead (Twomey et al., 2019). Because the conditional and marginal residuals are both usually reported simply as the adjusted phenotype, a brief methods description does not reveal which was used, so the distinction, and its effect on between-cow variation are easy to overlook; the fact that shrunken estimates are inappropriate as GWAS phenotypes is itself an established point (Garrick et al., 2009). Indeed, the conditional-residual mean is the natural default output of a standard mixed-model fit and was the definition initially adopted in this study before the artifact was recognized; it reflects a defensible common choice rather than an unusual mistake. The natural step is to use the residuals returned by a standard mixed-model fit. In lme4 and similar implementations, however, the default residuals() are conditional residuals: they automatically subtract the estimated, shrunken animal effect (the BLUP) in addition to the fixed effects. Because that shrinkage increases with the number of records, a cow with many records is pulled strongly toward the population mean, while a cow with few records retains most of her apparent deviation. The adjusted phenotype then encodes record count, a purely technical quantity, alongside biology. In the present data, this was not a subtle effect: the between-cow standard deviation of the per-cow conditional phenotype fell from approximately 2.76 mg/dL among cows with 2–3 records to approximately 0.17 mg/dL among cows with 11 or more records, and the between-cow SD of the per-cow phenotype was compressed roughly 2.6-fold relative to the corrected definition.

The more appropriate phenotype for this design is the fixed-effects-only marginal adjusted mean (Section 2.1), which preserves each cow’s full environment-corrected deviation regardless of record count. When no national genetic evaluation is available for the trait, as is the case for MUN here, de-regressed estimated breeding values cannot be constructed, and the marginal adjusted mean is the appropriate two-step phenotype; the standard derivation of de-regressed breeding values by Garrick et al. (2009) makes explicit why raw BLUP or shrunken estimates are inappropriate GWAS phenotypes, and worked examples in cattle GWAS confirm that adjusted records are used precisely when a national evaluation is unavailable (Twomey et al., 2019; Cai et al., 2023). Single-step approaches that blend pedigree and genomic information would in principle avoid the two-step problem, but in a cohort without a pedigree and without ungenotyped relatives, the single-step relationship matrix collapses to the genomic matrix and the method degenerates to standard genomic BLUP (Misztal et al., 2020), so the two-step adjusted-mean route is retained here. We do not claim the underlying statistical fact as new, since shrunken estimates are already known to be inappropriate GWAS phenotypes (Garrick et al., 2009); our contribution is to illustrate, quantify, and correct this easily overlooked pitfall in a concrete small two-step GWAS with uneven record counts.

This is a generalizable methodological point because its cause can be isolated and its conditions specified. The factorial experiment (Section 3.5, Table 2) shows that the phenotype switch alone drives the chromosome 18 collapse, acting through standard-error inflation rather than any change in effect size, while the covariate change on its own does not. The artifact is therefore a property of the phenotype rather than the covariate model, and it follows a predictable rule. It is most severe when three conditions coincide: record counts vary widely across animals, the trait is at least moderately repeatable (so the animal variance, and hence the BLUP shrinkage, is non-trivial), and the per-animal phenotype is taken as a mean of conditional residuals. It should be negligible under balanced recording, for weakly repeatable traits, or when de-regressed estimated breeding values from a national evaluation are available, since in each case the record-count-dependent shrinkage that drives the artifact is absent or uniform across animals. A simple safeguard is to inspect the correlation between record count and the absolute per-cow phenotype, computed within farm so that recording depth is not confounded with the herd effect. A strong negative correlation under the conditional definition versus no correlation under the marginal one, as here (–0.458 versus +0.025; Section 3.5), is a direct signature of the artifact. Stated this way, the lesson is portable to any small two-step livestock GWAS, including smallholder and local-breed datasets, where record depth is most uneven.

The marginal adjusted mean is not without cost. By averaging records without weighting, it gives a cow with two records the same influence as one with many, which adds residual noise and reduces power relative to a reliability-weighted phenotype. In the absence of a national evaluation from which such weights could be derived, we accept this trade-off as the lesser distortion, because unweighted averaging adds noise but does not generate spurious signals, whereas the conditional residual demonstrably does.

No robust locus survives under the corrected phenotype

Under the corrected phenotype, no SNP reached genome-wide significance, and none survived FDR control. The chromosome 18 signal seen under the conditional phenotype (Section 3.5) has two reinforcing causes. The dominant one is the conditional-residual artifact: it inflated between-cow contrasts in a record-count-dependent way that happened to align with the rare BTA 18 allele, and switching the phenotype alone, with the covariate model held fixed, produced essentially the entire collapse. The second, smaller cause is the irreducible breed-farm structure: the BTA 18 minor allele is concentrated in Brown Swiss and Holstein, each confined to a single farm, so no model can separate a between-farm allele-frequency difference from a breed or management effect, and only the within-Simmental analysis, in the breed spanning both farms, indicates that the residual signal is not a within-breed genetic effect (although at N = 324 and a minor-allele frequency of 0.064, this check is corroborative rather than definitive).

This pattern is consistent with what the genetics-of-discovery literature predicts. Effect-size estimates at the most significant loci of a genome scan are strongly upward-biased and can be nearly independent of the true effect, so apparent peaks in small samples frequently fail to replicate even when localization is nominally correct (Göring et al., 2001). The same winner’s-curse logic, that estimates conditioned on crossing a threshold are exaggerated relative to those obtained at follow-up, applies directly here (Kraft, 2008). Confident localization of a causal variant typically requires replication and refinement that a single small sample cannot provide (Ioannidis et al., 2009). The same reasoning governs the three suggestive leads of the corrected scan: none survived FDR control (minimum q = 0.723), their number matches the count expected by chance at the suggestive threshold, and at this power, the most significant leads have a low probability of marking a true association (Göring et al., 2001; Kraft, 2008; Ioannidis et al., 2009). We therefore treat them as sampling noise and do not assign them a biological interpretation. The appropriate conclusion is that this study identifies no MUN locus and that the apparent chromosome 18 candidate was an artifact of phenotype definition compounded by structural confounding.

Comparison with prior MUN GWAS and the power of this scan

To our knowledge, these are the first MUN descriptive statistics reported for Romanian cattle. The null genome-wide result is consistent with the broader MUN GWAS literature once power and population specificity are taken into account. Early scans reported only suggestive QTL of small effect, each explaining a few percent of variance (Bouwman et al., 2010). Even the largest meta-GWAS, with tens of thousands of cows, localized its strongest MUN signals, among 14 detected QTL, to the casein cluster, PAEP, and DGAT1 (van den Berg et al., 2022), none of which reached suggestive significance in the present sample. Studies of a size comparable to ours, with approximately a few hundred genotyped cows, have likewise reported only suggestive or window-level signals rather than single-SNP associations surviving genome-wide correction (Ariyarathne et al., 2021; Honerlagen et al., 2021), and single-step analyses in Holstein report low heritability and scattered, population-specific signals (Ma et al., 2023). We note that one such study did report an MUN-associated SNP on BTA 18 (Ma et al., 2023), but at approximately 13.6 Mb (ARS-UCD1.2), distinct from the artifactual signal discussed here, which lies at 18.16 Mb in the same ARS-UCD1.2 assembly (UMD3.1 18.24 Mb; liftover confirmed to the base pair, Section 2.10). The two positions are therefore approximately 4.5 Mb apart and are not the same locus, so our null result does not contradict that independent finding. Candidate-gene work in Brown Swiss similarly identified associations that have not since been replicated in other populations, reporting milk-urea-associated polymorphisms on BTA 7, 21, and 26 (ADRB2, PI, and SCD -1, respectively) (Cecchinato et al., 2014). Against this backdrop, a small, multi-breed, structurally confounded cohort would not be expected to detect the small-effect, population-specific QTL that characterize MUN; our result is best read as non-detection under limited power and population specificity, not as a contradiction of prior biology. This interpretation is borne out by a formal power analysis. For a single-SNP additive test (1 degree of freedom) with a non-centrality parameter of NCP = N q/(1 - q) (Sham and Purcell, 2014), where q is the proportion of phenotypic variance explained, the scan at N = 406 had 80% power only for a locus explaining approximately 6.9% of phenotypic variance at the genome-wide threshold and approximately 5.4% at the suggestive threshold (Figure 5). In per-allele terms, given the corrected phenotypic variance, this corresponds to an effect of roughly 1.4–2.3 mg/dL per copy of the minor allele, depending on its frequency. At the 2%–3% per-locus effects reported for MUN (Bouwman et al., 2010), power was only approximately 4%–14%, so only an unusually large-effect locus could have been detected; non-detection is therefore the expected outcome rather than evidence against a role for these regions.

The established MUN loci from larger studies showed no detectable transfer to the present cohort, consistent with non-detection under limited power, given that our cohort is approximately 50–85 times smaller than the discovery studies (van den Berg et al., 2022; Ma et al., 2023). The only nominal hint was the DGAT1 region on BTA 14; as the most breed-robust, large-effect MUN-associated milk-fat locus, it is the single region most likely a priori to transfer, yet even it did not survive multiple-testing correction. That the one locus most likely to appear still does not reach significance illustrates the scale of the power limitation rather than implicating DGAT1 in this cohort. The absence of a transferable signal therefore reflects the power and population specificity of this sample rather than a contradiction of the prior biology of MUN.

Moderate repeatability with heritability bounded rather than point-identified

MUN was moderately repeatable (ICC of 0.243 and 0.388; Section 3.3), in line with the within-lactation repeatability of 0.38 reported for grazing New Zealand cattle (Beatson et al., 2019), although lower values have been reported elsewhere, for example, 0.12 in Chinese Holsteins (Ma et al., 2023). Its narrow-sense heritability, however, is not separately identifiable in this cohort: with a near-identity GRM (mean off-diagonal approximately –0.0025), the additive-genetic and permanent-environment variances are confounded, so repeatability serves only as an upper bound. This ordering, with heritability well below repeatability once a permanent-environment term is included, is well recognized: dairy-cattle studies report low MUN heritabilities (approximately 0.06–0.28) and higher repeatabilities (Beatson et al., 2019; Bobbo et al., 2020; Ma et al., 2023), and the same pattern recurs across species, for example, a heritability near 0.10 and a repeatability near 0.42 for lactose in dairy sheep (Carta et al., 2023), with the same ordering for milk yield in goats (heritability approximately 0.15–0.28 versus repeatability approximately 0.68–0.73) (Rout et al., 2017). Resolving the genetic component for MUN here would require a more related and better-balanced sample, ideally with pedigree information.

Methodological recommendations for small multi-breed GWAS

In addition to the phenotype-definition lesson (Section 4.1), which is the principal methodological contribution of this study, three practical recommendations follow for small multi-breed GWAS; the first two make established good practice explicit, while the third is specific to confounded designs like ours. First, in a structured multi-breed sample, the GRM should carry breed structure, with residual axes accounted for by a few PCs and only genuinely non-genetic covariates (here, farm) fitted explicitly; this configuration was well calibrated (λ = 1.034), and a sensitivity analysis over 0, 2, 4, 6, and 10 PCs left the calibration near one (λ = 0.98–1.05) and the genome-wide null unchanged, with no SNP significant at any PC count (Supplementary Table 7) (Yang et al., 2011; Misztal et al., 2020). Second, multiple-testing thresholds derived from SimpleM depend on the declared variance cutoff, so that cutoff must be stated; we declared C = 0.995 (M eff = 13,988) and additionally report a parameter-free Bonferroni bound on all SNPs that is in fact more stringent than the SimpleM threshold, so the genome-wide null holds under every threshold considered (Gao et al., 2008, 2010). Third, and most importantly, structural breed-farm confounding cannot be removed by any model when a breed occupies a single farm; the only reliable mitigation is a within-breed analysis in the breed that spans farms (Sections 4.2 and 4.6). These considerations are directly transferable to other small, multi-breed livestock GWAS. The sweep also shows λ falling monotonically, from 1.055 at zero PCs through 1.034 at four and 1.011 at six to 0.983 at ten, while AX-124375815 (rs110407992) moves from P = 0.060 to P = 0.020 over the same range. The value at six components is nominally closer to 1 than the value at four that we adopt, which is precisely why the count was fixed a priori rather than tuned according to λ; beyond four components, the additional axes begin to absorb genuine between-animal variation, thereby deflating the test statistics without removing any residual structure. With three breeds, two principal components separate the breeds, and two more capture residual within-breed substructure, so four axes span the structure present in the sample; components beyond those no longer track stratification and begin to absorb genuine between-animal variation, thereby deflating the test statistics genome-wide. We therefore fixed the PC count at four a priori rather than selecting it based on calibration and retain that choice here, since tuning the number of components according to λ or the P-value of a particular SNP would constitute the post hoc selection the sweep is intended to guard against.

Limitations

Several limitations qualify the study. First and most importantly, the breed-farm structure is structural: Brown Swiss is confined to SCDCB Arad and Holstein to Wagner, so breed and farm effects are inseparable except through Simmental, the only breed present at both farms. Section 3.1 quantifies how much this matters: farm alone accounts for 84.7% of the variance in the corrected phenotype and 80.8% within Simmental, so the between-farm contrast dominates any genetic contrast this design can support. Second, the sample is small (N = 406) and mostly unrelated, which limits power to detect the small-effect, population-specific QTL typical of MUN and prevents separation of additive-genetic from permanent-environment variance. Third, the 65K array density may miss causal variants not in strong LD with genotyped markers; imputation to sequence would improve resolution. Fourth, no independent replication cohort was available, and a polygenic-score replication and allele-aligned directional validation against the larger discovery studies were precluded by the absence of openly released genome-wide per-SNP effect sizes with effect alleles for MUN (the prior studies report only top-hit or QTL summaries, and the largest datasets are industry-owned), despite our examination of their full texts and Supplementary Materials, PubMed, the GWAS Catalog, Animal QTLdb, and data repositories. Fifth, the marginal adjusted-mean phenotype weights all cows equally regardless of record count, which trades a modest loss of precision for the removal of the conditional-residual artifact (Section 4.1). Sixth, the two analytic choices that could most plausibly have been made differently, the four-PC count and the SimpleM variance cutoff, were tested directly and left the genome-wide null unchanged (Section 4.5). Finally, the study is limited to two farms, a state research and development station (SCDCB Arad) and a private commercial farm (Wagner), and does not represent the Romanian smallholder sector, which accounts for a large share of the national cattle population. Because smallholder milk recording can produce especially uneven record depth, the phenotype-definition pitfall documented here would be even more damaging in that setting, which strengthens rather than weakens the case for the marginal-mean approach in any future smallholder or other small, unevenly recorded analysis. A further limitation concerns the composition of the genotyped subset. The repeatability estimate is approximately 60% higher in the matched genotyped cows (ICC = 0.388, N = 406) than in the full recording cohort (ICC = 0.243, N = 954). Denser recording in the genotyped animals lowers their within-cow variance and accounts for the direction of the difference, but the reason those cows are more densely recorded is that they are not a random sample of the recording population: animals selected for genotyping tend to be those retained longer in the herd and monitored more intensively, which may correlate with productivity or deliberate selection for genomic evaluation. Both the repeatability estimate and the association results therefore describe a recording-enriched subset rather than the Romanian dual-purpose cattle population at large, which is why we consider 0.243, from the full recording cohort, as the more representative population value.

Implications for nitrogen management and for genomic selection

The corrected MUN of this cohort was low, with a record grand mean of approximately 8.5 mg/dL and a per-cow mean of 9.33 mg/dL, below the herd means of 13.6–17.3 mg/dL reported for North American dairy herds (Aguilar et al., 2012) and at or below the lower advisory band used in many milk-recording schemes. This reflects SCDCB Arad, which accounts for 347 of the 406 cows and uses a forage-based, moderate-crude-protein ration, whereas the 59 Wagner cows are at the opposite end of the range (17.67 mg/dL, Section 3.1), consistent with the brewer’s spent grains in their ration. The low cohort mean is therefore a property of SCDCB Arad rather than both farms. Because the urea/2.14 conversion is standard and the values are internally consistent, the low central value most plausibly reflects diet rather than an assay or conversion artifact. The variance decomposition points in the same direction: farm accounts for 84.7% of the variation in corrected MUN in this cohort (Section 3.1), so in this population, the lever with immediate effect on MUN is the ration rather than the genome.

The practical case for managing MUN rests on two distinct elements that should not be conflated. Because the present cohort already sits at the low end of the MUN range, the reduction argument that follows applies to MUN management generally rather than indicating an unmet need in this herd. The nutritional case is well supported and independent of any locus: lowering dietary protein surplus reduces MUN, and lower MUN is mechanistically linked to reduced urinary nitrogen excretion and, downstream, to lower ammonia emissions. This pathway is quantitatively material: dairy cows convert only approximately 25% of ingested nitrogen into milk nitrogen, so the majority is excreted (Lavery and Ferris, 2021), and because each 1-percentage-point change in dietary crude protein shifts MUN by approximately 1.1 mg/dL (Aguilar et al., 2012), modest ration adjustments move the indicator appreciably. The link from MUN to excreted nitrogen is itself quantitative: the allometric relationship UN (g/day) = 0.0259 × body weight (kg) × MUN (mg/dL), derived from Holstein and Jersey cows (Kauffman and St-Pierre, 2001), implies that a 1 mg/dL reduction in population-mean MUN corresponds to approximately 15.5 g/day per cow of avoided urinary nitrogen at a body weight of 600 kg. Physiologically, between-cow MUN differences appear to be governed substantially by renal urea handling rather than hepatic ureagenesis: cows with divergent MUN do not differ in hepatic urea-cycle gene expression, with the difference instead attributed to ruminal ammonia synthesis and renal clearance (Prahl et al., 2022), and lower renal AQP10 expression and kidney size limit blood-to-urine urea transfer in high-MUN cows (Prahl et al., 2023). This renal-transport mechanism is consistent with a moderately repeatable but polygenic trait, and a trait arising from many small renal-transport effects rather than one hepatic locus is precisely the kind of trait that requires very large samples to map, which aligns with both the null scan and the power analysis (Section 4.3). The nutritional MUN-reduction pathway is directly relevant to Romania’s commitments under the EU Nitrates Directive (91/676/EEC), which targets agricultural nitrate losses to water, and the National Emission Reduction Commitments Directive (EU) 2016/2284, which sets binding national ammonia-reduction targets, with livestock agriculture being the principal contributor (European Union, 1991, 2016).

The genomic-selection case, by contrast, is aspirational here. Because no locus replicated and narrow-sense heritability is under-identified, the present data do not support marker-assisted or locus-specific selection for MUN in this population. Concretely, separating additive-genetic from permanent-environment variance would require a substantially larger and more related reference population, on the order of several thousand cows with pedigree information or with deliberately sampled relatives; pooling records across Romanian and regional herds under a single shared marginal-mean phenotype and imputing the 65K genotypes to whole-genome sequence are the two highest-value next steps. That selecting for the trait can reduce nitrogen excretion has been demonstrated experimentally in grazing cows divergent for MUN breeding value, where lower breeding values were accompanied by measurably less urinary urea nitrogen (n = 48; Marshall et al., 2020); the present cohort is simply too small and too unrelated to support such selection. We therefore frame the environmental relevance of MUN around nutrition and management, with genomic selection positioned as a longer-term prospect contingent on replication, rather than as a conclusion of this study.

This study provides the first characterization of milk urea nitrogen repeatability and genetic architecture in Romanian dual-purpose cattle and the first MUN dataset reported for Romanian cattle. MUN was moderately repeatable (intraclass correlation of 0.243 in the full recording cohort and 0.388 in the matched genotyped subset), with a polygenic architecture in which no SNP reached genome-wide significance and none survived false discovery rate control. After correcting the phenotype from a record-count-biased conditional residual to a fixed-effects-only marginal adjusted mean, the chromosome 18 signal near ZNF423 collapsed (full-cohort P = 0.065 under the primary model; Simmental-only P = 0.795). A phenotype-by-model factorial experiment showed that the phenotype switch alone, with the covariate model held fixed, essentially produced the entire collapse, acting through an increase in the standard error rather than any change in effect size. Narrow-sense heritability could not be separately identified in this mostly unrelated cohort; we therefore report repeatability and, beyond noting that heritability is bounded above by the repeatability, we do not assign heritability a point value. Loci previously reported for MUN in larger populations were not detectable in this cohort, whether tested per locus, by region, or as an aggregate genomic partition, consistent with the limited power of a sample of this size. The value of the study is therefore concrete and threefold: it delivers the first MUN baseline and repeatability estimate for this Romanian dual-purpose cattle population, documents, isolates, and corrects a phenotype-definition pitfall of general relevance to small two-step GWAS, and provides a power benchmark for the cohort sizes that an adequately powered follow-up study would require. Replication in larger, more balanced, and more related cohorts, ideally with pedigree information and denser genotypes, is required before any locus or heritability estimate for MUN in this population can be considered established.