Multi-trait meta-analyses identify potential candidate genes for growth-related traits in Holstein heifers

Understanding the underlying pleiotropic relationships among growth and body size traits is important for refining breeding strategies in dairy cattle for optimal body size and growth rate. Therefore, we performed single-trait genome-wide association studies (GWAS) for monthly-recorded body weight (BW), hip height (HH), body length (BL) and chest girth (CG) from birth to 12 mo of age in Holstein animals, followed by stepwise multiple regression of independent or lowly-linked markers from GWAS loci using conditional and joint association analyses (COJO). Subse-quently, we conducted a multi-trait meta-analysis to detect pleiotropic markers. Based on the single-trait GWAS, we identified 170 significant single nucleotide polymorphisms (SNPs), in which 59 of them remained significant after the COJO analyses. The most significant SNP, located at BTA7: 3 ,676 ,741, explained 2.93% of the total phenotypic variance for BW6 (body weight at 6 mo of age). We identified 17 SNPs with potential pleiotropic effects based on the multi-trait meta-analyses, which resulted in 3 additional SNPs in comparison to those detected based on the single-trait GWAS. The identified quantitative trait loci (QTL) regions overlap with genes known to influence human growth-related traits. According to positional and functional analyses, we proposed HMGA2 , HNF4G , MED13L , BHLHE40 , FRZB , DMP1 , TRIB3 , and GATAD2A as important candidate genes influencing the studied traits. The combination of single-trait GWAS and meta-analyses of GWAS results improved the efficiency of detecting associated SNPs, and provided new insights into the genetic mechanisms of growth and development in Holstein cattle.


INTRODUCTION
Growth-related traits are directly linked to economic returns in beef cattle, and also more recently in dairy cattle herds employing "beef-on-dairy" crossbreeding schemes (Berry, 2021).Growth traits are also important in dairy cattle breeding schemes (Miglior et al., 2017, Yin andKonig, 2020) due to associations between body weight (BW) with feed efficiency (e.g., Pryce et al., 2015), methane emission (e.g., de Haas et al., 2017), and fertility (e.g., Handcock et al., 2020) traits.Relative economic weights for BW have increased in updated selection indexes in the Australian dairy industry (Byrne et al., 2016).Another important growth trait, height (stature) is also linked to productivity (Van De Stroet et al., 2016).Hip height (HH) is highly heritable with heritability estimates ranging from ~0.27 to 0.69 in Holstein cattle populations (Koenen and Groen, 1998;Wall et al., 2005;Haile-Mariam et al., 2013;Zhang et al., 2017), thus, being an important phenotype for studying the genetic background of growth and development.Moreover, previous research has indicated that many loci regulating growth are shared across mammalian species.For instance, a set of common genes related to stature was identified between human stature and cattle height (Pryce et al., 2011;Bouwman et al., 2018).Sahana et al. (2015) also reported 2 loci associated with human adult stature and calf size.A meta-analysis of daily weight gain in pigs also identified 23 reported human growth-related genes (Berndt et al., 2013).Thus, deciphering the genetic architecture of growth in Holstein not only benefits dairy cattle breeding but also provides insights into the mechanisms controlling stature and BW in other species.
Several genes related to growth and development have been reported, including ABCA12, FLRT2, LHX4, MAP3K5, NRAC, NTNG1, PIGN, and ZNF75A, which were identified based on a genome-wide association study (GWAS) for predicted BW in Holstein cattle (Cole et al., 2014).In addition, Zhang et al. (2017) identified multiple genetic loci for HH and chest girth (CG) at 6, 12, 18, and 24 mo of age in Holstein cattle.Yin and Konig (2019) detected candidate genes for direct genetic and maternal genetic effects influencing BW at birth, 2 to 3 mo of age, and 13 to 14 mo of age, although few significant single nucleotide polymorphisms (SNPs) were found associated with the maternal genetic effects.However, the aforementioned publications focused on few time point measurements or growth stages, such as birth weight, weaning weight, and yearling weight/size.Based on random regression models, a previous study reported heritability estimates for BW from birth to 2 years old ranging from 0.26 to 0.83 in Holstein cattle (Yin and Konig, 2018), suggesting that the genes regulating growth are likely not the same throughout the life of the animals.Yin and Konig (2019) reported that only 2 SNPs located on BTA5 significantly contributed to the 3 BW traits measured at different growth stages in 15,921 Holstein cattle.Hence, the records of growth traits with shorter time interval but longer growing period are more conducive to dissect the genetic architecture of growth.Moreover, automatic weighing systems installed in modern dairy cattle farming systems is a promising opportunity to further investigate growth traits.
Despite the success of GWAS in identifying many variants associated with complex traits, it is underpowered to find polymorphisms affecting multiple traits.Pleiotropy is one of the causes of genetic correlations between traits (Falconer and Mackay, 1996), so considering several correlated phenotypes together in the multi-trait meta-analysis can boost statistical power of the analyses (Evangelou and Ioannidis, 2013).In addition to weight and height measurements, body length (BL) and CG can also be recorded for better understanding growth mechanisms in cattle.Consequently, the main objectives of this study were to identify SNPs and genes associated with BW, HH, BL, and CG traits measured from birth to yearling ages in Holstein animals based on single-trait GWAS and multi-trait metaanalyses (GWAS summary statistics), and examine whether the candidate genes are known to be related to growth traits in humans.

Ethics Statement
The procedures for body measurements and blood sample collection were carried out in strict accordance with the protocol approved by the Animal Welfare Committee of the China Agricultural University (Protocol Number: DK996).

Animals and Phenotypes
Measurements from Holstein heifers from 2 large herds located in the Chinese provinces of Henan and Shandong were included in this study.Four growth traits were studied, including BW, HH, BL, and CG.Birth weight (BW0) records of 7,122 Holstein calves born between 2014 and 2020 were obtained from calving information extracted from the AfiFarm software (Afimilk, Kibbutz Afikim Israel).BW of 2,173 animals was measured repeatedly once a month during the period from 2 to 12 mo of age using the intelligent weighing scale (ID5000, Tru-Test Corp. Ltd., Auckland, New Zealand) and HH, BL, and CG of 2,123 animals were collected around 3, 6, 9, and 12 mo of age based on manual measurements.Based on the approach suggested by Motulsky and Brown (2006), phenotypic outliers were detected from nonlinear curve fitted with a false discovery rate (FDR) of 1%.The pedigree file consisted of 14,643 females and 1,924 males born from 1948 to 2020.Only animals with pedigree and/or genomic information were included in the analyses.The summary statistics of the data sets after data editing are presented in Table 1.

Genotypes and Quality Control
Among the cattle with phenotypic data, 1,002 animals were genotyped with the Illumina 150K Bovine Bead chip (Illumina, Inc., San Diego, CA, USA), which contains 139,376 SNPs.Imputation of missing SNPs was performed using the BEAGLE v5.1 software (Browning et al., 2018).The physical positions of the SNPs were updated from the UMD 3.1 assembly to the ARS-UCD 1.2 assembly of the bovine reference genome using the UCSC LiftOver tool (Kent et al., 2002).Genomic quality control was performed using the PLINK v1.90 software (Purcell and Chang, 2019).SNPs with call rate greater than 90%, minor allele frequency (MAF) higher than 0.05, no extreme departure from the Hardy-Weinberg equilibrium (P > 10 −6 ), known chromosome and genome position, and located in the autosomal chromosomes were retained for further analyses.Animals with a call rate less than 90% were removed.After the quality control, 971 animals and 114,658 SNPs remained for further analyses.

Pedigree-based (Co)Variance Components and Breeding Value Estimation
To obtain the pseudo-phenotypes for the GWAS, variance components and estimated breeding values for each trait were calculated based on the Restricted Maximum Likelihood (REML) method, as implemented in the DMU v6 software package (Madsen and Jensen, 2013).The statistical model for BW0 can be described as where y is a vector of phenotypic observations for BW0; b is a vector of fixed effects including herd, birth yearseason, and dam parity; a is a vector of direct additive genetic effects with ) where A is the ped- igree-based relationship matrix and σ a 2 is the direct additive genetic variance; m g is a vector of random maternal genetic effects with m where I is an identity matrix and σ me 2 is the maternal environmental variance; e is a vector of random residual effects with e I ~, , N e 0 2 σ ( ) where σ e 2 is the random residual variance; X, Z, W, and S are the incidence matrices for b, a, m g , and m e , respectively.The animal model fitted for the other BW traits can be described as: where y is a vector of phenotypic observations for BW traits from 2 to 12 mo of age (BW2, BW3, BW4, BW5, BW6, BW7, BW8, BW9, BW10, BW11, and BW12); b is a vector of fixed effects including herd, birth yearseason, dam parity, and age in days at the time of measurement as a linear covariate; a is a vector of additive genetic effects with ) where A is the pedigree-based relationship matrix and σ a 2 is the direct additive genetic variance; e is a vector of random residual effects with e I ~, , N e 0 2 σ ( ) where σ e 2 is the random residual variance; X and Z are the incidence matrices for b and a, respectively.We fitted bivariate models considering body size and weight traits measured at similar ages as it improved the reliability of the estimated breeding values for HH, BL, and CG (data not shown).The fixed effects were the same as defined in model [2].The estimated breeding values and residuals from models [1] and [2] were used to calculate the pseudo-phenotypes.For a specific genotyped animal, the pseudo-phenotype was the sum of the estimated direct breeding value, the maternal breeding value, the maternal environmental effect, and the residual from model [1] for BW0, and was the sum of estimated breeding value and the residual from model [2] for the other traits.

Single-trait GWAS
Individual records for the 24 growth traits were used to perform GWAS using the fixed and random model Circulating Probability Unification (FarmCPU; Liu et al., 2016) with the rMVP software (Yin et al., 2021), in which each SNP was tested for an association with the trait.The FarmCPU model uses the fixed and random effects iteratively to simultaneously control for false positives without compromising false negatives.The fixed effect model (FEM) can be written as: where y i is the pseudo-phenotype of the ith individual for each of the 24 growth traits; M ij are the genotypes of j pseudo quantitative trait nucleotides (QTNs); S ik is the genotype for the kth genetic marker of the ith individual; b j and d k are the corresponding effects; e i is the residual effect with distribution, e I ~, , where σ e 2 represents the residual variance.The random effect model (REM) was used to select the most appropriate pseudo-QTNs.The model can be written as follows: [4] where y i and e i are the same as in the FEM [3]; g i is the genetic effects of the ith individual and in which σ g 2 is the additive genetic variance and K is kinship derived from the pseudo-QTNs.The FEM and REM are iterated until no new pseudo QTNs are added, or the specified maximum number of iterations is reached.The obtained p-values were corrected for multiple testing by calculating the FDR at 5% genome-wise significance level.

Conditional and Joint Association (COJO) Analyses
To refine the single-trait GWAS within each QTL region and account for linkage disequilibrium (LD), we performed approximate conditional analyses (COJO; Yang et al., 2012), as implemented in the GCTA software (Yang et al., 2011), to the single-trait GWAS for each chromosome through a stepwise model selection procedure.As described by Yang et al. (2012), this method combines GWAS summary-level statistics and estimates the joint effects of multiple SNPs in a genomic region, which may be more powerful and/or interpretable than the standard/marginal single SNP-based analyses.The top markers from the association with P < 5 × 10 −6 were selected.The analyses were restricted to variants within 10 Mb from the conditional markers.The whole discovery sample of the GWAS analyses was used as the LD reference sample.
The meta-analysis used the SNP effect (β) and the corresponding standard error (SE) from single-trait GWAS to perform a chi-squared test statistic based on the following equation (Bolormaa et al., 2014): where n is the number of traits included in the respective meta-analysis; t i is a n × 1 matrix of the signed tvalues (β/SE) at the ith SNP; t i ' is the transpose vector of t i ; V −1 is the inverse of the n × n correlation matrix of the correlation between the signed t-values.The obtained P values were corrected for multiple testing.An FDR of 10% was used to determine the genome-wide significant threshold in the meta-analyses.

Phenotypic Variance Explained by Significant SNPs
The phenotypic variance explained (PVE) by the significant SNPs associated with growth traits was es-Ma et al.: META-ANALYSES FOR DAIRY CATTLE GROWTH timated using the GCTA software (Yang et al., 2011) based on the fastGWA model (Jiang et al., 2019): where y is a vector of pseudo-phenotypes; x snp is a vector of mean-centered genotype variables of a variant of interest with its effect β snp ; X c is the incidence matrix of the first 3 principal components with their corresponding coefficients β c ; g is a vector of the aggregate effects of all SNPs with ) where G is the genomic relationship matrix with all of the small offdiagonal elements (those < 0.05) set to 0 (Zaitlen et al., 2013); e is the vector of residuals with e The PVE was calculated as (Shim et al., 2015): where N is sample size and SE snp

Functional Annotation of Candidate Genes
Positional candidate genes located at up to 200 kb upstream and downstream of the significant SNPs were identified (Sanchez et al., 2017, Chen et al., 2021) based on the ARS-UCD 1.2 assembly of the latest cattle reference genome using the National Center for Biotechnology Information (NCBI; www .ncbi.nlm.nih.gov/gene/ ) and the Human Gene Database (www .genecards.org/).Furthermore, we searched for common genetic factors that underly human growthrelated traits.The catalog of human body mass index (EFO_0004340), obesity (EFO_0001073), body fat percentage (EFO_0007800), appendicular lean mass (EFO_0004980), body weight (EFO_0004338), body height (EFO_0004339), waist-hip ratio (EFO_0004343), hip circumference (EFO_0005093) and waist circumference (EFO_0004342) were downloaded from the NH-GRI-EBI Catalog of human genome-wide association studies (www .ebi.ac.uk/gwas/ home).

Genetic Parameters
Estimates of variance components and heritabilities for BW, HH, BL, and CG are shown in Table 2. Moderate direct heritabilities in the range from 0.31 to 0.54 were identified for BW at different age points.The maternal heritability for BW0 was 0.08, and the genetic correlations between direct and maternal genetic effects was −0.66.The heritabilities for body size traits were moderate to high ranging from 0.35 to 0.49 for HH, 0.49 to 0.66 for BL, and 0.19 to 0.76 for CG.Based on the trait categories of subsequent meta-analysis, genetic correlations for BW, HH, BL, and CG by age in months are presented in Figure 1.The genetic correlations between BW traits were close to 0.90 for adjacent months, but decreased with increasing time spans.Generally, genetic correlations of HH from different age points, ranging from 0.53 to 0.91, were larger than the estimates for BL and CG from different age points.Large and positive correlations between the 4 growth traits at 6 and 9 mo of age were observed, and a genetic correlation of 0.87 ± 0.08 was greatest between CG6 and BW67 (BW at either 6 or 7 mo of age, depending on which measurement was closer to the measurement time of CG6).However, the genetic correlation between BL3 and CG3 was negative (−0.54 ± 0.13).
The number of significant associations detected for each trait ranged from zero to 15, with the largest number for BW8, and no significant SNPs were identified for BL9 based on single-trait GWAS.For the BW traits, the most significant SNP (ARS-BFGL- ), located at BTA7: 3 ,676 ,741, explained 2.93% of the total phenotypic variance for BW6.The SNP Hapmap50112-BTA-80530, located at BTA7: 106 ,968 ,339, explained 3.52% of the total phenotypic variation for BW3, which was also the highest PVE by one significant SNP across the 24 studied traits.For body size traits, the SNP Bovine-HD0100024874 located at BTA1: 86 ,836 ,561 with the P-value of 1.85 × 10 −11 was the most significant marker and ARS-BFGL-NGS-108785 located at BTA11: 68 ,884 ,476 with the PVE of 2.85% was the top SNP for HH; the SNP related to CG with the highest significance (P = 2.05 × 10 −13 ) and the largest PVE (2.52%) was the BovineHD1600002786 located at BTA16: 9 ,745 ,951.The significance of the most significant SNP (Bovine-HD1000022893; BTA10: 79 ,999 ,486; P = 4.41 × 10 −8 ) for BL was weaker than that of the most significant SNP for the other traits, and the PVE by this SNP was only 1.16%.The SNP with the maximum PVE (1.80%) for BL was the BovineHD1400006012, which is located at BTA14: 19 ,363 ,068.

Conditional and Joint Association Analyses
The conditional analyses were performed for the 24 growth traits separately.In this study, a QTL region was defined as a chromosomal region at which adjacent pairs of associated SNPs are less than 10 Mb distant.We found 120 jointly associated SNPs (P < 5 × 10 −6 ) located in 42 independent QTL regions, including 59 leading and 62 additional SNPs (Supplemental Table S1; https: / / doi .org/ 10 .6084/m9 .figshare.22217719;Ma et al., 2023).In other words, 111 SNPs were not significant anymore after the conditional analyses.The reason why the sum of "59 leading and 62 additional SNPs" does not equal to 120 SNPs is because the SNP BovineHD1400006012, located at BTA14: 19 ,363 ,068, was significant only for BL12 in the single-trait GWAS but reached the significance level for both BL12 and BW0 when fitted jointly with the leading SNP(s) at its QTL region.We did not find any QTL region showing associated SNPs for BW7, BW9, and BL9, and no additional signals were detected for BW4, BW10, BL3, BL6, CG3, and CG12.Compared with the single-trait GWAS, the significance of 31 leading SNPs increased in the COJO analysis.There were 12 loci harboring more than 2 associated SNPs, with a maximum number of 17.The top 3 loci in terms of significance were centered on the variants BovineHD1600002786 located at BTA16: 9 ,745 ,951, BovineHD0600029059 located at BTA6: 102 ,489 ,992, and BovineHD0100013894 located at BTA1: 48 ,811 ,095.The leading SNP BovineHD1700017362, located at BTA17: 58 ,611 ,759, remained significant for both BW10 and CG9.

Candidate Genes and Functional Analyses
Taken together, based on the different analyses, 233 SNPs were found to be significantly associated with growth traits in Holstein cattle.Twelve SNPs were significant regardless of the analyses.Due to the reasonably small sample size and potential false-positive results, we focused on the 64 primary SNPs, the combination of the 59 leading SNPs identified by the COJO approach and the 17 pleiotropic SNPs identified by meta-analysis, to map positional candidate genes.Out of these, 12 primary SNPs overlapped between the 2 sets.A total of 191 protein-coding genes were detected based on the annotations of the ARS-UCD 1.2 assembly of the bovine reference genome.Furthermore, we searched for homologous genes of human and  e) growth traits at 3 to 4 mo of age, f) growth traits at 6 to 7 mo of age, g) growth traits at 9 to 10 mo of age, h) growth traits at 12 mo of age.cattle for growth-related trait from the GWAS catalog (www .ebi.ac.uk/gwas/ home).The corresponding traits of reported genes are listed in Supplemental Table S2 (https: / / doi .org/ 10 .6084/m9 .figshare.22217749;Ma et al., 2023).Out of the 191 genes, 84 overlapped with genes identified to influence growth and development traits in humans, with 29,10,7,14,20,56,16,4, and 2 associated with body mass index, obesity, body fat percentage, appendicular lean mass, body weight, body height, waist-hip ratio, hip circumference, and waist circumference, respectively.HMGA2, identified for BW12, has been previously associated with 6 human growth-related traits (www .ebi.ac.uk/gwas/ home).MED13L has been associated with body mass index, appendicular lean mass, body weight, and body height (www .ebi.ac.uk/gwas/ home), which were also found to be associated with BW and CG in Holstein cattle (this study).

DISCUSSION
The heritability estimates for BW, HH, BL, and CG from birth to 12 mo of age obtained in this study were of moderate to high magnitude.For BW0, a negative genetic correlation between direct and maternal effects was observed.Similar results were reported in previous Holstein cattle studies.Yin and Konig (2018) estimated a direct-maternal genetic correlation for calf birth weight of −0.39.In this regard, Bauman and Currie (1980) hypothesized that nutrient partitioning in mammals can lead to genetic antagonism between direct and maternal effects, as indicated by the negative genetic correlation between weight and milk yield observed by Berry et al. (2003).Besides, the number of progenies per dam and the reduced amount of information from recorded dams could have a large influence on the accuracy of genetic correlation estimates between direct and maternal genetic effects (Maniatis and Pollott, 2003).The heritability estimates for BW3 were lower compared with those of other BW traits.This finding can be attributed to specific factors within our experimental population.For instance, calves undergo challenging events such as weaning and group transition at around 60-70 d of age.Consequently, when measuring BW at 3 mo of age, the calves are in a post-transition adaptation phase, where there is usually greater environmental influence and, consequently, reduced heritability.
We first carried out single-trait GWAS for 24 growthrelated traits to identify trait-specific or age-specific loci in Holstein heifers, and subsequently more related traits were combined to detect markers with pleiotropic effects.Although a common strategy is to perform multivariate GWAS for multiple related phenotypes from one population, they are more computational demanding.Alternatively, multi-trait meta-analyses (or GWAS summary statistics) can be performed as an approximate method, as done in previous studies (e.g., Bolormaa et al., 2014;Fang et al., 2019;Chen et al., 2022).Our multi-trait meta-analyses resulted in the detection of 3 significant SNPs associated with growth traits in addition to those detected through single-trait GWAS.Despite the relatively high heritability for BW, HH, BL, and CG, few SNPs were detected.This might be due to the highly polygenic nature of these traits and relatively small sample size, as previous studies have shown that growth traits tend to be affected by many polymorphisms of small effects (Bouwman et al., 2018;Cai et al., 2022).Furthermore, clustering genetic  correlated traits may provide higher statistical power for detecting significant associations than combining less related traits (Cole et al., 2011;Fang and Pausch, 2019).Not all growth traits are genetically correlated with each other at a moderate to high level.Genetic correlations decreased with increasing measurement interval for BW, HH, BL, and CG, and a negative genetic correlation (−0.54) was also found between BL3 and CG3.This result is consistent with the findings of Baldi et al. (2010) and Yin and Konig (2018).Furthermore, based on the single-trait GWAS, 6 significant SNPs overlapped between the same or different traits at similar ages after 6 mo of age, such as BW8 and BW10, and HH6 and CG6.The maximum measurement interval between 2 traits with common significant SNPs was 4 mo.The genetic mechanism of growth might be that different genes are "switched on or off" with aging and development.The allele substitution effect of each SNP was also not constant for the same trait at different growth stages.This is evident in the case of the various effects of BTA7: 3 ,676 ,741 (ARS-BFGL-NGS-104759) on BW from birth to 12 mo of age, and it plays an important role in BW during the period between 4 and 6 mo of age only.However, the effect of significant SNPs associated with HH, BL, and CG showed no regular changes over time due to the greater measurement interval compared with BW, highlighting the need for continuous observations.Utsunomiya et al. ( 2013) also single-trait GWAS and multi-trait meta-analysis for body weight; GWAS_HH and Meta_HH: single-trait GWAS and multi-trait meta-analysis for hip height; GWAS_BL and Meta_BL: single-trait GWAS and multi-trait meta-analysis for body length; GWAS_CG and Meta_CG: singletrait GWAS and multi-trait meta-analysis for chest girth; Meta_M3: multi-trait meta-analysis for growth traits at 3 to 4 mo of age; Meta_M6: multi-trait meta-analysis for growth traits at 6 to 7 mo of age; Meta_M9: multi-trait meta-analysis for growth traits at 9 to 10 mo of age; Meta_M12: multi-trait meta-analysis for growth traits at 12 mo of age.
suggested that continuous measurement is important for the identification of genes and polymorphisms underlying growth traits.The inclusion of growth traits measured at different age points should be considered in Holstein selection strategies.
The multi-trait meta-analyses (or GWAS summary statistics) performed in this study have been widely applied in analyses of different traits of the same populations (Guo et al., 2017;Fang and Pausch, 2019;Tahir et al., 2021).One of the drawbacks of combining results from multiple studies is that the effect of all (significant and non-significant) SNPs and their standard error are usually unavailable.Therefore, we performed COJO using a single GWAS cohort with individual-level genotypes instead of the summary-level statistics from the meta-analysis.In the COJO analysis, 59 top SNPs remained significant, of which 52% were more powerful than the single-trait GWAS.To minimize the type-I error rate, the residual variance was considered constant at the same level of the phenotypic variance, even after fitting SNPs that cumulatively explain a substantial proportion of the total phenotypic variation based on the GCTA-COJO method (Yang et al., 2012).However, the problem of over-fitting still exists because of a small discovery set.The significance of the most associated SNP (BovineHD1600002786 located at BTA16: 9 ,745 ,951) increased from P GWAS = 2.05 × 10 −13 to P COJO = 3.29 × 10 −151 .In this case, to avoid false positives, we focused on these SNPs identified repeatedly in GWAS and COJO, though some important signals may have been ignored.Future studies using whole-genome sequence (WGS) data, larger data sets, and more comprehensive statistical models (e.g., random regression models) will be carried out.Some of the genes located nearby significant SNPs in our study have been described in other mammal species, thus confirming the relevance of these genomic regions for growth traits in dairy cattle.Among these genes, the high mobility group AT-hook 2 HMGA2 is a key candidate gene, suggested by the peak SNP on BTA5, which affects growth, fatness, and skeletal system and has been previously reported in humans (Yang et al., 2010;N'Diaye et al., 2011;Horikoshi et al., 2013), mice (Yuan et al., 2017;Lee et al., 2022;Negishi et al., 2022), horses (Makvandi-Nejad et al., 2012;Kader et al., 2016), dogs (Zapata et al., 2016), chicken (Song et al., 2011), pigs (Chung et al., 2018), and cattle (Bolormaa et al., 2014;Bouwman et al., 2018).Three independent association studies have revealed a significant association of the rs1042725 genotypes CT and CC in the 3′ UTR of HMGA2 with height across human populations (Weedon et al., 2007;Sanna et al., 2008;Guo et al., 2017).Buysse et al. (2009) and Lynch et al. (2011) provided strong evidence that a heterozy-gous deletion in HMGA2 is causing the abnormality of body height.It is consistent with that of study where the loss of one or both alleles of HMGA2 resulted in reduced body size in mice (Lee et al., 2022), and another gene knockout study of the mouse counterpart demonstrated that this gene is involved in diet-induced obesity (Yuan et al., 2017).HMGA2 can also lead to dwarfism in chicken (Ruyter-Spira et al., 1998) and pigs (Chung et al., 2018).In addition, 3 SNPs in HMGA2 were reported to be associated with body weight at the ages of 6, 7, 9, and 12 weeks in chickens (Song et al., 2011).Our findings agree with these literature reports, in which the significant variant BovineHD0500013922 located downstream of HMGA2 explained around 1.3% phenotypic variance for BW12 in Holstein cattle.Bouwman et al. (2018) identified 163 significant genomic regions regulating body size in cattle, including one missense variant in HMGA2.The frizzled related protein (FRZB), a member of the Wnt signaling pathway, near the significant SNP (ARS-BFGL-NGS-89201) at BTA2: 13 ,612 ,446 was also associated with BW12 in this study.FRZB is involved in the regulation of bone development (Hoang et al., 1996).FRZB gene knockout mice have lower body and muscle mass compared with wildtypes (Casas-Fraile et al., 2020).In pigs, FRZB is also a major candidate gene for growth traits.Wang et al. (2014) indicated that its expression has a negative association with muscle growth and a positive association with fat deposition.
The most significant SNP identified for BW was found based on the Meta_M6 and Meta_BW analyses.The closest genes to this top SNP include ZNF354B, ATP13A1, GMIP, LPAR2, PBX4, CILP2, YJEFN3, NDUFA13, TSSK6, and GATAD2A.Among the genes identified, GATAD2A (ATA zinc finger domain containing 2A) was previously reported to be associated with human height and appendicular lean mass, and it is an essential gene of early development in mice (Marino and Nusse, 2007).With strong support from functional activity, 2 genes associated with BW5 (HNF4G and TRIB3) and one gene associated with BW8 (BHLHE40) were considered to be key candidate genes.Ayari et al. (2022) demonstrated that HNF4G gene (hepatocyte nuclear factor 4 gamma) invalidation prevents diet-induced obesity via intestinal lipid malabsorption, however, this outcome is contrary to that of Gerdin et al. (2006) who found that HNF4G knockout mice had a higher body weight despite having reduced feed and water intake.In previous studies, genetic variants of HNF4G have been associated with height and obesity in humans (Berndt et al., 2013) and intramuscular fat deposition in beef cattle (Ramayo-Caldas et al., 2014).Another TRIB3 gene (tribbles pseudokinase 3) plays an important role in energy metabolism and Ma et al.: META-ANALYSES FOR DAIRY CATTLE GROWTH Ma et al.: META-ANALYSES FOR DAIRY CATTLE GROWTH NGS-104759; P = 5.06 × 10 −12 Figure 1.Estimates of genetic correlations for growth traits in Holstein heifers.The diagonals are the abbreviations of growth traits: BW, HH, BL, and CG represent the body weight, hip height, body length and chest girth, respectively, and the number suffix on the acronym represents the month of age.For example, HH3 = hip height at 3 mo of age; BW34 = body weight at either 3 or 4 mo of age (depending on which measurement was closer to the time point of body size).
Figure 2. Results of single-trait genome-wide association studies for 24 growth traits.a) Number of significant SNPs identified in each of the single-trait genome-wide association studies.b) The -log 10 (P) of 170 SNPs associated with growth traits in each of the single-trait genome-wide association studies.BW, HH, BL, and CG represents the body weight, hip height, body length, and chest girth, respectively, and the number suffix on the acronym represents the month of age (e.g., HH6 = hip height at 6 mo of age).The SNPs listed at the top of the plots are significantly associated with 2 traits.
Figure 3. Manhattan plots of the 8 multi-trait meta-analyses.In each plot, the red dots represent significant SNPs according to a false discovery rate (FDR) of 10% in the corresponding meta-analysis related to: a) body weight, b) hip height, c) body length, d) chest girth,e) growth traits at 3 to 4 mo of age, f) growth traits at 6 to 7 mo of age, g) growth traits at 9 to 10 mo of age, h) growth traits at 12 mo of age.
Ma et al.: META-ANALYSES FOR DAIRY CATTLE GROWTH

Figure 4 .
Figure 4. Number of significant SNPs single-trait GWAS and multi-trait meta-analysis.a) Overlapping SNPs between singletrait GWAS and multi-trait meta-analysis.b) Overlapping SNPs among methods and traits, denoted by connected black circles below the histogram.Orange bars represent that overlap among at least 3 categories, and horizontal bars show SNPs set size.GWAS_BW and Meta_BW:single-trait GWAS and multi-trait meta-analysis for body weight; GWAS_HH and Meta_HH: single-trait GWAS and multi-trait meta-analysis for hip height; GWAS_BL and Meta_BL: single-trait GWAS and multi-trait meta-analysis for body length; GWAS_CG and Meta_CG: singletrait GWAS and multi-trait meta-analysis for chest girth; Meta_M3: multi-trait meta-analysis for growth traits at 3 to 4 mo of age; Meta_M6: multi-trait meta-analysis for growth traits at 6 to 7 mo of age; Meta_M9: multi-trait meta-analysis for growth traits at 9 to 10 mo of age; Meta_M12: multi-trait meta-analysis for growth traits at 12 mo of age.
Ma et al.: META-ANALYSES FOR DAIRY CATTLE GROWTH

Table 1 .
Ma et al.: META-ANALYSES FOR DAIRY CATTLE GROWTH Numbers of animals, mean, standard deviation (SD), minimum, and maximum values for each growth-related trait in Holstein cattle ; BLi = body length (in cm) recorded at i month of age (i = 3, 6, 9, 12); CGi = chest girth (in cm) recorded at i month of age (i = 3, 6, 9, 12). 2 Number of animals used for estimating variance components and breeding values. 3Number of animals used for the single-trait GWAS analyses.

Table 3 .
Effects estimated in corresponding single-trait GWAS of genome-wide significant SNPs detected by the meta-analysis related to body weight (Meta_BW), hip height (Meta_HH), and chest girth(Meta_CG)