Dairy sheep breeding has growing importance worldwide mainly because of its processed dairy products. Dairy sheep are reared in several European countries, especially the southern regions surrounding the Mediterranean Sea (Carta et al., Reference Carta, Casu and Salaris2009). Italy produces over 4% of the world's ovine milk (FAOSTAT, 2016) mainly provided by ewes of Sarda breed, which is considered one of the most important Italian dairy breeds (Dettori et al., Reference Dettori, Pazzola, Pira, Paschino and Vacca2015). Sarda sheep milk is almost entirely used in the production of cheese, and three cheeses produced in Sardinia are recognized by the European Union (EU) as Protected Designation of Origin (PDO) (Cipolat-Gotet et al., Reference Cipolat-Gotet, Cecchinato, Pazzola, Dettori, Bittante and Vacca2016). Given the growing economic importance of the sheep cheese production sector, recent investigations have been devoted to better understand sheep milk coagulation properties (MCP), as a valid tool to be made available to the dairy industry.
MCP can be traditionally measured using a lactodynamographic instrument, which detects three single point parameters: rennet coagulation time (RCT, min), curd firming time (k 20, min) and curd firmness at 30 m of analysis (a 30, mm), first described by McMahon and Brown (Reference McMahon and Brown1982). In dairy cattle MCP have been proved to be independent from milk yield, and mainly influenced by the titratable acidity of milk (Bittante et al., Reference Bittante, Penasa and Cecchinato2012). In addition, RCT and a 30 are strongly and negatively correlated. In contrast, the correlation between RCT and a 30 was not evident in milk samples from Sarda breed ewes (Pazzola et al., Reference Pazzola, Dettori, Cipolat-Gotet, Cecchinato, Bittante and Vacca2014). The same authors extended the MCP analysis up to 60 m, obtaining also curd firmness at 45 (a 45) and 60 (a 60) min and observed that, in comparison with bovine milk, sheep milk has a very early gelation time (RCT = 8.6 vs. 10–20 min), a rapid increase in curd firming time (k 20 = less than 2 vs. 5–15 min), and a higher curd firmness at a 30 (a 30 = 50 vs. 35 mm) (Bittante et al., Reference Bittante, Penasa and Cecchinato2012; Pazzola et al., Reference Pazzola, Dettori, Cipolat-Gotet, Cecchinato, Bittante and Vacca2014). In addition to MCP, several research papers exploited all available lactodynamograph data to model curd firming over time (CFt) in milk from different species (Bittante et al., Reference Bittante, Pellattiero, Malchiodi, Cipolat-Gotet, Pazzola, Vacca, Schiavon and Cecchinato2014; Stocco et al., Reference Stocco, Cipolat-Gotet, Bobbo, Cecchinato and Bittante2017; Pazzola et al., Reference Pazzola, Stocco, Dettori, Cipolat-Gotet, Bittante and Vacca2018). These modeled parameters have proven more informative than the traditional MCP. A four parameter model was applied to cow milk coagulation and curd firming test prolonged from 30 to 90 min (Bittante et al., Reference Bittante, Contiero and Cecchinato2013). The four parameter model was modified and applied to test coagulation ability of sheep milk from Sarda sheep by Cipolat-Gotet et al. (Reference Cipolat-Gotet, Pazzola, Ferragina, Cecchinato, Dettori and Vacca2018).
MCP are clearly influenced by genetic factors such as species, breed and the individual animal (Bittante et al., Reference Bittante, Penasa and Cecchinato2012), and evidence has been given for the casein genotype (Ceriotti et al., Reference Ceriotti, Chiatti, Bolla, Martini and Caroli2005; Noce et al., Reference Noce, Pazzola, Dettori, Amills, Castelló, Cecchinato, Bittante and Vacca2016). The Growth Hormone Releasing Hormone Receptor (GHRHR), the Growth Hormone Receptor (GHR), and the Insulin-Like Growth Factor 1 (IGF1) genes have been considered as potential candidate genes for milk quality traits in cattle (Viitala et al., Reference Viitala, Szyda, Blott, Schulman, Lidauer, MakiTanila, Georges and Vilkki2006; Banos et al., Reference Banos, Woolliams, Woodward, Forbes and Coffey2008). These genes are involved in the growth hormone (GH) system. GH regulates many physiological functions, such as metabolism, growth, reproduction, feeding, osmoregulation and immune system function (Bergan-Roller and Sheridan, Reference Bergan-Roller and Sheridan2018), in addition to its effects on mammary development (mammogenesis) and milk production (galactopoiesis) (Akers, Reference Akers2017). The GHRHR protein, expressed in somatotropic cells, mediates the production and release of GH from the somatotropic cells, upon ligand binding with the hypothalamic factor GH releasing hormone (GHRH) (Pang and Chan, Reference Pang, Chan, Coleman and Tsongalis2010). The actions triggered by GH are mediated by its specific receptors (GHR) distributed among tissues, which in turn are regulated at the expression level by several factors reflecting the metabolic and nutritional status of the organism (Bergan-Roller and Sheridan, Reference Bergan-Roller and Sheridan2018). GHRs can be associated, within the cell, with different effectors, which in turn, can cause different responses upon GH activation. When GHR is associated with the Janus tyrosine kinase-signal transducer and activator of transcription (JAK/STAT), PI3K-protein kinase B (Akt) and extracellular signal regulated kinase (ERK), it causes the synthesis and secretion of IGF1 polypeptide hormone and therefore the pathway of cell growth (Herington and Lobie, Reference Herington and Lobie2012). In contrast, when GHRs are associated with intracellular effectors as cAMP/protein kinase C (PKC) pathways, they mediate lipolytic GH signaling by targeting expression and activation of lipases (Chaves et al., Reference Chaves, Frasson and Kawashita2011). Many of the growth-promoting effects of GH are mediated by IGF1. Circulating GH stimulates the synthesis and secretion of IGF1 from the liver, and IGF1, in turn, stimulates cell growth and differentiation in a variety of target tissues, through distinct IGF receptors (Laviola et al., Reference Laviola, Natalicchio and Giorgino2007).
Dettori et al. (Reference Dettori, Pazzola, Paschino, Amills and Vacca2018) investigated association between a panel of 36 SNPs within the GHR, GHRHR and IGF1 genes and milk production and quality traits in Sarda sheep, revealing that the GHR gene is associated with daily fat and protein yield. They also revealed the IGF1 gene is associated with milk protein and casein content. The present study aims to extend this work, exploring associations between the 36 SNP panel of the GHRHR, GHR and IGF1 genes and traditional and modeled MCP in Sarda sheep.
Materials and methods
No specific authorization from an animal ethics committee was required, because according to the EC Directive 86/609/EEC and Directive 2010/63/EU, none of the procedures met the criteria to be defined as an experiment or procedure. Blood samples for DNA isolation were collected by experienced veterinarians and milk samples were collected concurrently with official sampling procedures for performance controls of the flock book.
A total of 380 lactating ewes, in their first to seventh parity, were sampled from 19 farms (20 ewes per farm) located in Sardinia (Italy). The ewes were included in the selection scheme of the Sarda breed and registered in the flock book. Ewes were between 2 and 7 months after parturition. Detailed description of farms, animals and sampling is given in Pazzola et al. (Reference Pazzola, Dettori, Cipolat-Gotet, Cecchinato, Bittante and Vacca2014) and Vacca et al. (Reference Vacca, Pazzola, Dettori, Pira, Malchiodi, Cipolat-Gotet, Cecchinato and Bittante2015). Ewes from each flock were individually sampled in a single day (one sampling day for each flock). During the afternoon milking 200 ml of milk were collected from each ewe. Milk samples were maintained at 4 °C and were analyzed within 24 h. Individual blood samples were collected in K3EDTA vacuum tubes (BD Vacutainer, Becton Dickinson, Franklin Lakes, NJ) from each ewe for genomic DNA isolation, performed with the Puregene Blood Kit (Qiagen, Hilden, Germany). The concentration and purity of DNA were determined with an Eppendorf BioPhotometer instrument (Eppendorf, Hamburg, Germany).
MCP were measured with the Formagraph instrument (Foss Italia, Padova, Italy). Individual milk samples (10 ml × 2 replicates) were heated to 35 °C and they were added 200 µl of rennet solution (Hansen Naturen Plus 215, Pacovis Amrein AG, Bern, Switzerland), containing 80 ± 5% chymosin and 20 ± 5% pepsin (215 international milk clotting units per ml, IMCU/ml), which was diluted to 1.2% (wt/vol) in distilled water to achieve 0.0513 IMCU/ml milk. The traditional single point parameters RCT, k 20 and a 30 were recorded, and the analysis was extended to 60 m to obtain the values of curd firmness at 45 (a 45) and 60 (a 60) min. Six milk samples were excluded from the statistical analyses as did not coagulated. In addition to traditional single point parameters, we retrieved from the Formagraph instrument the specific file containing the complete record of curd firming values (expressed as the width of the oscillatory graph, in mm), detected every 15 s. This created a total of 240 CF values for each replicate, for a 60 min run. Data obtained were included in the four-parameter model described by Bittante et al. (Reference Bittante, Contiero and Cecchinato2013):
$${\rm C}{\rm F}_t = {\rm C}{\rm F}_{\rm P} \times \lsqb {1-{\rm e}^{-k_{{\rm CF}} \times (t-{\rm RC}{\rm T}_{{\rm eq}})}} \rsqb \times {\rm e}^{-k_{{\rm SR}} \times (t-{\rm RC}{\rm T}_{{\rm eq}})}$$where CFt is curd firmness at time t (mm); CFP is the asymptotical potential maximum value of CF at an infinite time (mm); k CF is the curd-firming instant rate constant (% × min−1); k SR is the curd syneresis instant rate constant (% × min−1), and RCTeq is the rennet coagulation time estimated by the model, on the basis of all data points (min). The CFP is conceptually independent from test duration and is not intrinsically dependent on RCT (unlike a 30). The parameter k CF describes the shape of the curve from the time of milk gelation to infinity, and is conceptually different from k 20 as it uses all available information. The parameter k CF is assumed to increase CF toward the asymptotic value of CFP, whereas k SR is assumed to decrease CF toward a null asymptotic value. In the initial phase of the test, the first rate constant prevails over the second, so that CFt increases to a point in time (t max) at which the effects of the 2 parameters are equal but opposite in sign; this is when CFt attains its maximum level (CFmax). Thereafter, CFt decreases, tending toward a null value due to the effect of curd syneresis and the resulting expulsion of whey.
The 36 SNP panel included 31 SNPs mapping to the sheep GHR gene, 2 SNPs of the GHRHR gene and 3 SNPs of the IGF1 gene, genotyped in the 380 Sarda sheep. Genotyping was carried out with a 12K Flex QuantStudio instrument (Thermo Fisher Scientific), based on a custom TaqMan Real-Time PCR assay (Thermo Fisher Scientific, Waltham, MA) as described in Dettori et al. (Reference Dettori, Pazzola, Paschino, Amills and Vacca2018).
The Haploview software package (Barrett et al., Reference Barrett, Fry, Maller and Daly2005) was used to estimate and plot pairwise linkage disequilibrium (LD) measures (D′ and r 2). The same tool was used to infer haplotype frequencies as well as to define LD blocks according to the Gabriel criteria (Gabriel et al., Reference Gabriel, Schaffner, Nguyen, Moore, Roy, Blumenstiel, Higgins, DeFelice, Lochner, Faggart, Liu-Cordero, Rotimi, Adeyemo, Cooper, Ward, Lander, Daly and Altshuler2002). Haplotype analysis revealed seven LD blocks within the GHR gene sequence, described in Dettori et al. (Reference Dettori, Pazzola, Paschino, Amills and Vacca2018).
The 240 CFt observations available for each sample were fitted with curvilinear regressions using the non-linear procedure (PROC NLIN) of SAS (version 9.4, SAS Institute Inc., Cary, NC). The Marquardt iterative method has been used according to Bittante (Reference Bittante2011).
Association analysis between GHRHR, GHR and IGF1 genotypes and experimental data regarding CFt modeling parameters was based on the following model (1):
$$Y_{ijklmn} = {\rm \mu} + G_i + F_j + P_k + {\rm DI}{\rm M}_l + {\rm SIRE}\lpar G \rpar _m + e_{ijklmn}$$where Y ijklmn is the observed trait (RCTeq, k CF, k SR, C FP, CFmax, and t max); μ is the general mean; G i is the fixed effect of the ith SNP genotype, one at a time (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 to 19 levels; the different farms where animals were reared); P k is the fixed effect of kth parity of the ewes (k = 1 to 4 levels; first to fourth or more parities); DIMl is the fixed effect of the lth days in milking (l = 4 levels; level 1: ≤100 d; 2: 101–140 d; 3: 141–160 d; level 4: ≥161 d); SIRE(G)m is the random effect of the mth sire (m = 108 different sires) nested within the genotype, and e ijklmn is the error random residual effect.
This model (1) was also used to investigate the association between both traditional and modeled MCP and each of the seven LD blocks, one at a time. In the single SNP and LD block analysis, we only considered SNPs with a MAF >0.05, to make sure that genotypic means are correctly estimated. The MIXED procedure of SAS (version 9.4, SAS Inst. Inc.) was used to carry out the association analysis and correction for multiple testing was implemented with the Bonferroni method (one milk trait for each SNP or LD block at a time).
Results and discussion
Descriptive statistics of traditional MCP and CFt model parameters of milk samples are shown in Table 1. All traits exhibited high variability, the coefficient of variation (CV) of traditional MCP traits was between 24.65% (for a 30) and 43.45% (for RCT), and k SR had the highest CV value (129.36%). Table 2 shows the F-values obtained from the analysis of variance of CFt model parameters, as a function of genotype of the GHR, GHRHR and IGF1 genes. The SNP genotypes exhibiting significant effects on phenotype variance are described in Table 3. Among the three genes analyzed, only the GHR polymorphism was significantly associated with the considered traits. The only physiological ligand of GHR is GH, and in the same breed, polymorphism of the GH gene was associated with milk yields (Vacca et al., Reference Vacca, Dettori, Balia, Luridiana, Mura, Carcangiu and Pazzola2013) and with lipid content, in addition to protein, casein and lactose contents (Dettori et al., Reference Dettori, Pazzola, Pira, Paschino and Vacca2015).
Descriptive statistics for traditional (MCP) and curd firming over time (CFt) coagulation properties of milk from Sarda sheep
CV, coefficient of variation; RCT, measured rennet coagulation time; k 20, time interval between coagulation and attainment of curd firmness of 20 mm; a 30, a 45 and a 60, curd firmness 30, 45 and 60 min after rennet addition; CFP, asymptotic potential curd firmness; k CF, curd firming instant rate constant; k SR, syneresis instant rate constant; CFmax, maximum curd firmness achieved within 45 min; t max, time at achievement of CFmax; RCTeq, RCT estimated according to curd firm change over time modeling (CFt).
Analysis of variance (ANOVA, F-values and significance) of curd firming over time (CFt) model parameters of milk from Sarda sheep
CFP, asymptotic potential curd firmness; k CF, curd firming instant rate constant; k SR, syneresis instant rate constant; CFmax, maximum curd firmness achieved within 45 min; t max, time at achievement of CFmax.
*P < 0.05; **P < 0.01.
Least square means of traditional (MCP) and curd firming over time (CFt) coagulation properties of milk from Sarda sheep according to the GHR gene SNP and genotypes (n = 380)
RCT, measured rennet coagulation time; k 20, time interval between coagulation and attainment of curd firmness of 20 mm; a 30, a 45 and a 60, curd firmness 30, 45 and 60 min after rennet addition; CFP, asymptotic potential curd firmness; k CF, curd firming instant rate constant; k SR, syneresis instant rate constant; CFmax, maximum curd firmness achieved within 45 min; t max, time at achievement of CFmax.
Means with different superscript capital or lower-case letters in each column differ significantly in genotype comparison at P < 0.01 and P < 0.05 respectively.
Statistical analysis highlighted a significant association of SNP rs404237321 with k SR. Figure 1 clearly depicted the effect of the SNPs on the pattern of coagulation, in particular, the rs404237321 CT genotype showed larger syneresis compared with CC genotype (Fig. 1a). The SNP rs404237321 was a missense variant of exon 5, causing the p.Gly147Asp variation in the GHR protein, and according to the SIFT (http://sift.jcvi.org/) prediction algorithm, it was not expected to affect protein function. As regards SNP rs426666828, the CC genotype showed higher k SR compared with CT and TT genotypes (Fig. 1b). This SNP was located in intron 3 of the GHR gene (13.8 kb from exon 4) and it was included in haplotype block 4, which was the largest in size, consisting of ten SNPs (Dettori et al., Reference Dettori, Pazzola, Paschino, Amills and Vacca2018). The SNP rs412881843 was significantly associated with both traditional RCT (data not shown) and the k SR value (P < 0.05); its effects on RCT, shorter for GG genotype, and k SR, are shown in Figure 1c. The SNP rs412881843 is localized in intron 3 (only 427 bp from exon 3) and linkage disequilibrium analysis revealed it was included in haplotype block 4, as was SNP rs426666828. In the resource population of the present paper the GG genotype of SNP rs412881843 was associated with an RCT value of 6.97 min, which is shorter than the average RCT value of 8.6 min found by Pazzola et al. (Reference Pazzola, Dettori, Cipolat-Gotet, Cecchinato, Bittante and Vacca2014). The heterozygote CG genotype of SNP rs412881843 was associated with a delayed value of RCT (9.09) with a similar mean value of RCT reported for the Brogna breed (Bittante et al., Reference Bittante, Pellattiero, Malchiodi, Cipolat-Gotet, Pazzola, Vacca, Schiavon and Cecchinato2014). Finally, the SNP rs402337124, located in the upstream region and included in haplotype block 7, was associated with k SR, lower for AA genotype (Fig. 1d).
Pattern of curd firming over time (CFt) of milk samples from SNP rs404237321 (a) rs426666828 (b), rs412881843 (c) and rs402337124 (d) according to their genotypes.
Although the literature is poor about this topic, Bittante et al. (Reference Bittante, Pellattiero, Malchiodi, Cipolat-Gotet, Pazzola, Vacca, Schiavon and Cecchinato2014) showed that an integration of the lipid fraction of the diet with rumen-protected conjugated linoleic acid, doubled the rate of whey expulsion (k SR trait) in Alpine sheep breeds. In a previous investigation on the same resource population (Dettori et al., Reference Dettori, Pazzola, Paschino, Amills and Vacca2018), the GHR gene has been shown to affect variation of the lipid content of milk, possibly indicating that the effect of GHR is not direct on coagulation, but mediated by the milk composition, in particular by the lipid content.
Linkage Disequilibrium (LD) analysis was performed from 29 informative SNPs in the GHR gene and seven regions of LD were identified (described in Dettori et al., Reference Dettori, Pazzola, Paschino, Amills and Vacca2018). Haplotype association analysis revealed a significant effect of block 1 on CFP (P < 0.05), with the lowest values recorded for haplotype H4 (CCG) (CFP of 24.24 mm vs. 63.32 of haplotype H1; Table 4). Haplotype H4 of block 1 was also associated with a reduction of lipid and casein contents and of milk energy in the same animals (Dettori et al., Reference Dettori, Pazzola, Paschino, Amills and Vacca2018). The molecular bases underlying the observed associations are unknown and need more investigation, especially in sheep. In fact, the GHR gene is characterized by high transcriptional complexity, due to the structural organization of the 5′ regulatory region of this gene (Adams, Reference Adams1995). Multiple forms of GHRs are known to exist in vertebrates, with specific tissue expression and differential expression in relation to the distinctive conditions of the organism (Bergan-Roller and Sheridan, Reference Bergan-Roller and Sheridan2018). In cattle there are multiple forms of GHR mRNA variants, with disparate tissue specific expression (Jiang and Lucy, Reference Jiang and Lucy2001), while two specific forms of GHR mRNA are currently known in the sheep: the ovine Pl promoter, with liver-specific expression in vivo (Adams et al., Reference Adams, Baker, Fiddes and Brandon1990) and the P2 promoter, with widespread tissue transcription (Adams, Reference Adams1995).
Least square means of milk traits according to the different haplotype blocks in Sarda sheep
SNPs in Block 1 are rs408890407 (C/T; exon 10, synonymous), rs55631463 (C/T; exon 10, missense) and rs405063669 (G/A; intron 8).
Means with different superscript capital or lower-case letters in each column differ significantly in haplotype comparison at P < 0.01 and P < 0.05 respectively.
In conclusion, this is the first research exploring the potential effects of the GHR, GHRHR and IGF1 genes on traditional MCP and CFt parameters, based on prolonged curd firmness recording. In particular, the study demonstrated that polymorphisms of the GHR gene are associated with milk rennet coagulation time and syneresis. These findings may be useful for the dairy industry, as well as for selection programs.
Acknowledgement
This research is a part of the MIGLIOVIGENSAR project, funded by Regione Autonoma della Sardegna, Italy (CUP: B82I13000580002).
