Novel Genes Associated with Dairy Traits in Sarda Sheep

Simple Summary In a population of 1112 individuals of Sarda sheep, we investigated the polymorphism of 45 single nucleotide polymorphisms (SNPs) in 15 different genes that have been previously investigated because of their relation to metabolism and innate immunity. For SNPs with a sufficient allele frequency, we also tested their influence on a wide panel of dairy traits. The genotyping evidenced some new genetic patterns, and seven SNPs had significant influences on the milk yield, composition and coagulation traits. The results may represent a starting point to acquire new knowledge in the field of dairy sheep sector and to improve the breeding schemes for the Sarda sheep breed. Abstract The aim of the present research was to analyze the variability of 45 SNPs from different genes involved in metabolism and innate immunity to perform an association analysis with the milk yield, composition and milk coagulation traits. A population of 1112 Sarda breed sheep was sampled. Genotyping was generated by a TaqMan Open ArrayTM. Thirty out of the 45 SNPs were polymorphic, and 12 displayed a minor allele frequency higher than 0.05. An association analysis showed that the variability at genes PRKAG3 and CD14 was significantly associated with the daily milk yield. The variability at PRKAG3 was also associated with the protein and casein content, somatic cell score and bacterial score. The variation at the PRKAA2 gene was associated with the milk lactose concentration. The SNPs at CD14 were also associated with the traditional milk coagulation properties, while the SNPs at GHR and GHRHR were associated with kSR, a derived coagulation parameter related to the rate of syneresis. The information provided here is new and increases our knowledge of genotype–phenotype interactions in sheep. Our findings might be useful in appropriate breeding schemes to be set up for the Sarda sheep breed, but these should be confirmed by further studies, possibly performed on independent populations.


Introduction
Sheep are farmed all over the world, because of their great adaptability and relatively high productions. Productive specialization is different among countries and geographical areas. In the Middle East and in the Mediterranean basin, sheep milk is of major importance [1]. In the European Union, sheep milk production has been markedly increased in the last decades, extending from 2 to about 3 million t in the period of 1960-2019, and dairy sheep farming is a relevant sector of agriculture in many areas, such as for Greece, the leading EU country with 944 thousand t of ewe milk produced in 2019, followed by Romania, Spain and Italy [2].
In dairy sheep farming, almost all milk is transformed into cheese; therefore, milk traits related to cheese-making, including milk coagulation properties (MCP), curd firming and syneresis, need to be evaluated, as they might impact the farm income [3]. The MCP are described as the rennet coagulation time (RCT), curd firming time (k 20 ) and curd firmness hormone system plays an important role in protein and lipid metabolism and affects the galactopoiesis in dairy species [23]. We investigated the polymorphism at the B4GALT1 (β1,4-galactosyltransferase 1) and ABCG2 (ATP-binding cassette subfamily G member 2) genes. The β1,4-galactosyltransferase 1 protein is involved in glucose metabolism, as it forms a heterodimer with alpha lactalbumin in mammary epithelial cells. That generates the lactose synthase enzyme, which synthesizes lactose by combining one molecule of UDP-galactose with one of glucose [24]. The ABCG2 protein belongs to a transporter family, including the proteins involved in the transport of molecules across cell membranes [25]. The ABCG2 gene is expressed in different tissues, including the mammary gland epithelial cell during lactation, and plays an important role in active secretion of molecules across the blood-milk barrier, including drugs and toxins [26].
As regards the genes involved in the immune response and defense against pathogens, CD14 (Cluster of Differentiation 14 Molecule), TLR2 (Toll-Like Receptor 2), TLR4 (Toll-Like Receptor 4), and PI (Serpin Peptidase Inhibitor), whose variability in dairy sheep and correlation with milk traits has not been investigated yet, the CD14 gene encodes a molecule that is fundamental for innate immunity and acts against a range of pathogens. The CD14 protein is a surface antigen that is preferentially expressed on monocytes and macrophages, which contributes to mediate the response to bacterial lipopolysaccharide [16]. The TLR4 gene encodes a protein, a member of the Toll-like receptor family, involved in pathogen recognition and innate immunity. As the TLR4 receptor is involved in the cellular response to the lipopolysaccharides found in most Gram-negative bacteria, it also plays a role in mastitis [27]. We also evaluated the variability at the Neuropeptide FF receptor 2 (NPFFR2) gene. NPFFR2 is a member of the G protein-coupled receptor superfamily, membrane proteins activated by neuropeptides NPAF and NPFF [28]. The NPFF neuropeptides were shown to display anti-inflammatory activity in a mouse model and were linked to T-lymphocyte proliferation and the immune response [29].
In addition to the above-mentioned metabolism and immune response-related genes, we investigated the variability of the PLCE1 and TYRP1 genes. The PLCE1 (Phospholipase C Epsilon 1) gene encodes a phospholipase enzyme, which catalyzes the hydrolysis of phosphatidylinositol-4,5-bisphosphate, and produces the second messengers inositol 1,4,5triphosphate (IP3) and diacylglycerol (DAG). The two second messengers can affect cell growth, differentiation, and gene expression [16]. The TYRP1 gene (Tyrosinase Related Protein 1) plays an important role in the synthesis of eumelanin and is associated with the coat color in domestic animals, including sheep [30].
The present study aimed to explore, in a population of 1112 dairy ewes of the Sarda breed, the variability of 45 SNPs located on fifteen genes responsible for metabolism and immune defense reactions in order to perform an association analysis between the genotyped SNPs and a wide range of dairy traits, including milk yield and composition, as well as traditional and model coagulation traits.

Farms, Animals, and Samples
For the present study, 1112 ewes in 23 different farms were sampled. Detailed features of the animals and farms were reported by Pazzola et al. [5]. The main aspects were the following: ewes were of Sarda breed, regularly recorded in the flock book, and daughters of 120 different sires (minimum 4, maximum 40 daughters per sire); parity order and stage of lactation of ewes were, respectively, between the first and ninth, and from the second to the seventh month after parturition; farms were of commercial type, located in the regional territory of Sardinia (Italy) and managed under semi-extensive farming systems.
At each farm, a variable group of ewes were sampled (minimum 22, maximum 89 ewes per farm) in a single sampling day. Individual milk samples were collected in 200-mL sterile plastic containers during the evening milking, which was performed by using manually operated milking machines. The daily milk yield (dMY) was calculated as the sum of the morning and evening milking of the same sampling day. Blood samples were individually collected in K 3 EDTA vacuum tubes (BD Vacutainer, Becton Dickinson, Milan, Italy) from the jugular vein. Milk and blood samples were transported at 4 • C to the laboratories to be stored and analyzed.

DNA Extracton and Genotyping
Blood samples were processed for DNA extraction within 24 h after sampling. Genomic DNA was extracted using the Gentra Puregene blood kit (Qiagen, Hilden, Germany). DNA concentration and purity were assessed using an Eppendorf BioPhotometer (Eppendorf, Hamburg, Germany).
We choose each SNP to be analyzed in the present experiment, respecting the technical requirements of the TaqMan 25 and a further analysis of the individual genotypes using SVS software version 7 (Golden Helix Inc., Bozeman, MT, USA). Samples or SNPs with call rates < 0.9 were removed from the dataset. The minor allele frequencies (MAF) observed and expected heterozygosity, consistency of the genotype distributions with the Hardy-Weinberg equilibrium, definition of the linkage disequilibrium (LD), and LD blocks were estimated by the Haploview software package [31].

Milk Analysis
Analysis of the milk samples was performed within 24 h after collection, as reported by Dettori et al. [32]. As regards the fat, protein, casein, lactose content, and pH, these were achieved by using the standard of the International Organization for Standardization and International Dairy Federation standard [33] and the MilkoScan FT6000 instrument (Foss Electric A/S, Hillerød, Denmark); the total bacterial count (TBC) was measured according to the ISO-IDF standard [34] by using a BactoScan FC150 instrument (Foss Electric); the somatic cell count (SCC) was measured according to the ISO-IDF method [35] by using a Fossomatic 5000 instrument (Foss Electric). The TBC and SCC were transformed into respective logarithmic scores to normalize the distribution: LBC (log-bacterial count = log 10 (total bacterial count/1000), according to the ISO-IDF standard [36])) and SCS (somatic cell score = log 2 (SCC × 10 −5 ) + 3, according to Shook [37]).
As regards the coagulation traits, both the traditional milk coagulation properties (MCP) and curd firmness over time traits (CF t ) were measured. The following MCPs were measured according to the method by McMahon and Brown [4] and the revisions for sheep species by Pazzola et al. [5] by using a Formagraph instrument (Foss Italia, Padova, Italy): RCT (rennet coagulation time in min); k 20 (curd firming time in min); and a 30 , a 45 , and a 60 (curd firmness 30, 45, and 60 min after rennet addition in mm). The following CF t traits were measured according to a method used for ovine species by Vacca et al. [6] by using the curd-firming individual point observations recorded by a Formagraph instrument: CF P (the maximum potential curd firmness at an infinite time in mm), k CF (curd-firming rate constant in % × min −1 ), k SR (syneresis rate constant in % × min −1 ), CF max (maximum curd firmness in mm), and t max (time to attain CFmax in min).

Statistical Analysis
As regards the editing of the MCP and CF t traits from the 1112 milk samples, the values of RCT higher than 60 min, also called noncoagulating samples, were labeled as missing (n = 26); k 20 was labeled as missing when higher than 5 min and for noncoagulating samples (n = 36); curd firmness was labeled as missing at 0 mm (a 30 : n = 29, a 45 : n = 24, a 60 : n = 26, CF P : n = 19, and CF max : n = 19).
The association analysis between the genotypes of the polymorphic SNPs with the milk composition, MCP, and CF t traits was achieved using the MIXED procedure of SAS (version 9.4, SAS Inst. Inc., Cary, NC, USA), which was derived from our previous study by Dettori et al. [32]. The following model (1) was used to investigate one milk trait for each SNP at a time only for SNPs with a MAF higher than 0.05 and genotypes with a frequency higher than 0.01: where Y ijklm is the observed traits of the milk composition, MCP, and CF t ; µ is the general mean; G i is the fixed effect of the ith SNP genotype (i = 2 to 3 levels: the two homozygotes and the heterozygote); F j is the fixed effect of the jth farm, which also includes animal management and feeding (j = 1-23 levels, i.e., the different farms where animals were reared); P k is the fixed effect of the kth parity of the ewes (k = 4 levels: first, 455 sheep; second, 215 sheep; third, 212 sheep; and fourth or more parities, 230 sheep); DIM l is the fixed effect of the lth days in milking (l = 4 levels: level 1, ≤100 days, 210 sheep; level 2, 101-140 days, 381 sheep; level 3, 141-160 days, 133 sheep; and level 4, ≥161 days, 388 sheep); SIRE(G) m is the random effect of the mth sire (m = 124 different sires) nested within the genotype; and e ijklm is the error random residual effect. The association analysis between the milk traits and the LD block was performed by a further model (2) derived from model (1), with LD i (i = 2 to 3 levels) instead of G i and SIRE(LD) m instead of SIRE(G) m , analyzing one milk trait for each LD block at a time.
The effects of the models were declared significant at p < 0.05, and multiple testing, for both models (1) and (2) and for the factors with more than two levels was performed using the Bonferroni method at α = 0.05.

Allele Frequencies and Linkage Disequilibrium
The results of 45 SNPs' genotyping are reported in Table 1. Thirty (33.3%) SNPs were polymorphic and, among these, 12 showed a MAF higher than 0.05 and were submitted to the association analysis. A linkage disequilibrium (LD) analysis was performed for all the SNPs mapping to the same chromosome. LD was highlighted in chromosome OAR:16 ( Figure 1) within the GHR gene. The SNPs analyzed ranged from the 5 UTR region to exon 10, the last exon of the GHR gene. One LD block was located at exon 10 and included two htSNPs: rs408890407, synonymous mutation, and rs55631463, missense mutation. Three haplotypes were recorded with CT showing the highest frequency, 0.493, while CC and TT had 0.277 and 0.230, respectively.

Association Analysis
The descriptive statistics of the milk composition and coagulation traits achieved from the 1112 sampled ewes is reported in Table 2. The fixed effects included in models (1) and (2), i.e., farm, parity, and stage of lactation, showed high levels of significance, at p < 0.001 and p < 0.01, for almost all the milk traits (data not shown in the tables). The results about these effects have been widely discussed in previous papers by the same authors [5,13] using datasets connected to the present one. The proportion of variance explained by the random effect of the sire computed in both the models was always lower than 10% for almost all the SNPs and the LD block. The lowest values were recorded for the CFp, kCF, kSR, Tmax, and LBC, at about 0%, and the highest for the fat, protein, casein, and lactose, spanning from 5% to 9% (data not shown in the tables).
Supplementary Table S2 displays the results of the statistical analysis computed by model (1) regarding the genotype effects of the 12 polymorphic SNPs, with the MAF > 0.05, on the milk traits. The least square means of the milk yield and composition, MCP, and CFt, according to the genotypes with a frequency > 0.01, are reported in Table 3. SNP rs119102735 at PRKAA2 was the only investigated polymorphic locus affecting the lactose

Association Analysis
The descriptive statistics of the milk composition and coagulation traits achieved from the 1112 sampled ewes is reported in Table 2. The fixed effects included in models (1) and (2), i.e., farm, parity, and stage of lactation, showed high levels of significance, at p < 0.001 and p < 0.01, for almost all the milk traits (data not shown in the tables). The results about these effects have been widely discussed in previous papers by the same authors [5,13] using datasets connected to the present one. The proportion of variance explained by the random effect of the sire computed in both the models was always lower than 10% for almost all the SNPs and the LD block. The lowest values were recorded for the CF p , k CF , k SR , T max , and LBC, at about 0%, and the highest for the fat, protein, casein, and lactose, spanning from 5% to 9% (data not shown in the tables). Supplementary Table S2 displays the results of the statistical analysis computed by model (1) regarding the genotype effects of the 12 polymorphic SNPs, with the MAF > 0.05, on the milk traits. The least square means of the milk yield and composition, MCP, and CF t , according to the genotypes with a frequency > 0.01, are reported in Table 3. SNP rs119102735 at PRKAA2 was the only investigated polymorphic locus affecting the lactose content, and the milk samples from ewes carrying genotype TT were characterized by a lower concentration than CC and CT. The SNP with the most conspicuous effects on the milk traits was rs159573167 at PRKAG3, with AA homozygous ewes producing the highest daily milk yield, lower concentrations of protein and casein, the lowest value of SCS and LBC, and the shortest rennet coagulation time. SNP rs160202315 at the TLR4 (Toll-Like Receptor 4) gene was associated with the milk composition and cell traits, milk samples from heterozygous TC ewes characterized by the best contents, i.e., higher, of protein and casein, and the highest scores of SCS and LBC. The effects of rs409504706 at GHRHR and rs55631463 at GHR were directed on the k SR trait, with TT and CC ewes, respectively, showing the fastest syneresis process. SNP rs428862267 at GHR affected the curd firmness, and AA ewes produced milk with the highest a 30 . SNP rs160087383 at CD14 had a highly significant effect on the milk yield (p < 0.001) and coagulation, and CC homozygous ewes showed both higher daily yields and better coagulation traits than TT, i.e., faster rennet coagulation times and larger curd firmness (p < 0.05): dMY 2012 vs. 1654 g/day, RCT 7.60 vs. 9.28 min, and a 30 52.52 vs. 48.03 mm, respectively. Table 3. Least square means of the milk yield and composition, the milk coagulation properties (MCP), and curd firmness over time traits (CF t ), according to the genotypes of the SNPs investigated in the Sarda sheep (n = 1112).  The results of the statistical analysis computed by model (2) regarding the genotype effect of the haplotype blocks at chromosome 16 on the milk traits are reported in Supplementary Table S3. No significant effect on the investigated milk traits was found for the LD block at GHR (SNPs rs408890407 and rs55631463).

Genes Related to Sheep Metabolism
Important metabolic variations occur during pregnancy and lactation, which involve nutrient partitioning in order to guarantee productive functions while maintaining the animal's well-being [38]. In the present study, we evaluated the variability of genes responsible for metabolic activities and their association with milk traits in a population of Sarda sheep.
AMPK, or 5 AMP-activated protein kinase, is an enzyme (EC 2.7.11.31) formed by three protein subunits (α, β, and γ) conserved from yeast to humans. The AMPK enzyme helps to maintain the energy homeostasis of the cell, activating glucose and fatty acid uptake and their oxidation when the AMP:ATP ratio varies [39]. The alpha subunit is encoded by the PRKAA2 gene and the gamma subunit by the PRKAG3 gene. In the present study, we highlighted that SNP rs119102735 of the PRKAA2 gene was associated with the milk lactose concentration. Based on the available literature, this is the first study to report a significant association between the polymorphisms at this gene and a trait regarding milk composition in dairy sheep. In accordance with our study, a recent finding by Huang et al. [40] elucidated that, in bovine mammary epithelial cells, the activation of AMP-activated protein kinase systems is linked with the reduction of lactose synthesis. The AMPK protein, composed of three subunits in humans, is expressed in different tissues, mainly in the skeletal muscle, heart, liver, and hypothalamus [41]. Hao et al. [42] did not reveal it by RNA-seq among the mRNAs expressed in the lactating mammary gland of sheep. Bionaz-Loor et al. [43] revealed, by using Real-Time PCR, that PRKAA2 is expressed in the mammary glands of lactating cattle, and the same authors did not reveal variations in the expression of this gene during the analyzed period. They observed that, although high-protein synthesis persisted during lactation, the level of mRNA synthesis for the AMPK subunits, or other genes corresponding to related biosynthetic pathways, did not vary. In sheep, Crisà et al. [44] analyzed a polymorphism in exon 4 of PRKAA2, which was not associated with changes in the synthesis of fatty acids in milk. AMPK is a molecule that functions as a metabolic switch, which, in the absence of ATP, blocks consumption and activates the catabolic pathways to produce ATP, and turns off the anabolic pathways, which consume ATP [45]. Since AMPK is involved in linking energy sensing to the metabolic pathways, it plays an important role at the whole-body level by translating endocrine communications into adapted metabolic responses [46]. For this reason, it is possible that the observed association between PRKAA2 and lactose might not to be directly linked to a variation of the AMPK expression or activity in the secretory cell of the mammary gland but, rather, to the regulation of the animal's metabolism in the physiological phase of lactation.
We also genotyped three SNPs of the PRKAG3 gene and evidenced that SNP rs159573167 AA was associated with an increase of dMY in the sampled population (+17% in comparison with genotype GG; data not shown in the tables) and a decrease of the protein and casein contents. The variability at the PRKAG3 gene was investigated in New Zealand Suffolk sheep [47], and the variability and expression profile in skeletal muscle were investigated in five local Spanish breeds, where PRKAG3 was found to be highly expressed [48]. Here, for the first time, we highlighted an association between PRKGA3 and dMY and the protein and casein contents. Milk from PRKAG3 rs159573167 AA ewes also showed the lowest SCS and LBC scores, which was not directly explained by its function, and the rennet coagulation time (RCT). In the present study, the occurrence of lower concentrations of protein, somatic, and bacterial cells and higher daily yields found in PRKAG3 rs159573167AA ewes might be attributable, in addition to the genetic effect, to the so-called dilution effect, as the higher the milk yield, the lower the level of total solids and cells and vice versa [1].
We also genotyped twelve SNPs from genes IGF1, GHRHR, and GHR, which are part of the growth hormone pathway and play an important role in lipid, protein, and glucose metabolism [27]. We highlighted an association with curd firmness a 30 (GHR rs428862267) and the syneresis instant rate constant, k SR (GHR rs55631463 and GHRHR rs409504706). Previous studies, carried out in a subset of the resource sheep population, have revealed an association of several SNPs at the GHR gene with the daily fat and protein yields [49] and with the k SR [50]. The k SR is a derived parameter depending on the speed of syneresis, which, in turn, is related to the lipid content of the milk [7].

Genes Related to Innate Immunity
The role of specific genes in the innate immune response to some pathogens in ruminants has been investigated [51], due to its potential positive impact in sustainable agriculture systems. Here, we evaluated the association of some genes involved in innate immunity with milk traits in dairy Sarda sheep. Genotype CC at SNP rs160087383 of the CD14 gene was associated with an increase of dMY. Li et al. [52] investigated the role of the CD14 gene in relation to the Gram-negative bacteria-induced innate immune response and revealed that the variability of this gene is associated with the morbidity of Chinese Holsteins. Thus, the CD14 gene, which is expressed in macrophages and plays an important role in a number of pathways related to the immune response, in our set of animals was found to be associated with the daily milk yield. The CD14 rs160087383 CC ewes were characterized by the highest daily milk yield and the concurrence of favorable coagulation traits, i.e., a faster RCT and a larger curd firmness a 30 . As observed also for the PRKAG3 gene, the milk samples from high-yielding ewes were characterized by the shortest, and consequently better, values of the rennet coagulation time. This positive feature is in accordance with a study using a dataset interconnected with the present one [5] that has evidenced this peculiar result regarding sheep species, because cow milk yields and RCT commonly show a positive (i.e., unfavorable) correlation [53].
In our investigation, SNP rs160202315 at the TLR4 gene was associated with both the protein and cellular traits. The significant effect of the TLR4 gene on the protein and casein contents is analogous to that recorded for cattle species [54], while effects on milk coagulation have been found for TLR2, encoding for a protein of the same family, which influences the curd-firming time in dairy cows [55]. Milk samples from TLR4 rs160202315 CC homozygous ewes were characterized by the lowest scores of SCS and LBC; this association can be explained in relation to the gene function, pathogen recognition, and innate immunity, as described above.
We have to keep in mind that the observed associations between the genotyped SNPs and dairy traits may be due to linkages with unknown causal variations.

Conclusions
The designed open array revealed some new genetic patterns and the related influences on the milk yield, composition, and coagulation traits. The results obtained in the present study from the association analysis evidenced several significant effects of the investigated SNPs on the milk traits, and these may represent an advancement for the dairy sheep sector. In particular, for the Sarda sheep breed, this information might be useful in appropriate breeding schemes, but these should be confirmed by further studies, possibly performed on independent populations. Supplementary Materials: The following are available online at https://www.mdpi.com/article/ 10.3390/ani11082207/s1: Table S1: Gene names, SNP ID, and context sequences of the 45 SNPs genotyped in the sampled population of Sarda sheep (n = 1112). Table S2: p-value and significance for the milk yield and composition, milk coagulation properties (MCP), and curd firmness over time traits (CF t ), according to the effects of 12 polymorphic SNPs out of 45 investigated in Sarda sheep (n = 1112). Table S3: p-value and significance for the milk yield and composition, milk coagulation properties (MCP), and curd firmness over time traits (CF t ), according to the effects of the block evidenced among the 45 SNPs investigated in Sarda sheep (n = 1112).  Institutional Review Board Statement: Ethical review and approval were waived for this study. No specific authorization from an animal ethics committee was required. DNA was isolated from blood samples that were collected by official veterinarians of the local health authorities (ASL) during scheduled control or eradication sanitary programs not linked with the present study. The milk samples were collected concurrently with official and performance controls of the flock book not directly linked with the present study. The sheep belonged to commercial private flocks and were included in this trial with the agreement of the farmers on a voluntary basis.
Data Availability Statement: Not applicable.