A sequence variant in the diacylglycerol O-acyltransferase 2 gene influences palmitoleic acid content in pig muscle

The bulk of body fat in mammals is in the form of triacylglycerol. Diacylglycerol O-acyltransferase 2 (DGAT2) catalyses the terminal step in triacylglycerol synthesis. The proximity of DGAT2 with stearoyl-CoA desaturase (SCD) in the endoplasmic reticulum may facilitate provision of de novo SCD-mediated fatty acids as substrate for DGAT2. Here, we first searched for sequence variants in the DGAT2 gene to then validate their effect on fat content and fatty acid composition in muscle, subcutaneous fat and liver of 1129 Duroc pigs. A single nucleotide polymorphism in exon 9 (ss7315407085 G > A) was selected as a tag variant for the 33 sequence variants identified in the DGAT2 region. The DGAT2-G allele increased DGAT2 expression in muscle and had a positive impact on muscular C14 and C16 fatty acids at the expense of C18 fatty acids. Although there was no evidence for an interaction of DGAT2 with functional SCD genotypes, pigs carrying the DGAT2-G allele had proportionally more palmitoleic acid relative to palmitic acid. Our findings indicate that DGAT2 preferentially uptakes shorter rather than longer-chain fatty acids as substrate, especially if they are monounsaturated, and confirm that fatty acid metabolism in pigs is subjected to subtle tissue-specific genetic regulatory mechanisms.

. The role of DGAT2 in the synthesis pathway of triacylglycerol. Triacylglycerol is (1) synthesized de novo by sequential addition of fatty acyl moieties to a glycerol 3-phosphate (the G3P pathway) or (2) formed through the re-esterification of partial glycerides that results from either the hidrolysis of pre-existing triacylglycerol or via the monoacylglycerol (MAG) pathway. De novo transformation of diacylglycerol to triacylglycerol is catalyzed by DGAT2 (1) but re-esterification of partial glycerides is done by either DGAT2 (2a) or DGAT1 (2b). However, while DGAT1 specifically incorporates pre-formed fatty acids, DGAT2, which is able to form protein complexes with MGAT2, FATP1 and SCD, can use both endogenous and pre-formed fatty acids as substrates. DGAT1 diacylglycerol O-acyltransferase 1, DGAT2 diacylglycerol O-acyltransferase 2, ELOVL fatty acid elongase, FATP1 fatty acid transport protein 1, MGAT2 monoacylglycerol acyltransferase 2, SCD stearoyl-CoA desaturase.

Animals and phenotypes.
A total of 1129 pigs from 158 sires and 559 dams of the same Duroc line were used in the experiment. Pigs were raised in 15 fattening batches between 2006 and 2016 following a similar protocol for data recording and tissue sampling 25 . In each batch, pigs were raised from 75 days of age until slaughter in the same farm under identical conditions. During this period, pigs had ad libitum access to commercial feed (Esporc, Riudarenes, Girona, Spain). At 206 (SD 8) days of age, all pigs were weighted and backfat (BT) and loin (LT) thickness were ultrasonically measured at 5 cm off the midline between the third and fourth last ribs using the portable equipment Piglog 105 (Frontmatec, Kolding, Denmark). Pigs were slaughtered at 211 (SD 9) days in a commercial slaughterhouse equipped with a carbon dioxide stunning system. Carcass composition traits were recorded, including carcass weight, BT and LT. Both BT and LT were measured at 6 cm off the midline between the third and fourth last ribs using an ultrasound automatic scanner (AutoFOM, SFK-Technology, Denmark). Immediately after slaughter, a sample of semimembranosus muscle (SM, n = 40) and subcutaneous fat (SF, n = 226) were collected, snap-frozen and stored at -80ºC. After chilling for about 24 h at 2 °C, a large sample of the muscles gluteus medius (GM, n = 1093) and longissimus thoracis (LM, n = 526) were collected, vacuum packaged, and stored at − 20 °C until required. The IMF content in GM and LM, as well as the fatty acid composition in GM, LM and SF were determined in duplicate by quantitative gas chromatography 26 . The amount of each fatty acid was expressed as the percentage of each individual fatty acid relative to total fatty acid. The proportion of SFA (C14:0, C16:0, C18:0 and C20:0), MUFA (C16:1n-7, C18:1n-9, C18:1n-7 and C20:1n-9) and PUFA (C18:2n-6, C18:3n-3, C20:2n-6 and C20:4n-6) were calculated.
Genome-wide association study. In order to identify candidate genomic regions associated with fatrelated traits, GWAS for each individual fatty acid were performed using a subset of 254 pigs. The genomic DNA of these pigs was isolated as described in Ref. 24 and used for single nucleotide polymorphism (SNP) genotyping with either the PorcineSNP60 v2 Genotyping BeadChip (n = 138) or the GGP Porcine HD Array (n = 116) porcine arrays (Illumina, San Diego, CA, USA). Those SNPs that displayed a minor allele frequency below 0.10, a call rate below 0.95, or that could not be mapped to Sus scrofa reference genome (Sscrofa11.1) 27 were filtered out. A total of 36,000 SNPs remained after data quality control. For each trait, a GWAS was performed by fitting a linear mixed model using GEMMA 28 , where phenotypes were adjusted for batch (11 levels) and IMF content as a covariate. The association of each SNP was tested using the Wald statistic considering the Bonferroni correction for multiple testing. Significance was set to a level of P ≤ 1.4 × 10 -6 . Two regions were found to be associated with C16:1n-7, one on SSC9 and another on SSC14, the latter corresponding to a reported polymorphism in the SCD gene that is known to affect MUFA and particularly C16:1n-7 24 . Candidate genes mapping within the SSC9 region were explored with Ensembl Genes Database using BioMart (https:// www. ensem bl. org/ bioma rt/ martv iew). Functional analyses like Gene Ontology and Reactome Pathway Enrichment Analysis were performed using Enrichr 29 . DGAT2 was retrieved as the most promising candidate gene for C16:1n-7 on SSC9.
Sequence variation in DGAT2. Variant discovery in porcine DGAT2 was examined by retrieving all sequence variants from the coding region, the 3′-untranslated region and 500 bp upstream on the proximal promoter of the gene (SSC9 10,031,627 to 10,068,464 bp) in a subset of 199 pigs of the same line with whole-genome sequencing data available. Sequenced pigs covered all representative sire families used in the experiment. DNA samples were submitted to Centre Nacional d' Anàlisi Genòmica (CNAG-CRG, Barcelona, Spain) for sequencing. Libraries were prepared and sequenced with paired-end reads with a NovaSeq 6000 instrument (Illumina) according to the manufacturer's protocol. Libraries were aligned to the Sscrofa11.1 27 using the BWA-MEM algorithm 30 . The average realized sequencing coverage was 6.8x (SD = 1.2×; min = 4.4×; max = 12.2×). Variants were identified following GATK HaplotypeCaller 3.8.0 software 31,32 . The SNP ss7315407085 G > A in exon 9 (SNP5, 9:10,065,826; Table 1) was selected for further validation with the whole set of pigs.
Genotyping DGAT2. All pigs (n = 1129) used in the experiment were genotyped for SNP ss7315407085 in exon 9 of DGAT2 using the primers described in Supplementary Table S1. Amplifications were performed by real-time PCR (QuantStudio3, Applied Biosystems, Thermo Scientific, Waltham, MA, USA) with High-Resolution Melt analysis (Luminaris Colour HRM Master Mix, Thermo Scientific) using 20 ng of genomic DNA and 0.4 µM of each primer in 5 µL final volume reaction. Thermocycling conditions were 50 °C 2 min, 95 °C 10 min, and 40 cycles of 95 °C 15 s, 60 °C 1 min, followed by a high-resolution melting curve starting with a denaturation at 95 °C for 15 s, annealing at 60 °C for 1 min and a slow ramp at 0.015 °C/s up to 95 °C. High Resolution Melt software v3.1 (Applied Biosystems, Thermo Scientific) was used for melting data analysis and sample genotyping. All pigs were also genotyped for the SCD (rs80912566 T > C; on SSC14) and leptin receptor (rs709596309 C > T; on SSC6) SNPs following the protocols described in Refs. 24,33 , respectively. DGAT2 expression. DGAT2 expression was measured by quantitative real-time PCR (qPCR) in the SM of 40 pigs from a single batch. RNA was isolated with TRI-Reagent according to the manufacturer's protocol and purity was assessed by spectrophotometry with a Nanodrop-1000. Total RNA (1.5 µg) was reverse-transcribed using SuperScript IV Reverse Transcriptase (Invitrogen, Thermo Scientific) with 100 µM random hexamers at 23 °C for 10 min, 50 °C for 20 min and 80 °C for 10 min. The cDNA was diluted 1:30 in water. Real-time PCR assays were carried out in triplicate using SYBR Green Supermix (Bio-Rad, Hercules, CA, USA), 0.2 µM of each www.nature.com/scientificreports/ primer and 3 µl of diluted cDNA. The B2M and RPL32 genes were used as reference genes to quantify and normalize the DGAT2 expression data. Primers used for DGAT2 expression (qDGAT2) are given in Supplementary  Table S1.
Estimation of DGAT2 tag-SNP effects. The effect of the DGAT2 ss7315407085 G > A SNP genotypes on production traits (body weight, carcass weight, BT and LT) and IMF content and fatty acid composition was estimated using a mixed model that included the batch (15 levels) and the genotype for DGAT2 (GG, AG and AA), SCD (TT, CT and CC) and leptin receptor (TT, CT and CC) as fixed effects and the sire and the dam as random effects. The leptin receptor polymorphism was included in the model because of its impact on fat content and fatty acid composition 33 . The slaughter age, for production traits, and IMF content, for fatty acids, were added as covariates. The additive and dominant effects of the DGAT2 SNP were estimated by replacing the genotype effect with the covariates (1, 0, − 1) and (0, 1, 0) for the GG, AG and AA genotypes, respectively. Gene expression was analysed only as a function of the DGAT2 genotype. The effects of the genotype and covariates were tested using the F-statistic. Multiple pairwise comparisons among DGAT2 genotypes were tested with the Tukey HSD test.
Results are presented as least-square means ± standard error. All the analyses were performed using the statistical package JMP Pro 14 (SAS Institute Inc., Cary, NC).

Consent for publication.
All the authors read and agree to the content of this paper and its publication. Table 1. Position, alelles and location of the variants identified in the genomic region from SSC9 10,031,627 to 10,068,464 bp, which includes coding region, 3′-untranslated region and 500 bp upstream on the proximal promoter of the DGAT2 gene. a n: Number of animals genotyped. b MAF: Minor allele frequency.

Results
The preliminary GWAS revealed two regions associated with C16:1n-7, one on SSC9 and another on SSC14 (Fig. 2). The region on SSC9 (9.8 Mb to 12.8 Mb) contained 18 significant SNPs (P < 1.4 × 10 -6 ). Candidate gene DGAT2 mapped within this region, specifically at 10.0-10.1 Mb. This association was not identified in previous analyses in this Duroc line 32 and evidence for potential implications of sequence variation in DGAT2 on fat content and composition in pigs is still very scarce 14,21 . On the other hand, the region on SSC14 matched the variant in SCD described above 24 . No evidence of association was found between the genomic region of DGAT1 (on SSC4) and fat content and fatty acid composition.
Sequence variation in DGAT2. Using whole-genome sequencing data from pigs of the same line, we identified a total of 33 sequence variants in DGAT2. Three of them were exonic, although synonymous (Table 1) and, from these, only SNP5 in exon 9 (ss7315407085) had a minor allele frequency greater than 0.25 (0.33).
With 1000 samples, which are around the current sample size of our biorepository, 0.25 is the minimum minor allele frequency required to detect (P < 0.05) a difference of around 5% (0.5 SD) between genotypes for C16:1n-7 with a power of at least 80%. This SNP has been previously reported in other pig breeds and crossbreds 14,21 and, in line with our results, the A allele was the minor allele in all genetic types, with a frequency that ranged from 0 to 0.32. The other 30 variants were located in the promoter region (2 SNPs) and in the 3′-UTR (21 SNPs and 7 indels). Using only the sequenced pigs, we confirmed that DGAT2 exerted an additive influence on C16:1n-7 (Supplementary Table S2), with SNP5 and the sequence variants in the interval between INDEL13 and SNP22 presenting the most relevant associations. Linkage disequilibrium of all variants was analysed with Haploview 4.2 software ( Supplementary Fig. S1), showing that SNP5 was in linkage disequilibrium with both INDEL13 (r 2 = 0.51) and SNP22 (r 2 = 0.47). In view of these results, we selected SNP5 (ss7315407085) as a tag variant of this haplotype for further validation.
Validation of DGAT2 tag SNP. The expression of DGAT2 by ss7315407085 genotype is shown in Fig. 3.
The DGAT2 SNP did not show dominance effects nor evidence for interaction with the SCD genotype for any of the fatty acids. The additive behaviour of the DGAT2-G allele was mantained across all SCD genotypes (Fig. 4). Thus, the maximum difference in C16:1n-7, which accounted for around 25% of the mean, occurred between the two extreme DGAT2/SCD haplotypes (from 4.05% ± 0.05, for GG/TT, to 3.23% ± 0.08, for AA/CC, P < 0.0001). Despite this, the difference between DGAT2-GG and AA genotypes was greater for the SCD-CC genotype (+ 0.32% ± 0.09, P < 0.001) than for the SCD-TT genotype (+ 0.15% ± 0.10, P = 0.13).
We did not find evidence of consistent effects of DGAT2 on body weight or carcass fat content (Table 4) or on the fatty acid composition of SF (Supplementary Table S3) or liver (Supplementary Table S4). Only a minor  (Table 4). Although DGAT2-GG pigs weighted 2.0 kg ± 0.7 (P < 0.05) less than DGAT2-AG pigs, the G allele did not show a clear additive behaviour (− 1.0 kg ± 0.6, P = 0.08).

Discussion
Understanding the regulation of fat metabolism is a milestone in the prevention and treatment of human lipid diseases, such as obesity, and in animal science. In this context, DGAT2 has received attention as the enzyme that catalyzes the final step in the synthesis of TG, the most common type of body fat. Here, we validated the effects of the ss7315407085 G > A SNP on fat content and fatty acid composition as a tag SNP for the haplotype in DGAT2 that was identified following a preliminary GWAS and after mRNA expression and sequence variant association analyses. Using a larger sample, we confirmed that the minor allele A segregates in the studied Duroc line at a moderate frequency (0.33), and that it has a negative impact on C16:1n-7 content in muscle. The Table 3. Least square means and additive values (± standard errors) for fatty acids content and fatty acid ratios in longissimus dorsi by DGAT2 genotype. A SFA: saturated fatty acids (C14:0 + C16:0 + C18:0 + C20:0); MUFA: monounsaturated fatty acids (C16:1n-7 + C18:1n-9 + C18:1n-7 + C20:1n-9); and PUFA: polyunsaturated fatty acids (C18:2n-6 + C18:3n-3 + C20:2n-6 + C20:4n-6). B Additive allele substitution of G for A. a,b,c Within trait, means with different superscripts differ significantly (P < 0.05). Bold font indicates statistical significance. Our validation scheme evidenced a consistent effect of DGAT2 SNP on C16:1n-7 in muscle that is independent of IMF content. A genetic variant in SCD with relevant effects on MUFA also segregates in this line 24 . The DGAT2 SNP had a lower impact on C16:1n-7 than the SCD SNP variant, particularly in terms of genetic variance. While the additive value of DGAT2 for C16:1n-7 was around half of that of SCD, the additive variance only reached 22% of that attributed to SCD. It has been shown that SCD co-expresses with DGAT2 34 and that their respective enzymes locate very close in the endoplasmic reticulum 11,17 , two findings that support the hypothesis that SCD is involved in the TG synthesis by providing a more accessible pool of MUFA through substrate channelling 12,17 . This also would explain why cells overexpressing DGAT2 show a greater proportion of MUFA and very particularly of C16:1n-7 34 . Our results are in line with this hypothesis, because C16:1n-7 was the fatty acid most affected by the DGAT2 SNP and the G allele, the one with higher DGAT2 expression, was the one that led to accumulate more C16:1n-7. However, we did not find evidence for a positive interaction between SCD and DGAT2 SNPs for C16:1n-7. Using RNAseq data from 40 other pigs of the same line, we were able to confirm that DGAT2 genotype (GG: n = 19, AG: n = 14, AA: n = 5) affected muscular DGAT2 expression, but we did not detect differential expression of SCD by DGAT2 genotype. In absence of interaction, the two genes mostly behaved in an allele-dose manner (Fig. 3).
The positive effect of the DGAT2 G allele on C16:1n-7 was not observed on C18:1n-9. Zhang et al. 34 evidenced that the impact of DGAT2 on C16:1n-7 was proportionally greater as fatty acid composition turned more dependent on de novo fatty acids. For C18:1n-9, the same effect was only seen in mature adipocytes, especially when they were cultured in a fatty acid-rich medium. In line with these results, a plausible explanation to the differential trend between C16:1n-7 and C18:1n-9 could be that the most immediate (short-term) effect of DGAT2 is to alter fatty acid composition by selectively capturing as substrate the first available de novo fatty acids rather than to enhance the endogenous biosynthesis of longer-chain fatty acids. This may be the case here, since differences in DGAT2 expression across genotypes were relatively small and no correlated responses in IMF content and SCD expression were observed. The favourable effect of the G allele on C14 and C16 fatty acids would reinforce this hypothesis. On the other hand, the impact of de novo fatty acids on final fatty acid composition is more apparent in C16:1n-7 and C18:1n-7 than in C18:1n-9, which is proportionally more abundant in pig feed than C16:1n-7 and thus more likely to be influenced by the diet. Similarly, the effect of DGAT2 on fatty acid composition was only observed in IMF, which is precisely the adipose tissue that develops later 35 and is less sensitive to dietary fat 36 .
Changes in DGAT2 expression are related to lipid accumulation. Results in human and murine cells indicate that DGAT2 promotes TG synthesis and storage in cytosolic lipid droplets 37 and accordingly DGAT2 is expressed more abundantly in cells with greater de novo fatty acid synthesis 19 . This is in line with findings in pigs, in which DGAT2 was most expressed in the liver of fatter breeds, where TG are expected to be more actively synthesized 38 . The expression of DGAT2 was twofold higher in fatter than in leaner breeds and around fivefold and tenfold higher in liver than in subcutaneous fat and muscle, respectively. However, the mRNA expression of DGAT2 was positively correlated with IMF content but not with BT. In view of these results, it would appear that the pigs carrying the DGAT2 A allele, which down-regulates the expression of DGAT2, should be less fatty and, specifically, have less IMF content. Even though this association was reported in other studies 14,21 , we did not find evidence that DGAT2 A entails a reduction in IMF content, BT or even fat content in liver. This may contradict biological expectations, but, in comparison with previous reports, validation here was performed using a much larger set of pigs from a single line so as to avoid breed or family biases. A 13-bp deletion in DGAT2 3′-UTR region has also been associated with increased DGAT2 mRNA expression and BT 39 . This deletion allele showed higher transcriptional activity, most likely owing to a less stable 3′-UTR secondary structure. This indel was located at 905 bp downstream from the stop codon and matches INDEL14 (Table 1). In line with the results observed for www.nature.com/scientificreports/ SNP5, none of these polymorphisms including INDEL14 (r 2 = 0.57, Supplementary Fig. S1) affected BT in our Duroc line (Supplementary Table S2).
The results obtained are more appealing for lipid research than for pig breeding. In recent years, C16:1n-7 has received a lot of attention for its potential role as a lipokine, i.e. as a lipid hormone that acts in distant organs. Despite this, the effects of C16:1n-7 are still under debate. On the one hand, preclinical experiments with cell and rodent models show that C16:1n-7 supplementation has anti-inflammatory properties 40,41 that protect against metabolic disorders. But, on the other hand, research in humans reported elevated blood levels of C16:1n-7 in patients with obesity and metabolic syndrome 42,43 . Findings in pigs may contribute to disentangle the biological implications of C16:1n-7 as well as to understand the genetic regulation of fat metabolism.

Conclusions
A sequence variant in the porcine DGAT2 influences muscular gene expression and fatty acid composition, with changes mostly limited to swapping C18 fatty acids for C14 and C16 fatty acids of the same degree of saturation. Differences across DGAT2 genotypes are quantitatively small, but the variant that enhances DGAT2 expression is associated with increased content of shorter-chain fatty acids, thereby indicating that DGAT2 preferentially uptakes the first available de novo fatty acids as substrate. Although DGAT2 and SCD act additively, the fact that DGAT2 overexpression impacts C16:1n-7 more than C16:0 indicates that DGAT2 and SCD are closely interlinked. Further enzyme characterization is needed, but our findings provide evidence that C16:1n-7 is the immediate substrate for DGAT2 and corroborate that fatty acid metabolism in pigs is subjected to subtle tissue-specific genetic regulatory mechanisms.

Data availability
Data that support the findings of this study are available within the article and Supplementary Information, or from the authors upon reasonable request.