A Candidate Gene Association Study for Economically Important Traits in Czech Dairy Goat Breeds

Simple Summary The milk production traits of goats are economically important. In the Czech Republic, goat milk is processed directly on farms and distributed as cheese. Although goat breeding is not a main focus of animal production in the Czech Republic, it is essential for the agricultural sector. A group of 14 SNPs (single-nucleotide polymorphisms) within four candidate genes (ACACA, BTN1A1, LPL, and SCD) were analysed in two Czech dairy goat breeds, White Shorthaired (WSH) goats and Brown Shorthaired (BSH) goats. The SNPs were significantly associated with milk-production traits (daily milk yield, protein, and fat percentage) and somatic cell count. This information may be useful for marker-assisted selection or related techniques to increase the accuracy of selection. Abstract Milk production is influenced by many factors, including genetic and environmental factors and their interactions. Animal health, especially udder health, is usually evaluated by the number of somatic cells. The present study described the effect of polymorphisms in the ACACA, BTN1A1, LPL, and SCD genes on the daily milk yield, fat, and protein percentages and somatic cell count. In this study, 590 White Shorthaired (WSH) and Brown Shorthaired (BSH) goats were included. SNP genotyping was performed by PCR-RFLP and multiplex PCR followed by SNaPshot minisequencing analysis. The linear mixed model with repeated measurement was used to identify the genetic associations between the studied genes/SNPs and chosen traits. All selected genes were polymorphic in the tested goat populations and showed significant associations with milk traits. Only BTN1A1 (SNP g.599 A > G) showed a significant association with the somatic cell score. After Bonferroni correction, a significant effect of LPL g.300G > A on daily milk yield and fat percentage, LPL g.185G > T on protein percentage, and LPL G50C, SCD EX3_15G > A, and SCD EX3_68A > G on fat percentage was found. The importance of environmental factors, such as the herd-year effect, month of milking, and lactation order on all milk performance indicators was confirmed.


Introduction
Milk and dairy products are essential for human nutrition. Goat milk is rich in minerals, vitamins, and bioactive components, is easily digestible, and contains fewer allergic proteins than cow milk. These characteristics also suggest the possible use of goat milk for therapeutic purposes [1,2]. In Europe, the small ruminant milk industry is not widespread because of the low number of animals and insufficient milk volume in goats compared to cows. However, the number of goats used for milk production is growing due at chromosome 19, consists of 52 exons, the transcription length is 6990 bps, and codes for 2329 AA residues (NCBI gene 100861224). In all, 3555 variant alleles were found in the transcript.
The SCD gene, which is located on goat chromosome 26, has an important role in the cellular biosynthesis of monounsaturated fatty acids (MUFAs), because most of the conjugated linoleic acids are synthesised in the mammary gland by the action of SCD in circularizing vaccenic acid [21,22]. The SCD mRNA has been identified by Bernhard et al., 2001 [23] and Yahyaoui et al., 2003 [24]. There are many SNPs that have been described and identified in exon 3, intron 3, intron 4, exon 5 and 6, and a deletion of a nucleotide triplet in the 3 UTR [23][24][25][26]. The LPL gene is involved in the hydrolysation of triglycerides to glycerol and free fatty acids and in lipoprotein transportation [27,28]. It is synthesised in the mammary gland's epithelial cells and influences the release of fatty acids in the mammary tissue [29]. LPL enzymes are encoded by the LPL gene, which consists of nine exons and eight introns, for a total of 3555 nucleotides. This gene encodes a protein containing 478 amino acids, JQ670882. A few SNPs have been described in goats: a missense mutation responsible for a Ser17Thr amino acid substitution at position 17 of the signal peptide (DQ370053:c.G50C), a C2094T (DQ370053) substitution in the 3 UTR of the gene, and a substitution ss522928251:C > T in intron 7 [21,30].
Our aim was to perform an association analysis of SNPs in the ACACA, BTN1A1, LPL, and SCD genes with milk production traits: daily milk yield (DMY), protein, and fat percentage (PP, FP), and SCC of goats on organic farms in the Czech Republic.

Ethical Approval
The experiment was carried out under Directive 2010/63/EU of the European Parliament and European Council of 22 September 2010, on the protection of animals used for scientific purposes.

Animals
In this study, a total of 590 animals were included. All individuals belonged to the two Czech national goat breeds: White Shorthaired (WSH) and Brown Shorthaired (BSH) goats or their crosses with other breeds. Thus, three breed groups were determined: purebred WSH, purebred BSH, and crossbreds of both breeds, where the proportion of WSH or BSH was 50% or higher. The experiment was performed on two farms that were located in the Czech Republic. Both farms specialize in organic goat milk production, which means that goats grazed and were fed only organic feed; no hormones, antibiotics, or similar substances were applied (except a form of veterinary treatment for individuals), no genetically modified organisms were included, and animals were kept under welfare conditions according to legislation in the European Commission Regulation 889/2008, and the European Commission Regulation 834/2007. The winter feed ration consisted of haylage at approximately 2 kg a day, hay ad libitum, and a grain mix, which was dosed during milking in the milking parlour with a total amount of 300 g a day. The summer feed consisted of grass at approximately 2 kg a day (loaded in the stable), hay ad libitum, and grain mix at 300 g a day. Goats were machine-milked twice a day.

Performance-Testing Database
Phenotype data of genotyped animals were obtained from the performance-testing database of the Czech-Moravian Breeders' Association. In this study, the daily milk yield (DMY) in kg, milk protein (PP) and fat (FP) per cent, and somatic cell count (SCC) were considered. The analysis of DMY, PP, and FP was performed in the group of 590 goats belonging to the White Shorthair breed (n = 490), Brown Shorthair breed (n = 76), and crossbreeds between WSH, BSH, and Saanen (n = 24). Purebred individuals were approximately 75-100% pure. The goats were sired by 37 bucks; the average number of daughters per buck was 8 (minimum 1, maximum 35). They were on the first to eleventh lactation. In all, 8640 milking records from two farms were analysed. At farm A, 2241 repeated-milk records from 279 dairy goats in 2013, 2014, and 2016 were collected. Milk records from farm B included 311 dairy goats with 6399 milk records between 2010 and 2017. Milk samples were collected repeatedly during the milking seasons. There were 2 to 49 repeated records per goat, on average, there were 14 records (approximately four milk controls per year). Milk samples for DMY, PP and FP analysis were collected throughout the whole year as follows: January (n = 571), February (n = 533), March (n = 735), April (n = 1061), May (n = 1161), June (n = 1203), July (n = 1035), August (n = 1069), September (n = 950), October (n = 74), November (n = 57), and December (n = 191).
Samples from only one farm were analysed for somatic cell count, n = 146 goats, that is, White Shorthair goats n = 100, Brown Shorthair goats n = 38, and crossbreeds n = 8. The goats were sired by 27 bucks, and the average number of daughters was 5 (minimum 1, maximum 19). Milk samples were taken during the third to eleventh lactation. The number of repeated records for SCC was 857. Data were collected during 2016 (n = 728) and 2017 (n = 129). The analysis was conducted repeatedly throughout lactation. There were 2-11 repeated records per goat, on average 5 samples per goat, with approximately three controls per year. Milk samples for SCC analysis were collected throughout the milking season as follows: March (n = 26), April (n = 148), May (n = 148), June (n = 128), July (n = 163), August (n = 124), and September (n = 120).

DNA Extraction and SNP Genotyping
Blood samples from 590 animals (5 mL of each) were collected from the jugular vein and preserved in 0.5 mM EDTA (pH 8.0). Genomic DNA was extracted from blood using GeneAll, Exgene TM , and a Clinic SV mini isolation kit (GeneAll Biotechnology cp., Ltd., Seoul, Korea; Bohemia Genetic Ltd., Prague, Czech Republic) according to the manufacturer's recommendations. The SNPs analysed in our study are described in Table 1. SNPs in the BTN1A1 and LPL genes were genotyped according to previously described methods [30][31][32].  Figure S1), AA-Amino Acid.
Individual SNPs of the ACACA and SCD loci were detected by primer extension analysis with the SNaPshot Commercial Kit (Applied Biosystem, Foster City, CA, USA). Primers used for the PCR, extension analysis and electropherogram of the SNaPshot product along with the GeneScan TM 120LIZTM size standard are given in Supplementary  Table S1.
For the PCR reaction (512 bp fragment) of 3 SNPs in the 5 UTR [33], of ACACA locus, we used the set of primers, designed on the basis of the ovine sequence AJ292286 [34]. The PCR assay was performed in a 10 µL reaction mixture, consisting of 2 µL genomic DNA (10-100 ng), 1× PPP Master Mix (Top Bio Ltd., Prague, Czech Republic), and 0.5 µL (10 pmol/µL; 0.01 mM) of each primer (GENERI Biotech Ltd., Prague, Czech Republic), and sterile water up to volume. Thermal cycling conditions are presented in Supplementary  Table S2.
A 536bp fragment of the SCD locus, at region exon3 and intron 3, was amplified by PCR with the following set of primers: F: 5 -TCCTAAgCTTATTCCAgCCCC-3 and R: 5 -gCCAgTCACTCAgAAgTACCC-3 , designed on the basis of the GenBank goat sequence (GenBank AH011188.2; AF422168.1) using Primer 3 software [35]. PCR assay was performed in a 20 µL reaction mixture consisting of 2 µL genomic DNA (10-100 ng), 1× PPP Master Mix (Top Bio Ltd., Prague, Czech Republic), of each primer (GENERI Biotech Ltd., Prague, Czech Republic) and sterile water up to volume. Thermal cycling conditions are presented in Supplementary Table S2.
The presence of fragments obtained in this phase was confirmed by gel-electrophoresis stained with ethidium bromide. The obtained PCR products were purified by using 1 unit of FastAP Thermosensitive Alkaline Phosphatase and Exonuclease I (Fermentas, Ltd., Prague, Czech Republic) to remove unwanted subproducts and incubated at 37 • C for 60 min, followed by 15 min at 85 • C.
The PEA assay utilises internal unlabelled primers which bind to a complementary PCR-generated template in the presence of AmpliTaq DNA Polymerase and fluorescently labelled ddNTPs. The polymerase extends the primer one nucleotide, adding a single ddNTP to its 3 end. Primers were designed to allow size and colour discrimination between the different alleles (Table S1) and were optimised to be used simultaneously.
The single-base extension (SBE) reaction for ACACA locus was performed in a reaction mixture with final volume of 5.0 µL, containing 1.5 µL of purified multiplex PCR product, For electrophoretic detection, 0.5 µL of purified multiplex SBE reaction was mixed with 0.5 µL of GeneScan-120 LIZ size standard (Applied Biosystems, Foster City, CA, USA) and 9.0 µL of Hi-DiTMFormamide (Applied Biosystems, Foster City, CA, USA), following denaturation step at 95 • C for 5 min and analysed by capillary electrophoresis using the Applied Biosystems ® 3130 Genetic Analyzer, an E5-Matrix Standard Set DS-02, a 36 cm capillary, and POP7 polymer (Applied Biosystems, Foster City, CA, USA). The results of genotyping were analysed and evaluated using GeneMapper v 3.5 software (Applied Biosystems, Foster City, CA, USA). The dye colour of the fragment was used to identify the nucleotide of interest ( Figure S1).

Statistical Analysis
The dataset was edited, and unreliable data were removed. SCCs less than 13 and more than 9998 were removed from the analysis. To achieve an approximately normal distribution, the SCC was log-transformed into somatic cell score (SCS). The transformation was performed as follows: where SCC is somatic cell count, which is expressed in thousands per 1 mL of milk.
Hardy-Weinberg equilibrium was tested by the χ 2 test. The effect of gene polymorphisms on milk performance traits and SCS was analysed using the PROC MIXED procedure of SAS with repeated measurements [36]. Tested effects were considered statistically significant at p < 0.05, but biological importance was also considered. The following linear model was used for all traits (DMY, PP, FP, SCS): where Y ijklmno = DMY, FP, PP, SCS; G i = fixed effect of the genotype (class effect i = 1, 2, 3); HY j = combined fixed effect of herd-year (class effect j = 1, . . . , 11); month k = fixed effect of the month of the year of milking (class effect l = 1, . . . , 12); lac l = fixed effect of the lactation order (class effect l = 1, . . . , 9 for DMY, FP and PP or l =1, . . . , 9 for SCS); breed m = fixed effect of the breed (class effect m = 1, . . . , 3); sire n = random effect of the father of the goat; goat o = permanent environment of the goat (repeated measurement); and e ijklmno = random residual effect.
The post hoc comparison was performed by Scheffe s method. A Bonferroni correction for multiple comparisons was applied to all significant associations. The correction factor was derived from the number of SNPs tested. The significance threshold (p < 0.05 and p < 0.01) was divided by the number of tests. Thus, Bonferroni-corrected significance levels of 0.05/13 = 0.004 and 0.01/13 = 0.0008 were applied.

Descriptive Statistics and Phenotypic Correlations
Genotype and allelic frequencies and the number of animals and records included in the analysis are shown in Table 2.  The descriptive statistics of the analysed traits are shown in Table 3. Phenotypic correlations between milk traits (DMY, PP, FP) and SCS are shown in Table 4. The relationships of DMY with PP and FP were negative, which corresponds to the well-known dilution effect reflecting the reductions in fat and protein contents as milk yield increases [37]. Among milk components and SCC, positive correlations were found. Some authors found a negative and significant phenotypic correlation between the logarithm of somatic cell count and milk yield [38]. Milk protein content consistently showed a significant positive correlation to the logarithm of SCC. Their study showed a similar correlation between SCC and milk yield, or milk protein content of dairy goats' milk as found in dairy cows' milk. However, they stated the impossibility of employing commonly used physiological parameters for dairy cows in evaluating the mammary health status of dairy goats. According to other results, the goat milk composition did not change when milk SCC varied among three groups from 214,000 to 1450,000 cells/mL [39]. In ovine milk, the components can be expected to vary independently of milk SCC [40].

Environmental Factors
Environmental factors play an essential role in milk production and mastitis occurrence, so controlling environmental factors could permit or prevent animals from expressing their genetic potential [15]. In our analysis, we considered the following fixed effects: the combined effect of the herd-year (HY), month of milking, lactation order, and breed. All of the mentioned factors were significant for daily milk yield, protein, and fat percentages. The exception was an effect of the breed on DMY, where significance was observed only in models with a few SNPs LPL g.103G > A, g.185G > T, g.257C > T, and g.300G > A. The combined effect of HY explained 7% of the total of 26% variability explained by all tested factors. The HY effect included farm management, milking routines, milking frequency, and feed quality. The Czech Republic is situated in the middle of Europe. The climate is mild with four seasons (spring, summer, autumn, winter). Seasonal kidding usually starts in January in the winter season. After this, goat milk production increases with the growing needs of kids. The month of milking reflects changes during time as well as the nutritional condition of pasture which is rising with the increasing temperature, and also the variation of climate and phase of lactation throughout the year. Milk yield is also negatively correlated with milk fat and protein contents, so with increasing milk production, the milk fat and protein levels decrease [41]. The fixed effect of a month of milking explained approximately 10% of a total of 26% variability. The highest milk production was observed between the lactations. This is in agreement with many authors [41][42][43]. During this period, the highest milk yield is probably caused by the physical appearance, size, and quality of the udder. These morphological characteristics are affected by the breed and genotype of goats as well [44]. Contrary, younger goats tend to have a higher milk fat content than older goats [44]. Milk components were significantly affected by breed; for daily milk yield, the impact of breed varied. Differences between breeds were also confirmed by many authors [42,44].
The somatic cell count could be affected by numerous factors, such as milking routine, stage of lactation, lactation order, breed, feed quality, and health status [3,5,45]. In our analysis of the SCS, only the HY effect was significant, but not the effect of month of milking, lactation order, and breed. HY comprises farm management, feed quality, milking frequency, milking hygiene, pasture quality throughout the year, etc. This analysis was performed on only one farm over multiple years. Thus, the significance of the HY effect probably indicates the difference between the tested years. A trend of increasing somatic cell counts was observed throughout the year (from March to September); the differences among succeeding months were significant or near the significance threshold. Bergonier et al., 2003, claimed that higher rates of mastitis occurrence are observed at the beginning of machine milking and during the first third of lactation, but mastitis is rarely observed during drying-off or parturition [46]. The stage of lactation was not included in this study because of a lack of data. Nevertheless, the mastitis incidence in goats does not vary with the lactation stage, in contrast with cows [46]. Goat milk contains naturally higher SCCs than cow milk due to the apocrine secretory process in goats [47]. There is no consensus on whether the SCC is related to milk production (MY, FP, PP). Several authors claimed that a significantly lower fat content was found in goat milk infected with S. aureus than in noninfected milk. Even an SCC of approximately 3,300,000 cells/mL might not be connected with mastitis or pathological differences in the goat mammary gland [5].
Analysed farms were specialized in organic goat milk production, so the environmental conditions and herd management might differ in feeding routine, health management, and welfare issues compared to conventional goat farms. There is an assumption that there should be differences in SCC between organic and conventional farming. Goats from conventional farming could be exposed to various chemical agents such as disinfections, mycotoxins, pesticides, or residues of antibiotics. These agents affect the development and activity of the microbiological profile of milk. However, many studies on small ruminants and dairy cattle did not confirm differences between organic and conventional farms in milk quality parameters [48,49].

Associations between SNP, Milk Production Traits, and SCS
In total, 14 SNPs were included in the association analysis. SNP g.1255A > G of the ACACA gene was excluded because of monomorphism. The associations between 13 SNPs, milk production traits (MY, FP, PP), and SCS are shown in Table 5, and the differences between genotypes are shown in Table 6.
Polymorphisms other than those in the ACACA gene were associated with fat content in milk. Seven SNPs, namely, BTN1A1 g.599A > G; LPL g.185G > T, g.300G > A, and G50C, and SCD EX3_15G > A; EX3_68A > G, and IVS3+46 C > T, were detected as significant (Table 5). Significant differences between genotypes were found for five SNPs: BTN1A1 g.599A > G; LPL G50C; and SCD EX3_15G > A, EX3_68A > G, and IVS3+46 C > T ( Table 6). After Bonferroni correction, only SNPs in the SCD and LPL genes showed a significant association with fat content in milk. The non-significance of ACACA polymorphisms is surprising, as other authors found an effect of the gene, especially on the milk fat content, but not for milk yield and protein content in goat and sheep milk [33,50], which keeps the gene promising for the future research [51]. The role of LPL in fat synthesis was also stated by other authors [52].
All analysed genes showed an association with protein content. A strong association (p < 0.01) was found between BTN1A1 g.599A > G and LPL g.185G > T. Other SNPs, such as ACACA g.1322T > C and g.257C > T and SCD IVS3+46 C > T, reached significant levels of p < 0.05 (Table 5). Significant differences between genotypes were found for ACACA g.1322T > C, BTN1A1 g.599A > G, LPL 257C > T, and SCD IVS3+46 C > T (Table 6). After Bonferroni correction, the number of statistically significant SNPs decreased from 5 SNPs to just LPL g.185G > T. However, for this SNP, the differences between genotypes were not significant. Only two SNPs were found to be associated with daily milk yield, ACACA g.1322T > C (p < 0.05) and LPL g.300G > A (p < 0.01) ( Table 5). The differences between genotypes showed only slightly significant differences (Table 5). For the SCS, the only association was detected with the SNP BTN1A1 g.599A > G (p < 0.05). In this polymorphism, genotype AA was connected with the highest somatic cell score in milk, and significant differences were found between the AG and GG genotypes (p < 0.05). Unfortunately, this association was not significant after Bonferroni correction. No other significant associations were found for the SCS.
When comparing the effect of SNPs retained after Bonferroni correction and the differences among genotypes, only the LPL G50C and SCD EX_15G > A and EX3_68A > G polymorphisms showed a significant effect at the SNP and genotype levels. Additionally, other gene polymorphisms in goats are studied. For example, the AGPAT6 gene was found to significantly influence both the fat and protein contents and milk yield [53]. Additionally, the variability of the haplotypic sequences at the loci of the casein complex was studied [15]. The authors point out that a complete definition of the haplotypes in the casein complex in goats is difficult given the high genetic variability.
LPL, ACACA, and SCD ovine genes expression were found to be influenceable by diets [54]. The ability of nutrigenomic regulation of the transcription confirmed that these genes play a critical role in the regulation of lipid metabolism processes in sheep and could be associated with fatty acid profiles in milk and meat. Ovine LPL gene should also be studied due to its expression to microRNA-148a [55]. SCD and BTN1A1 genes were not found to be differentially expressed when comparing two Spanish sheep breeds [56].

Conclusions
Concisely, some SNPs in the included genes showed association with milk traits but not with SCS. After more accurate Bonferroni correction, the significant associations of SNPs were only rare. Thus, our results support the generally accepted fact that environmental factors are more important than genetic for milk-production traits. However, along with population genetic analyses, the study of major genes in goats helps to better understand the genetic background of the milk-performance complex. This may contribute to future more effective selection in dairy goat populations.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/ 10.3390/ani11061796/s1, Table S1: Primers used for the amplification PCR product and extension primers of the goat Stearoyl-CoA desaturase, Acetyl-CoA carboxylase, cDNA, and annealing temperature, Table S2: Thermal cycling conditions of Stearoyl-CoA desaturase and Acetyl-CoA carboxylase. Figure S1: Electropherogram of SNPs in the Stearoyl-coenzyme A desaturase and Acetyl-CoA carboxylase α gene analysed with Gene Mapper software.