All animal experimental procedures were inspected and reviewed by the Animal Welfare and Ethic Committee of the Institute of Animal Science, Chinese Academy of Agriculture Sciences (approval number: IAS 2021-73) and conduced according to the guidelines established by the Institutional Animal Care and Use Committee of the Chinese Academy of Agriculture Sciences. Animal and sample collection work was performed at the National Tan Sheep Conservation Farm. A total of 15 female Tan sheep (fat-tail sheep) and five female Suffolk sheep (thin-tail sheep) aged from three to 4 years old were maintained under the same feeding regimen, management and environment.

All sheep were treated with estrus synchronization and artificial insemination. In brief, this involved the following: 1) a controlled internal drug releasing (CIDR) progesterone implant for 11 days; 2) removal of the progesterone implant and injection of 1 ml cloprostenol (0.1 mg/ml); 3) estrus identification after 48 h; 4) only sheep that were in estrus at the same time were used for artificial insemination; and 5) A B-ultrasonic machine (Gandaofo, Henan, China) was used to evaluate pregnancy status after approximately 30 days. When the Tan sheep embryos developed to the 60-day (E60), 70-day (E70) and 80-day (E80) stages individually, and the Suffolk sheep embryo reached the 70-day (E70_SFK) stage, the pregnant ewes were transferred to the slaughterhouse. Samples of embryonic tail, skin, heart, muscle, liver, and kidney were collected from at least three individual animals, and fat tissues from the tail of every adult sheep were also harvested.

Samples were immersed in 4% paraformaldehyde solution for fixation, dehydrated in a graded ethanol series, cleaned with xylene and embedded in paraffin wax. Using a RM2255 Automated Rotary Microtome (Leica, Wetzlar, Germany), the paraffin-embedded tissue blocks were sectioned into 5 μm-thick slices, which were then stained using hematoxylin and eosin (HE staining). For oil red staining, tissue was frozen with embedding medium swiftly after fixation and cut into sections as described above; it was then stained with oil red O solution. A BX51 microscope (Olympus, Shinjuku, Japan) and a DP72 digital imaging system (Olympus, Shinjuku, Japan) were used to visualize stained sections.

The tail tissues of Tan sheep at four time points (E60, n = 4; E70, n = 3; E80, n = 4; Adult, n = 3) and the tail tissues of Suffolk sheep at E70 (E70_SFK, n = 4) were used for RNA extraction with the RNeasy lipid tissue Kit (Qiagen, Hilden, Germany). The concentration and quality of RNA samples were evaluated using an Agilent 2,100 Bioanalyzer (Agilent, California, United States). Samples with RIN values higher than 7.0 and concentrations greater than 40 ng/ul were used for RNA-seq. The mRNA selection, library preparation and sequencing were performed by the BerryGenomics Company (Beijing, China) and sequencing was performed on the HiSeq X Ten high-throughput sequencing (HTS) platform (Illumina, California, United States).

The NGS QC Toolkit (v2.3.3) (Patel and Jain, 2012) was used to perform quality control of the raw data, including removal of reads with overall quality scores less than 20 according to the Phred+33 scale for at least 70% of the bases and adapter sequences. The remaining clean reads were mapped to the sheep reference genome (Oar v4.0) (Jiang et al., 2014) using the TopHat2 software tool (v2.1.1) (Kim et al., 2013) with default parameters.

Assembly of the mapped reads and transcriptome quantification (in Fragments Per kilobase of transcript per Million mapped reads; FPKM) were performed using the Cufflinks package (v2.2.1) (Trapnell et al., 2012) with the sheep reference genome and annotation file. The Cuffmerge program was then used to generate a new genome annotation file with unified annotation information. Following this, the Cuffdiff program was used to identify DEGs between experimental groups. We removed genes that, at least in one group, exhibited mean FPKM value less than 0.5. Genes that exhibited absolute log₂ fold change values ≥ 1 and FDR-adjusted q values <0.05 were considered to be differentially expressed.

The gmodels package was used to perform principal component analysis (PCA) in R (version 3.6.1) with the FPKM values for all annotated transcripts from the eighteen transcriptomes. The heatmap and Venn diagrams were generated using the FPKM values and gene symbols of corresponding DEGs with the Pheatmap and VennDiagram packages in R (version 3.6.1).

The standardized and centralized FPKM values were used for hierarchical analysis, with cluster numbers determined by the Calinsky criterion. Using the “timeclust” parameter of the TCseq package (Gu, 2021) to perform the hierarchical analysis with corrected FPKM. The cluster dendrogram was divided using the “complete” parameter to classify the genes based on expression trend.

Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses was implemented using the clusterProfiler package (Wu et al., 2021). The GSEA software (version 4.0.2) (Subramanian et al., 2005), which uses predefined gene sets from the Molecular Signatures Database, was applied to perform gene set enrichment analysis based on FPKM values for the E70 and E70_SFK experimental groups. In this study, we used H: hallmark gene sets, CP:KEGG: KEGG gene sets and C5: GO gene sets as molecular signatures databases and gene sets with FDR q values <0.05 were considered to represent significant pathways. The minimum and maximum criteria for selection of gene sets from the collection were 10 and 500 genes, respectively.

Using previously described methods for adipose tissue digestion in mouse (Sun et al., 2020), we isolated the primary adipose derived stem cells (ADSCs) from embryonic fat tail tissues at E80 by collagenase digestion. Following this, the ADSCs were cultured in DMEM/F12 medium containing 10% FBS and 1% penicillin/streptomycin (growth medium). The differentiation process to generate adipocyte from ADSCs required 6 days in total. Isolated primary cells were cultured in induction medium (growth medium supplemented with 10 ng/ml insulin, 1 μM dexamethasone, 1 μM rosiglitazone and 0.5 mM IBMX) for 2 days, and following this, the induction medium was replaced with keep medium (growth medium supplemented with 10 ng/ml insulin and 1 μM rosiglitazone) for 4 days and the keep medium was refreshed every 2 days. At d6, a large number of droplets were observed by oil red staining. The differentiating cells at d0, d2, d4 and d6 were then used to extract RNA with a standard Trizol method (Invitrogen, Waltham, United States) and examined for expression of MTFP1 using RT-qPCR.

Small interfering RNA (siRNA) transfections were conducted in isolated primary ADSCs. To minimize potential off-targets effect of siRNA, three different siRNAs (siMTFP1-1, siMTFP1-2 and siMTFP1-3) targeting the MTFP1 gene and a negative control siRNA (siNC) were designed (Supplementary Table S9). Lipofectamine 3,000 (Invitrogen, California, United States) was used with these siRNAs for downregulation of MTFP1 gene expression according to the manufacturer’s protocol. First, we measured the knockdown efficiency for the three siRNAs and the most significant of these (siMTFP1-1) was used to perform the downstream experiments. The cells culture in growth medium were treated with siMTFP1 to examine its effects on the expression of mitochondrion- and adipogenesis-related genes, and proliferation and migration ability. The RT-qPCR experiment was carried out as previously described (Dong et al., 2020). The sequence of the siRNAs and the primer sequences are provided in Supplementary Table S9. Also, the cells cultured in induction medium and keep medium were treated with siMTFP1 for 6 days to examine its effect on cellular adipogenesis at d6.

Three methods were used to test cell proliferation ability depending on whether MTFP1 is expressed or knocked down. The first method used a Cell Counting Kit-8 (CCK-8, Sigma-Aldrich, Germany) (Yang et al., 2021). This method involved addition of 10% CCK-8 solution to the cell culture medium, followed by incubation for 1 h at 37°C and colorimetric changes was measured using optical density at 450 nm with a microplate reader (Tecan Infinite 200 Pro, Männedorf, Switzerland). The second method used a different colorimetric assay based on 3-(4,5-dimethylthiazol-2-yl)-2,5-diphenyltetrazolium bromide (MTT, Sigma-Aldrich) (Prabst et al., 2017). MTT solution (10%) was added to the cell culture medium, which was then incubated for 3 h at 37°C. The mixture was discarded, and the same volume of Formazan solution was added for 10 min and the colorimetric change was measured as optical density at 490 nm. The third method involved the use of alamarBlue (Sigma-Aldrich) (Prabst et al., 2017) as a 10% solution added to the cell culture medium, followed by incubation for 2–6 h at 37°C and detection of relative fluorescence units according to the manufacturer’s instructions.

The wound healing assay was performed as follows: the cells were allowed to grow to 100% coverage and a scratch was introduced into the middle of plate; the growth medium was replaced with Opti-MEM (Gibco, Waltham, United States); and photographic images of the scratch at 0 and 12 h were analysed using the ImageJ software (Schneider et al., 2012) for area quantification.

Oxygen consumption rate (OCR) assay was carried out using a BBoxiProbe™ R01 kit (Bestbio, China) to evaluate cell metabolism and mitochondrial state. In brief, this involved the following: 1) cells were cultured in a 96-well black plate with transparent bottom and treated with siMTFP1 for 48 h; 2) refreshing the growth medium and adding 5 μl oxygen fluorescent probe; 3) oxygen blocking buffer was added to prevent external oxygen; 4) relative fluorescence unit (RFU) was detected using a microplate reader every 5 min for 1.5 h.

To examine the morphogenesis of the fat-tail trait by external observation, HE and oil red staining were performed for embryonic tail tissue at three time points (E60, E70 and E80). The external morphology of tail tissue at E60 and E70 did not show any change except that the thickness increased gradually with development (Figure 1A). In the sections collected from the E60 and E70 stages, circular preadipocytes accumulated to form regular leaf-like cell masses separated by fibers without fat deposition (Figure 1B and Supplementary Figure S1A). Obvious morphogenesis occurred at E80, including increased thickness, and altered shape characteristics. The sections obtained from E80 showed visible lipids, highly dense adipocytes, and an irregular cell mass shape (Figure 1B and Supplementary Figure S1A). These results demonstrate that morphological changes were caused by fat deposition.

The higher proportion of unilocular adipocytes are mainly distributed in the central area of the cell masses and multilocular adipocytes at the edge (Supplementary Figure S1A), which indicate that the lipid droplets gradually fuse to form a single larger droplet during the morphogenesis of fat tail. In addition, the edge areas of cell masses at E80 were surrounded by fibers developed from veined fiber networks at E60 and E70 (Supplementary Figure S1A). To determine whether specific fat deposition occurred in the tail tissue at E80, we also investigated fat deposition in five additional tissues (muscle, heart, skin, liver and kidney); these results showed that fat droplets stained red were only observed in tail tissue (Supplementary Figure S1B).

To investigate the genomic regulatory network that promote tail development and to obtain a comprehensive view of the transcriptional changes involved in fat-tail morphogenesis, we performed RNA-seq for tail samples from Tan sheep (fat-tail) at four time points, including E60 (n = 4), E70 (n = 3), E80 (n = 4) and adult (Fat) sheep (n = 3), as well as Suffolk sheep (thin-tail) at E70 (E70_SFK) (n = 4) (Supplementary Figure S2A). Between 32 and 63 million 150-bp paired-end reads were generated for each sample with a mean mapping rate of 78.8% (Supplementary Table S1). After filtering for lowly expressed genes, there were 15,427 genes remaining for subsequent analyses that were expressed in all samples. The expression values of these genes were transformed into FPKM values. The PCA of the expression data partitioned the samples into five broad clusters with the E70_SFK samples dispersed between the E60 and E70 groups (Supplementary Figure S2A). The heatmap generated using all expressed genes showed that E60 and E70_SFK samples clustered on the same branch, with E70 and E80 also sharing a similar transcriptome profile, and with the fat-tail tissue of adult sheep emerging in a single branch (Supplementary Figure S2B).

Subsequently, differentially expressed gene (DEG) analysis among different groups was performed to identify the key regulators involved in tail fat deposition. The comparison between the E60 and E70 stages was termed transition 1 (Trans1), that between E70 and E80 as transition 2 (Trans2), and that between E70 and E70_SFK as transition 3 (Trans3) (Figure 1C). A total of 2,929, 1702 and 1909 DEGs were identified across Trans1, Trans2 and Trans3, respectively, and there were 4,328 unique genes in the union set of DEGs (Supplementary Table S2). There were 1,293 DEGs shared between Trans1 and Trans3 with a totally different expression trend (Figure 1C), including 1,285 DEGs with the opposite direction of expression and eight DEGs with the same direction of expression, an observation consistent with the fact that Trans1 (from E60 to E70) was a process of positive fat tail morphogenesis, while Trans3 (from E70 to E70_SFK) was a comparison between fat and thin-tailed sheep, which represents a negative process of fat tail morphogenesis. There is a significant time gap between the E80 and adult stages; therefore, we examined three different contrasts among the three embryonic time points and the adult time point (Comp1: E60-Fat; Comp2: E70-Fat; and Comp3: E80-Fat). These analyses detected 5,697, 4,474 and 3,858 DEGs for Comp1, Comp2 and Comp3, respectively. The total number of DEGs decreased across the Comp1, Comp2, and Comp3 comparisons, which was reflected in the numbers of DEGs exhibiting increased or decreased expression (Table 1).

Since Trans3 was the comparison between fat-tail and thin-tail sheep at the same time point, we targeted the Trans3 DEGs that overlapped with the other two transition groups (Trans1 and Trans2) as well as the three comparison groups (Comp1, Comp2, and Comp3) (Figure 1C). A total of 1,058 DEGs were detected (Supplementary Table S2) with different expression patterns across four stages of fat tail development (E60, E70, E80, and Fat), which can be divided into six clusters (K1–K6) by the Calinski-Harabasz Criterion (Calinski and Harabasz, 1974) (Figure 1D, Supplementary Figures S2C, S3). The top two clusters, K2 and K4, containing 271 and 386 genes and accounted for 25.6 and 36.5% of the DEGs, respectively. Cluster K2 represents a cluster of periodic genes that were increased in expression during E60 to E70 and decreased in expression during E70 to E80. Cluster K4, on the other hand, consisted of the periodic genes that exhibited a pattern of expression opposite to cluster K2 (Figure 2).

Through KEGG pathway enrichment analysis for each cluster (Supplementary Figure S3 and Supplementary Table S3), five pathways were specifically enriched in K2, including Oxidative phosphorylation and Thermogenesis (Figure 2), which are closely related to energy metabolism. The top 10 most significant pathways in K2, in addition to the Proteasome and Ribosome, shared a subset of DEGs (Supplementary Figure S2D). The pathways overrepresented in K4 are highly associated with cell activity, including Focal adhesion, ECM-receptor interaction, and PI3K-Akt signaling pathway (Figure 2). Furthermore, gene ontology (GO) enrichment analysis was conducted for these two clusters and the top ten most significant terms in each GO category were determined (Supplementary Figure S4 and Supplementary Table S4). Almost all the significant enriched GO terms in K2 were associated with the mitochondrion, such as Respiratory chain, Mitochondrial part, which was the most significant term, and Electron transfer activity (Supplementary Figure S4).

In order to further understand the function of these genes, we performed gene set enrichment analysis (GSEA) (Subramanian et al., 2005) using the gene expression data between E70 and E70_SFK. GSEA is an enrichment method that evaluates transcriptome profile data at the level of gene sets. As a result, 37 gene sets were specifically enriched at E70 (q < 0.05), which can be mainly divided into two groups according to their biological function (Supplementary Table S5). One group related to different components of the mitochondrion as well as energy metabolism, and another group is involved in ribosome function and translation. The top 10 most significant gene sets in each part are displayed in Figure 3A. However, the 52 gene sets enriched at E70_SFK were difficult to assign to a primary functional category, and the GO term DNA binding transcription factor activity exhibited the highest NES score (Supplementary Figure S5B).

Furthermore, 183 genes have been identified in the 37 gene sets that were enriched at E70 by leading edge gene analysis compared with E70_SFK (Subramanian et al., 2005) (Figure 3B and Supplementary Figure S5C). Sixty-nine DEGs were associated with various components of the mitochondrion (Supplementary Table S7), including TIMM50, MPC1, MPC2, and TOMM5 genes; other genes were involved in membrane transport in mitochondrion and subunits of NADH dehydrogenase (NDUFA, NDUFB, NDUFC, and NDUFS). There are also some genes encoding proteins with a functional role in regulating gene expression and cellular activities. The GNG11 and GNG5 genes encode proteins belonging to the G protein gamma family that are involved in various transmembrane signaling systems (Gautam et al., 1989). The POLR2G, POLR2I, and POLR3K genes encode subunits of the DNA-dependent RNA polymerase that synthesizes mRNA precursors and many functional non-coding RNAs (Ujvari and Luse, 2006). The NME2 and NME4 genes (Hu et al., 2019) encode proteins within the NDK (nucleoside diphosphate kinase) family and have multiple functions in cellular energetics, signaling, proliferation, differentiation, and tumor invasion (Boissan et al., 2009). In addition, there were four genes that were most overexpressed in E70 compared with E70_SFK, including GMPR, GSTP1, SLC25A4, EIF4EBP genes, which may encode important regulators of nucleotide metabolism (Figure 3B).

Among these 183 enriched genes, 17 genes have been reported in previous tail phenotype related studies in sheep (Wang et al., 2014; Li et al., 2018; Ma et al., 2018; Bakhtiarizadeh et al., 2019; Kalds et al., 2021; Wang et al., 2021; Zhang et al., 2021). Expression of these genes fluctuated from the E60 to Fat stages and was significantly higher at E70 than E70_SFK (Supplementary Figures S6A,S6B). The reported genes can be divided into two groups according to the pathways they were assigned to (Supplementary Table S8): energy metabolism related genes (MTFP1, ACAA1, GPX4, MRPL28, MSRB1, NDUFA6, PPP1CA, SLC25A39, and TOMM5) and regulation of transcription related genes (AP2S1, DBI, EEF1D, MVD, NUDT16L1, PDGFD, PSMB5 and RPLP2). The expression of MTFP1 and EEF1D genes was highly positively correlated to other previously reported genes (Figure 3C). It has been reported that mitochondrial fission process protein 1 (encoded by MTFP1) is one of the key regulators of mitochondria fission, and involved in cell apoptosis, carcinogenesis, and tumor progression. In particular, the expression of MTFP1 was most significantly different between E70 and E70_SFK with the highest enrichment score (Figures 3B,C); this warrants further research for its functional role in adipogenesis in vitro.

Firstly, we found that the expression of MTFP1 fluctuated during the process of ADSCs differentiation (d0, d2, d4 and d6) and the expression at the late stage is higher than that at the early stages, which is consistent with the expression trend in fat-tail tissues from the E60 to Fat stages (Figure 4A). Small interfering RNA (siRNA) was used to downregulate the expression MTFP1 and to investigate its role in proliferation, migration and adipogenesis. Three siRNAs were designed to specifically downregulate expression of MTFP1 and the siRNA with the most effective knockdown ability (siMTFP1-1) was selected to perform downstream assays and minimize potential off-target effects (Supplementary Figure S6C). After treatment with siMTFP1, the gene expression decreased to 20% of normal expression (Figure 4B). Three methods, including CCK-8, alamarBlue and MTT, jointly showed that knockdown of MTFP1 inhibited cell proliferation (Figure 4C), which was in concordance with the results of the wound healing experiment showing that downregulation of this gene significantly reduced cell migration ability (Figure 4D). Furthermore, we examined the expression of other mitochondrion-related genes that also has been identified in this study (Supplementary Table S7) and adipogenesis-related genes.

Knockdown of MTFP1 resulted in decreased expression of DRP1, NME2, MRPL18, TIMM50, and TOMM5, which are genes involved in the mitochondrial membrane, ribosome and membrane protein (Smirnov et al., 2011; Morita et al., 2017; Chaudhuri et al., 2021). The genes (NDUFA6, NDUFB8, and NDUFS8) that encode subunits of NADH dehydrogenase (Stroud et al., 2016) retained their original expression levels. Besides, expression of ARL2 and NME4 was upregulated (Figure 5A). Although it has been reported that expression of NME4 can stimulate respiratory ATP regeneration (Tokarska-Schlattner et al., 2008), the oxygen consumption ability within these cells was still inhibited after silencing MTFP1, which may caused by decreased cell number (Figure 5B). For adipogenesis-related genes, decreased expression of MTFP1 induced upregulation of adipocyte-related genes (PPARG and LPL) and downregulation of progenitor-associated genes (TOP2A and BIRC5), reflecting a positive effect on adipogenesis after silencing MTFP1 (Figure 5C) We also treated the ADSCs with siMTFP1 during cell differentiation, which significantly promoted fat deposition (Figure 5D).