## Materials and Methods
### Genome-Wide Associations  
Single-trait GWAS with interaction effects between farming systems and SNP effects were conducted using the self-written R program \"GWAInter.R\". All relevant data and the program source code are available at JLUdata at the University of Giessen ([[https://jlupub.ub.uni-giessen.de/](https://jlupub.ub.uni-giessen.de/)]([https://jlupub.ub.uni-giessen.de/](https://jlupub.ub.uni-giessen.de/)), accessed on 7 January 2021). The application of GWAIner.R focuses on the detection of possible GxE for SNP markers in a one-step approach. The approach simultaneously enables the estimation of main SNP effects, as well as of SNP by environmental interaction effects, for the two different farming systems--compost and cubicle. Compared with other approaches, it is possible to estimate main SNP effects and interaction effects in one model without splitting the data set. The corresponding statistical model was:  
\[\mathrm{y}=\mathrm{Xb}+\mathrm{x}_{\mathrm{snpi}}\;\mathrm{b}_{\mathrm{snpi}} +\mathrm{x}_{\mathrm{interi}}\;\mathrm{b}_{\mathrm{interi}}+\mathrm{Zg}+ \mathrm{e} \tag{1}\]  
where y was a vector of observations for the traits LYM, MAC, PMN, sN, bN, MAJOR, MINOR, CNS, COR, AER, and AESC; b was a vector of fixed effects including parity, the general health status of the cow based on somatic cell count levels (with >200,000 cells indicating inferior cow health), farming system (compost or cubicle), farm, herd-test day for trait recording, and the person analyzing the milk samples in the laboratory; xspin was a vector of genotypes for marker i; bsmpi was a regression coefficient for the i_th_ SNP (the main SNP effect); xniteri was a vector of genotypes for cows; binteri was a regression coefficient of i_th_ SNP (the interaction effect); and g was a vector of genetic effects following N(0, G\(\mathrm{\sigma}^{2}\)g) and contained the polygenic effects, where G was a relationship matrix [30] and \(\mathrm{\sigma}^{2}\)g was the genetic variance. The G matrix was constructed considering all SNPs, but we excluded the SNP markers from the same chromosome where the candidate SNP was located, i.e., the \"leave one chromosome out strategy\"; e was the vector of random residual effects following N(0, I\(\sigma^{2}\)e), where I was an identity matrix and \(\sigma^{2}\)e was the residual variance. X and Z were incidence matrices for b and g, respectively. The \"leave one chromosome out approach\", as also implemented in most of the commercial GWAS software packages, was realized by connecting GWAlInter.R with the R package \"gaston\" [31], implying the permanent estimation of variance components per trait and chromosome. That is to say, we gave 29 G matrices (G matrix without chromosome 1, chromosome 2... chromosome 29) to \"gaston\", separately for each trait. For the estimation of main SNP and SNP interaction effects, we implemented a generalized least squares approach in GWAlInter.R. The \(p\)-value for each SNP was calculated by applying a test statistic which follows a chi-square distribution (Wald-test). The test statistic for the \(i\)_th_ main SNP effect was: \(\chi^{2}_{spi}=\frac{\widehat{\mathbf{b}_{spi}}^{2}}{\mathrm{var}\left( \widehat{\mathbf{b}_{spi}}\right)}\), with 1 degree of freedom (df). The test statistic for the \(i\)_th_ interaction effect was: \(\chi^{2}_{interi}=\frac{\widehat{\mathbf{b}_{spi}}^{2}}{\mathrm{var}\left( \widehat{\mathbf{b}_{spi}}\right)}\), with 1 df. Significantly associated SNPs were detected according to the Bonferroni corrected significance threshold \(p_{Bn\prime}=0.05/43\),095 SNP markers. Additionally, a less conservative threshold was defined as \(p_{sugg}=0.05/\)(number of independent SNPs). The number of independent SNPs was 4479 when setting the linkage disequilibrium as \(\leq\)0.15 (squared correlation between two SNPs) for consecutive genomic windows including 500 SNPs along the bovine genome.