Determination of genetic associations between indels in 11 candidate genes and milk composition traits in Chinese Holstein population

Background We have previously identified 11 promising candidate genes for milk composition traits by resequencing the whole genomes of 8 Holstein bulls with extremely high and low estimated breeding values for milk protein and fat percentages (high and low groups), including FCGR2B, CENPE, RETSAT, ACSBG2, NFKB2, TBC1D1, NLK, MAP3K1, SLC30A2, ANGPT1 and UGDH those contained 25 indels between high and low groups. In this study, the purpose was to further examine whether these candidates have significant genetic effects on milk protein and fat traits. Results With PCR product sequencing, 13 indels identified by whole genome resequencing were successfully genotyped. With association analysis in 769 Chinese Holstein cows, we found that the indel in FCGR2B was significantly associated with milk yield, protein yield and protein percentage (P = 0.0041 to 0.0297); five indels in CENPE and one indel in MAP3K1 were markedly relevant to milk yield, fat yield and protein yield (P < 0.0001 to 0.0073); polymorphism in RETSAT was evidently associated with milk yield, fat yield, protein yield and protein percentage (P = 0.0001 to 0.0237); variant in ACSBG2 affected fat yield and protein percentage (P = 0.0088 and 0.0052); one indel in TBC1D1 was with respect to fat percentage and protein percentage (P = 0.0224 and 0.0209). Significant associations were shown between indels in NLK and protein yield and protein percentage (P = 0.0012 to 0.0257); variant in UGDH was related to the milk yield (P = 0.0312). The two exonic indels in FCGR2B and CENPE were predicted to change the mRNA and protein secondary structures, and resulted in the corresponding protein dysfunction. Conclusion Our findings presented here provide the first evidence for the associations of eight functional genes with milk yield and composition traits in dairy cattle. Electronic supplementary material The online version of this article (10.1186/s12863-019-0751-y) contains supplementary material, which is available to authorized users.

investigations on genetic and phenotypic diversities in human, chicken, pig and dairy cattle [9][10][11][12][13]. A previous study found that 2-18 base pairs (bp) indel located upstream of TAL bHLH transcription factor 1 (TAL1) was responsible for the T-cell acute lymphoblastic leukemia (T-ALL) [14]. In chicken, the 9-15 bp indel of premelanosome protein (PMEL17) gene was confirmed to be the causative mutation for the plumage color (Dominant white, Dun and Smoky) [10]. In pig, an intronic inserted retrotransposon of sperm flagellar 2 (SPEF2) led to the immotile short-tail sperm defect [15]. In Belgian blue cattle, a 11-bp indel in myostatin (MSTN) gene resulted in double-muscled phenotype [9], and an exonic 15-bp insertion in coagulation factor XI (F11) gene caused the factor XI deficiency in Japanese black cattle [11]. However, up to now, limited research of indel polymorphisms associated with milk production traits in dairy cattle has been reported [16].
With the rapidly emergence of next-generation sequencing (NGS), whole genome resequencing has been an important tool in the efforts to detect polymorphsims which were contributed to the complex traits or economic traits in human and domestic animals [17][18][19][20]. In our previous whole genome resequencing study, we identified over 0.9 million short indels and 3625 common differential indels with the same allelic distribution directions based on the 8 Holstein bulls with extremely high or low estimated breeding values (EBVs) of milk protein and fat percentages (high and low groups) [21]. Based on this, 11 genes were identified as the promising candidates affecting milk compositions traits in dairy cattle, including FCGR2B, CENPE, RETSAT, ACSBG2, NFKB2, TBC1D1, NLK, MAP3K1, SLC30A2, ANGPT1 and UGDH, which contained 25 differential indels [21]. Thus, the aim of this study was to further validate whether these identified indels in the 11 genes significantly impact on milk yield and compositions traits in Chinese Holstein population.

Indel verification and genotyping
Based on two DNA pools from 40 Holstein sires, with PCR product sequencing, 22 of 25 indels identified by whole genome resequencing [21] were confirmed as true ones (Additional file 2), among them, four indels were identified for the first time (Table 1). Subsequently, 13 indels in 8 genes were successfully genotyped and performed for association analysis. Of the 13 indels, two indels, including rs381714237 in FCGR2B, ss2137349053 in CENPE were located in the exons, whilst, the remaining 11 indels were located in the intronic regions. Chi-squared test showed that all the 13 indels were in Hardy-Weinberg equilibrium (P > 0.05). The genotype frequencies and allele frequencies of the 13 indels were summarized in Table 2.

Associations between indels and five milk production traits
The results of associations between the 13 indels and five milk production traits were shown in Table 3. It was observed that all these indels were significantly associated with at least one of the milk traits (P < 0.0001 to P = 0.0312) as described below.

Intronic indels
The four intronic indels (rs385060942, ss2137349051, rs453960300 and rs378415122) in CENPE were significantly associated with milk yield (P < 0.0001), fat yield (P = 0.0004 to 0.0073) and protein yield (P < 0.0001 to 0.0002). Additionally, the five indels (four intronic indels and one exonic indel above) of CENPE gene were found to be highly linked (r 2 > 0.98), and one haplotype block was inferred as presented in Fig. 1. Haplotype-based association analysis showed that the haplotype combination was evidently associated with milk yield, fat yield and protein yield as well (P < 0.0001 to P = 0.0076) ( Table 4).
Indel rs134985825 in the intron 6 of RETSAT showed remarkable effects on milk yield, protein yield, fat yield and protein percentage (P = 0.0001 to 0.0237). For ACSBG2, indel rs377943075 in the intron 7 was significantly associated with fat yield (P = 0.0088) and protein percentage (P = 0.0052). Variant rs136639319 in the intron 3 of TBC1D1was significantly associated with fat percentage (P = 0.0224) and protein percentage (P = 0.0209).
For the indel rs379188781 in the intron 1 of NLK, it was found to be associated with protein percentage (P = 0.0047), the other indel rs134444531 in the intron 3 was associated with protein yield (P = 0.0012) and protein percentage (P = 0.0257). While, no LD was observed between such two indels (r 2 = 0.14).
For MAP3K1, indel ss2137349058 in the intron 16 was markedly associated with milk yield, fat yield and protein yield (P < 0.0001).
For the intronic indel of UGDH gene, the indel ss2019489562 located in intron 2 was significantly associated with milk yield (P = 0.0312).
Additionally, the significant additive, dominant and allele substitution effects of the 13 indels on the five milk traits were observed as well (Table 5).

Prediction the mRNA and protein structures
Using a statistical folding algorithm, the alteration of the most stable mRNA secondary structures caused by the two exonic indels for FCGR2B and CENPE were observed for both the ins/ins and del/del genotypes. As illustrated in Fig. 2, obvious structural differences spanning the position 971-980 between the ins/ins and del/del genotypes of the indel rs381714237 in FCGR2B gene were observed. The free energy (ΔG) of the ins allele was predicted to be higher (ΔG = − 468.70 kcal/mol) than the del allele (ΔG = − 470.30 kcal/mol). Correspondingly, the ins allele was deduced to form one larger single loop structure, which potentially decreasing the stability of mRNA (ΔΔG = + 1.6 kcal/mol). It is worth mentioning that previous studies have evidenced that the ΔΔG ranged from − 3.9 kcal/mol to + 4.0 kcal/mol could affect the mRNA stability [22][23][24][25][26][27][28]. In addition, indel rs381714237 of FCGR2B was predicted to decrease the number of amino acid by 38, which might change protein structure and function. As a result, differences of the protein secondary structures were predicted between the FCGR2B proteins corresponding to alleles del and ins with regard to alpha helix (21.64% vs. 16.45%), extended strand (23.10% vs. 24.67%), beta turn (7.02% vs. 7.89%) and random coil (48.25% vs. 50.99%) using the SOPMA program.
For the non-frameshiting indel, subtle change of mRNA secondary structures between the two homozygous genotypes of indel ss2137349053 in CENPE was occurred (data not shown). The free energy was altered from − 1816.10 kcal/mol for the del allele to − 1818.90 kcal/mol for the ins allele. While, slight difference was predicted for the CENPE protein in accordance between the del/del and ins/ins genotypes, alpha helix (72.57% vs. 72.69%), and random coil (13.94% vs. 13.82%). There was no change of extended strand and beta turn for CENPE protein. Note: 1 indels were genotyped successfully; 2 indels were failed to genotype using MALDI-TOF MS; 3 indels were not polymorphic in current population; 4 primers of indel were failed to design

Discussion
In the present work, we confirmed that 13 indels belonging to 8 candidate genes (FCGR2B, CENPE, RETSAT, ACSBG2, TBC1D1, NLK, MAP3K1 and UGDH) for milk compositions identified by our previous whole genome resequencing study [21] showed significant genetic effects on at least one of milk traits in dairy cattle. As far as our knowledge, this is the first report to connect these genes to milk production traits of dairy cattle.    Among the total 25 differential indels with the same allelic distribution directions between the bulls in high and low groups identified by our previous whole genome resequencing study [21], indel rs383700527 (3 N ins) located upstream of ACSBG2 gene was found to contribute to milk fat in a cis-regulatory manner (unpublished data). Thus, we investigated another 24 indels in the present study. Among them, one intronic indel (4 N ins in CENPE) was failed to be verified by Sanger sequencing due to the special characteristic of the flanking sequence with lower GC% and repetitive DNA sequences. Two indels (3 N del in CENPE and 2 N ins in NFKB2) didn't show polymorphic in this study. Eight indels (1 N del and 21 N ins in CENPE, 2N del in RETSAT, 1 N ins and 1 N del in NLK, 6 N del in SLC30A2, 2 N ins in ANGPT1 and 1N ins in UGDH) were failed to be genotyped by using MALDI-TOF MS. The possible reason may be that MALDI-TOF MS for multiplex genotyping was relied on multiplex-PCR primers and extended primers to genotype multiple loci [29], simultaneously, the primer design was depended on sequence composition, molecular weight, annealing temperature and reaction efficiencies of each locus [29]. Hence, a total of 13 polymorphic indels were successfully genotyped and performed for association analysis.

Significant associations between candidate genes and milk production traits Six indels in FCGR2B and CENPE
For indel rs381714237 in FCGR2B, we demonstrated that ins/ins genotype had higher protein percentage. As a regulator, FCGR2B was contributed to immune response [30]. Additionally, bovine mammary gland is a product of the innate immune system and active during lactation. Thus, these evidences indicated that FCGR2B might affect milk protein percentage through impacting the cows on immune response during lactation.
For the five indels in CENPE, the association results revealed that ins/ins genotypes were dominant compared with del/del genotypes for milk yield, fat yield and protein yield. Previous report has found that CENPE acted as a monitor protein and was necessary for cell cycle [31]. Thus, it appeared that the CENPE might affect these traits through modulating bovine mammary gland development.

Seven indels located in six genes
Our association analysis confirmed that the ins/ins genotype of the indel rs134985825 in RETSAT gene increased milk yield, fat yield, and protein yield. RETSAT was considered as a regulator for liver metabolism, and was critical for lipid accumulation and adipogenesis promotion [32]. Previous research has investigated that the polymorphisms of RETSAT gene were associated with premium cut yields and backfat thickness in pig [33]. Taken together, we speculated that RETSAT might affect milk traits through influencing the lipid metabolism.
Herein, we found that individuals with ins/ins genotype of indel rs377943075 in ACSBG2 showed higher fat yield than those with del/del genotype. The ACSBG2 gene encodes the protein that belongs to a member of the acyl-CoA synthetase family and participated in PPAR signaling pathway and involved in lipid metabolism and lipid droplet formation [34,35]. Previous researchers have found that polymorphisms of ACSBG2 showed positive effects on yolk development and abdominal fat weight [36].
In current study, our results also showed a significant relationship between the indel rs136639319 in TBC1D1 and fat percentage as well as protein percentage. It was worth mentioning that TBC1D1, as a member of Rab GTPase-activating proteins (GAPs), was involved in translocation of GLUT4 to the plasma membrane. Polymorphisms in TBC1D1 have been observed to show significant effects on severe obesity or carcass in human [37] and chicken [38], respectively, suggesting exhibiting functions related to lipid and energy homeostasis as reported by Hargett et al. [39,40].
Two intronic indel (rs379188781 and rs134444531) in NLK showed strong associations with protein yield and protein percentage. Interestingly, Cole et al. reported that one single nucleotide polymorphism (SNP) (ARS-BFGL-NGS-106227) significantly associated with protein percentage (P = 5.59 × 10 − 8 ) was merely 90 kb away from the NLK gene [41].  Furthermore, NLK, as a member of MAPK subfamily, had an essential role in mediating the mTORC1 signaling pathway which was involved in milk protein synthesis [42,43]. Together, these data suggested that significant variation of protein yield and protein percentage might be regulated by NLK.
As for MAP3K1, individuals with del/del genotype of indel polymorphism ss2137349058 had higher milk yield, fat yield and protein yield. MAP3K1 was known to be involved in the MAPK signaling pathway, and was considered to be a metabolic stimuli inducing cell proliferation [44,45]. Meanwhile, it functioned as a candidate gene for type 2 diabetes (T2D) by interacting with insulin signaling pathway [46]. Thus, we concluded that MAP3K1 might regulate milk composition traits by modulating bovine mammary gland development.
The intronic indel ss2019489562 in UGDH showed significant effect on milk yield. UGDH encodes the protein that was implicated with biosynthesis of glycosaminoglycans, hyaluronan, chondroitin sulfate, and heparan sulfate. Previously, Xu et al. have demonstrated that two exonic SNPs in UGDH showed significant associations with milk production traits in Chinese Holstein population [47]. In particular, UGDH was close to the peak location of two reported QTLs for fat yield, fat percentage and protein yield [48][49][50][51]. Further, two previously reported significant SNPs for fat yield, protein yield, fat percentage and protein percentage [41] were near to the UGDH gene. Moreover, expression pattern in InnateDB showed that UGDH have the highest expression in liver which plays an indispensable role in metabolism of carbohydrates, fats and proteins in dairy cattle. Hence, these data demonstrated that UGDH gene might be a vital regulator for milk traits by affecting liver metabolism.

Conclusion
In the present study, we performed association analysis for the 13 short indels within 8 candidate genes for milk compositions identified by our previous whole genome resequencing study, including FCGR2B, CENPE, RETSAT, ACSBG2, TBC1D1, NLK, MAP3K1 and UGDH. As a result, the 13 indels were shown to have significant genetic effects on at least one of milk yield and composition traits. These results not only validated the candidate genes and indels from the previous whole genome resequencing work, but also provided novel molecular information for genetic improvement program of dairy cattle.

Ethics statement
All the procedures for sample collections and phenotypic observations of experimental individuals were

Animals
The animals used for association analysis included a total of 769 Chinese Holstein cows those were daugters of 40 sire families. These daughters were collected from 22 herds of Beijing Sanyuanlvhe Dairy Farming Center, a leading dairy company in China. Phenotypic data of the five milk production traits including 305-day milk yield (MY), fat yield (FY), protein yield (PY), fat percentage (FP) and protein percentage (PP) those were calculated based on at least 6 test-day records in each lactation using a multiple trait random regression test-day model by the Dairy Data Center of China, Dairy Association of China (http://www.holstein.org.cn/).
Genomic DNA was isolated from whole blood of cows and frozen semen of sires as previously described by Yang et al. [16].

Indels selection, PCR amplification, sequencing and genotyping
Of the 25 short indels that identified by our previous whole-genome resequencing study, 24 indels were investigated the associations with the five milk production traits except for a three-nucleotide insertion (3 N ins) in ACSBG2 gene.
A total of 23 pairs of PCR primers were designed with Primer Premier 5.0 and Oligo 7.0 softwares based on the genomic sequences of the 11 candidate genes in Bos_ taurus_UMD3.1 assembly (Additional file 1). To identify the twenty-four potential indel polymorphisms, two DNA pools for the above 40 sires were constructed with equal concentration of 50 ng/μl of each bull (20 individuals/pool). PCR products basd on the pooled DNA were purified with an EasyPure PCR Purification Kit (Trans-Gen Biotech, Beijing, China) and then bi-directionally sequenced using ABI3730xl DNA Analyzer (Applied Biosystems, Foster City, CA, USA).
To further confirm the position and sequence of the insertions and deletions, the purified PCR products were cloned into the pClone007 vector with a pClone 007 Vector Kit (TsingKe Biological Technology, Beijing, China). Positive clones including target indels were sequenced to search potential indels. The BLAST software (https://blast.ncbi.nlm.nih.gov/Blast.cgi) and Chromas 2.0 (Technelysium, Australia) were applied for sequence alignment to the reference sequence of the corresponding gene referring to Bos_taurus_UMD_3.1 assembly. Finally, genotyping for the identified indels in 769 chinese cows was performed by using the Sequenom MassArray matrix-assisted laser desorption ionization time-of-flight mass spectrometry (MALDI-TOF MS).

Bioinformatics analysis
To further explore the potential impact of the exonic indels in FCGR2B and CENPE on the mRNA secondary structures as well as the second structures of corresponding proteins, the online RNA FOLDING FORM (version 2.3) software [52] and SPOMA program (http:// npsa-pbil.ibcp.fr/) [53] were used, respectively.

Association analysis
Allele frequencies and genotype frequencies between the insertion and deletion genotypes, as well as the Hardy-Weinberg equilibrium were determined through a chisquare test. Associations between the 13 investigated indels and the five milk production traits were carried out by applying the mixed procedure in SAS 9.2 [54] based on the following linear mixed regression model: Where Y is the phenotypic record for the analyzed trait of the cows, μ is the overall mean of the phenotypic record for each trait, hys is a fixed effect of herd, year and season, b is linear regression coefficient on calving month (M), M is effect of calving month, G is a fixed effect of indel genotype or haplotype, a is a random polygenic effect account for all known pedigree relationships, and e is a random residual.

2
; d ¼ AB-AA þ BB . 2 and α ¼ a þ dðq-pÞ where AA, BB and AB were the least square means of the phenotypic values for corresponding genotypes, and p and q indicates the allele frequencies of the corresponding alleles. Multiple t-tests with Bonferroni correction were used to compare the effects of the genotypes on each indel.
The linkage disequilibrium (LD) extent among the genotyped indels (five indels in CENPE gene and two indels in NLK gene) and haplotype blocks was estimated using