A healthy 1-day-old calf was born and raised at the animal farm of Gansu Agricultural University (LanZhou, China) and humanely slaughtered in a slaughterhouse, followed by the collection of the perirenal adipose tissue along with heart, spleen, kidney, liver, lung, and leg muscle tissues. The study protocol was approved by the Gansu Agricultural University, China, as well as the Ministry of Science and Technology of the People's Republic of China (Approval number 2006–398). Animal experimentation, including sample collection, was performed in agreement with the ethical commission's guidelines of Gansu Agricultural University. Furthermore, the experimental protocol complied with the local animal welfare guidelines.

The perirenal adipose tissue of 1-day-old calf (from the animal farm of Gansu Agricultural University) was collected and shredded in PBS, followed by the addition of type I collagenase (200 U/mL). The solution was shaken at 37°C for 1 h, and enzymatic digestion was stopped by the addition of an equal volume of complete medium. The digested product was serially filtered through 100 μm and 40 μm mesh filters, followed by centrifugation (1,500 rpm, 10 min). The erythrocyte lysate was mixed with the pellet and centrifuged (1,500 rpm, 10 min). After discarding the supernatant, a complete medium was added to the pellet, mixed, and centrifuged twice (1,500 rpm, 10 min). The pellet was resuspended, and the cells were seeded in a 10 cm cell culture dish and cultured in 5% CO₂ at 37°C. The medium was refreshed after 12 h to remove the non-adherent cells and was refreshed every 48 h. At 90% confluence, the cells were removed with trypsin and frozen in liquid nitrogen.

Bovine pre-adipocytes were grown and passaged to the F3 generation, and cells were isolated on differentiation days 0, 2, 4, 6, and 8. We performed qRT-PCR to assess ACSL1 expression at different stages of differentiation.

The complete sequence of bovine ACSL1 (GenBank No. 537161) was used to design and synthesize three siRNA sequences (Guangzhou Ruibo Biological Company). Table 1 shows the sequence information.

Next, the synthesized siRNA-ACSL1 and the corresponding siRNA-Non-specific Control (siRNA-NC) were transfected into bovine adipocytes. After 48 h of culture, we performed qRT-PCR to determine ACSL1 mRNA expression to screen the siRNA-ACSL1 with the highest interference efficiency.

Next, we sequenced the bovine adipocytes in which ACSL1 expression was silenced and those that exhibited normal expression.

We used TRIzol reagent to isolate total RNA from siRNA-ACSL1 or siRNA-NC bovine adipocytes, followed by strict quality control of the RNA samples. First, 1% agarose gel electrophoresis (AGE) was used to analyze RNA integrity and the presence of DNA contamination. RNA purity was measured using an IMPLEN NanoPhotometer® spectrophotometer and quantified using an RNA analysis kit and a Qubit®2.0 Fluorometer. Finally, the RNA Nano 6000 Assay Kit was used to assess RNA integrity using the Agilent Bioanalyzer 2100 system.

The cDNA library was prepared using the NEBNext Ultra Directional RNA Library Prep Kit for Illumina sequencing using a combination of equal amounts of total RNA of two parallel bovine adipocytes cultures in equal ratios of siRNA-ACSL1/siRNA-NC. The total RNA was from 6 dishes of bovine adipocytes, of which 3 dishes were treated with ACSL1 gene interference, and the other 3 dishes were the control group. The Epicenter Ribo-zero rRNA Removal kit was used for the removal of the rRNA, followed by the construction of six libraries, which were labeled si_ACSL1_1, si_ACSL1_2, si_ACSL1_3, NC_ACSL1_1, NC_ACSL1_2, and NC_ACSL1_3, respectively. Briefly, poly-T oligo-attached magnetic beads were used to obtain mRNA, followed by fragmentation using divalent cations at high temperatures in NEBNext First Strand Synthesis Reaction Buffer (5X). First-strand cDNA was prepared using a random hexamer primer and M-MuLV Reverse Transcriptase (RNaseH-), followed by the application of DNAPI and RNase H. These activities facilitated the conversion of the remaining overhangs into blunt ends. The NEBNext adaptor with hairpin loop structure was ligated to prepare for hybridization after adenylation of the 3′ ends of the DNA fragments. The Beckman Coulter AMPure XP system was used to obtain cDNA fragments of 150 ~ 200 bp in length. These fragments were treated with 3 μL of USER Enzyme at 37°C for 15 min, followed by 5 min at 95°C before PCR. Next, Phusion High-Fidelity DNA polymerase, Universal PCR primers, and Index (X) Primer were used for PCR analysis. Finally, the Agilent Bioanalyzer 2100 system was used to assess the quality of the library generated. The index-coded samples were clustered using a cBot Cluster Generation System for TruSeq PE Cluster Kit v3-cBot-HS (Illumina). An Illumina HiSeq PE150 platform was used to sequence the library at Novogene Bioinformatics Institute (Beijing, China), resulting in the generation of data on the 150 bp paired-end reads.

In-house Perl scripts were used to analyze the raw reads. Clean reads were obtained after filtering the raw data for low-quality reads/adaptors/polyNs, and at the same time, we obtained the Q20, Q30, and GC information. The STAR (v2.5.1b) software (32) was used to create a list of the reference genome, and the corresponding paired-end clean reads were aligned to the cattle reference genome (Bos_taurus_Ensembl_94). We used HTSeq v0.6.0 to assess the read counts mapped to each gene. The gene lengths and read counts aligned to the genes were calculated using FPKM. Normalization of the expression of lncRNAs/mRNAs removed the impact of sequencing depth, gene length, and sample difference on gene expression. The R package DESeq2 was used to identify DE-lncRNAs (DELs) and protein-coding genes (33) following a negative binomial distribution. Genes and lncRNAs with an adjusted P < 0.05 (P-adjust) were labeled as differentially expressed.

Cuffcompare software was used to screen the lncRNA transcripts. The assembled transcripts were filtered to determine probable lncRNAs to reduce the false-positive rate using the following steps: (1) New transcripts were identified by comparison with reference annotations; (2) Transcripts with exon numbers ≥ 2 and lengths > 200 nt were selected for the exon number screening and transcript length screening of the transcripts; (3) Transcripts were assessed for coding potential by integration using coded potential calculator 2 (CPC2) (34), Pfam (35), and Coding-Non-Coding-Index (CNCI) (36) to predict the coding potential and the intersection of the non-coding transcripts, allowing selection of the candidate lncRNA dataset; (4) The candidate novel lncRNAs were finally screened using the HUGO Gene Nomenclature Committee (HGNC) naming guidelines for long non-coding RNA (37) and named based on their positional relationships with coding genes.

The target genes for the lncRNAs in the RNA-seq results were predicted based on cis-and trans-regulation. We investigated the coding genes within 100 kb of the lncRNAs to determine their target genes (38). The trans role implied an expression-based association between the lncRNA and target genes. We determined the correlation between the lncRNAs and coding genes with custom scripts, and the genes with Pearson correlation coefficient > 0.95 were identified as the target genes of the lncRNAs (39).

The GOseq R package was used to perform Gene Ontology (GO) analysis of the target genes of differentially expressed mRNAs (DEMs)/lncRNA target genes, with corrected gene length bias (40). The KEGG database was used to evaluate the functional significance of the biological systems (cells, organisms, and the ecosystem) based on genome sequencing and other high-throughput experimental data. The cluster Profiler R package was used for KEGG pathway analysis of the DEMs/lncRNA target genes. A P < 0.05 indicated a significant relationship between the terms and the DE-lncRNA target genes.

Next, we obtained the differentially expressed lncRNAs (DELs) and their target genes from the cis and trans prediction results, respectively, to further explore the interactions between the DEMs and DELs. These were incorporated with the DEMs from the sequencing results. We screened the lipogenic-related mRNAs based on literature searches to improve the visibility of the results and finally constructed a DEL-DEM co-expression network using Cytoscape (v3.6.0) (41).

We randomly selected 12 significant DEMs and 7 DELs and analyzed their expression using qRT-PCR. The Evo-MLV kit was used to synthesize the RNA extracted from bovine adipocytes. First, we used Primer 5.0 to design specific primers for mRNAs and lncRNAs. Tables 2, 3 show the primers of DEMs and DELs used, respectively. The PCR products were analyzed using 1.5% AGE and then sequenced using the Sanger method. Then, we verified the mRNAs/lncRNAs from bovine adipocytes by comparing the amplified sequences of the PCR products with the sequences of the specific mRNAs and lncRNAs obtained from RNA-Seq using MEGA 5.0. Second, the Applied Biosystems QuantStudio®6 Flex was used to perform qRT-PCR using the TB Green® Premix Ex Taq™ II, and the bovine GAPDH gene as the internal reference. The 2^−ΔΔCT method was used to calculate relative expression with GAPDH as the internal control (42). Also, the expression of mRNAs and lncRNAs from both the techniques were expressed as Log₂ (fold change) for the si-treated group in comparison with the NC-treated group. Paired t-tests using GraphPad Prism v6 were used for data analysis. Data were expressed as the mean ± SE deviation of triplicate experiments. A P < 0.05 indicated a significant difference.

Next, the lncRNA levels were determined in various bovine tissues to further verify the role of lncRNA in the production of bovine UFAs. Briefly, we used TRIzol to isolate total RNA from the collected heart, spleen, lung, liver, kidney, leg muscle, and perinephric fat, followed by reverse transcription to generate cDNA. Finally, the expression of lncRNAs was determined using qRT-PCR. The 2^−ΔΔCT method was used to determine the relative expression using GAPDH as the internal control. The paired t-test using GraphPad Prism v6 was done for data analysis. Data were expressed as the mean ± SE deviation of triplicate experiments.

The cultured F3 generation bovine pre-adipocytes were collected on differentiation days 0, 2, 4, 6, and 8, and qRT-PCR was performed to determine the relative expression of ACSL1. The rise in ACSL1 levels was followed by a reduction, with the maximum expression on day 4 (P < 0.01). Therefore, day-4 bovine adipocytes were used for subsequent experiments (13).

We transfected three pairs of siRNAs into the bovine adipocytes and determined their expression after 48 h to screen the most efficient siRNA fragments for silencing ACSL1. There was a >70% reduction in the mRNA expression of ACSL1, with si3-ACSL1 showed the highest efficiency of interference (86%; P < 0.01). Therefore, si3-ACSL1 was selected for subsequent experiments (13). Next, we transfected differentiation day 4 adipocytes with si3-ACSL1 (si-treated group) (13), followed by sequencing of both the si-treated cells and NC-treated cells.

We constructed six cDNA libraries for the NC-treated group (n = 3) and si-treated group (n = 3) to explore the regulation of UFA synthesis by mRNAs and lncRNAs. After sequencing, the raw reads were acquired using the Illumina HiSeq PE150 Platform and deposited in GenBank (accession numbers: SRR13358757-13358762). Next, we tested the correlation between the expression of different samples (Figure 1A). In total, we obtained 98,462,108–117,507,866 and 93,238,078–110,700, 276 raw reads from the libraries from the NC-treated group and si-treated group, respectively, which resulted in 96,857,744–115,542,046 and 91,106,662–109,206,440 clean reads after removing low-quality/adaptor sequences, resulting in ~15 GB data/sample. The GC contents were 51.54–56.29%. Q20 was >96%, and Q30 was >90% (Table 4). More than 87.76% of clean reads aligned to the reference genome, resulting in 173,090 assembled transcripts (Supplementary File 1). Our study identified 5,432 putative lncRNAs, which included 65.2% lincRNAs, 20.7% antisense lncRNAs, and 14.1% sense_overlapping lncRNAs (Figures 1B,C). At the transcriptional level, 15,547 mRNAs were obtained, of which 454 were novel mRNAs.

We analyzed the gene structure and expression to determine the dissimilarities between the lncRNAs and mRNAs obtained by RNA-seq. The results showed that mRNAs were expressed at a relatively higher degree with the lncRNAs, based on the FPKM values (Figure 2A). The mRNAs (5,321 bp on average) were substantially longer than the lncRNAs (1,813 bp on average) (Figure 2B). The mRNAs had a higher number of exons (6.04 on average) than those of the lncRNAs (2.91 on average). Also, 61.23% of mRNAs had four or more exons, while 81.41% of lncRNAs had three or fewer exons (Figure 2C). Furthermore, most of the mRNAs (402 bp on average) had longer ORFs than the lncRNAs (147 bp on average) (Figure 2D).

We detected 75 DELs and 107 DEMs in the si-treated group compared with the NC-treated group (P-adjusted <0.05). We performed differential mRNA cluster analysis to examine the clustering pattern of DEMs under varying experimental conditions. We found that the same group of DEMs was clustered together (Figure 2E), supporting the accuracy and reliability of the samples. Furthermore, we found that in the si-treated group, 31 lncRNAs were significantly upregulated, and 44 lncRNAs were significantly downregulated (P-adjusted <0.05). Figure 2F shows the clustering results. Similarly, we observed significantly upregulated expression of 47 mRNAs and downregulated expression of 60 mRNAs in the si-treated group (P-adjusted <0.05). The volcano plots (Figures 2G,H) and the Supplementary Files 2, 3 show the DEL and DEMs, respectively.

The results of the GO analysis revealed that the significant DEMs were classified into 386 functional groups (P < 0.05; Supplementary File 4). Of these, 64 terms were significantly enriched in molecular function (GO-MF), and metalloendopeptidase inhibitor activity (GO: 0008191; P = 0.0017) was the most significantly enriched GO term, followed by EP4 subtype prostaglandin E2 receptor binding (GO: 0031867; P = 0.0019) and oxidoreductase activity (GO: 0016671; P = 0.0034). For the cellular component (GO-CC), 33 terms were significantly enriched, with the extracellular region (GO: 0005576; P = 0.0001) the most significantly enriched, followed by the invadopodium membrane (GO: 0071438; P = 0.0004). For the biological processes (GO-BP), 289 GO terms were related to various processes or regulation, such as glycosaminoglycan metabolic process (GO: 0030203; P = 0.0009), aminoglycan metabolic process (GO: 0006022; P = 0.0014), and positive regulation of MAPK cascade (GO: 0043410; P = 0.0055). Figure 3A shows the top 50 functional GO annotations for the DEMs. The KEGG pathway analysis showed that the DEMs were involved in 116 signaling pathways (Supplementary File 5). Figure 3B shows the top 20 signaling pathways, and the PPAR signaling pathway showed the highest level of significance (P < 0.05), followed by the Rap1 signaling pathway and PI3K-Akt signaling pathway. Our results indicated that these pathways were probably involved in UFA synthesis.

We identified 20,724 protein-coding genes as the nearest neighbors of 5,423 out of 5,432 lncRNAs using 100 kb as the cutoff (Supplementary File 6). Subsequent GO analysis revealed 435 significantly enriched GO terms (P < 0.05; Supplementary File 7). Figure 4A shows the top 50 GO functional annotations of the DEL target genes. The top five GO terms were embryonic skeletal system development, anterior/posterior pattern specification, sequence-specific DNA binding, regionalization, and embryonic skeletal system morphogenesis. The results of KEGG analysis indicated the involvement of seven pathways (P < 0.05; Supplementary File 8), and Figure 4B shows the top 20 KEGG pathways, such as MAPK signaling pathway, PPAR signaling pathway, Ribosome biogenesis in eukaryotes, and ether lipid metabolism.

Next, we identified 5,432 lncRNAs and 698,719 genes to investigate the trans target genes of lncRNAs using the Pearson correlation ≥ 0.95 as the cutoff (Supplementary File 9). Figure 4C shows the top 50 GO functional annotations of the DEL target genes. GO analysis identified 1,154 significantly enriched GO terms (Supplementary File 10), and the top five GO terms were enzyme binding, kinase regulator activity, rRNA binding, G-protein coupled nucleotide receptor activity, and G-protein coupled purinergic nucleotide receptor activity. We also found that 274 KEGG pathways were involved in the trans-target genes for lncRNAs (Supplementary File 11), such as the Ribosome, Citrate cycle (TCA cycle), protein digestion and absorption, FA metabolism, and FA degradation (Figure 4D).

We obtained the DELs and their targets from the cis and trans prediction results, respectively, and incorporated them with the DEMs from the sequencing results (Supplementary File 12). Then, we screened lipogenic-related mRNAs based on literature searches and finally constructed a DEL-DEM co-expression network (Supplementary File 13). The DEL-DEM co-expression network involved 14 DELs and 22 trans-targets (Figure 5A), as well as 2 DELs and 2 cis-targets (Figure 5B). The co-expression network results showed that some lncRNAs interacted with multiple DEMs; for example, 14 DEMs co-expressed with TCONS_00069661, and 2 DEMs co-expressed with TCONS_00142457, TCONS_00040771, and TCONS_00074138, respectively, indicating that these lncRNAs and mRNAs belonged to the core non-coding/coding RNAs and had important regulatory effects on the synthesis of UFAs.

We randomly chose 12 and 7 genes from DEMs and DELs, respectively, to analyze the validity of the RNA-seq data by qRT-PCR. The results were comparable to the results of mRNA and lncRNA expression derived from the RNA-seq data (Figures 6A,B). Furthermore, the sequence of the qRT-PCR-amplified fragments of the specific primers were consistent with the sequences obtained by RNA-seq (Supplementary File 14).

We used qRT-PCR to study the expression of lncRNA TCONS_00069661, located at the center of the co-expression network and expressed in the kidney, spleen, heart, liver, lung, leg muscle, and perirenal fat of cattle. We found that TCONS_00069661 was highly expressed in the adipose tissue of cattle (Figure 6C). Moreover, our previous studies showed that ACSL1 gene was also highly expressed in bovine adipose tissue (12). Thus, these results substantiated the reliability of the co-expression network, and suggested the regulatory roles played by DELs and DEMs in UFA synthesis.