GPS Coordinates for Modelling Correlated Herd Effects in Genomic Prediction Models Applied to Hanwoo Beef Cattle

Simple Summary It is widely known that the environment influences phenotypic expression and that its effects must be accounted for in genetic evaluation programs. The most used method to account for environmental effects is to add herd and the contemporary group to the model. Although generally informative, the herd effect treats different farms as independent units. However, if two farms are located physically close to each other, they potentially share correlated environmental factors. We introduce a method to model herd effects using physical distances between farms based on GPS coordinates as a proxy for the correlation matrix of these effects, aiming to account for similarities and differences between farms due to environmental factors. A population of beef cattle was used to evaluate the impact on the variance components and on the genomic prediction, of modelling herd effects as correlated, in comparison to assuming the farms as completely independent units. The main result was an increase in the reliabilities of the predicted genomic breeding values compared to reliabilities obtained with traditional models, a finding of practical relevance for genetic evaluation programs. Abstract It is widely known that the environment influences phenotypic expression and that its effects must be accounted for in genetic evaluation programs. The most used method to account for environmental effects is to add herd and contemporary group to the model. Although generally informative, the herd effect treats different farms as independent units. However, if two farms are located physically close to each other, they potentially share correlated environmental factors. We introduce a method to model herd effects that uses the physical distances between farms based on the Global Positioning System (GPS) coordinates as a proxy for the correlation matrix of these effects that aims to account for similarities and differences between farms due to environmental factors. A population of Hanwoo Korean cattle was used to evaluate the impact of modelling herd effects as correlated, in comparison to assuming the farms as completely independent units, on the variance components and genomic prediction. The main result was an increase in the reliabilities of the predicted genomic breeding values compared to reliabilities obtained with traditional models (across four traits evaluated, reliabilities of prediction presented increases that ranged from 0.05 ± 0.01 to 0.33 ± 0.03), suggesting that these models may overestimate heritabilities. Although little to no significant gain was obtained in phenotypic prediction, the increased reliability of the predicted genomic breeding values is of practical relevance for genetic evaluation programs.


Introduction
The concept of estimated breeding values and the development of best linear unbiased prediction (BLUP) methodology [1,2] ushered in the modern approach to genetic evaluation in animal breeding through the integration of pedigree information to calculate genetic relationships between animals [3]. This allowed the inclusion of performance measures of relatives to estimate the breeding value of an individual-the animal model. However, it was only in the late 1980s that computers became powerful enough to enable the use of animal models for genetic evaluation in practice. From there on, BLUP became widely adopted as a tool to estimate breeding values of animals for livestock production and grew into the cornerstone of successful genetic improvement programs of livestock populations [4].
The accuracy of genetic evaluations using BLUP relies on the pedigree relationships between individuals [5,6] and on the connections between contemporary groups [7]. With the development of genomic technologies and models that incorporate dense marker information [8], these genetic evaluations moved from relying solely on the pedigree relationships between individuals, to relying on the realized genetic relationships between individuals [9] estimated from their genetic markers. These genomic relationship matrices (GRM) can then be used instead of the traditional pedigree relationship matrix in a BLUP extension termed genomic best linear unbiased prediction-GBLUP [9]. In practice, industry largely uses both pedigree and genomic information simultaneously in a combined model known as single step GBLUP [10].
Genetic evaluation models largely assume that economically important traits are quantitative, that is, that the traits are influenced by a large number of genes. Many of these models also try to account for the various environmental effects to which animals are exposed. The common method to consider the effects of the different environments is to include herd and the contemporary group in the model [6,11]. Although there is a case to fit the effects of the contemporary groups (CG) as fixed [12], there are strong arguments to fit them as random effects in models for genomic prediction [13]. Henderson did conceptualize CG effects as random, however opted to treat CG as fixed under the argument that it eliminated bias from the sire model [12]. When sire models were widely applied, the choice of fixed CG effects was sensible. Nowadays, most genetic evaluation systems no longer consider the sire model, preferring the animal model instead. Nonetheless, CG effects often continue to be treated as fixed, even though in the animal model CG could be treated as random without compromising BLUP's unbiasedness property [13]. Moreover, when CG are small, or when animals have relatively low genetic connections, fitting the CG as random can improve the accuracy of prediction, compared to fitting CG as fixed [14][15][16][17][18]. Traditionally, the herd effect is modeled by assuming different farms as completely independent units. However, if two farms are physically close to each other, it is not unrealistic to expect that some of the climate, geographical factors, management strategies, and even social aspects of production would be more similar between these two farms than between two farms farther apart from each other. Treating farms as non-independent units in genomic prediction is not yet a usual practice. In a study aiming to estimate genotype by environment (GxE) interactions, management strategies, information about the CG, and climate data have been used to create correlated environmental information to be included in the model for genomic prediction [19]. A more recent study that accounted for the geographical location of farms indicated that adding a correlation structure between herd effects benefits the genetic evaluation in smallholder breeding programs [20]. Spatial modelling is widely used in geostatistics studies [21][22][23]; however, the use of spatial models is still limited in the context of animal breeding. Most spatial models are designed with the objective of making inferences for regions where no observations have been previously collected. For such scenarios, a Gaussian random field [24] is a popular choice, with the Matérn covariance function [25] being the most commonly used structure to model spatial covariance.
Here we introduce a method for genomic prediction that models the herd effects through a covariance matrix derived directly from the Euclidian distances between farms based on their Global Positioning System (GPS) coordinates, which serves as a proxy for the underlying unknown environmental correlation matrix between herds. Using the physical distances between individual farms, we defined a herd covariance matrix that enabled the model to account for their similarities or differences due to climate and other environmental factors, when these factors are not explicitly measured. The objective of this study was to evaluate how modelling the herd effects as correlated impacts the variance components and the genomic prediction, in comparison to assuming the farms as completely independent units.

Material and Methods
The Animal Care and Use Committee of the National Institute of Animal Science (NIAS), Rural Development Administration (RDA), South Korea, approved the experimental procedures (data used is part of project with approval no. 2018-293), and appropriate animal health and welfare guidelines were followed.

Genotypes, Phenotypes, and GPS Data
We used industry production trait data from a population of Korean Hanwoo cattle to compare our proposed model with models that did not include the GPS information. Genomic and phenotypic data were available for a total of 4168 commercial animals sampled throughout South Korea from a population of the Hanwoo beef cattle breed. Individuals were genotyped using the 50 k Illumina (San Diego, CA, USA) Bovine SNP50V2 BeadChip array, and after quality control filtering (SNPs with minor allele frequency <0.01 and individuals with >2% missing genotypes were excluded), 43,749 SNP genotypes were included in the analysis. These were commercial and genetically unrelated animals (no pedigree was available to these animals, and relationships in the genomic relationship matrix were low), spread across 124 finishing farms in South Korea. Table 1 summarizes the number of animals sampled across the farms. All 124 farms had their GPS coordinates described by their latitude and longitude, and Figure 1 shows their dispersion in South Korea along with a summary about the locations' altitude, average temperature, and average cumulative precipitation (historical means on measurements between the years 2000 and 2019).  All animals were born between 2014 and 2015 and slaughtered in 2017, at ages that ranged from 25 to 35 months. The data set comprised only male animals that were classified as bulls or steers. Beef production traits were recorded, and we performed our analyses on four of these traits: backfat thickness (BFT), eye muscle area (EMA), carcass weight (CWT), and marbling score (MS). BFT was measured in millimeters at the cross-sectional slice between the thoracic and the 1st lumber vertebrae perpendicular to the vertebral column, EMA was also measured from the same region in squared centimeters using a dotgrid; carcass weight was measured in kilograms at the end of refrigeration for 24h after slaughter; marbling scores were visually classified according to standard grading guidelines into nine ordinal levels [28,29]. Table 2 presents a summary of the traits and number of animals per sex and age group.  All animals were born between 2014 and 2015 and slaughtered in 2017, at ages that ranged from 25 to 35 months. The data set comprised only male animals that were classified as bulls or steers. Beef production traits were recorded, and we performed our analyses on four of these traits: backfat thickness (BFT), eye muscle area (EMA), carcass weight (CWT), and marbling score (MS). BFT was measured in millimeters at the cross-sectional slice between the thoracic and the 1st lumber vertebrae perpendicular to the vertebral column, EMA was also measured from the same region in squared centimeters using a dot-grid; carcass weight was measured in kilograms at the end of refrigeration for 24 h after slaughter; marbling scores were visually classified according to standard grading guidelines into nine ordinal levels [28,29]. Table 2 presents a summary of the traits and number of animals per sex and age group.

Prediction Models
Five prediction models were compared, with their equations fully described in Table 3. Each model consisted of the fixed effects (denoted along with their design matrix as Xb), and of a varying combination of the breeding values (g) and the herd effects (η), in our analyses represented by the farms, either assuming the farms as completely independent units (η FARM ) or correlated by their geographical distances (η GPS ). We assumed herd effects as random in all models that accounted for these effects. Models GRM, FARM, and GPS fitted each of these three random components alone, and models GRM + FARM and GRM + GPS fitted the breeding values with the herd effects simultaneously, with the two components (breeding values and herd effects) fitted as independent from one another.  Table 3. Models used for prediction.

Model Equation
A preliminary study indicated that age (fitted as a continuous covariate), sex (bull/steer), and size of the contemporary group of the herd (fitted as a continuous covariate, with size being the number of animals from each herd) were statistically significant enough to be included as fixed effects in all the models. Year, month, and season of birth were not found to be significant effects in our models. Our data consisted of animals slaughtered in the same year (2017), therefore year of slaughter was not included in the model. Although animals were slaughtered from March 2017 to November 2017, the vast majority (~65%) were slaughtered in the summer months (July-September), and neither season nor month were found to be significant effects.
Breeding values were assumed to be distributed g ∼ N 0, Gσ 2 g , where G denoted the genomic relationship matrix (GRM) as per [9]. Herd effects were considered random in all models that accounted for these effects, with Z being the design matrix indicating the herd to which the phenotypic observations y belonged to. Models FARM and GRM + FARM assumed the herd effects as independent with distribution η FARM ∼ N 0, I F σ 2 η , where I F is an identity matrix of order F, the number of herds. Models GPS and GRM + GPS assumed the herd effects as correlated with distribution η GPS ∼ N 0, Eσ 2 η , in which E is the herd covariance matrix, obtained with standardized distances between the farms on which each individual was located derived from their GPS coordinates. The following steps were used to obtain the E matrix: First, a matrix D was calculated with the raw pairwise distances between farms. This was performed using the function distm from the R [30] package geosphere [31], setting the argument fun = distGeo. The second step was to rescale the distances and ensure that all values were between zero and one by dividing them by the largest distance observed, D * = D/max D ij : i, j = 1, . . . , F . Finally, each element of E was defined as E ij = 1 − D * ij , with values 0 < E ij < 1, for every i, j = 1, . . . , F. The closer the farms, the closer to one E ij will be. Likewise, the farther the farms, the closer to zero E ij will be. Random residuals (ε) were assumed to be distributed ε ∼ N 0, I n σ 2 ε , where I n is an identity matrix of order n, the number of individuals. In addition, σ 2 g , σ 2 η , and σ 2 ε denoted the total genomic, herd and residual variances, respectively. In our study, we did not aim to infer herd effects to herds that were not in the reference population. Therefore, kernels commonly used in spatial analysis and geostatistics to infer unobserved data were not the target kernels to correlate herd effects in our study. That is because we accounted all herds as observed in the reference population, and therefore the herd effects estimated using, for example, the Matérn covariance are mathematically equivalent to the herd effects estimated using the herd covariance matrix based on the standardized distances between the farms (Supplementary Material A). Nonetheless, we ran a set of analyses considering the Matérn covariance as the kernel to correlate herd effects, to confirm empirically the theoretical equivalence of the kernels, and to verify that results would match those obtained using the herd covariance matrix based on the standardized distances between the farms (Supplementary Material B; results presented in Figures S1-S4).

Methods for Variance Components Estimation and Model Assessment
Variance components (σ 2 g , σ 2 η , and σ 2 ε ) were estimated on a training group using the restricted maximum likelihood (REML) [32], and genomic prediction was performed using GBLUP [9]. All analyses were conducted using the R [30] programming language, with functions implemented by the authors in the GenEval package, available online at https://github.com/bcuyabano/GenEval.
The performance of the different models for genomic prediction was assessed considering the reliability of the predicted genomic breeding values (PGBV) ( g test ) and of the predicted herd effects ( η test ) on a test group (here, η test refers to either η FARM,test or η GPS,test , depending on whether the herd model matrix I F or E was used). The reliabilities were calculated as R 2 PGBV = r 2 PGBV /h 2 and R 2 η = r 2 η /e 2 , a proxy for the average individual reliabilities of each breeding value, such that r PGBV = cor g test , y test − X testb and r η = cor Z η test , y test − X testb are the accuracies of prediction of the genomic breeding values and herd effects, respectively. The values of R 2 PGBV and R 2 η are constrained to the (0,1) interval, and desired to be as high as possible. We also compared the accuracy of the predicted phenotypes, calculated as r y = cor( y test , y test ). The predicted phenotypes were y test = X testb + g test for the GRM model, y test = X testb + Z η test for the FARM and the GPS models, and y test = X testb + g test + Z η test for the GRM + FARM and GRM + GPS models. Each model was replicated 100 times, with each replicate consisting of reassigning 80% of the individuals to the training group and the other 20% of the individuals to the test group. Tukey's multiple comparison test [33] was used to compare the results obtained with the different models within each trait, at a significance level of 0.05. To perform the contrast between models within each trait, the prediction accuracies and reliabilities (after taking their square root) were normalized using Fisher's z-transformation [34]. Since the heritabilities were obtained using REML, they are therefore assumed to be normally distributed, thus no transformation was done to their values. To compare the PGBVs obtained with models GRM + FARM and GRM + GPS to the PGBVs obtained with model GRM, Pearson and Spearman correlations were calculated, as well as the mean squared differences between the PGBVs, defined as MSD model 1, To rule out any potential confounding effects between G and W η FARM = ZZ , and between G and W η GPS = ZEZ , we compared their variance structures using the kappa statistics (κ η:GRM ) plotted against the eigen-values (λ η ) of W η FARM and W η GPS , as proposed by [5]. If κ η:GRM and λ η are very aligned, there is confounding between the matrices. If κ η:GRM ≈ 1, there is no confounding between the matrices [35]. Each element of κ η:GRM was computed as κ η:GRM,i = ∑ n j=1 U η,i U GRM,j 2 λ GRM,j for every i = 1, . . . , n (n being the number of observations), where U η,i and U GRM,j are the i th and j th eigen-vectors of W η (FARM or GPS) and G, respectively, and λ GRM,j is the j th eigen-value of G. Figure 2 presents the principal components plot of the genomic relationship matrix G, and this plot shows a grouping of individuals into two main genomic groups; however, these groups are not associated with any particular herd.  Figure 2 presents the principal components plot of the g , and this plot shows a grouping of individuals into two main these groups are not associated with any particular herd.  Table 4 presents the values obtained for heritability, env accuracies, and reliabilities for each of the traits and differen ferences between models within a trait were tested using Tu test, at a significance level of 0.05. Results will be presented based on the groups obtained with the multiple comparison te that values are equal, lower, or higher, we mean statistically significance level of 0.05.  Table 4 presents the values obtained for heritability, environmentability, prediction accuracies, and reliabilities for each of the traits and different models. The statistical differences between models within a trait were tested using Tukey's multiple comparison test, at a significance level of 0.05. Results will be presented, compared, and discussed based on the groups obtained with the multiple comparison test. Therefore, when we state that values are equal, lower, or higher, we mean statistically equal, lower, or higher at a significance level of 0.05.

Results
Heritability estimates obtained for EMA, CWT, and MS obtained with the GRM + GPS model were lower than those obtained with the GRM and GRM + FARM models. Heritability estimates obtained for BFT were equal in all three models that included breeding values. The environmentability estimates obtained with the GPS and GRM + GPS models were higher than those obtained with the FARM and GRM + FARM models for all traits. Environmentability estimates obtained with the FARM and GRM + FARM models were equal for BFT and EMA, while environmentabilities estimated with the FARM model were higher than those estimated with the GRM + FARM model for CWT and MS. Except for CWT, which had a higher environmentability estimate with the GPS model than with the GRM + GPS model, environmentabilities estimated with the GPS and GRM + GPS models were equal.
Prediction accuracies of the PGBV obtained with all the models were equal for all traits. Table 5 presents the Pearson and Spearman correlations and the mean squared differences used to compare the PGBVs obtained with models GRM + FARM and GRM + GPS to the PGBVs obtained with model GRM. Although both the Pearson and the Spearman correlations were consistently higher for model GRM + GPS, their values were between 0.97 and 1 for both the GRM + FARM and GRM + GPS models, indicating that including the herd effects to the genomic prediction did not produce significant changes in the ranking of the PGBVs, when compared to the GRM model that did not account for the herd effects. Comparing the PGBVs by the mean squared differences, we observed that the GRM + GPS model obtained PGBVs more similar to those obtained with the GRM model than did the GRM + FARM model. This indicates that in our data, the inclusion of correlated herd effects interfered less with the PGBVs than did the inclusion of independent herd effects. Table 4. Heritability and environmentability estimates (ĥ 2 = σ 2 g /σ 2 y andê 2 = σ 2 η /σ 2 y ), prediction accuracy of genomic breeding values and herd effects (r PGBV = cor( g test , y test ) and r η = cor( η test , y test )), reliability of predicted genomic breeding values and herd effects (R 2 PGBV = r 2 PGBV /ĥ 2 and R 2 η = r 2 η /ê 2 ), for all four traits and all the five models evaluated. Values in this table represent the mean over 100 cross-validation replicates of each model.  Table 5. Comparisons of the predicted genomic breeding values (PGBV) between the model that modelled the PGBV alone (GRM) and the models that accounted for the herd effects, either as independent (GRM + FARM), or as correlated (GRM + GPS), for all four traits and all the five models evaluated. The comparisons were performed based on the Pearson (cor P ) and Spearman (cor S ) correlations between the PGBVs, and on the mean squared differences between the PGBVs Prediction accuracies of the herd effects obtained with the GPS model were higher than those obtained with the FARM model for all traits. Except for CWT that had equal herd effect prediction accuracy with the GRM + GPS and GRM + FARM models, prediction accuracies of the herd effects obtained with the GRM + GPS model were higher than those obtained with the GRM + FARM model. Prediction accuracies of the herd effects obtained with the FARM and GRM + FARM models were equal for BFT and CWT, while prediction accuracies of herd effects obtained with the FARM model were lower than those obtained with the GRM + FARM model for EMA and MS. Prediction accuracies of the herd effects obtained with the GPS and GRM + GPS models were equal for all traits.
Reliabilities of PGBV and environmental effects are shown in Figure 3. All PGBV reliabilities from the GRM + GPS model were higher than those from the GRM and GRM + FARM models, and the difference is particularly large for EMA, CWT, and MS. The PGBV reliabilities between the GRM and the GRM + FARM models are equal except for MS. Reliabilities of the predicted herd effects were higher in the FARM and GRM + FARM models than in the GPS and GRM + GPS models, except for MS that was the opposite.  (blue letters) and of 2 (red letters) obtained with the different models, within each trait. Figure 4 shows the accuracy of phenotypic prediction. The inclusion of a herd effect component in the genomic model provided a significant phenotypic prediction gain for BFT, EMA, and MS. For BFT, this gain was observed only with the GRM + GPS model. Additionally, BFT was the only trait in which we observed a significant phenotypic prediction gain by using the GRM + GPS model instead of the GRM or GRM + FARM models; however, it was a small gain: the mean observed with the GRM, GRM + FARM, and GRM + GPS models were respectively 0.372, 0.363, and 0.383. From the boxplots of the accuracy of phenotypic prediction we could observe that, for BFT and EMA, the prediction accuracies obtained over the 100 replicates with the model GRM + GPS varied less than the prediction accuracies obtained with the model GRM + FARM.   Figure 4 shows the accuracy of phenotypic prediction. The inclusion of a herd effect component in the genomic model provided a significant phenotypic prediction gain for BFT, EMA, and MS. For BFT, this gain was observed only with the GRM + GPS model. Additionally, BFT was the only trait in which we observed a significant phenotypic prediction gain by using the GRM + GPS model instead of the GRM or GRM + FARM models; however, it was a small gain: the mean observed r y with the GRM, GRM + FARM, and GRM + GPS models were respectively 0.372, 0.363, and 0.383. From the boxplots of the accuracy of phenotypic prediction we could observe that, for BFT and EMA, the prediction accuracies obtained over the 100 replicates with the model GRM + GPS varied less than the prediction accuracies obtained with the model GRM + FARM.   ′, and between and = ′, using the kappa statistics [35]. We observed that : and are not aligned at all, and that : ≈ 1, indicating that there is no confounding between the two relationship matrices.  . Accuracy of phenotypic prediction (r y = cor( y test , y test )), for all four traits and all the five models considered. Boxplots based on 100 cross validation replicates of each model. The predicted phenotypes were y test = X testb + g test for the GRM model, y test = X testb + η test for the FARM and the GPS models, and y test = X testb + g test + η test for the GRM + FARM and GRM + GPS models. Letters A, B, C, D indicate the groups obtained with Tukey's multiple comparison test at a significance level of 0.05, comparing the values of r y obtained with the different models, within each trait. Figure 5 presents the comparison of the variance structures between G and W η FARM = ZZ , and between G and W η GPS = ZEZ , using the kappa statistics [35]. We observed that κ η:GRM and λ η are not aligned at all, and that κ η:GRM ≈ 1, indicating that there is no confounding between the two relationship matrices.

Discussion
This work evaluated how modelling the herd effects as correlated impacts the variance components and the genomic prediction, in comparison to assuming the herds as completely independent units. To estimate and predict the correlated herd effects, we presented a method that used the physical distances (based on GPS coordinates) to create a herd covariance matrix, which was used as a proxy to model the unknown underlying correlation matrix of these effects. Conceptually, this herd covariance matrix should enable the model to partially account for similarities or differences between herds due to climate and other geographical factors that were not explicitly measured and included in the model. We tested this assumption on four production traits in a commercial population of Hanwoo beef cattle sampled from across South Korea by modelling herd effects using GPS coordinates. The main result was an increase of the reliabilities of PGBV for all traits in comparison to the reliabilities obtained with a model that assumed the herd effects to be independent or with a model that did not account for the herd effects at all.

Discussion
This work evaluated how modelling the herd effects as correlated impacts the variance components and the genomic prediction, in comparison to assuming the herds as completely independent units. To estimate and predict the correlated herd effects, we presented a method that used the physical distances (based on GPS coordinates) to create a herd covariance matrix, which was used as a proxy to model the unknown underlying correlation matrix of these effects. Conceptually, this herd covariance matrix should enable the model to partially account for similarities or differences between herds due to climate and other geographical factors that were not explicitly measured and included in the model. We tested this assumption on four production traits in a commercial population of Hanwoo beef cattle sampled from across South Korea by modelling herd effects using GPS coordinates. The main result was an increase of the reliabilities of PGBV for all traits in comparison to the reliabilities obtained with a model that assumed the herd effects to be independent or with a model that did not account for the herd effects at all.
The source of this increase in the reliabilities of PGBV is of particular interest. It is basically because the GRM + GPS model reallocated to the herd effect a proportion of the phenotypic variance that the other models assigned to the breeding values. In all four traits, the heritability estimates of the GRM + GPS model were significantly lower than the estimates of the GRM and GRM + FARM models. The decrease in the heritability estimates led to an increase in the PGBV reliabilities since R 2 PGBV = r 2 PGBV /h 2 . At this point, our results suggest that the heritabilities of the traits estimated with the GRM + GPS model are more accurate than the higher values estimated by both the GRM and GRM + FARM models. If the GRM + GPS model was underestimating the heritabilities in comparison to GRM and GRM + FARM models, we would expect changes in the PGBV variances or a reduction in the prediction accuracies. However, across all four traits, the PGBV variances and the prediction accuracies were statistically equal for the GRM, GRM + FARM, and GRM + GPS models. The reduction in heritability estimates with the GRM + GPS model did not affect the PGBV accuracies, which suggests that the GRM and GRM + FARM models probably overestimated the variance of the breeding values. If the PGBV reliabilities are in fact higher than what traditional approaches indicate, there are relevant real-world implications for breeding programs, since accurate estimates are of practical importance and are broadly used to design breeding programs, predict response to selection, and determine selection strategies [36][37][38].
The covariance matrix built from GPS coordinates to model herd effects was expected to account for some of the similarities or differences between herds due to climate and other environmental factors. This was confirmed by the much higher environmentability estimates from the GPS and GRM + GPS models compared to the environmentabilities estimated with the FARM and GRM + FARM models. Even though the GRM + GPS reallocated genetic variance into environmental variance, the total variance explained was much higher than with the other models (Table 4), which would intuitively suggest a better predictive value as well; however, as predictors, the GPS information was not sufficiently informative to improve phenotypic prediction; little to no significant gain in phenotypic prediction was achieved with the GRM + GPS model. In part, this could be due to the data used in this study being relatively small and to South Korea being also a small country with a relatively uniform climate and production environment. It is possible that prediction accuracies would increase with a larger dataset from a country with more diverse conditions. It is, however, more likely that to substantially improve phenotypic prediction, direct and detailed environmental/management data modelling will be necessary.
Although the results from our study suggest that herd effects should be modelled as correlated, we understand that the herd covariance matrix based on GPS coordinates is just an approximation that tries to simplistically capture a very wide range of environmental and farm conditions. While modelling herds as correlated effects clearly accounted for significant amounts of phenotypic variation, and it did not negatively impact genomic prediction accuracies, it was still not descriptive enough to increase prediction accuracies in this study. Different methods to correlate herd effects based on GPS coordinates should be further explored in follow-up studies. Moreover, recent studies using different approaches for phenotypic prediction have incorporated climate information measured in more diverse environments than what we observed in South Korea, and have reported positive gains in prediction accuracy in livestock [39] and in crops [40,41]. We believe that the inclusion of environmental parameters to quantify their influence on phenotypic variance and prediction accuracies and to understand how they correlate with the GPS herd covariance matrix can further improve models for genomic and phenotypic prediction.
Finally, although the models with a herd effect that we explored in this study provided little improvement in phenotypic prediction, the environmentability estimates and the reliabilities of the PGBV achieved with the models that used the GPS coordinates are relevant for genetic evaluations. Modelling herd effects as correlated instead of independent indicated that results from traditional models might overinflate the heritability estimates.
These overestimates can lead to wrong reliabilities of the PGBV and, consequently, to misleading values of expected genetic gains. In this study, we modelled the correlations between herds through their physical distances using GPS coordinates, but other methods and direct environmental factors such as differences in altitude, temperature, precipitation, or even management or nutrition could be considered instead.

Conclusions
Environmental effects have a large influence on phenotypic expression and must be accounted for in genetic evaluations. This is usually achieved through the inclusion of a herd effect term in the model. Accurately modelling these herd effects is fundamental to adequately assign phenotypic variance to breeding values and random herd effects, in order to maximize prediction accuracy and reliability of the PGBV. Our results showed that modelling herd effects by considering the herds as being correlated with each other was more successful in capturing environmental variance (environmentability) and the reliability of the PGBV was greater than when herds were considered as completely independent units. The increased reliabilities of the PGBV are particularly relevant for decision making in breeding programs.

Data Availability Statement:
The data that support the findings of this study were available from the Rural Development Administration, Republic of Korea. Restrictions apply to the availability of these data, which were used under license for this study.