All experimental procedures were performed following the Guide for Animal Care and Use of Laboratory Animals from the Institutional Animal Care and Use Committee at China Agricultural University. The experimental protocol was approved by the Department Animal Ethics Committee at China Agricultural University (Permit No. DK996).

Three Landrace piglets from a pig breeding farm in Ninghe (Tianjin, China) were obtained for the present study. The piglets were euthanized at the age of 7 days through an intraperitoneal injection of pentobarbital sodium (50 mg/kg body weight) and exsanguination. Subcutaneous adipose tissue samples were carefully resected under sterile conditions. The adipose tissue was digested with 0.1% type I collagenase (Sigma, Beijing, China) for approximately 2 h at 37°C and then centrifuged at 1,000 × g for 8 min (Zhou et al., 2007). The resulting digestion mixture was successively filtered through 100-µm mesh filters by centrifuging for 8–10 min to obtain preadipocyte pellets. The cell pellets were then resuspended in Dulbecco’s modified Eagle’s medium/nutrient mixture F-12 (DMEM/F12) containing 10% fetal bovine serum (FBS) and plated in glass culture dishes (Gibco, NY, United States). The culture medium was replaced every other day.

When the cells reached 90% confluence, they were transferred to six-well plates (Corning Costar, NY, United States) at a density of 3 × 10⁶ cells/ml in 2 ml per well, incubated at 37°C with 5% CO₂ and 95% O₂, and cultured until 90% confluence. The standard culture medium was then replaced with adipogenic induction medium (DMEM containing 10% FBS, 0.5 mM 3-isobutyl-1-methylxanthine, 1 µM dexamethasone, and 10 μg/ml insulin; Sigma) and cultured for 2 days. The differentiation medium was then replaced with maintenance medium (DMEM containing 10% FBS and 10 μg/ml insulin). Cell morphology was observed under a microscope. Twelve cell culture samples were obtained from each pig. The samples were collected from each pig at 0, 2, 4, and 8 days in triplicate and stored in liquid nitrogen until RNA isolation.

Oil Red O staining was performed to identify lipid droplets. The cell culture plates were gently washed with phosphate-buffered saline (PBS) three times and fixed in 10% formaldehyde for 15 min. The cells were then washed with PBS three times and stained with Oil Red O for 20 min. Finally, the cells were washed three times with PBS and photographed using an inverted microscope (Leica, Wetzlar, Germany). An equal volume of 100% isopropanol solution was added to the wells of each culture plate, and the absorbance at 500 nm was measured after thoroughly homogenizing the samples. Each experiment was repeated three times.

The total RNA was purified from the 12 samples using TRIzol reagent (Invitrogen, Carlsbad, CA, United States) according to the manufacturer’s instructions. RNA quantity and quality were assessed using the Agilent 2100 Bioanalyzer (Agilent, Santa Clara, CA, United States). All RNA samples with RNA integrity numbers >8.0 and an absorbance 260:280 ratio >1.9 were selected for library construction and deep sequencing.

From each of the 12 samples, 10 μg of RNA was used for RNA sequencing (RNA-seq) library preparation using the TruSeq® Stranded Total RNA Sample Preparation Kit (Illumina, CA, United States) according to the manufacturer’s instructions. The ligation products were size-selected by agarose gel electrophoresis and amplified by PCR. After purification and enrichment, 12 libraries were sequenced using the Illumina HiSeq 2500 system by Gene Denovo Biotechnology Co. (Guangzhou, Guangdong, China).

Raw reads were cleaned by removing adapter and primer sequences, reads containing more than 10% of unknown nucleotides (N), and more than 50% of low-quality (Q-value ≤ 20) bases. The clean reads were mapped to the pig reference genome (Scrofa 11.1, ftp://ftp.ensembl.org/pub/release-68/fasta/sus_scrofa/cdna) (Langmead et al., 2009) using Bowtie 2 (Langmead and Salzberg, 2012) and TopHat2 (version 2.0.3.12) (Trapnell et al., 2012) with default parameters. HT-seq (version 0.6.1) was used to calculate the read number mapped to each gene in each sample (Glaser et al., 2010).

Differential gene expression, based on the normalized read count of expressed genes, was analyzed in three time-point contrasts (0 vs. 2 days, 0 vs. 4 days, and 0 vs. 8 days) using the edgeR package (Robinson et al., 2009) in R (version 4.0.2). Genes with a false discovery rate (FDR) ≤0.05 and absolute log₂ fold change (|log₂FC|) >1 were identified as differentially expressed genes (DEGs) between two groups (Zhao et al., 2019).

As mentioned previously, fat deposition is affected by preadipocyte differentiation. Thus, to identify candidate genes influencing preadipocyte differentiation, we applied WGCNA (Langfelder and Horvath, 2008) to fragments per kilobase of exon model per million mapped fragments (FPKM) values obtained from the RNA-seq of the 12 samples to identify genes associated with the degree of preadipocyte differentiation. Based on the expressed coding genes with FPKM >0, hierarchical clustering was performed. Outliers above a height threshold of 2,000 were filtered out. The remaining samples were used to establish an unsigned co-expression network.

We then used a one-step method to construct the network and determine the gene module. According to Zhang and Horvath (2005), gene co-expression networks should have scale-free characteristics and follow a power-law distribution. Hence, a weighted adjacency matrix was created, which can be defined as follows: A i j   =   | c o r   ×   ( x i , x j ) | β , (1) where xi and xj are the ith and jth genes, respectively.

Adjacent to the adjacent network is the combination of the soft-thresholding power parameter β, which is required to improve the co-expression similarity for computing the adjacency. To keep the network consistent with scale-free topology, the pickSoftThreshold() function in R was used to analyze its topology and select an appropriate soft-thresholding power value (β) to build it, and the mean connectivity of all genes in the module was evaluated. A power of 18 was selected as the soft threshold to ensure a scale-free topology network (R ² = 0.85).

Subsequently, the topology overlap matrix (TOM) (Zhang and Horvath, 2005), which measures the connectivity of a pair of genes, was calculated from the adjacency matrix. Once the power value was determined, the TOM and dissTOM = 1−TOM were obtained. Based on hierarchical average linkage clustering using a dynamic tree-cutting algorithm, genes with similar expression patterns were clustered into the same modules and assigned a color (Zhang and Horvath, 2005).

After construction, the module eigengene (ME) was calculated using the first principal component of the expression profiles, which represented the weighted average expression profile. To identify biologically significant modules and select potentially critical modules for downstream analysis, the WGCNA approach was used to define the module–trait relationships (MTRs) (Kogelman et al., 2014) and gene significance (GS) of each module (Liu et al., 2016). We considered the time point of cell differentiation as a dichotomous variable (Langfelder and Horvath, 2008). A Pearson’s correlation test was run between the ME and time point of cell differentiation and between the expression profiles and time point of cell differentiation to estimate the MTRs and GS, respectively. The module significance (MS) was the mean GS value of the module genes. According to the selection criteria for critical modules reported in a previous study, modules with MTR > 0.30 and MS > 0.25 were considered candidate modules for functional enrichment analysis (Howard et al., 2017). Highly connected genes, also known as hub genes, may play important roles in a module (Wilhelm et al., 2004). Hub genes are relatively conserved at the core of the gene co-expression network, which can act as a genetic buffer to reduce the effect of other gene mutations (Kadioglu et al., 2021). Here, genes with GS > 0.3, module membership (MM) > 0.85, and intramodular connectivity > 5 were considered hub genes. Cytoscape (Shannon et al., 2003) was used to map the gene–gene interaction network for visualizing gene relationships.

Gene Ontology (GO; http://www.geneontology.org/) is widely used in the field of bioinformatics to classify genes into terms from three different biological categories: cellular components, molecular functions, and biological processes. The default values were adapted for the parameters of the phenotype-related module, and three GO term enrichment analyses were performed on the genes in the module. Adjusted p-values <0.05 were considered significant, and the 10 most prominent entries for each analysis were retained.

Kyoto Encyclopedia of Genes and Genomes (KEGG, http://www.genome.jp/kegg/) is a database for systematic analysis of gene function and genomic information that facilitates the study of genes and gene expression as a part of an entire network. “ClusterProfiler” (Yu et al., 2012) and “ggplot2” packages were used to analyze and visualize such genetic information, respectively. The R package BioMart (http://www.biomarbiomart.org/) (Haider et al., 2009) was used to annotate genes in the module with Sscrofa11.1 as a reference genome. We selected a subset of modules based on their functional annotation and selected genes related to fat development. Based on the above information, the candidate genes affecting fat deposition were determined in this experiment.

In our previous study (Zhao et al., 2019), ANGPTL4 was identified as a DEG in all time-point contrast groups and was related to the PPAR signaling pathway and cell differentiation process. Similarly, in the present study, ANGPTL4 was identified in a key module; therefore, we hypothesized that ANGPTL4 has some regulatory effects on preadipocyte differentiation. According to the sequence of pig ANGPTL4 (ID: 397,628) in GenBank, three pairs of siRNAs targeting and corresponding negative controls were designed and synthesized by GenePharma (Suzhou, Jiangsu, China) (Table 1). The siRNAs were centrifuged at 10,000 rpm for 2 min and dissolved in 125 μl of Nuclease-Free Tubes, Tips, and Buffers (DEPC) water (Gibco, NY, United States) to a final concentration of 20 μmol/L.

Once the cell confluence reached 70%–80%, Lipofectamine 2000 (Invitrogen) was used for transfection. According to the manufacturer’s instructions, the medium was replaced with 2.5 ml of DMEM/F12 medium without serum and antibody, 30 min before transfection. Thereafter, 5 µl of Lipofectamine 2000 was diluted with 250 µl of serum and double-antibody-free DMEM/F12 medium per well, mixed gently, and incubated at room temperature (21°C–25°C) for 5 min; this was mixture I. Next, 10 µl of siRNA was collected from each well, diluted with 250 µl of serum-free and double-resistant DMEM/F12 medium, mixed gently, and incubated at RT for 5 min; this was mixture II. Mixtures I and II were gently mixed and allowed to stand for 20 min at RT. Thereafter, mixtures I and II were added to each well to a total volume of 500 µl and mixed thoroughly. After incubation at 37°C for 6 h, the medium was replaced with a complete medium. These steps were followed to induce cells and were repeated three times for each experiment.

Based on our previous study (Zhao et al., 2019), ANGPTL4 has the same expression pattern as Acetyl CoA acyltransferase 2 (ACAA2), aldehyde dehydrogenase two family member (ALDH2), and solute carrier family 27 member 1 (SLC27A1). Furthermore, lipoprotein lipase (LPL), stearoyl-CoA desaturase (SCD), and fatty acid synthase (FASN) are markers of fat deposition. Thus, we observed whether the expression of those genes changed before and after the interference of ANGPTL4.

The total RNA was extracted using a TRIzol reagent (Invitrogen) and reverse-transcribed into cDNA according to the manufacturer’s instructions. Reverse-transcription quantitative polymerase chain reaction (RT-qPCR) was performed using the Light Cycler 480 Real-Time PCR system (Roche, Hercules, CA, United States). The primers used for RT-qPCR detection of the selected genes are listed in Supplementary Table S1. All RT-qPCR experiments were performed using three biological replicates with three technical replicates for each sample. The 2^−ΔΔCt method was used to measure the change in mRNA abundance, and glyceraldehyde 3-phosphate dehydrogenase (GAPDH) was used as the internal control.

All qRT-PCR results are presented as mean ± standard deviation (SD). Statistical analyses were conducted using SPSS software (version 20.0; IBM Corp., Armonk, NY, United States); results with p < 0.01 were considered extremely significant, and those with p < 0.05 were considered significant.

Compared to cell shapes observed during the initial phase (day 0), preadipocytes gradually changed from fibrous to spherical on day 2. Subsequently, lipid droplets became visible on day 4, gradually increasing in number until day 8 (Figure 1). These results indicated that the preadipocytes differentiated successfully and could be used for further analysis.

Genes with a sum of FPKM values >0 (n = 12,816) were selected for a co-expression network analysis (Supplementary Table S2). The WGCNA package in R was used to construct co-expression modules. No outlier samples were found in the hierarchical clustering of samples using the flashClust package in R (Supplementary Figure S1).

According to the standard of a scale-free network, we selected an appropriate weighted parameter of the adjacency function, namely, the soft threshold, which was 18 in this study (Supplementary Figure S2). We then calculated the correlation and adjacency matrices and combined them into the topology matrix. We finally identified 27 gene modules based on genetic similarity after merging modules with dissimilarities less than 0.25 and a minimum size of 30 (Figure 2; Supplementary Table S3).

The Pearson correlation coefficient between the eigengenes of modules and the corresponding variables represents the linear correlation between the module and phenotypic information (Figure 3). We found that the blue and brown modules were significantly correlated with day 8 (r = 0.59, p < 0.05; Figure 3), which was the time point with the most preadipocytes, suggesting that genes in these two modules promote preadipocyte differentiation. In addition, the magenta module was significantly correlated with day 0 when there was little differentiation (r = 0.78, p < 0.01), suggesting that these genes inhibit preadipocyte differentiation.

As a final module assessment, scatter plots of GS for preadipocyte differentiation versus MM for the blue, brown, and magenta modules are presented in Figure 4. We found a significant correlation between GS and MM, suggesting that genes in the preadipocyte differentiation-related modules tended to highly correlate with fat deposition.

For the genes in the blue, brown, and magenta modules, GO term and KEGG pathway enrichment analyses were performed (Table 2, Table 3). We found that genes in the blue module play a role in fatty acid degradation, those in the brown module play a role in mitogen-activated protein kinase (MAPK) signaling pathways, and genes in both modules participate in fatty acid beta-oxidation. Genes in the magenta module are related to fatty acid beta-oxidation using acyl-CoA dehydrogenase, the Wnt/PI3K-Akt/TGF-beta signaling pathway, and the regulation of stem cell pluripotency. These results suggested that genes in the blue, brown, and magenta modules are related to lipid metabolism.

Genes with the highest degree of connectivity were selected from the three candidate modules, and Cystoscape software was used to draw the gene interaction network diagrams (Figure 5). We found ANGPTL4 in the blue module (Supplementary Table S4). Thus, we knocked down ANGPTL4 to observe its effect on preadipocyte differentiation.

Three pairs of siRNAs were synthesized. The expression level of ANGPTL4 was determined 0, 2, 4, and 8 days after differentiation following transfection (Figure 6). The transfection efficiency of siRNA-752 was the highest (>90%); thus, it was used in the subsequent experiments.

As shown in Figure 7A, the mRNA levels of ANGPTL4, ACAA2, SLC27A1, ALDH2, PPAR, SCD, FASN, and LPL increased with preadipocyte differentiation. After transfection, we observed changes in the expression levels of these genes. The expression level of ANGPTL4 significantly decreased, confirming successful transfection. Except for LPL, whose expression level was opposite to that of ANGPTL4, the expression of the other genes initially increased and then decreased (Figure 7B).

Two days after transfecting preadipocytes with siRNA-752, the inducer was added to promote differentiation. The cells were then stained with Oil Red O and observed under a microscope. Compared with the negative control group, the siRNA-752 transfection group showed a significant increase in lipid droplet production (Figures 8A,B) and OD_500nm values (p < 0.05) (Figure 8C), indicating that the knockdown of ANGPTL4 promoted the differentiation of porcine preadipocytes.