Variation in the stearoyl-CoA desaturase gene (SCD) and its influence on milk fatty acid composition in late-lactation dairy cattle grazed on pasture

Abstract Gene markers have become useful tools for improving animal genetics and breeding since they improve the accuracy of selection for superior breeding stock. In this study, the stearoyl-CoA desaturase (Δ-9-desaturase) gene (SCD) was investigated in New Zealand pasture-grazed Holstein–Friesian × Jersey cows. Three nucleotide substitutions were identified in exon 5 of the gene (c.702A/G, c.762T/C and c.878C/T), and a single nucleotide substitution was identified in intron 5 (c.880+105A/G). The c.878C/T substitution would, if expressed, result in the amino acid substitution p.A293V. Four nucleotide substitutions (c.*1783A/G, c.*1883C/T, c.*1984G/A and c.*2066T/C/G) were identified in the 3′-untranslated region (3′-UTR), and these resulted in three nucleotide sequence variants (named a, b and c). The sequence that would encode valine (V) at position 293 of SCD was linked to 3′-UTR variant a, and the sequence that would encode alanine (A) was linked to variants b and c. The frequency of the genotypes was as follows: VV (equivalent to aa: 15.1 %), VA (equivalent to ab+ac: 50.0 %) and AA (equivalent to bb+cc+bc: 34.9 %). The cows with the V variant produced less C10:1, C12:1 and C14:1 fatty acid (FA) but more C10:0, C11:0, C14:0, C16:1 and C18:2 FA than the A variant cows (P<0.001). Effects of c.*1783A/G and c.*2066T/C/G on milk fat composition were also found for the AA cows. The presence of c was associated with decreased levels of C16:1 (P<0.001), C17:1 (P=0.001), C18:2 cis-9, trans-13 (P=0.045), C18:2 cis-9, trans-12 (P=0.018) FA and C16:1 FA index (P<0.001). The presence of b was associated with increased levels of C13:0 iso FAs (P<0.001), monounsaturated FA (MUFA; P=0.002) and C12:1 (P<0.001).

The gene (SCD) encoding for SCD is located on bovine chromosome 26, and it is expressed in a variety of tissues (Chung et al., 2000). In lactating ruminants, the expression of SCD occurs at high levels (Bernard et al., 2005;Bionaz et al., 2008;McDonald et al., 1973), and it is considered to be important with respect to milk fat composition in these animals (Gautier et al., 2006). Ntambi et al. (2004) reported that dietary factors could regulate the expression of SCD, and hence the expression patterns of SCD in cows that are grazing outdoor on pasture might therefore be different to cows that are fed indoors on supplements. This is supported by the observation of Elgersma (2015) that milk from grazing-based production systems has less SFA and more polyunsaturated FA (PUFA), which is considered beneficial for health.
Other than the effect of dietary factors, nucleotide sequence variation in SCD is reported to be another factor that can change milk traits. For example, a nonsynonymous nucleotide substitution in SCD exon 5 (c.878C/T), which causes the substitution of valine (V) with alanine (A) at position 293 of the protein (p.A293V), has been associated with some milk traits.
In Japanese beef cattle, Taniguchi et al. (2004) identified more variants in the 3 -untranslated region (3 -UTR) of SCD, including c.*829C/T, c.*2066T/C/G, c.*2273G/A, c.*2458G/A and c.*3649A/T (these were labelled as 1905, 3143, 3351, 3537 and 4736 in their study), but the effect of these variations on intra-muscular fat composition was not significant, and their effect on milk FA is unknown. Hence, the objective of this study was to investigate the relationships between SCD variation in two regions -the first spanning p.A293V and the second spanning the nucleotide variation identified by Taniguchi et al. (2004) -and milk traits in Ki-wicross™ cows during late lactation and in a wholly pasturebased outdoor dairy production system.

Cows studied and milk sample collection
This research was approved by the Lincoln University Animal Ethics Committee (AEC Number 521) under the provisions of the Animal Welfare Act 1999 (NZ Government).
In total, 450 Holstein-Friesian × Jersey (HF × J) crossbred (Kiwicross™) dairy cows from two herds (124 cows in herd 1 and 326 cows in herd 2) were studied. The cows were between 3 to 10 years of age, and they were grazed solely outdoors on pasture (a mixture of perennial ryegrass and white clover) on the Lincoln University Dairy Farm (LUDF; Canterbury, NZ). The cows calved over the period August-September, and they were then milked twice a day until the end of May in the following year.
The milk yield in litres per day was recorded using Trutest milk meters (Tru-test Ltd, Auckland, NZ). Milk samples for analysis were collected once a month from September to February, and these were analysed for fat percentage (%) and protein percentage (%) using Fourier-transform infrared spectroscopy (MilkoScan FT 120,Foss,Hillerød,Denmark). The milk samples for individual FA analysis were collected at one afternoon milking in mid-January (days in milk (DIM) = 148 ±19 d) and frozen at −20 • C. After freezing, the samples were freeze-dried and then individually ground to a fine powder for component analysis.

Gas chromatography of the fatty acids in the milk sample
Fatty acids in the milk samples were analysed by gas chromatography as FA methyl esters (FAME), as described in Li et al. (2019).

PCR-SSCP analysis and genotyping
A blood sample from each of the cows was collected onto FTA cards. These were air-dried and stored for analysis. For the genetic analysis, leukocyte DNA was purified from a 1.2 mm punch of the dried blood spot using the method described by Zhou et al. (2006). Two pairs of PCR primers were designed based on the cattle reference sequence (NM_173959.4) to amplify two target regions of SCD. The first pair (forward: 5 -AATCAGGTAGGTCTCAGCG-3 and reverse: 5 -TTCTAATACTGTCCCTTAG-3 ) amplified a fragment of 436 bp (Region 1) spanning part of intron 4, exon 5 and part of intron 5. The second pair (forward: 5 -GAACCACTGTTTCTCTTTAC-3 and reverse: 5 -CACTTTGGAACCTGCCTTTG-3 ) amplified a fragment of 397 bp (Region 2) containing part of the 3 -UTR. The primers were synthesized by Integrated DNA Technologies (Coralville, IA, USA).
The PCR amplifications were performed as 15 µL reactions containing the purified genomic DNA on a 1.2 mm punch of the purified FTA paper, 0.25 µM of each designed primer, 150 µM of each dNTP (Bioline, London, UK), 2.5 mM of Mg 2+ , 0.5 U of Taq DNA polymerase (Qiagen, Hilden, Germany) and 1× the reaction buffer supplied with the DNA polymerase enzyme.
The amplifications were undertaken using S1000 thermal cyclers (Bio-Red, Hercules, CA, USA). The thermal profile included an initial denaturation for 2 min at 94 • C, which was followed by 35 cycles of 30 s at 94 • C, 30 s at 50 • C (for Region 1) or 56 • C (for Region 2) and 30 s at 72 • C. A final extension was undertaken for 5 min at 72 • C.

Sequencing of the SCD variants and sequence analysis
Homozygous PCR amplicons identified using PCR-SSCP were sequenced at the Lincoln University DNA Sequencing Facility. If homozygous samples were not found upon PCR-SSCP analysis, individual variants were isolated from heterozygous cattle and sequenced using an approach described previously (Gong et al., 2011). Briefly, single bands of interest from the heterozygous cattle were recovered directly from the SSCP gels as a gel slice. This was macerated and the DNA was eluted into 50 µL TE buffer by incubating at 70 • C for 20 min. The original primers and 1 µL of the eluted solution (as a template) were used for a second round of PCR amplification to produce a simple SSCP gel pattern which could be directly compared to or found in, the pattern derived from the original heterozygous cow. When banding patterns could be matched and identified, then the second homozygous PCR amplicons were directly sequenced at the Lincoln University DNA Sequencing Facility. The computer program DNAMAN (version 5.2.10, Lynnon BioSoft, Canada) was used for sequence alignment and comparisons. The BLAST algorithm was used to search the NCBI GenBank database (https://blast.ncbi.nlm.nih.gov, last access: 23 January 2020) for homologous sequences.

Statistical analyses
Hardy-Weinberg equilibrium (HWE) for the SCD genotypes was analysed using an online chi-square calculator (http://www.husdyr.kvl.dk/htm/kc/popgen/genetik/applets/ kitest.htm, last access: 2 April 2020). Other statistical analyses were carried out using IBM SPSS version 22 (IBM, NY, USA). Associations between variation in SCD and variation in gross milk traits and milk FA traits were tested using general linear mixed-effects models (GLMMs).
Single-variant presence/absence models (fixed effects: DIM, age and herd) were used to ascertain which variants should be analysed in subsequent multi-variant models. The multi-variant models included any variant that had a variant -gross milk trait or variant -FA trait association in the single-variant presence/absence analysis with a P value of less than 0.200 (and thus were potentially affecting the trait). The multi-variant models were again corrected for the fixed effects of (DIM, age and herd) and with other variants fitted as random effects.
A GLMM (fixed effect: genotype, DIM, age and herd) and multiple pair-wise comparisons with Bonferroni corrections were used to ascertain the effect of the different genotypes with a frequency greater than 5 % (thus insuring adequate sample size), on gross milk production traits (milk yield, fat percentage and protein percentage). Next, a GLMM (fixed effect: genotype, DIM, age and herd) and multiple pair-wise comparisons with Bonferroni corrections were used to ascertain the effect of genotypes on milk FA component levels.
The effect of cow sire could not be included in the GLMMs. Some semen straws (sire genetics) used in NZ dairy cattle artificial -insemination-based breeding approaches contain mixed-sire semen purchased from commercial semen producers. In these cases, individual sire identity is impossible to ascertain, but because the straws were mixed-semen straws and because different sires are used for different inseminations in different years, it is unlikely that sire was a strongly confounding effect. Cow age and herd might also be confounded with sire, but this cannot be confirmed.

Variation in SCD
Eight nucleotide substitutions were found in the two regions of SCD investigated (Fig. 1).
Three nucleotide substitutions (c.702A/G, c.762T/C and c.878G/A) in exon 5 and one nucleotide substitution (c.880+105A/G) in intron 5 (Region 1) and four nucleotide substitutions (c.*1783A/G, c.*1883C/T, c.*1984G/A, and c.*2066T/C/G) in the 3 -UTR (Region 2) were found in the Kiwicross™ cows studied. The variants A and V (from p.A293V), described in previous studies, were identified by these nucleotide variations (Fig. 1a), and three variants (a, b and c) could be identified when the 3 -UTR variations were analysed (Fig. 1b). A linkage was observed between c.702A/G, c.762T/C, c.878C/T, c.880+105A/G, c.*1883C/T and c.*1984G/A. The V variant in Region 1 corresponded to 3 -UTR variant a in Region 2, and the A variant in Region 1 corresponded to variants b and c in Region 2.
The overall frequency of the genotypes was VV (equivalent to aa: 15.1 %), VA (equivalent to ab + ac: 50.0 %) and AA (equivalent to bb + cc + bc: 34.9 %) across both herds, with the individual herds having frequencies of A (58.13 %) and V (41.87 %) for Herd 1 and A (60.35 %) and V (39.66 %) for Herd 2. Overall, the dominant variant was A (59.9 %) and the frequency of V was 40.1 %. The P value for the chi-square for deviation from HWE for p.A293V was 0.388, suggesting that across both herds the population was at equilibrium.
Six genotypes aa, ab, ac, bb, bc and cc were identified in the 3 -UTR region amplified, with overall frequencies of 15.1 %, 10.7 %, 39.3 %, 1.8 %, 10.7 % and 22.4 % respectively. The frequencies for the individual SNPs in both herds are illustrated in Table 1. The most common variant was c (47.4 %), and the frequencies of a and b were 40.1 % and 12.4 %. The P value for the chi-square for deviation from HWE was 0.724, suggesting that across both herds the population was at equilibrium.

Milk traits, milk fat compositions and SCD variation
The average phenotypic measures for the gross milk traits of milk yield, milk fat percentage and milk protein percentage were 20.86 ± 0.389 L d −1 , 5.07 ± 0.044 % and 4.22 ± 0.024 % respectively for Herd 1. In Herd 2 they were 22.58± 0.211 L d −1 , 5.06±0.031 % and 4.11±0.017 % respectively, but no associations were revealed between the SCD variation and these gross milk traits. In contrast, Table 2 summarizes the associations revealed between the composition of milk fat and the SCD variants defined by p.A293V.
Tables 3 and 4 summarize the associations revealed between the SCD variants a, b and c and the composition of individual and grouped FAs respectively. The results are not presented if no association was found with the variants. Overall, the presence of variant a was associated with lower C10:1 index, C12:1 index, and C14:1 index values but elevated C16:1 levels. Variant b or c appeared to have an opposite effect on the milk FAs when compared to variant a. There were higher C10:1 index, C12:1 index, C14:1 index and C18:1 index values when variant b was present. Table 5 summarizes the associations revealed between the composition of milk fat and the SCD genotypes from the 3 -UTR variants. Results are only presented if an association was found. The frequency of the bb genotype was low (1.8 %). Thus bb cows were not included in the analyses. The effects of c.*2066T/C/G and c.*1783A/G on milk fat composition were observed as the difference between bc and cc genotype cows, with there being less C16:1, C17:1, C18:2 cis-9, trans-13 and C18:2 cis-9, trans-12 in the milk from cc cows.

Discussion
The enzyme SCD is a rate-limiting enzyme in the biosynthesis of many monounsaturated FAs (Ntambi et al., 2004), and sequence variation in SCD has been revealed to affect milk fat composition. The SCD gene is therefore considered to be a useful candidate gene for use in breeding programmes targeted at improving the nutritional value of milk. The most commonly described variation in SCD is the variant c.878C/T located in exon 5. It underpins the substitution of the amino acid valine (V) with alanine (A) at amino acid position 293 in the SCD protein. The function of SCD is likely to be affected by p.A293V because it is located in the third histidine-rich region of the enzyme. This histidine-rich region has been revealed to be important for the catalytic activity of the enzyme (Shanklin et al., 1994).
There was linkage between c.878C/T and other variants. Taniguchi et al. (2004) identified linkage between three nucleotide substitutions c.702A/G, c.762T/C and c.878C/T in exon 5 of Japanese black cattle. Based on these three variants, they described the haplotype p.293V (GCT) and p.293A (ATC). This linkage was also found by Baeza et al. (2013) with the variant g.10153A> G (equivalent to c.702A/G here) being in complete linkage disequilibrium with c.878C/T in the beef cattle they studied. In this study, three nucleotide variations in exon 5 (c.702A/G, c.762T/C and c.878C/T), one variant in intron 5 (c.880+105A/G) and four variants in the 3 -UTR (c.*1883C/T, c.*1984G/A, c.*1783A/G and c.*2066T/C/G) were revealed. There was linkage between In an in vitro study, Enoch et al. (1976) revealed that the acyl-CoA derivatives with 12 to 19 carbon atoms were required as substrates for SCD enzyme activity. Schennink et al. (2008), Kgwatalala et al. (2009b) and this research (Table 1) all suggested that SCD p.A293V has effects on individual FA and FA index levels, and especially in this study on the level of C10:0, C10:1, C12:1, C14:0, C14:1 and C16:1, as well as the C10:1 index, C12:1 index, C14:1 and C16:1 index levels. In addition, Enoch et al. (1976) found that the SCD enzyme has substrate specificity, with a preference for longer-chain FAs. However, the effect of p.A293V on C18 FA levels (such as the C18:1 trans-11, C18:1 cis-9, C18 index) was not confirmed here.
The effect of p.A293V on C18 FA levels could be influenced by different factors, including seasonal variation, the impact of other genes, variation between breeds and the stage of lactation. Duchemin et al. (2013) revealed that the p.293V allele was negatively associated with C18:1 trans-11 and that this negative effect was larger in summer than in winter. Schennink et al. (2008) investigated the joint effect of SCD and DGAT1 variation and suggested that the genetic variation explained by DGAT1 p.232K and by SCD p.293A is additive with respect to their effect on C16, C18 and CLA levels. Moioli et al. (2007) revealed that p.A293V affected C10:1, C14:1, C16:1, C10 index and C14 index levels in Pied-montese (n = 81), Jersey (n = 75) and Valdostana (n = 730) cows, but the effect of p.293A on C18 levels was only observed when the variant frequency was high (i.e. when at a frequency of 0.94 in Jersey cattle and 0.65 in Valdostana cattle). With a lower frequency of occurrence of p.293A (i.e. 0.42 in Piedmontese cattle), the effect of SCD variation on C18 levels was not significant. Mele et al. (2007) investigated the effect of p.293A and lactation stage on C18:1 cis-9 FA levels. A negative effect was observed when the DIM increased (P < 0.01), and that p.293A occurrence (the frequency was 0.57) had a positive effect on the C18:1 cis-9 FA level (P < 0.05). Finally Carvajal et al. (2016) revealed associations between SCD p.A293V genotypes and C6:0, C8:0, C10:0, C11:0, C14:0, C16:0, C18:0, C20:0, C16:1, C18:1 cis-9, linoleic acid, α-linolenic acid and CLA cis-9 trans-11 levels.
In our study, a higher frequency of the p.293A allele (59.9 %) was observed in the Kiwicross™ cows. This is consistent with the findings of Schennink et al. (2008), Kgwatalala et al. (2009b) and Mele et al. (2007), who reported similar results, with frequencies of 73 % in Dutch HF cows, 69 % in Canadian HF cows and 57 % in Italian HF cows, respectively.
The 3 -UTR can play a critical role in translation termination and post-transcriptional gene expression (Barrett et al., 2012). Baeza et al. (2013) reported that the 3 -UTR variant g.15001A>G (equivalent to c.*1783A/G here) affected meat    1 Predicted means and standard error of those means derived from GLMM. "Cow age", "days in milk (DIM)" and "herd" were fitted to the models as fixed effects. 2 P < 0.05 in bold; 005 < P < 0.2 in italics.
fat composition (C14:1 level, P < 0.05) in their Argentinian Brangus beef cattle. In the results reported here, after differentiating the p.293A variant to b and c (based on 3 -UTR typing), the different effects of the variants on long-chain FA levels were revealed. The decline in C16:1, C17:1, C18:2 cis-9, trans-13 and C18:2 cis-9, trans-12 FA levels and C16:1 index level was associated with the presence of variant c. The presence of variant b was associated with an increase in C18:1 cis-9 and total C18:1 FA levels, MUFA levels, C18:1 index level and MUFA index level. At the level of genotype (Table 5), there was a significant difference between aa and cc cows in C16:1, C17:1, C18:2 cis-9, trans-13 and C18:2 cis-9, trans-12 FA levels. In addition, there was a significant difference between aa and bb cows for C20:0 levels and MUFA levels. This suggests the 3 -UTR (specifically the substitutions c.*1783A/G and c.*2066T/C/G) is better able to resolve the effect of SCD on milk FA levels. In this respect, Kgwatalala et al. (2009b) regarded the effect of different lactation stages on milk MUFA was more important than p.A293V, and this along with the results presented here suggests more research into SCD variation and lactation stage is needed to gain clarity into what is driving variation in FA levels. The influence of 3 -UTR variation was also described by Kgwatalala et al. (2009a) in 46 Holstein and 35 Jersey cows in Canada. In their study, three haplotypes H1, H2 and H3 (equivalent to the variants c, b and a here) were identified with the frequency 67.1 %, 2.3 % and 30.6 % respectively. Significant differences were only found in two milk FA levels, with their H1H1 cows producing more C10:1 and C12:1 than the H3H3 cows (similar results were observed when comparing cc and aa cows here). Moreover, Kgwatalala et al. (2009a) reported that an internal ribosome entry site (IRES) could be found in the c variant only. The presence of an IRES motif may ultimately affect SCD protein turnover or the quantity of the SCD enzyme produced, because it could enhance translation of the constituent mRNA. Kgwatalala et al. (2009a) suggested that the nucleotide vari- Table 5. Association between milk fatty acid levels and 3 -UTR genotypes 1 .

FAME
Mean FAME level ± SE 2 (g 100 g −1 milk FA) P ation in 3 -UTR region might lead to the absence of the IRES in the a and b variants. With a low frequency of b in their Holstein (0.054) and Jersey (0.000) cattle, Kgwatalala et al. (2009a) suggested that milk fat composition was not affected by variant b significantly. Therefore, the variation in C10:1 and C12:1 FA levels in their study could be due to either the variation in the 3 -UTR (c.*1783A/G) or the variation in exon 5 (c.878C/T). If this was true, then the effect of c on milk fat composition should be the opposite of the effects of a, b or a + b. In the Kiwicross™ cows studied here the frequency of b was 0.124, and opposite effects were observed between variant c and the other two variants for C16:1, C18:1 cis-9, C18:2 trans-9, 12, C18:2 cis-9, trans-12, C18:2 cis-9, trans-13, C19:0 FA levels and C16:1 index levels (Tables 3 and 4). At the genotype level (Table 5), the cc cows produced more C10:1 and C12:1 FA than the aa cows. This is similar to what was reported by Kgwatalala et al. (2009a). Moreover, the ab cows produced less C10:1 and C14:1 cis-9 FA but more C18:2 cis-9, trans-12 and C16:1 FA than the cc cows.
The effect of SCD on gross milk traits (milk yield, fat and protein percentage) is still disputed. Macciotta et al. (2008) investigated 313 Italian Holstein cows and found that their VV cows (referring to p.A293V) had higher milk yields and protein yields than their VA and AA cows. The associations they found appear to be consistent across different stages of lactation. However, Mao et al. (2012) and Signorelli et al. (2009) reported a significant negative effect of the V allele on milk yield in Chinese Holstein, Piedmontese and Valdostana breeds. In addition, Schennink et al. (2008) did not find any significant associations between p.A293V and milk traits in Dutch Holstein, a result that was consistent with our findings here in and Kiwicross™ cows. Data availability. The original data are available upon request to the corresponding author.
Author contributions. YL, HZ, LC and JH conceived the research idea. YH performed the experiments, analysed the data and drafted the manuscript. HZ participated in statistic analysing and supplied the technical guidance. LC participated in the sampling. JZ performed the GC analysing. JH provided the experimental environment and participated in the manuscript writing.