Association Analysis between SPP1, POFUT1 and PRLR Gene Variation and Milk Yield, Composition and Coagulation Traits in Sarda Sheep

Simple Summary The purpose of our research was to analyze the association between the three candidate genes secreted phosphoprotein 1 (SPP1), protein O-fucosyltransferase 1 (POFUT1) and prolactin receptor (PRLR) with milk production, quality and coagulation properties in 380 Sarda dairy sheep. Results revealed an association between SPP1 and somatic cells count, in line with the function of this gene and with its genomic position. We revealed an association of POFUT1 variation with milk coagulation properties, and PRLR with quality. This information can be useful for future breeding schemes in sheep. Abstract Many studies focus on the identification of genomic regions that undergo selective processes, where evidence of selection is revealed and positional candidate genes are identified. The aim of the research was to evaluate the association between positional candidate genes, namely secreted phosphoprotein 1 (SPP1, sheep chromosome Ovis aries OAR6, 36.651–36.658 Mb), protein O-fucosyltransferase 1 (POFUT1, OAR13, 61.006–61.027 Mb) and prolactin receptor (PRLR, OAR16, 38.969–39.028 Mb) with milk yield, composition and coagulation traits. Eight single nucleotide polymorphisms (SNPs) mapping to the three genes were genotyped in 380 Sarda dairy sheep. Statistical analysis revealed an association between SNP rs161844011 at SPP1 (chromosome position Oar_v3 OAR6:36651870, gene region exon 7) and somatic cell score, while POFUT1 SNP rs424501869 (OAR13:61007495, intron 1) was associated with curd firmness both 45 and 60 min after rennet addition (p = 0.015 and p = 0.007, respectively). SNP rs400874750 at PRLR gene (OAR16:39004070, intron 2) had a significant association with lactose content (p = 0.020), somatic cell score (p = 0.038), rennet coagulation time (p = 0.018) and curd firming time (p = 0.047). The outcome of this research confirmed predictions based on genomic studies, producing new information regarding the SPP1, POFUT1 and PRLR genes, which may be useful for future breeding schemes.


Introduction
In the last decades, the increase of sheep milk and cheese production has been recorded worldwide [1]. The enhancement of cheese and other dairy products is often pursued by official labelling, such as in European Union countries where protected designation of origin (PDO) and protected geographical indication (PGI) marks have a significant positive impact on the economic value of products [2]. Nevertheless, the starting point for the improvement of dairy products is based on the study of milk quality [3]. Phenotypic traits of milk, i.e., fat and protein content, have a key role in influencing cheese yield and the quality of dairy products, and several studies have also Animals 2020, 10, 1216 3 of 13 veterinarians, concurrently with official performance controls of the flock book, which were not directly linked with the present trial. Sheep were included in this trial with the agreement of the farmers on a voluntary basis.

Animals, Farms and Sampling
A population of 380 ewes from 19 farms (20 ewes at each farm) located in Sardinia (Italy) was used for the study. A detailed description of the farms, animals and sampling procedures are reported in Pazzola et al. [25] and Vacca et al. [26]. In brief, all the ewes and farms were officially registered in the flock book of the Sarda sheep breed; farms were mainly managed with semi-extensive methods; multiparous sheep sampled for the present study lambed in November-December, whereas primiparous sheep lambed in the successive February-March, as normally happens in semi-extensive sheep farming in Sardinia. Ewes were between the first and ninth parity, and two to seven months after lambing. Individual milk samples were taken at each farm in a single day (one day sampling for each farm) throughout a single lactation period from spring to summer in 2012. During the afternoon milking, milk samples were collected in 200 mL sterile plastic containers and kept at 4 • C until analyses. Daily milk yield was recorded as the morning plus evening milking of the same sampling date.
A blood sample was taken from each ewe in K 3 EDTA vacuum tubes (BD Vacutainer, Becton Dickinson, Franklin Lakes, NJ, USA) for DNA extraction.

Milk Analysis
Milk samples were analyzed within 24 h after collection for composition and milk coagulation properties. Fat, protein, casein, lactose content and pH were measured using a MilkoScan FT6000 device (Foss Electric A/S, Hillerød, Denmark), according to the International Organization for Standardization and International Dairy Federation (ISO-IDF) standard [27]. Daily fat and protein yield in g/day (dFPY) was calculated as the sum of fat protein content multiplied by daily milk yield. Total bacterial count (TBC) was measured using a BactoScan FC150 instrument (Foss Electric) according to the ISO-IDF standard [28] and later, in order to normalize the distribution, transformed into the log-bacterial count (LBC = log 10 (total bacterial count/1000)), according to the ISO-IDF standard [29]. Somatic cell count (SCC) was measured using a Fossomatic 5000 equipment (Foss Electric) according to the ISO-IDF method [30] and transformed into the somatic cell score (SCS = log 2 (SCC × 10 −5 ) + 3) according to Shook [31]. Milk coagulation properties (MCPs) were measured using a Formagraph instrument (Foss Italia, Padova, Italy) and the method first described in McMahon and Brown [32] and modification reported for the sheep species in Pazzola et al. [25]. Briefly, a volume of 10 mL of each milk sample was heated to 35 • C and mixed with rennet enzyme (200 µL of a solution with final milk clotting units (IMCUs) per mL of 0.0513/mL of milk, obtained with the dilution in distilled water 1.2% (wt/vol) of Hansen Naturen Plus 215 (Pacovis Amrein AG, Bern, Switzerland; 80 ± 5% chymosin and 20 ± 5% pepsin; 215 international milk clotting units per mL, IMCU/mL)) and analyzed for 60 min after rennet addition. The following MCPs were measured: 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). Curd firmness over time (CF t ) traits were calculated using the records of curd firmness downloaded from the Formagraph and the method first described in Bittante [33] and modification reported for the sheep species in Vacca et al. [26]. The following CF t were measured: RCT eq (rennet coagulation time estimated by the CF t equation, in min); 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 CF max , in min).

DNA and Haplotype Analyses
Genomic DNA was extracted using the Gentra Puregene blood kit (Qiagen, Hilden, Germany), and purity and concentration were measured with an Eppendorf BioPhotometer (Eppendorf, Hamburg, Germany). A custom open array, based on the TaqMan real-time PCR assay, was designed for single Animals 2020, 10, 1216 4 of 13 nucleotide polymorphism (SNP) genotyping. It included a total of 8 SNPs: 2 from the SPP1 gene, 3 from POFUT1 and 3 from PRLR (Table 1). Context sequences are given in Supplementary Information Table S1. Genotyping was carried out using a 12 K Flex QuantStudio instrument (Thermo Fisher Scientific, Waltham, MA, USA). Genotypes were visualized with Taqman Genotyper v.1.3 software (Applied Biosystems, Waltham, MA, USA). The Haploview software package [34] was used to estimate and plot pairwise linkage disequilibrium (LD) measures (D and r2) and to infer haplotype frequencies as well as to define LD blocks according to the Gabriel criteria [35]. Haploview was also used to estimate minor allele frequencies (MAF) and observed and expected heterozygosites and to identify significant departures from the Hardy-Weinberg equilibrium at each polymorphic locus.

Statistical Analyses
Association analysis between the genotypes of the polymorphic SNPs at SPP1, POFUT1 and PRLR, with milk composition, MCP and CF t traits was performed using the MIXED procedure of SAS (version 9.4, SAS Inst. Inc., Cary, NC, USA) and the following model (1): where Y ijklm is the observed trait of milk composition, MCP and CF t ; µ is the general mean; G i is the fixed effect of the i th SNP genotype (i = 2 to 3 levels: the two homozygotes and the heterozygote); F j is the fixed effect of the j th farm, which also includes animal management and feeding (j = 1 to 19 levels; i.e., the different farms where animals were reared); P k is the fixed effect of k th parity of the ewes (k = 1 to 4 levels; first to fourth or more parities); DIM l is the fixed effect of the l th days in milking (l = 4 levels; level 1: ≤100 days; 2: 101-140 days; 3: 141-160 days; level 4: ≥161 days); SIRE(G) m is the random effect of the m th sire (m = 108 different sires) nested within the genotype; and e ijklm is the error random residual effect. We analyzed one milk trait for each SNP at a time. We only considered SNPs with MAF > 0.05. In order to investigate the association between milk traits and each of the LD blocks, the same model (1) was slightly modified into model (2), with LD i (i = 3 levels) instead of G i , one milk trait for each LD block at a time. Correction for multiple testing for both models (1) and (2) were performed using the Bonferroni method at α = 0.05.

Results
Descriptive statistics regarding all the milk traits recorded from the sampled population of 380 Sarda ewes are summarized in Supplementary Information Table S2. Values of curd firmness at 0 mm, both for MCP and CF t , were labelled as missing and were not computed in the statistical analysis (a 30 : n = 3; a 45 : n = 2; a 60 : n = 3; CF P : n = 7; CF max : n = 4).

Allele Frequencies at SPP1, POFUT1 and PRLR, and Association Analysis with Milk Traits
The SNPs analyzed in the present experiment were chosen taking into account the technical requirements of a TaqMan ® assay, and by using the Ensembl Genome Browser (https://www.ensembl. org/index.html) to localize and choose each SNP. For the SPP1 gene we analyzed one coding SNP (rs161844011) causing the amino acid substitution p.Gln235Arg, and one exonic SNP (rs426249393) in the 5' untranslated region (UTR); the other SNPs analyzed were intronic. Among the eight SNPs genotyped, only rs421284407 (POFUT1 intron 2) was monomorphic; all the others were polymorphic, with minor allele frequencies (MAF) higher than 0.1. The population analyzed was in Hardy-Weinberg equilibrium at all loci ( Table 1).
Results of the statistical analysis (F-values and significance) from model (1), performed to investigate the influence between each of the seven polymorphic SNPs and milk traits, are reported in Supplementary Information Table S3. The results regarding the fixed effects of farm, parity and stage of lactation included in model (1) showed high levels of significance for almost all the milk traits (data not reported in tables), as already discussed in previous papers using datasets connected to that of the present study [7,25].
In the SPP1 gene, the rs161844011 SNP marker was associated with SCS (p = 0.003), and ewes carrying the homozygote TT genotype displayed the highest value ( Figure 1a). The SPP1 SNP rs426249393 was associated with t max (p = 0.046), and homozygous ewes AA and GG showed delayed times to attain the maximum curd firmness, about six to seven minutes later than heterozygous AG (Figure 1b). The SNPs investigated at POFUT1 were associated with curd firmness traits. Milk samples from rs424501869 AA ewes were characterized by the highest value of a 45 (48.10 min; p = 0.015) and a 60 (45.15 min; p = 0.007) (Figure 1c,d), whereas milk samples from homozygous rs408068827 AA ewes were associated with the longest time to attain the maximum curd firmness, 38.20 min (p = 0.028; Figure 1e).
The SNP rs400874750 at PRLR gene was significantly associated with milk composition and coagulation times: milk from CC homozygous ewes displayed the lowest lactose concentration at 4.66 mg/100 mL (p = 0.020; Figure 1f) and the highest SCS at 5.69 (p = 0.038; Figure 1g). In addition, milk produced by rs400874750 CC ewes had delayed times both for RCT (p = 0.018) of about 2.5 min of delay in comparison with TC and TT, and k 20 (p = 0.047), about 0.30 min delay (Figure 1h,i).

Linkage Disequilibrium and Association Analysis between Haplotype Blocks and Milk Traits
The LD analysis performed on the three genes revealed only one LD block (Figure 2). The haplotype tagging SNPs were evidenced at the POFUT1 gene, rs424501869 and rs408068827, which were enclosed in Block 1 (Figure 2b

Linkage Disequilibrium and Association Analysis between Haplotype Blocks and Milk Traits
The LD analysis performed on the three genes revealed only one LD block (Figure 2). The haplotype tagging SNPs were evidenced at the POFUT1 gene, rs424501869 and rs408068827, which were enclosed in Block 1 (Figure 2b) showing three haplotypes: AA (frequency 0.610); GA, (frequency 0.042); GC (frequency 0.348). Haplotypes at Block 1 and milk traits were submitted to statistical analysis by using model (2). Table S4. Similarly to the results regarding the investigation of single SNPs, fixed effects of farm, parity and stage of lactation included in model (2) showed high levels of significance, but these are not reported in tables.

Results of the statistical analysis (F-values and significance) are reported in Supplementary Information
Haplotypes at Block 1 were significantly associated with milk lactose content (p = 0.048), with the lowest concentration for GA ewes (Figure 3a). Regarding RCT and k20, milk samples from ewes carrying the GA haplotypes showed rennet coagulation time about 7 minutes (p = 0.008; Figure 3b), and curd firming time about 1.70 minutes (p = 0.001; Figure 3c) longer than AA and GC. Finally, in GA ewes curd firmness measured by a30 was smaller (p = 0.003; Figure 3d) and tmax delayed (p = 0.003; Figure 3e). Haplotypes at Block 1 and milk traits were submitted to statistical analysis by using model (2). Table S4. Similarly to the results regarding the investigation of single SNPs, fixed effects of farm, parity and stage of lactation included in model (2) showed high levels of significance, but these are not reported in tables.

Results of the statistical analysis (F-values and significance) are reported in Supplementary Information
Haplotypes at Block 1 were significantly associated with milk lactose content (p = 0.048), with the lowest concentration for GA ewes (Figure 3a). Regarding RCT and k 20 , milk samples from ewes carrying the GA haplotypes showed rennet coagulation time about 7 min (p = 0.008; Figure 3b), and curd firming time about 1.70 min (p = 0.001; Figure 3c) longer than AA and GC. Finally, in GA ewes curd firmness measured by a 30 was smaller (p = 0.003; Figure 3d) and t max delayed (p = 0.003; Figure 3e).

Discussion
We evaluated the association between SPP1, POFUT1 and PRLR genes and milk production, composition and coagulation properties in 380 Sarda breed ewes. We considered the abovementioned genes because they had been reported to be positional candidates in relation to specific dairy traits [11], and they have been poorly studied for their effects in many sheep populations.
Osteopontin (encoded by the SPP1 gene) is involved in immune regulation and tissue remodeling [36], and it plays an important role in mammary gland development and local immunity,

Discussion
We evaluated the association between SPP1, POFUT1 and PRLR genes and milk production, composition and coagulation properties in 380 Sarda breed ewes. We considered the abovementioned genes because they had been reported to be positional candidates in relation to specific dairy traits [11], and they have been poorly studied for their effects in many sheep populations.
Osteopontin (encoded by the SPP1 gene) is involved in immune regulation and tissue remodeling [36], and it plays an important role in mammary gland development and local immunity, as well as in milk production [16]. Gutierrez-Gil et al. [11] found that SPP1 is located in a genomic region on OAR6, a candidate for milk production traits and lactation regulation in sheep. Ruiz-Larrañaga et al. [37] reported that SPP1 is in a highly selected genomic region displaying selection signatures in Latxa sheep, a breed which has undergone long dairy selection pressure. Yurchenko et al. [38] performed high-density genotyping of 15 local sheep breeds from Russia and found that SPP1 is in a candidate region for milk and lactation traits. However, other studies carried out in cattle and sheep populations are not in accordance with that last hypothesis. The SPP1 SNP c.8514C>T was investigated in Holstein-Friesian cattle [39] and in the dairy cattle breed Girolando [40], and no significant associations with disease resistance or milk parameters were evidenced. García-Fernández et al. [41] did not find any significant association between SPP1 polymorphisms and milk production traits in a commercial population of Churra ewes. In the present study, one exonic SNP of the SPP1 gene, rs161844011, was significantly associated with SCS, with the TT genotype displaying the highest least mean square. The SNP rs161844011 is a missense variant located on exon 7, causing the amino acid variation p.Gln235Arg, which has a score of 0.24 according to the SIFT (scale-invariant feature transform) algorithm (http://sift.jcvi.org/), indicating it is tolerated by the protein, based on sequence homology. The influence of the SPP1 SNP rs161844011 on SCS confirmed predictions based on genomic studies, searching for regions under selection. Alain et al. [16] studied three SNPs in the promoter region and one SNP in the 3'UTR of SPP1 and evidenced the occurrence of a link with SCS and resistance to mastitis in dairy cows. They also found evidence of the presence of haplotype blocks among the different SNPs, which were not evidenced in the present study (Figure 2a). Even if in sheep SCC is still far from being accepted as a universal tool to reveal mastitis because of the different cut-off limits based on the management or the breed [42], it is known that in the Sarda breed appropriate values of heritability can be exploitable to improve somatic cells values in breeding schemes at the field level [43]. In addition, the identification of novel favorable SPP1 genotypes in sheep is a good perspective to investigate the key role of osteopontin in improving cell-mediated immunity in infected udders.
The real implication of the influence on t max by the rs426249393 SNP at SPP1 is still to be clarified. The parameter t max was much shorter in heterozygous AG, and we can speculate this was attributable to one of the many secondary and still unknown roles exerted by osteopontin, a multifaceted protein [17]. In fact, Sheehy et al. [15] reported that during lactation of dairy cows, the gene product of SPP1 regulates the level of expression of two casein genes codifying for β and κ-casein, CSN2 and CSN3, the latter being strictly involved in rennet coagulation of milk [44].
The POFUT1 gene emerges as a positional candidate gene in relation to mammary gland development traits [11], and it also appears in the database of Ogorevc et al. [45] of cattle candidate genes for milk production. The expression of this gene is fundamental for the correct formation of the epithelial and myoepithelial cells that form the mammary alveolus [46]. We genotyped three intronic SNPs at POFUT1, one of which was monomorphic while two others were heterozygous and in high LD, forming a haplotype block. POFUT1 SNP rs424501869 was associated with two different traits of curd firmness, a 45 and a 60 . These two curd firmness traits have been measured in addition to the traditional value of a 30 , which is the curd firmness in mm that is observed 30 min after rennet addition, while a 45 and a 60 showed decreasing values and are linked to the phases of curd syneresis and whey expulsion [25]. The association of POFUT1 with this particular phase of milk coagulation was also found for SNP rs408068827, showing association with t max . The t max parameter is derived from the modeling of milk firmness [47] and represents the point at which CF t attains its maximum level and the point at which the effects of the two parameters k CF and k SR are equal but opposite in sign [26]. To the best of our knowledge, an association between POFUT1 polymorphisms and milk coagulation properties was highlighted in this investigation for the first time, although the underlying mechanisms still have to be elucidated. Genome-wide and pathway-based association, which are more informative than the analysis of a single or a limited number of SNPs, have revealed that milk coagulation and cheese-making traits in the bovine species are affected by the combined additive effect of clustered genes [48,49].
The PRLR gene encoding the prolactin receptor has emerged from genome wide studies because of its position in a genome region where positive selection was highlighted in sheep [11,20]. Here we investigated three intronic SNPs. Two of them had a high LD value, but none of the three were involved in a haplotype block. Only SNP rs400874750 had a significant association with the lactose and SCS content. The prolactin receptor has been mainly investigated as a regulator of reproduction, lactation and the association with milk traits in dairy cattle [20]. A stimulating finding from the present study showed by the polymorphism at rs400874750 was the unfavorable values recorded for concentration of lactose (lower) and SCS (higher) in milk samples from CC homozygote ewes (Figure 1f,h). Among the several processes occurring in the mammary epithelium, it is worth noting that macrophage activation and the synthesis of milk lactose are mediated by the PRL/PRLR system [21].
The effect of rs400874750 on milk coagulation traits, RCT and k 20 (Figure 1h,i), may be explained by the link between the PRL/PRLR system and the synthesis of milk proteins. Indeed, some negative isoforms of PRLR expressed in vitro in bovine mammary cells inhibit the transcription of milk protein genes [21], and those may be consequently involved in the complex stages from milk to gel formation and cheese-making [44]. The results recorded in the present study are in accordance with what has been recorded for the sheep species, since the PRLR short form has a negative effect on the activation of milk protein gene transcription [50], and lactose and SCS are indirect markers of udder inflammation and inversely correlated in milk (the higher the SCS the lower the lactose, and vice versa), causing the delay of rennet coagulation and curd firming times [43].

Conclusions
The polymorphism of SNPs at SPP1, POFUT1 and PRLR were investigated for the first time in the Sarda sheep, a specialized dairy breed. The present study evidenced a significant association between SNPs at the candidate genes and many milk traits in the sheep species. Despite the possibility that other SNPs that were not investigated in the present study might be responsible for the observed effects, the results represent encouragement to conduct further research and achieve the improvement in production from dairy sheep farms and cheesemaking plants.
Supplementary Materials: The following are available online at http://www.mdpi.com/2076-2615/10/7/1216/s1, Table S1: Context sequences of the 8 SNPs investigated at SPP1, POFUT1 and PRLR in the population of the Sarda sheep breed. Table S2. Descriptive statistics of milk yield and composition, milk coagulation properties (MCP) and curd firmness over time traits (CFt) from the sampled population of Sarda sheep (n = 380); Table S3. F-value and significance for milk yield and composition, milk coagulation properties (MCP) and curd firmness over time traits (CFt) according to the effect of each of the 7 polymorphic SNPs out of the 8 investigated at ovine SPP1, POFUT1 and PRLR genes in Sarda sheep (n = 380); Table S4. F-value and significance for milk yield and composition, milk coagulation properties (MCP) and curd firmness over time traits (CFt) according to the effect of LD Block1 at POFUT1 gene in Sarda sheep (n = 380).