All experimental procedures in this study were consistent with the National Research Council Guide for the Care and Use of Laboratory Animals. All protocols were approved by the Institutional Animal Care and Use Committee at Inner Mongolia University.

Three fat-1 transgenic cattle fetuses (FT_1, FT_2 and FT_3) and three wild-type Luxi cattle fetuses (WT_1, WT_2 and WT_3) at 50–60 days were obtained from Shengmu Experimental Base (Holinger District, Hohhot). The fetuses were rinsed with normal saline, placed in phosphate-buffered saline (PBS) containing dual antibodies (15140-122, Gibco, United States), placed in an icebox, and brought back to the laboratory for analysis.

The head, limbs, internal organs, and bones of the fetus were removed and discarded. The remaining tissues were rinsed three times with PBS and then once with 75% ethanol for 15 s and transferred back to PBS. Then, the samples were sectioned into 1–3-mm slices and placed in Petri dishes. collagenase IV at a concentration of 1 mg/ml was added; cultures were maintained at 38.5°C with 5% CO₂ for 2–3 h. After the tissues were dissolved, the mixtures were centrifuged at 1,500 rpm/min for 5 min. The cell pellet was collected, resuspended in DMEM F12 (11330107, Gibco, United States) supplemented with 20% FBS (10099141, Gibco, United States), and cultured in 60-mm dishes at 38.5°C with 5% CO₂. bEFs were passaged when they reached 90% confluence.

Fatty acids were analyzed from bEFs when they reached 80% confluency. Briefly, bEFs were collected using 0.05% trypsin (25300062, Gibco, United States). Digested cells were placed into PBS-containing centrifuge tubes and centrifuged at 1,500 rpm for 5 min. Next, 1 ml 2.5% H₂SO₄/methanol solution was added to the isolated supernatant, which was then incubated in an 80°C water bath for 90 min. The solution was allowed to cool to room temperature before 1.5 ml 0.9% NaCl solution and 1 ml n-hexane was added. The solution was vigorously mixed for 5 min and then centrifuged at 1,500 rpm for 5 min. The organic phase was transferred to a new centrifuge tube. Saturated KOH-methanol solution was added, vigorously mixed for 5 min, and centrifuged at 1,500 rpm for 10 min. The upper liquid was collected into gas chromatography–mass spectrometry (GC-MS) sample vials for fatty acid analysis. The gas chromatographic conditions were as follows: carrier gas helium, chromatographic column HP-88, constant linear velocity 20.0 cm/s, separation ratio 20.0%, injection volume 1 μL, heating program 60°C for 1 min, 40°C/min heating to 140°C, maintained at 140°C for 10 min, increased to 240°C at 4°C/min, and maintained at 240°C for 15 min.

Total RNA of bEFs was extracted using Trizol reagent (15596026, Invitrogen, United States). Poly (A) RNA was purified from total RNA (5 μg) by poly-T oligo-attached magnetic beads and then broken into small fragments. The spliced RNA fragments were reverse-transcribed to construct the final cDNA library. Paired-end sequencing was performed in LC Science using an HiSeq 4000 (Illumina, United States). The read values of the samples were aligned with the UCSC (http://genome.ucsc.edu/) reference genome using the HISAT package. This package initially removed a portion of the reads based on the quality information accompanying each read value, and then mapped the reads to the reference genome. All the transcripts in the samples were combined to reconstruct a comprehensive transcriptome using Perl Scripts. StringTie and edgeR were used to estimate expression levels of all transcripts. The differentially expressed mRNAs and genes with log2 (fold change) > 1 or log2 (fold change) < −1 and with statistical significance (p value <0.05) were screened using R.

The activity of key restriction enzymes and substrate contents in glycolysis, TCA cycle, and OXPHOS were measured using commercial kits following the manufacturer’s instructions. The following commercial kits, all obtained from Comibio Biotechnology, China: the content of ATP (ATP-1-Y), glucose (PT-1-Y), FDP (FDP-1-G), PA (PA-1-Y), CA (CA-1-Y), NADH (NAD-1-Y) and the activity of HK (HK-1-Y), PFK (PFK-1-Y), PK (PK-1-Y), ICDHm (ICDHm-1-Y), α-KGDH (KGDH-1-Y), SDH (SDH-1-Y), MDHm (MDHm-1-Y), mitochondrial complex I (FHTA-1-Y), complex III (FHTC-1-Y), complex IV (FHTD-1-Y), complex V (FHTE-1-Y). Briefly, 10⁷ cells were placed in a homogenizing tube containing 1 ml kit extract with ∼15 ceramic beads. The solution was homogenized at a low temperature with a homogenizer (Bertin, France). The supernatant was then added to a 96-well plate along with other reagents, as instructed. The absorbance was recorded with a microplate spectrophotometer (Thermo, United States) to calculate the enzyme activity or metabolite concentration.

Total RNA was isolated from tissues using the RNAiso Plus kit (9109, Takara, Japan) according to the manufacturer’s instructions. Briefly, RNA was then reverse-transcribed into cDNA using a cDNA reverse transcription kit (RR036A, Takara, Japan). RT-qPCR was performed using a real-time PCR detection system (7500, ABI, United States). The reaction mixture (20 μL) contained 1 μL cDNA, 0.5 μL of each primer, and 10 μL TB Green Supermix (RR820A, Takara, Japan). The primers used are listed in Supplementary Table S1. The protocol was as follows: initial denaturation step at 95°C for 30 s and then 40 cycles of [95°C for 15 s, 60°C for 30 s]. Relative abundance was quantified using the 2^−ΔΔCt method. RPLP0 was used as a housekeeping gene.

Total DNA was extracted using a QIAamp Fast DNA Tissue Kit (51404, Qiagen, Germany) following the manufacturer’s protocol and quantified using a spectrophotometer and the A260/A280 ratio. DNA samples were fragmented using sonication and subjected to bisulfite conversion. The Accel-NGS Methyl-Seq DNA Library Kit (Swift, MI, United States) was utilized for attaching adapters to single-stranded DNA fragments. Briefly, the Adaptase step is a highly efficient, proprietary reaction that simultaneously performs end repair, tailing of 3′ ends, and ligation of the first truncated adapter complement to 3′ ends. The Extension step incorporates the truncated adapter 1 by a primer extension reaction. The Ligation step is used to add the second truncated adapter to the bottom strand only. The Indexing PCR step increases the yield and incorporates full-length adapters. Bead-based Solid Phase Reverse Immobilization (SPRI) clean-ups are used to remove both oligonucleotides and small fragments, as well as to change enzymatic buffer composition. Finally, the pair-end 2 150 bp sequencing was performed in LC Science using an HiSeq 4000 (Illumina, United States).

DNA methylation and hydroxymethylation were quantified using the MethylFlash methylated and hydroxymethylated colorimetric DNA quantification kits (P-1034; P-1036; Epigentek, United States), respectively, according to the manufacturer’s instructions. The percentage of methylated DNA was proportional to the optical density intensity measured. All samples were obtained in triplicates.

BEFs were rinsed with PBS and lysed in 150 μL of ice-cold Radio Immunoprecipitation Assay buffer. Proteins were extracted from tissues with cell lysis buffer, boiled for 5 min, and stored at −80°C. Samples were separated using a 10% sodium dodecyl sulfate-polyacrylamide gel electrophoresis (SDS-PAGE), and then electrically transferred to polyvinylidene fluoride membranes blocked with 5% nonfat dry milk in TBS-Tween 20 (blocking buffer) for 1 h. The membrane was subsequently incubated with the following antibodies, all obtained from Santa Cruz Biotechnology, United States unless otherwise stated: anti-DNMT1 (1:500; sc-271729), anti-DNMT3a (1:500; sc-373905), anti-DNMT3b (1:500; sc-393279), anti-TET1 (1:500; sc-293186) and anti-TET2 (1:1,000; ab94580; Abcam, United States) and anti-TET3 (1:1,000; ab139311; Abcam, United States) or with a polyclonal antibody against α-Tubulin (11224-AP, proteinch, China) at 1:1,000 in TBST at 4°C overnight. Membranes were then washed three times and incubated with a horseradish peroxidase-conjugated goat anti-mouse and anti-rabbit secondary antibody at 1:5,000 in blocking buffer.

The CpG islands and CpG sites in the promoters of PKM, PFKM, IDH2, MDH2, NDUFS1 and ATP5F1C were examined using MethPrimer software online (http://www.urogene.org/cgi-bin/methprimer/methprimer.cgi).

The MeDIP analysis was performed according to the manufacturer’s instructions (ab117133; Abcam, United States). Sonication of total genomic DNA (three pulses for 10–12 s at level 2 with 30–40 s intervals between pulses while resting on ice) produced DNA fragments between 200 bp and 1,000 bp in length. Single-stranded DNA was produced by heating and denaturing 1 µg of DNA. Then, 1 μL of anti-5mC antibody (or normal mouse IgG for a negative control) was added to induce immunoprecipitation at room temperature for 2 h. Proteinase K was added and the mixture incubated at 65°C for 1 h to release mDNA from the DNA/antibody complex. Finally, the mDNA was captured, eluted, and detected by real-time quantitative PCR. The qPCR primers are listed in Supplementary Table S2.

GraphPad Prism 8.0 software was used for statistical analyses. The fatty acid contents of each tissue and the enzymatic activities and mRNA expressions in each group of cells were expressed as mean ± SEM. The mean of fat-1 transgenic and wild-type cattle samples was compared by one-way analysis of variance (ANOVA). p < 0.05 was considered statistically significant.

Both DNA and mRNA of exogenous fat-1 gene were significantly expressed in the transgenic bEFs, but not detected in WT cells (Supplementary Figure S1). n-3 PUFA were quantified using gas chromatography; % total n-3 PUFA (sum of C18:3n3, C20:3n3, C20:5n3 and C22:6n3) were significantly increased in FT groups in comparison to WT groups (17.149 ± 0.070 vs. 20.764 ± 0.761, p = 0.009, Figure 1A), while % total n-6 PUFA (sum of C18:2n6t, C18:2n6c, C18:3n6, C20:3n6 and C20:4n6) were significantly reduced in FT cells (20.976 ± 1.377 vs. 16.591 ± 0.102, p = 0.034, Figure 1B). Subsequently, the n-6/n-3 ratio was significantly lower in the FT group (Figure 1C). These results suggest the exogenous fat-1 gene impacts n-6 and n-3 PUFA synthesis and subsequently alters the balance of n-6/n-3 PUFA.

RNA-seq analysis was performed on FT and WT cells. Out of a total of 27,386 differential genes, 1,157 genes were significantly up-regulated, and 1,200 genes were significantly down-regulated (Figure 1D). The KEGG pathway identified ECM-receptor interaction, axon guidance, and MAPK signaling pathways were significantly upregulated (Figure 1E). Among the pathways identified by the KEGG pathway that were significantly downregulated included OXPHOS, thermogenesis, and the TCA cycle (Figure 1F). GO enrichment analysis suggested the upregulated genes were mainly related to cell migration, positive regulation of cell migration, and substrate adhesion-dependent cell spreading (Figure 1G). Downregulated genes identified included enriched cytoplasm, oxidation-reduction process, single-stranded DNA binding. However, many genes related to ATP synthesis, oxidation-reduction process, and the proton-transporting ATP synthase complex were also downregulated (Figure 1H). These results suggest that the transgenic fat-1 gene not only changes the fatty acid composition to decreased the rato of n-6: n-3 PUFA, but also regulates and influences energy metabolism, possibly through the ATP synthesis pathway.

ATP levels in FT and WT cells decreased by approximately 39% (30.101 ± 0.599 vs. 18.318 ± 0.762 nmol/min/10⁷ cells, p < 0.0001, Figure 2A), while glucose content increased by 47% (10.265 ± 1.633 vs. 19.347 ± 3.207 μmol/10⁷ cells, p = 0.032, Figure 2A). The activity of key glycolytic rate-limiting enzymes hexokinase (HK), phosphofructokinase (PFK) and pyruvate kinase (PK) in FT cells were significantly lower than those in WT cells (Figure 2B). The catalytic product of PFK, fructose-1,6-diphosphate (FDP) was also significantly reduced (0.700 ± 0.011 vs. 0.528 ± 0.021 mg/10⁷ cells, p < 0.0001, Figure 2A). Pyruvate (PA), a product of PK, was also significantly lower in FT cells (1.897 ± 0.039 vs. 1.673 ± 0.046 μg/10⁷ cell, p = 0.004, Figure 2A). RT-qPCR further confirmed that the mRNA expression levels of HK1, PFKM, and PKM subunits were significantly down-regulated (Figure 2C).

Detection of citric acid (CA) and NADH, products of the TCA cycle, were significantly lower in FT cells (0.425 ± 0.008 vs. 0.387 ± 0.002 nmol/10⁷ cells, p = 0.0009 and 1.151 ± 0.104 vs. 0.620 ± 0.056 nmol/10⁷ cells, p = 0.001, Figure 2A). In FT cells, the activity of mitochondrial isocitrate dehydrogenase (ICDHm) and mitochondrial malate dehydrogenase (MDHm), key rate-limiting enzymes in the TCA cycle, were significantly reduced in FT cells when compared to WT (Figure 2D). RT-qPCR results also confirmed that IDH2 of ICDHm and MDH2 were significantly downregulated (Figure 2E). There was no significant difference between the activity of α-ketoglutarate (α-KGDH) (Figure 2D) and the expression of OGDH encoding it (Figure 2E). Although there was no significant change in succinate dehydrogenase (SDH) activity (Figure 2D), the expressions of SDHA, SDHB and SDHD were downregulated to varying degrees (Figure 2E).

NADH dehydrogenase (complex I), pigment reductase (complex III), cytochrome oxidase (complex IV), and ATP synthase (complex V) participate in OXPHOS. Subsequently, these enzymes were examined and found (with the exception of complex IV) to be significantly downregulated (Figure 2F), which was additionally confirmed by RT-qPCR (Figure 2G).

Previous studies have shown exogenous addition of n-3 PUFA increased the overall degree of DNA methylation in adipocytes (Perfilyev et al., 2017). Therefore, we performed global methylation on FT cells compared to WT cells. Higher levels of % total DNA methylation was observed in FT cells in comparison to WT (Figures 3A,B), although levels of 5-hmc methylation was not statistically significant. Genome-wide DNA-methylation analysis also identified differences between the two groups (Figure 3C). KEGG analysis revealed that approximately 80% of the differential genes were enriched in metabolic pathways, and significant changes occurred in DNA methylation of genes related to thermogenesis, fatty acid elongation, TCA cycle, and OXPHOS between FT and WT cells (Figures 3D,E).

Methylation of genes associated with the key limiting enzymes in glycolysis, TCA cycle, and OXPHOS were all significantly increased in FT cells compared with WT cells (Figure 4A). PFKM and PKM in glycolysis, IDH2 and MDH2 in the TCA cycle, and NDUFS1 and ATP5F1 in OXPHOS were all hypermethylated at their promoters’ CpG island in FT cells (Figures 4A–C), correlating to fat-1 transgene. This suggests balance of n-6: n-3 PUFA ratios could play a role in energy metabolism-associated gene expression. We then analyzed activity of the methylases DNMT1, DNMT3A, DNMT3B, TET1, TET2, and TET3. As expected, expression of DNMT1 and DNMT3B increased, expression of TET2 and TET3 significantly decreased, and expression DNMT3A and TET1 did not change (Figures 4D–F). These results further suggest that exogenous fat-1 transgene in the FT group can influence hypermethylation of genes related to energy metabolic rate-limiting enzymes by upregulating the methyltransferases DNMT1 and DNMT3B and downregulating TET2 and TET3 to reduce the demethylation of genes encoding these enzymes.