Genome-wide Association Analysis of Body Conformation Traits in Chinese Holstein Cattle

doi:10.21203/rs.3.rs-5024087/v1

Background

The body conformation traits of dairy cattle are closely related to their production performance and health. The present study aimed to identify gene variants associated with body conformation traits in Chinese Holstein cattle and provide marker loci for genomic selection in dairy cattle breeding. The study findings could offer robust theoretical support to optimize the health of dairy cattle and enhance their production performance.

Results

This study involved 586 Chinese Holstein cows, using the predicted transmitting abilities (PTAs) of 17 body conformation traits evaluated by the Council on Dairy Cattle Breeding in the USA as phenotypic values. These traits were categorized into body size traits, rump traits, feet/legs traits, udder traits, and dairy characteristic traits. Based on the genomic profiling results from the Genomic Profiler Bovine 100K SNP chip, genotype data were quality-controlled using PLINK software, retaining 586 individuals and 80,713 SNPs for further analysis. Genome-wide association studies (GWAS) were conducted using the GEMMA software, employing both univariate linear mixed models (LMM) and multivariate linear mixed models (mvLMM). The Bonferroni method was used to determine the significance threshold, identifying gene variants significantly associated with body conformation traits in Chinese Holstein cows. The single-trait GWAS identified 24 SNPs significantly associated with body conformation traits (P < 0.01), with annotation leading to the identification of 21 candidate genes. The multivariate GWAS identified 54 SNPs, which were annotated to 57 candidate genes, including 39 new SNPs not identified in the single-trait GWAS. Additionally, 14 SNPs in the 86.84–87.41 Mb region of chromosome 6 were significantly associated with multiple traits such as body size, udder, and dairy characteristics. Four genes—SLC4A4, GC, NPFFR2, and ADAMTS3—were annotated in this region.

Conclusions

A total of 63 SNPs were identified as significantly associated with the 17 body conformation traits in Chinese Holstein cows through both single-trait and multivariate GWAS analyses. Sixty-six candidate genes were annotated, with 12 genes identified by both methods, including SLC4A4, GC, NPFFR2, and ADAMTS3, which are involved in biological processes such as active glucose transport, adipogenesis, and neural development. Thus, the study findings provided potential genetic marker information related to body conformation traits for the breeding of Chinese Holstein cattle.

Chinese Holstein Cattle

Body Conformation Traits

Single-trait GWAS

Multi-trait GWAS

The Chinese Holstein cattle is the first dairy breed developed in China and is the dominant breed in the country’s dairy cattle population, accounting for over 85% of the national herd^[1]. Although the body conformation traits of dairy cattle do not directly translate into economic benefits, these traits are closely related to the milk production capacity and overall health of cows^{[2, 3]}. These traits are an essential component for evaluating the overall performance of dairy cattle. Identifying candidate genes associated with body conformation traits in Chinese Holstein cattle is therefore crucial to provide effective molecular markers for genomic selection in breeding programs, further optimize breeding strategies, and improve the production performance of dairy cattle. The main body conformation traits in dairy cattle include (1) udder traits (such as fore udder attachment, front teat placement, teat length, rear udder height, rear udder width, and rear teat placement), (2) feet/legs traits (such as foot angle, heel depth, bone quality, rear legs side view, and rear legs rear view), (3) body size traits (such as stature, body depth, chest width, and strength), (4) rump traits (such as rump width and rump angle), and (5) dairy characteristics (such as angularity)^[4]. Since 1990, many countries have incorporated body conformation traits into dairy cattle breeding programs^[5] to improve the overall performance of dairy cows by optimizing these traits.

To gain a deeper understanding of the genetic basis of body conformation traits, genome-wide association studies (GWAS) are commonly used to identify candidate genes for economically important traits in dairy cattle. The mixed linear model proposed by Yu et al.^[6] has been widely recognized as the best GWAS analysis model currently available, as it effectively accounts for population structure and complex relationships within populations. The application of this model provides robust support for understanding the genetic mechanisms underlying body conformation traits in dairy cattle. Recently, Nazar et al.^[7] conducted a GWAS analysis on Chinese Holstein cattle using the mixed linear model and identified 18 SNPs significantly associated with five udder traits. Haque et al.^[8] were the first research team to report GWAS results for body conformation traits in Korean Holstein cattle, wherein they identified 24 SNPs significantly associated with 24 body conformation traits. Nazar et al.^[9] also conducted a GWAS analysis on three udder traits in Chinese Holstein cattle and detected nine SNPs significantly associated with udder traits after Bonferroni correction. Čítek et al.^[10] identified 32 SNPs significantly or nearly significantly associated with body conformation traits in the GWAS results of 25 body conformation traits in Czech Holstein cattle.

To date, most GWAS on body conformation traits in dairy cattle have been based on independent analyses of single traits. When multiple measured traits in an individual show a correlation, the lack of linkage between causal loci (linkage disequilibrium [LD]) is often overlooked. However, multivariate GWAS can jointly analyze genetically correlated traits by simultaneously considering within-trait and between-trait variations, thereby improving detection efficiency and accuracy^[11]. Numerous studies have shown that different body conformation traits are genetically correlated^{[4, 12, 13]}. To more comprehensively reveal the genetic variation and underlying mechanisms of body conformation traits in dairy cattle, multivariate GWAS methods are more advantageous. The present study conducted single-trait and multivariate GWAS analyses to determine the genetic variants associated with body conformation traits in Chinese Holstein cattle and provide marker information for genomic selection in dairy cattle; the findings of this study could offer strong theoretical support to further optimize the health of dairy cattle and enhance their production performance.

Genotype data and quality control

Genomic DNA was extracted from 586 Chinese Holstein cows by using the Tiangen Blood Genome Extraction Kit. Genotyping was performed using the GeneSeek Genomic Profiler Bovine 100K single nucleotide polymorphism (SNP) chip. Quality control of the individual and SNP data was conducted using PLINK (v1.9) software^[14], according to the following criteria: (1) individuals with an SNP missing rate of > 10% were excluded; (2) SNPs with a call rate of < 90% were excluded; (3) SNPs with a minor allele frequency (MAF) of < 5% were excluded; and (4) SNPs with a P-value of < 1.0 × 10^− 7 were excluded. LD analysis was performed using Haploview software^[15]. A total of 80,713 high-quality SNP markers from 586 individuals were ultimately selected. These SNP markers were evenly distributed across the chromosomes and were suitable for the subsequent GWAS analysis of Chinese Holstein cattle (Fig. 1).

Phenotype data collection

According to the updated standard methods of the Dairy Cattle Breeding Committee^{[16, 17]}, the predicted transmitting abilities (PTAs) of 17 body conformation traits, as assessed by the Council on Dairy Cattle Breeding, USA, were used as the phenotypic data for this study on 586 Chinese Holstein cows. These traits were categorized into five groups: (1) body size traits: stature (STA), strength (STR), and body depth (BDE); (2) rump traits: rump angle (RPA) and rump-thurl width (RTW); (3) feet and legs traits: rear legs side view (RLS), rear legs rear view (RLR), foot angle (FTA), feet/legs score (FLS), and feet/legs composite (FLC) index; (4) udder traits: fore udder attachment (FUA), rear udder height (RUH), rear udder width (RUW), udder cleft (UCL), udder depth (UDP), and udder composite (UDC) index; and (5) dairy form traits: dairy form (DFM).

Estimation of genetic parameters for body conformation traits

Genetic correlations were determined using the “-reml-bivar” parameter in GCTA software^[18]. Descriptive statistical analysis for the 17 body conformation traits, including maximum, minimum, mean, variance, and standard deviation, was conducted using SPSS 19 software. The frequency distribution histograms for each trait were plotted using the R package.

Single-trait and multi-trait GWAS

Before conducting a GWAS, principal component analysis (PCA) was performed using PLINK (v1.9) software to correct for potential false positives caused by population stratification. GWAS analysis between single conformation traits and genome-wide SNPs was performed using the univariate linear mixed model (LMM) in GEMMA software[19] by using the following model:

y = Χβ + Ζκγκ + ξ + ε

where y is the PTA vector, Χβ represents age and population structure effects (the first five principal components), Ζκγκ represents the effect of the marker to be tested, ξ ~ N (0,Kφ²) represents the polygenic effect, and ε ~ N (0,Iσ²) represents the residual effect. In the polygenic effect, K is the kinship matrix inferred from the markers.

Multi-trait GWAS analysis between multiple traits and genome-wide SNPs was conducted using the multivariate linear mixed model in GEMMA software^[20] by using the following model:

Y = Χβ + Ζκγκ + ξ + Ε

where Y is the n×d PTA matrix, n is the number of individuals in the population, and d is the number of traits analyzed; Χβ represents age and population structure effects (the first five principal components); Ζκγκ represents the effect of the marker to be tested; ξ ~ N (0,K,Vg) represents the polygenic effect; and Ε ~ N (0,In×n,Ve) represents the residual effect. In the polygenic effect, K is the kinship matrix inferred from the markers, Vg is the d×d polygenic variance-covariance matrix, and Ve is the d×d residual variance-covariance matrix. The number of independent SNPs was calculated using PLINK software with the “--indep pairs 50 5 0.2” command. The significance threshold was determined using Bonferroni correction (P_value < 0.05/number of independent SNPs). Manhattan plots and QQ plots were generated using the CMplot function in the R package.

Gene annotation

The SNP position information in the GeneSeek Genomic Profiler Bovine 100K SNP chip was based on Bos taurus UCD 1.2 version. The ARS-UCD 1.2 bovine reference genome information was downloaded from the ENSEMBL website. Significant SNPs within a 50-kb upstream and downstream range were annotated using ANNOVAR software^[21], and potential candidate genes related to the traits were identified.

Descriptive statistical analysis of the phenotypic data

Descriptive statistical analysis was performed on the 17 body conformation traits of 586 Chinese Holstein cows (Table 1). The phenotypic values of each trait generally followed a normal distribution(Figure S1).

Table 1. Descriptive statistics of body conformation traits

Traits		Min	Max	Mean	SD	Var
Body size	STA	-2.79	1.86	-0.3634	0.81616	0.666
	STR	-2.17	1.76	-0.3168	0.64178	0.412
	BDE	-2.03	1.57	-0.3865	0.67475	0.455
Rump	RPA	-2.83	2.01	-0.0458	0.74377	0.553
Rump	RTW	-2.53	1.94	-0.3065	0.79765	0.636
Feet/Legs	FLC	-2.06	0.99	-0.2341	0.50674	0.257
	RLS	-1.78	2.75	0.1781	0.73851	0.545
	RLR	-2.47	1.26	-0.44	0.67612	0.457
	FTA	-2.46	1.62	-0.335	0.63803	0.407
	FLS	-2.16	1.05	-0.2268	0.52009	0.27
Udder	UDC	-3.13	1.71	0.1534	0.64248	0.413
	FUA	-4.13	2.12	0.0933	0.82106	0.674
	RUH	-3.54	2.48	0.2525	0.85083	0.724
	RUW	-3.43	2.84	0.0274	0.92465	0.855
	UCL	-2.12	1.79	0.0703	0.64524	0.416
	UDP	-3.64	2.18	-0.0003	0.89008	0.792
	FTP	-2.33	2.45	0.1246	0.76152	0.58
	RTP	-2.3	2.68	0.1985	0.82223	0.676
	TLG	-2.53	1.95	-0.4231	0.7033	0.495
Dairy characteristics	DFM	-3.01	2.46	-0.124	0.91252	0.833

Determination of genetic correlations

The results of the determination of genetic correlations between the phenotypic traits indicated a strong positive correlation (r > 0.72) among the body size traits (STA, BDE, and STR; Table 2). In the rump traits, a weak negative correlation was observed between RTA and RPA (r = -0.02) (Table 3). Among the feet/legs traits, RLS showed varying degrees of negative correlation with other traits (-0.5 < r < -0.2), while strong positive correlations (r > 0.6) were found between RLR, FTA, FLS, and FLC index, except for RLS (Table 4). In the udder traits, varying degrees of positive correlation (0.1 < r < 1) were noted among UDC index, FUA, RUH, RUW, central ligament, and UDP (Table 5).

Table 2. Genetic correlation between body size traits

	STA	STR	BDE
STA	1.000
STR	0.723	1.000
BDE	0.810	0.854	1.000

Table 3. Genetic correlation between rump traits

	RPA	RTW
RPA	1.000
RTW	-0.024	1.000

Table 4. Genetic correlation between feet/legs traits

	FLC	RLS	RLR	FTA	FLS
FLC	1.000
RLS	-0.286	1.000
RLR	0.931	-0.347	1.000
FTA	0.661	-0.472	0.729	1.000
FLS	0.930	-0.220	0.912	0.760	1.000

Table 5. Genetic correlation between udder traits

	UDC	FUA	RUH	RUW	UCL	UDP
	UDC	FUA	RUH	RUW	UCL	UDP
UDC	1.000
FUA	1.000	1.000
RUH	0.893	0.686	1.000
RUW	0.756	0.510	0.871	1.000
UCL	0.659	0.438	0.567	0.635	1.000
UDP	0.697	0.864	0.472	0.221	0.353	1.000

Note: r ≤ 0.3 indicates weak correlation; 0.3 < r ≤ 0.4 indicates moderate correlation; r > 0.4 indicates strong correlation.^[22]

PCA

As shown in Figure 2, population stratification was observed. The explained variance percentage of the first 10 principal components (PCs) was calculated, and the results indicated that the first 5 PCs accounted for 80% of the variance. Therefore, in this study, the first 5 PCs were selected as covariates and included in the LMM for GWAS analysis.

Results of single-trait and multi-trait GWAS analyses

A single-trait GWAS was performed on the 17 body conformation traits, which identified 24 significant SNPs across 12 traits (Table 6). Figure 3 shows the QQ plots and Manhattan plots. In body size traits, one significant SNP was identified on chromosome 11, with one candidate gene, LDAH, annotated within the 50-kb upstream and downstream region of the SNP. In rump traits, four significant SNPs were identified on chromosome 7, with six candidate genes annotated, including OR2T4_1, ELAVL1, and LMAN2. In feet/legs traits, two significant SNPs were identified on chromosomes 8 and 29, with two candidate genes (PIP5K1B and NTM) annotated. In udder traits, eight significant SNPs were identified on chromosomes 5, 6, 7, and 19, with annotation of seven candidate genes, including ADGRE5, CCND2, and ARAP2. In milk production-related traits, nine significant SNPs were identified on chromosome 6, with annotation of four candidate genes, including SLC4A4, GC, NPFFR2, and ADAMTS3. One SNP on chromosome 19 (BovineHD1900013254) was significantly associated with UDC index, RUH, and RUW traits in the single-trait GWAS, thus indicating that it may be a pleiotropic locus.

Multi-trait GWAS was used to identify new significant SNPs, and 54 significant SNPs were identified. Figure 4 shows the QQ plots and Manhattan plots. Compared to the single-trait GWAS, 39 new SNPs were identified in the multi-trait GWAS (Table 7). In body size traits, 19 significant SNPs were identified across three combinations (STR-BDE, STA-STR, and STA-BDE), which were located on chromosomes 5, 6, 7, 8, 11, 19, and 29, and 14 candidate genes were annotated, including SLC4A4, GC, NPFFR2, and ADAMTS3. Of these 19 SNPs, 18 were newly identified SNPs. In rump traits, two significant SNPs were identified in the RPA-RTW combination on chromosome 7, and four candidate genes were annotated, including OR2T4_1, ATP8B3, and KLF16. Both SNPs were newly identified. In feet/legs traits, 11 significant SNPs were identified across five combinations (FLC-RLR, RLS-FTA-FLS, and FLC-RLS-FTA), which were located on chromosomes 8, 13, 21, 28, and 29, and 14 candidate genes were annotated, including GNAQ, SAMHD1, and SOGA1. Of these 11 SNPs, 9 were newly identified SNPs. In udder traits, 33 significant SNPs were identified across 41 combinations (FUA-RUH, FUA-RUH-RUW, and FUA-RUH-RUW-UCL), which were located on chromosomes 4, 6, 7, 11, 16, 17, 19, and 28, and 30 candidate genes were annotated, including BMT2, IGFBP1, and IGFBP3. Of these 33 SNPs, 28 were newly identified SNPs. In total, 32 SNPs were repeatedly detected across multiple trait combinations in the multi-trait GWAS, with 9 SNPs associated with body size and udder traits and 2 SNPs associated with feet/legs traits and udder traits, thus suggesting these SNPs may be pleiotropic loci.

A summary of the single-trait and multi-trait GWAS results for the 17 body conformation traits revealed 63 SNPs showing a significant association with body conformation traits in Chinese Holstein cows and 66 candidate genes annotated within the 50-kb region upstream and downstream of the significant loci. Additionally, the results of multi-trait GWAS showed that a region on chromosome 6 (86.84–87.41 Mb) contained 14 significant SNPs associated with udder and body size traits, and six haplotype blocks composed of 3, 2, 2, 9, 2, and 5 SNPs were observed (Figure 5). Eight of these SNPs exhibiting a significant association with milk production-related traits were also detected in the single-trait GWAS.

Table 6. Significant loci and genes identified by single-trait GWAS

SNP	Chr	Position	Traits	P_value	Nearby genes	Function area
BovineHD0600024228	6	86847656	DFM	1.97E-08	SLC4A4,GC	intergenic
ARS-BFGL-NGS-118182	6	86860291		2.27E-08	SLC4A4,GC	intergenic
BovineHD0600024243	6	86877334		2.27E-08	SLC4A4,GC	intergenic
BovineHD0600024355	6	87184768		7.02E-07	GC,NPFFR2	intergenic
BovineHD0600024357	6	87187812		7.02E-07	GC,NPFFR2	intergenic
chr6_89051385	6	87316810		7.29E-07	NPFFR2	exonic
DB-443-seq-rs110326785	6	87324678		7.29E-07	NPFFR2	exonic
DB-2033-seq-rs110186820	6	87368855		8.32E-07	NPFFR2	exonic
Hapmap60852-rs29024026	6	87725832		1.97E-06	ADAMTS3	exonic
ARS-BFGL-NGS-38413	11	78083413	STA	5.97E-06	LDAH	exonic
BovineHD0700011653	7	38822927	RPA	4.15E-06	LMAN2	exonic
BovineHD0700012421	7	41261062		6.06E-06	LOC526765,OR2W3,TRIM58	exonic
ARS-BFGL-NGS-13798	7	42322361		8.36E-07	OR2T4_1	exonic
BovineHD0700005042	7	16748554	RTW	7.32E-06	ELAVL1	exonic
BovineHD2900010725	29	34925638	RLS	3.21E-06	NTM	intronic
BTB-00344991	8	45049204	FLS	6.49E-06	PIP5K1B	exonic
Hapmap30832-BTA-144704	7	11394736	UDP	6.95E-06	ADGRE5	exonic
BovineHD0500013405	5	46349963	RUH	4.38E-06	LOC112446700,CAND1	intergenic
BovineHD0500013407	5	46351738	RUH	7.13E-06	LOC112446700,CAND1	intergenic
DB-364-seq-rs378727865	5	105784987	UCL	6.80E-06	CCND2	exonic
BovineHD0600015583	6	55313091		2.15E-06	ARAP2	intergenic
BTA-26162-no-rs	6	55327944		1.83E-06	ARAP2	intergenic
BovineHD1900012330	19	42905488	FUA	1.25E-06	LOC100138645,SAO	exonic
BovineHD1900013254	19	46942067	UDC	1.40E-06	TLK2	exonic
			RUH	6.46E-06
			RUW	2.70E-06

Table 7. Significant loci and genes identified by multi-trait GWAS

SNP	Chr	Position	Nearby genes	Function area
BTA-16397-no-rs	4	55471733	BMT2	exonic
BovineHD0400021231	4	76139944	IGFBP1,IGFBP3,LOC112446404	exonic
BovineHD0600015583	6	55313091	ARAP2	intergenic
BTA-26162-no-rs	6	55327944	ARAP2	intergenic
BovineHD0600024228	6	86847656	SLC4A4,GC	intergenic
ARS-BFGL-NGS-118182	6	86860291	SLC4A4,GC	intergenic
BovineHD0600024243	6	86877334	SLC4A4,GC	intergenic
BovineHD0600024315	6	87068809	GC,NPFFR2	intergenic
BovineHD0600024338	6	87130864	GC,NPFFR2	intergenic
BovineHD0600024345	6	87143505	GC,NPFFR2	intergenic
BovineHD0600024350	6	87153414	GC,NPFFR2	intergenic
BovineHD0600024355	6	87184768	GC,NPFFR2	intergenic
BovineHD0600024357	6	87187812	GC,NPFFR2	intergenic
BovineHD0600024365	6	87213962	GC,NPFFR2	intergenic
DB-442-seq-rs110392219	6	87314427	NPFFR2	exonic
chr6_89051385	6	87316810	NPFFR2	exonic
DB-443-seq-rs110326785	6	87324678	NPFFR2	exonic
DB-2033-seq-rs110186820	6	87368855	NPFFR2,ADAMTS3	intergenic
BovineHD0600028629	6	101050942	MAPK10	exonic
Hapmap30832-BTA-144704	7	11394736	ADGRE5	exonic
ARS-BFGL-NGS-13798	7	42322361	OR2T4_1	exonic
BovineHD0700013141	7	44152142	ATP8B3,KLF16,REXO1	exonic
ARS-BFGL-NGS-30237	8	53819922	GNAQ	intronic
BovineHD1100011070	11	37533658	EML6	exonic
Hapmap51861-BTA-86131	11	38522047	EFEMP1	exonic
chr11_38494447	11	38640411	MIR216A,MIR216B,MIR217	ncRNA_exonic
BovineHD1100011351	11	38650894	MIR216A,MIR216B,MIR217	ncRNA_exonic
BovineHD1100011736	11	39896854	TRNAY-AUA_2,TRNAC-GCA_141	intergenic
ARS-BFGL-NGS-38413	11	78083413	LDAH	exonic
BovineHD1100022676	11	79037363	TTC32,LOC104973438	intergenic
BovineHD1100022680	11	79049547	TTC32,LOC104973438	intergenic
BovineHD1100022725	11	79217388	LOC112448820,OSR1	intergenic
Hapmap47169-BTA-107308	11	80349202	LOC112448882,KCNS3	intergenic
BovineHD1100029423	11	101175170	FIBCD1	exonic
BovineHD1300018848	13	65929413	SAMHD1,SOGA1,TLDC2	exonic
ARS-BFGL-NGS-108133	16	52327847	TMEM51	UTR5
BovineHD1700014238	17	48992546	TMEM132C,TRNAC-GCA_192	intergenic
BovineHD1900012330	19	42905488	LOC100138645,SAO	exonic
BovineHD1900013254	19	46942067	TLK2	exonic
ARS-BFGL-NGS-115719	19	48149532	CD79B,SCN4A	exonic
BovineHD1900013592	19	48203740	LOC616254,PRR29	exonic
BovineHD1900013860	19	49052596	BPTF	exonic
BovineHD1900013865	19	49067122	BPTF	exonic
Hapmap47630-BTA-45710	19	49085782	BPTF	exonic
BovineHD2100019874	21	66165706	MEG9,LOC112443172	intergenic
BovineHD2100019929	21	66415334	LOC101907771,LOC112443172	ncRNA_exonic
ARS-BFGL-NGS-38270	21	66496288	DIO3	exonic
ARS-BFGL-NGS-86477	21	66743529	PPP2R5C	exonic
BovineHD2400000350	24	1265748	LOC104975719,TRNAK-UUU_41	intergenic
BovineHD2800005128	28	18663865	ZNF365,LOC112444734	intergenic
BovineHD2800005144	28	18784655	LOC101905431	exonic
BovineHD2800013715	28	19401828	JMJD1C	exonic
BovineHD2800005280	28	19436656	JMJD1C	intronic
UA-IFASA-6129	29	34835983	NTM	intronic

Genetic correlations among body conformation traits in Chinese Holstein cows

A very strong positive genetic correlation (r > 0.72) was observed between STA, BDE, and STR in body size traits, which is consistent with the findings of Ning^[13] and Degroot et al.^[23]. In feet/legs traits, negative genetic correlations were observed between RLS and other leg and hoof traits (-0.47 < r < -0.22). This was similar to the genetic correlation (-0.34) between RLS and RLR reported by Huang et al.^[24] in their estimation of genetic parameters for body conformation traits of dairy cows, although RLS showed a positive correlation with other leg and hoof traits in their study. Because leg and hoof traits are easily influenced by farm management and external environmental factors, the leg and hoof structure may vary between different cattle populations. In rump traits, the genetic correlation between RTA and RPA showed a weak negative correlation (r = -0.02), which is similar to the findings of Peng^[25] on the genetic correlation between RTA and RPA in Holstein cows in Hebei Province. However, the genetic correlations reported by Huang et al.^[24] and An et al.^[26] differed significantly (0.22 < r < 0.38). In udder traits, the genetic correlations ranged from 0.22 (RUW and UDP) to 1 (UDC index and FUA); this finding is similar to the results reported by Degroot et al.^[23]. The genetic correlations observed among body conformation traits suggest that the selection of one trait can indirectly influence other traits. Additionally, the positive correlations between traits indicate that these traits may share some common genetic basis, thus implying that certain genes or gene combinations may simultaneously affect multiple body conformation traits. Therefore, this genetic correlation can be used to develop more effective selection strategies.

Advantages of multi-trait GWAS

Body conformation traits are often controlled by multiple genes. Multi-trait GWAS can leverage the correlation between traits and combine weak genetic effects to enhance the statistical power of GWAS and improve the ability to detect new SNP loci^{[27, 28]}. In the present study, the body conformation traits showed genetic correlations. Compared to single-trait GWAS, 39 new SNPs were identified in the multi-trait GWAS. By using a similar strategy, Li et al.^[28] conducted a multi-trait GWAS on weaning weight and yearling weight in sheep and identified 93 new SNPs. Gao et al.^[22] discovered three new SNPs in carcass weight, carcass length, and chest depth in Huaxi cattle. These results suggest that when traits show genetic correlations among them, multi-trait GWAS can complement the findings of single-trait GWAS, thereby increasing the statistical power of GWAS.

Critical candidate genes

By using a combination of single-trait and multi-trait GWAS analyses, we detected 63 significant SNP loci and annotated 66 candidate genes. Among these, 12 genes were identified by both single-trait and multi-trait GWAS, including SLC4A4, GC, NPFFR2, ADAMTS3, LDAH, OR2T4_1, NTM, ADGRE5, ARAP2, TLK2, SAO, and LOC100138645. Four genes were located in the 86.84–87.41 Mb region of chromosome 6. Among these genes, SLC4A4 is a solute transporter and a member of a major transporter superfamily involved in active glucose transport^[29]. GC is a gene encoding the vitamin D-binding protein, which is specifically expressed in tissues such as the liver^{[30, 31]}. NPFFR2 is a member of the G-protein-coupled neuropeptide receptor subfamily activated by neuropeptides A-18-amide and F-8-amide^[32]. The ADAMTS3 gene activates vascular endothelial growth factors and promotes lymphangiogenesis^[33]. SNPs in this region were associated with dairy traits in single-trait GWAS and with body size and udder traits in multi-trait GWAS. Jiang et al.^[34] also found that the four genes in this region were related to milk yield and milk protein content in Holstein cows. Liang et al.^[35] conducted a GWAS of over one million US Holstein cows and found that the SLC4A4, GC, and NPFFR2 genes were related to fertility traits. Wu et al.^[36] reported that SLC4A4 and NPFFR2 were candidate genes for mastitis susceptibility in Danish Holstein cows. These results suggest that the four genes in this region may exhibit pleiotropy.

LDAH, a lipid droplet-associated hydrolase, is highly expressed in tissues that primarily store triacylglycerol and plays a key role in lipogenesis^[37]. Previous studies have shown that it is associated with hoof and leg diseases in Danish Holstein cows^[38]. NTM is a neurotrimin protein with an important role in neurodevelopment. Xu et al.^[39] analyzed the imprinting status of the NTM gene in cattle and identified an SNP (rs42185569) within the NTM gene through direct sequencing of PCR products. RT-PCR amplification showed monoallelic expression of the NTM gene in bovine placenta and adult tissues, thus suggesting that NTM is an imprinted gene in cattle. ADGRE5 primarily functions in cell adhesion and transport proteins and is associated with angiogenesis in cattle^[40]. ARAP2, involved in the endocytosis pathway, could be a candidate gene affecting loin strength in Chinese Holstein cows^[41]; we also speculate that it could be a candidate gene influencing udder traits in dairy cows. TLK2 is a Tousled-like kinase, and its main function involves phosphorylation of histone chaperones ASF1a and ASF1b and promotion of DNA replication-coupled nucleosome assembly, which is crucial for genome maintenance and proper cell division in plants and animals^[42]. SAO encodes a copper-containing amine oxidase that oxidizes spermine and plays an important role in polyamine metabolism in cattle^[43]. LOC100138645 is a primary amine oxidase and a liver isoenzyme that functions in amine metabolism. Thus, the candidate genes identified in the GWAS are involved in various important biological functions such as active glucose transport and lipogenesis. Notably, the four genes (SLC4A4, GC, NPFFR2, and ADAMTS3) in the 86.84–87.41 Mb region of chromosome 6 exhibit significant pleiotropy and may play roles in multiple economically important traits in dairy cattle, including milk production, body conformation, and reproductive health. These findings provide valuable genetic markers for further research on molecular breeding and functional validation in dairy cows.

In the present study, individual genotyping was performed using the Genomic Profiler Bovine 100K SNP chip, with a focus on the expected transmitting abilities of 17 body conformation traits in 586 Chinese Holstein cows. By performing single-trait and multi-trait GWAS analyses, 63 significantly associated SNP loci were identified, and 66 candidate genes were annotated, with detection of 12 genes by both methods. These genes are widely involved in various biological processes such as active glucose transport, lipogenesis, and neurodevelopment. Additionally, a genomic region significantly associated with body conformation traits was identified in the 86.84–87.41 Mb region on chromosome 6; this region included four candidate genes (SLC4A4, GC, NPFFR2, and ADAMTS3), which may be significantly related to body conformation traits in Chinese Holstein cows. The results of this study provide potential genetic markers for genomic selection breeding and related analyses in dairy cattle.

PTAs	predicted transmitting abilities
GWAS	Genome-wide association study
SNP	Single Nucleotide Polymorphism
FLC	Feet/Legs Composite
UDC	Udder Composite
STA	Stature
STR	Strength
BDE	Body Depth
DFM	Dairy Form
RPA	Rump Angle
RTW	Rump-Thurl Width
RLS	Rear Legs Side View
RLR	Rear Legs Rear View
FTA	Foot Angle
FLS	Feet/Legs Score
FUA	Fore Udder Attachment
RUH	Rear Udder Height
RUW	Rear udder width
UCL	Udder Cleft
UDP	Udder Depth
PCA	Principal Component Analysis
LD	Linkage disequilibrium

Acknowledgments

Thanks to all the authors for their contributions to the study.

Author contributions

Conceptualization: SL, YC and YM; Data curation: SL, LC and YL; Formal analysis and visualization: SL and LC; Writing the paper: SL and LC; Critical revision of the manuscript: SL, LC, YL, FG, HJ, HW, YC and YM.

Funding

This work was supported by the Breeding Industry Special Project of Tianjin Academy of Agricultural Sciences (2023ZYCX011 and 2024ZYCX012), Tianjin Seed Industry Special Project (22ZXZYSN00020), and a special financial aid from the Xizang Autonomous Region.

Data availability

The data were confidential and not deposited in an official repository.

Ethics approval and consent to participate

The animal study protocol was approved by the Science Research Department of the Institute of Animal Science, Tianjin Academy of Agricultural Science.

Consent for publication

Not applicable.

Competing interests

The authors declare no competing interests.

ZHANG S L, SUN D X. The past, present, and future of the dairy cattle breeding industry [J]. China Animal Husbandry Industry, 2021(15):22-26. (in Chinese)
ABO-ISMAIL M K, BRITO L F, MILLER S P, et al. Genome-wide association studies and genomic prediction of breeding values for calving performance and body conformation traits in Holstein cattle[J]. Genet Sel Evol, 2017,49(1):82.
KOCK A, LEDINEK M, GRUBER L, et al. Genetic analysis of efficiency traits in Austrian dairy cattle and their relationships with body condition score and lameness[J]. J Dairy Sci, 2018,101(1):445-455.
XUE X, HU H, ZHANG J, et al. Estimation of Genetic Parameters for Conformation Traits and Milk Production Traits in Chinese Holsteins[J]. Animals (Basel), 2022,13(1).
WU X, FANG M, LIU L, et al. Genome wide association studies for body conformation traits in the Chinese Holstein cattle population[J]. BMC Genomics, 2013,14:897.
YU J, PRESSOIR G, BRIGGS W H, et al. A unified mixed-model method for association mapping that accounts for multiple levels of relatedness[J]. Nat Genet, 2006,38(2):203-208.
NAZAR M, ABDALLA I M, CHEN Z, et al. Genome-Wide Association Study for Udder Conformation Traits in Chinese Holstein Cattle[J]. Animals (Basel), 2022,12(19).
HAQUE M A, ALAM M Z, IQBAL A, et al. Genome-Wide Association Studies for Body Conformation Traits in Korean Holstein Population[J]. Animals (Basel), 2023,13(18).
NAZAR M, LU X, ABDALLA I M, et al. Genome-Wide Association Study Candidate Genes on Mammary System-Related Teat-Shape Conformation Traits in Chinese Holstein Cattle[J]. Genes (Basel), 2021,12(12).
CITEK J, BRZAKOVA M, BAUER J, et al. Genome-Wide Association Study for Body Conformation Traits and Fitness in Czech Holsteins[J]. Animals (Basel), 2022,12(24).
KORTE A, VILHJALMSSON B J, SEGURA V, et al. A mixed-model approach for genome-wide association studies of correlated traits in structured populations[J]. Nat Genet, 2012,44(9):1066-1071.
SONG Y T, ZHANG R M, LI Y Q, et al. Estimation of genetic parameters for body conformation traits and the influence of pedigree generation in Shandong Holstein cattle [J]. Acta Veterinaria et Zootechnica Sinica, 2022, 53(05):1384-1395. (in Chinese)
NING J H, TIAN J Y. Correlation analysis of body conformation traits in Beijing Holstein cattle [J]. Shaanxi Agricultural Science, 2023, 69(11):91-94. (in Chinese)
PURCELL S, NEALE B, TODD-BROWN K, et al. PLINK: a tool set for whole-genome association and population-based linkage analyses[J]. Am J Hum Genet, 2007,81(3):559-575.
BARRETT J C, FRY B, MALLER J, et al. Haploview: analysis and visualization of LD and haplotype maps[J]. Bioinformatics, 2005,21(2):263-265.
COLE J B, WAURICH B, WENSCH-DORENDORF M, et al. A genome-wide association study of calf birth weight in Holstein cattle using single nucleotide polymorphisms and phenotypes predicted from auxiliary traits[J]. J Dairy Sci, 2014,97(5):3156-3172.
LI B, VANRADEN P M, NULL D J, et al. Major quantitative trait loci influencing milk production and conformation traits in Guernsey dairy cattle detected on Bos taurus autosome 19[J]. J Dairy Sci, 2021,104(1):550-560.
YANG J, LEE S H, GODDARD M E, et al. GCTA: a tool for genome-wide complex trait analysis[J]. Am J Hum Genet, 2011,88(1):76-82.
ZHOU X, STEPHENS M. Genome-wide efficient mixed-model analysis for association studies[J]. Nat Genet, 2012,44(7):821-824.
ZHOU X, STEPHENS M. Efficient multivariate linear mixed model algorithms for genome-wide association studies[J]. Nat Methods, 2014,11(4):407-409.
WANG K, LI M, HAKONARSON H. ANNOVAR: functional annotation of genetic variants from high-throughput sequencing data[J]. Nucleic Acids Res, 2010,38(16):e164.
GAO H, GE F, WANG Z Z, et al. GWAS analysis of carcass traits in Huaxi cattle [J]. China Animal Husbandry Journal, 2022, 58(11):92-99. (in Chinese)
DEGROOT B J, KEOWN J F, Van VLECK L D, et al. Genetic parameters and responses of linear type, yield traits, and somatic cell scores to divergent selection for predicted transmitting ability for type in Holsteins[J]. J Dairy Sci, 2002,85(6):1578-1585.
HUANG Y C, ZHANG H L, XU W, et al. Estimation of genetic parameters for body conformation traits in dairy cattle in the Ningxia region [J]. China Animal Husbandry & Veterinary Medicine, 2024, 51(07):2908-2922. (in Chinese)
PENG P, JIANG G E, YANG C D, et al. Genetic parameter analysis of body conformation traits in Holstein cattle in Hebei Province [J]. China Dairy Cattle, 2021(08):23-27. (in Chinese)
AN Y F, LI P J, WANG J H, et al. Parameter analysis of body conformation traits in Hebei Holstein cattle [J]. Chinese Journal of Animal Science, 2011, 37(02):6-10. (in Chinese)
ZHOU X, XIANG X, CAO D, et al. Selective sweep and GWAS provide insights into adaptive variation of Populus cathayana leaves[J]. Forestry Research, 2024,4(1).
LI Y, YANG H, GUO J, et al. Uncovering the candidate genes related to sheep body weight using multi-trait genome-wide association analysis[J]. Front Vet Sci, 2023,10:1206383.
PEDROSA V B, SCHENKEL F S, CHEN S Y, et al. Genomewide Association Analyses of Lactation Persistency and Milk Production Traits in Holstein Cattle Based on Imputed Whole-Genome Sequence Data[J]. Genes (Basel), 2021,12(11).
FREEBERN E, SANTOS D, FANG L, et al. GWAS and fine-mapping of livability and six disease traits in Holstein cattle[J]. BMC Genomics, 2020,21(1):41.
LEE Y L, TAKEDA H, COSTA M M G, et al. A 12 kb multi-allelic copy number variation encompassing a GC gene enhancer is associated with mastitis resistance in dairy cattle[J]. PLoS Genet, 2021,17(7):e1009331.
SUN Y L, ZHANG X Y, SUN T, et al. The anti-inflammatory potential of neuropeptide FF in vitro and in vivo[J]. Peptides, 2013,47:124-132.
JANSSEN L, DUPONT L, BEKHOUCHE M, et al. ADAMTS3 activity is mandatory for embryonic lymphangiogenesis and regulates placental angiogenesis[J]. Angiogenesis, 2016,19(1):53-65.
JIANG J, MA L, PRAKAPENKA D, et al. A Large-Scale Genome-Wide Association Study in U.S. Holstein Cattle[J]. Front Genet, 2019,10:412.
LIANG Z, PRAKAPENKA D, VANRADEN P M, et al. A Million-Cow Genome-Wide Association Study of Three Fertility Traits in U.S. Holstein Cows[J]. Int J Mol Sci, 2023,24(13).
WU X, LUND M S, SAHANA G, et al. Association analysis for udder health based on SNP-panel and sequence data in Danish Holsteins[J]. Genet Sel Evol, 2015,47(1):50.
GOO Y H, SON S H, PAUL A. Lipid Droplet-Associated Hydrolase Promotes Lipid Droplet Fusion and Enhances ATGL Degradation and Triglyceride Accumulation[J]. Sci Rep, 2017,7(1):2743.
WU X, GULDBRANDTSEN B, LUND M S, et al. Association analysis for feet and legs disorders with whole-genome sequence variants in 3 dairy cattle breeds[J]. J Dairy Sci, 2016,99(9):7221-7231.
DA X, JUNLIANG L, CUI Z, et al. The analysis of splice variants and genomic imprinting status of NTM gene in cattle (Bos taurus).[J]. Journal of Agricultural Biotechnology, 2018,26(10):1707-1713.
TALKER S C, BARUT G T, LISCHER H, et al. Monocyte biology conserved across species: Functional insights from cattle[J]. Front Immunol, 2022,13:889175.
LU X, ABDALLA I M, NAZAR M, et al. Genome-Wide Association Study on Reproduction-Related Body-Shape Traits of Chinese Holstein Cows[J]. Animals (Basel), 2021,11(7).
SIMON B, LOU H J, HUET-CALDERWOOD C, et al. Tousled-like kinase 2 targets ASF1 histone chaperones through client mimicry[J]. Nat Commun, 2022,13(1):749.
CERVELLI M, LEONETTI A, CERVONI L, et al. Stability of spermine oxidase to thermal and chemical denaturation: comparison with bovine serum amine oxidase[J]. Amino Acids, 2016,48(10):2283-2291.

No competing interests reported.

FigureS1Distributionoffrequencyof17bodyconformationtraits..pdf
Figure S1 Distribution of frequency of 17 body conformation traits.FLC, feet/legs composite index; UDC, udder composite index; STA, stature; STR, strength; BDE, body depth; DFM, dairy form; RPA, rump angle; RTW, rump−thurl width; RLS, rear legs side view; RLR, rear legs rear view; FTA, foot angle; FLS, feet/legs score; FUA, fore udder attachment; RUH, rear udder height; RUW, rear udder width; UCL, udder cleft; UDP, udder depth.
TableS1SignificantlociandgenesidentifiedbydifferentcombinationsofmultitraitGWAS.xlsx
Table S1 Significant loci and genes identified by different combinations of multi-trait GWAS.

Genome-wide Association Analysis of Body Conformation Traits in Chinese Holstein Cattle

Status:

Version 1

Abstract

Background

Results

Conclusions

Figures

Background

Materials and Methods

Genotype data and quality control

Phenotype data collection

Estimation of genetic parameters for body conformation traits

Single-trait and multi-trait GWAS

Gene annotation

Results

Discussion

Genetic correlations among body conformation traits in Chinese Holstein cows

Advantages of multi-trait GWAS

Critical candidate genes

Conclusion

Abbreviations

Declarations

References

Additional Declarations

Supplementary Files

Status:

Version 1