Multi-omics integrated analysis reveals the molecular mechanism of tail fat deposition differences in sheep with different tail types

Abstract

Background

The accumulation of tail fat in sheep is a manifestation of adaptive evolution to the environment. Sheep with different tail types show significant differences in physiological functions and tail fat deposition. Although these differences reflect the developmental mechanism of tail fat under different gene regulation, the situation of sheep tail fat tissue at the single cell level has not been explored, and its molecular mechanism still needs to be further elucidated.

Results

Here, we characterized the genomic features of sheep with different tail types, detected the transcriptomic differences in tail adipose tissue between fat-tailed and thin-tailed sheep, established a single-cell atlas of sheep tail adipose tissue, and screened potential molecular markers (SESN1, RPRD1A and RASGEF1B) that regulate differences in sheep tail fat deposition through multi-omics integrated analysis. We found that the differential mechanism of sheep tail fat deposition not only involves adipocyte differentiation and proliferation, but is also closely related to cell-specific communication networks (When adipocytes act as signal outputters, LAMININ and other signal pathways are strongly expressed in guangling large tailed sheep and hu sheep), including interactions with immune cells and tissue remodeling to drive the typing of tail fat. In addition, we revealed the differentiation trajectory of sheep tail adipocytes through pseudo-time analysis and constructed the cell communication network of sheep tail adipose tissue.

Conclusions

Our results provide insights into the molecular mechanisms of tail fat deposition in sheep with different tail types, and provide a deeper explanation for the development and functional regulation of adipocytes.

Background

The accumulation of sheep tail fat in specific breeds is a manifestation of evolutionary adaptability to the environment and has important biological and economic value [1].According to the amount and distribution of fat deposition, sheep are divided into thin tail type, short fat tail type and long fat tail type [2]. Sheep with different tail types show significant differences in physiological functions and gene expressions. These differences not only reflect the developmental mechanism of tail fat under the regulation of different genes, but also reveal its unique role in cell signaling and metabolic pathways [3, 4]. Sheep with long fat tails generally show higher fat storage capacity than sheep with short fat tails, while sheep with thin tails have less fat deposits in their tails (Fat deposition capacity: thin-tailed sheep < short fat-tailed sheep < long fat-tailed sheep). This phenotypic difference is mainly regulated by hub genes (such as PPARγ and FAS) in the process of fat synthesis and decomposition [5, 6]. The differentiation of adipocytes and their functional changes are closely related to a variety of transcription factors, hormone signals and molecular pathways related to the extracellular matrix [7]. Therefore, further analysis of the molecular mechanism of tail fat deposition in sheep with different tail types will help understand the biological laws of lipid metabolism and provide theoretical support for animal husbandry production.

As a cutting-edge method in modern biological research, multi-omics analysis is promoting a comprehensive and systematic understanding of complex biological processes [8]. Single-cell transcriptomics is a high-resolution technology that can capture cell-specific gene expression patterns compared to traditional population-level analysis, which helps to study different cell types, cell states, and functional characteristics. By analyzing the transcriptional characteristics of different cell populations, single-cell transcriptomics can reconstruct the cellular composition of tissues and reveal cell-to-cell interactions and key regulatory networks [9]. Genomics and transcriptomics play complementary roles in multi-level biological research. Genomics mainly analyzes the DNA sequence information of individuals, including the arrangement and mutation of genes, to explore the impact of genetic variation on biological traits [10]. Transcriptomics focuses on the dynamic changes of gene expression, and studies the activity of genes and their regulation by detecting the expression profile at the RNA level [11]. Combined with single-cell technology for integrated analysis, it has greatly promoted the development of precision medicine and complex trait research.

Adipose tissue, once considered morphologically and functionally bland, is now recognized to be dynamic, plastic, and heterogeneous and to participate in a wide range of biological processes, including energy homeostasis [12], glucose and lipid handling [13], blood pressure control [14] and host defense [15]. The single-cell atlas of human and mouse subcutaneous fat includes cell types such as adipocytes, adipose stem cells, progenitor cells and immune cells. Through further analysis, a new adipocyte subpopulation CYP2E1 + ALDH1A1 + was identified. This subpopulation can control the thermogenic function of other adipocytes in a paracrine manner by regulating the level of short-chain fatty acid acetate [16]. Comparing the immune cells in visceral adipose tissue of obese and healthy mice, it was found that obesity can change the composition of immune cells in adipose tissue, mainly including an increase in the proportion of macrophages, and an increase in the proportion of regulatory T cells and type 2 innate lymphoid cells ratio decreases [17]. The identification of primary adipocytes in the brown adipose tissue of mice revealed that in addition to the classic high-thermogenic brown adipocytes, there is also a type of low-thermogenic brown adipocytes; compared with the former, the latter have low expression of thermogenic and lipolysis genes and high expression of fatty acid uptake genes. When the ambient temperature changes, the two types of brown adipocytes can transform into each other, thereby changing the thermogenic capacity of the entire brown adipose tissue [18]. In livestock animals, subcutaneous and visceral adipose tissue of dairy cows is mainly composed of mature adipocytes, three stem cell subtypes, and a variety of immune cells. Among them, the ASPC subtype is defined as adipogenic (PPARG+), while the other two have fibroadipogenic characteristics (PDGFRA+) [19]. However, the expression of sheep tail adipose tissue at the single-cell level has not been explored, and the genetic and molecular mechanisms underlying the differences in tail fat deposition in sheep with different tail types remain to be elucidated.

In this study, we characterized the genomic characteristics of sheep with different tail types, explored the transcriptomic differences in tail adipose tissue between fat-tailed and thin-tailed sheep, mapped the single-cell atlas of sheep tail adipose tissue and identified its cell types. We compared the differences between Guangling Large-tailed Sheep (GLT; long fat-tailed) and Hu sheep (Hu; short fat-tailed) tail adipocytes at the cellular level to further explore potential molecular markers that regulate fat deposition in sheep tail. In addition, this study explored the differentiation trends of sheep tail fat cells and analyzed the communication patterns between tail fat cells and other cells. Through the integrated analysis of multi-omics data and the step-by-step comparison method (fat-tailed sheep vs. thin-tailed sheep, long fat-tailed sheep vs. short fat-tailed sheep), we aim to explain the molecular regulatory mechanism of tail fat deposition differences in sheep with different tail types. This study provides theoretical and practical basis for genetic improvement and molecular design breeding of sheep tail types, and contributes to the improvement of the single-cell atlas of ruminants.

Fig. 1

Research ideas of this study

Fig. 2

Schematic view of the procedures for data collection and analyses in the present study

Result

Screening genome-wide selection signals between fat-tailed and thin-tailed sheep by three methods

Fig. 1 shows the research ideas of this study. Fig. 2 shows the technical route of this study. After quality control, we obtained 536,156 SNPs from 911 sheep from all over the world (Table S1), and the geographical distribution of the sampling is shown in Fig. 3A. Subsequently, we used the genome-wide pairwise F_ST, π ratio and cross-population composite likelihood ratio (XP-CLR) score to calculate the selection signal between fat-tailed and thin-tailed sheep populations (Table S2) and annotated the candidate regions (Table S3). Specifically, under the top 5% threshold, the F_ST method screened out 5,153 candidate genomic regions and annotated 3,010 genes. The three candidate regions with the highest F_ST values were located on chromosomes 6, 10, and 14, respectively (Fig. 3B). The π ratio method screened out 5,151 candidate genomic regions and annotated a total of 3,238 genes. The three candidate regions with the highest π ratios were located on chromosomes 26, 6, and 12, respectively (Fig. 3C). XP-CLR screened out 5,129 candidate genomic regions and annotated a total of 4,268 genes. The three candidate regions with the highest XP-CLR scores were located on chromosomes 27, 1, and 2, respectively (Fig. 3D).

We intersected the genes annotated by the three selection signal methods and found that there were 1,899 genes identified by two methods and 550 genes identified by all three methods (Fig. 3E). Functional annotation of these genes revealed 28 candidate genes associated with adipogenesis and lipid droplet dynamics (e.g., GSK3B, FOXO4, NDEL1, GHR, SOCS2, GABRA1, CHMP4B, LOC101108715, GABRG3, SLC4A7, CHRM3, LC4A2, SLC16A3, TSHR, GNAQ, CREB3L2, ITPR1, JAK2, BMP2, PDGFD, GLIS1, ALG3, NRIP1, AR, VEGFA, PSMD1, ENPP2, and ACSL3), 20 genes associated with ECM remodeling (SVIL, ABLIM3, TNFAIP3, FGD2, LOC101111980, MYH11, MYH10, SH3YL1, MPRIP, GJB6, CTNNA2, ADD1, P2RX4, FAT2, FAT3, EPS15L1, BMP2, ACSL3, SOCS2 and VEGFA), 31 candidate genes related to cytoskeleton and cell migration (PRKCE, PXN, KLHL17, DYRK1A, MICAL3, ANK2, ARHGAP26, ADD1, FGD2, PPP2R2B, JAK2, SNTA1, PTPN3, TNFAIP3, GPC5, PIK3C2G, CTNNA2, GPC6, NDEL1, LOC101111980, MYH11, MYH10, TMOD1, MYBPC2, COL24A1, ANK2, SGCZ, CSRP1, TNNT2, TNNI1 and SNTB1), and 21 candidate genes related to immunity (ARHGEF11, RB1, GSK3B, PDGFRA, PPP3CA, SP1, GNAQ, CREB3L2, PXN, ITPR1, TSC1, TRAF2, DIABLO, CAMK2D, PRKCH, GNAQ, PRKCE, ASIC3, SOCS2, ACSL3and ETV3)(Table S4). GO and KEGG enrichment analysis of the 550 genes revealed their important roles in responding to insulin receptor (GSK3B, FOXO4, and NDEL1), growth hormone synthesis, secretion, and action-related (GHR, SOCS2, GSK3B, GNAQ, CREB3L2, ITPR1, JAK2), and parathyroid hormone synthesis, secretion, and action (ARHGEF11, SP1, GNAQ, CREB3L2, ITPR1, OXSR1) (Fig. S1A).

Fig. 3

Genome-wide selection signals in fat-tailed and thin-tailed sheep. (A) Geographical distribution of sampling for whole genome and transcriptome sequencing. A total of 663 fat-tailed sheep and 248 thin-tailed sheep whole genome data were collected, and transcriptome data of 13 fat-tailed sheep and 47 thin-tailed sheep were collected. (B) Whole-genome selective signals between fat-tailed sheep and thin-tailed sheep based on the pairwise F_ST selection test. (C) Whole-genome candidate selective regions between fat-tailed sheep and thin-tailed sheep by the log₂(π ratio) selection test. (D) Genome-wide selective signals between fat-tailed sheep and thin-tailed sheep by the XP-CLR test. The horizontal red dashed line corresponds to the genome-wide significance threshold (top 5%: F_ST = 0.1608, log₂(π ratio) = 0.5215, and XP-CLR = 0.8870). (E) The upset plot shows the number of genes and their intersection obtained from the three genome-wide selection signal analyses when the significance threshold is top 5%

Differential analysis of gene expression

First, we corrected and merged the gene expression profile data of 60 samples, and the geographical distribution of the samples is shown in Fig. 3A. Next, we compared the 13 fat-tailed samples and 47 thin-tailed samples in the merged data (Table S5), and obtained a total of 1,311 differentially expressed genes (DEGs) (Table S6), of which 119 genes were upregulated and 1,193 genes were downregulated (Fig. 4A).

Identification of module genes associated with fat-tail trait by weighted gene co-expression network analysis

To identify the module genes associated with the fat tail trait in sheep, we performed weighted gene co-expression network analysis (WGCNA) on 27,054 genes in the sheep samples. We selected 9 as suitable soft intensities (Fig. 4B), and then clustered them using the sample dendrogram and fat tail trait heat map (Fig. 4C), showing the overall situation of the fat tail group and the thin tail group in the sample (Fig. 4D). The selection of the blue module was determined by the module gene clustering (Fig. 4E). Finally, we identified 359 genes associated with the fat tail trait through WGCNA (Table S7). By crossing the DEGs with the module genes, we found that the 359 module genes were all included in the 1,311 DEGs (Fig. 4F). After functional annotation and classification of these genes, it was found that 83 genes related to adipogenesis, 44 candidate genes related to ECM remodeling, 33 candidate genes related to splicing, mRNA processing and RNA binding, 25 candidate genes related to cell adhesion and motility, and 47 candidate genes related to signal transduction were revealed (Table S8). In addition, GO and KEGG enrichment analysis of 359 genes revealed their important roles in responding to apoptosis (such as SH3GLB1, IGFBP3, PMP22), RNA binding (such as KHDRBS1, CPSF7, RBM25) and cell adhesion (such as LAMA5, POSTN, F11R) (Fig. S1B).

Fig. 4

Screening for differentially expressed genes regulating sheep tail fat deposition using multiple methods. (A) DEGs regulating tail fat deposition in sheep; red represents upregulation, blue represents downregulation. (B) Soft threshold power versus scale-free topological model fit index and average connectivity, with 9 selected as the appropriate soft power. (C) Gene clustering based on dissimilarity metrics. (D) Clustering of samples and indicator traits. (E) Association analysis between gene modules and different traits, the blue module is the most specific module for the fat-tail trait. (F) Venn plot showing the intersection of DEGs with blue module genes

Single cell atlas of sheep tail adipose tissue

Single-cell transcriptome sequencing was performed on the tail fat tissue samples of GLT and Hu, and a total of 19,471 cells were obtained, including 10,191 cells from GLT and 9,280 cells from Hu. These cells were divided into 16 clusters by uniform manifold approximation and projection (UMAP) analysis after removing the batch effect (Fig. 5A). Considering the expression pattern of cell type-specific marker genes, this study assigned an identity to each cluster, that is, the cells in the tail fat of sheep were divided into 10 cell types (Fig. 5B): Adipocyte (9059), B cells (1027), CD8 + T cells (1430), Dendritic cells (256), Hematopoietic stem cells (HSC) CD34+ (2365), Monocytes (84), NK cells (1417), Macrophages (300), Pre B cell CD34- (2491) and T cell (1042).

This study used the function “FindAllMarkers” to screen the differential genes of each cluster for cell annotation, and displayed the top 2 marker genes of each cell type (Fig. 5C). Specifically, the top five marker genes of adipocytes are FGF10, LAMA2, UST, LOC101117891 and ABI3BP; the remaining cells are mostly different types of immune cells, including B cells (ADGRL3 and MYO1B) and T cells (CCDC102B and ARHGAP152) (Table S9). In addition, to ensure the accuracy of the results, the bubble plot shows the top 20 marker genes for each cell type (Fig. S2A).

Differential expression analysis of sheep adipocytes with different tail types

In order to compare the differences between GLT and Hu tail fat at the single-cell level, we extracted their adipocytes separately for comparison. Specifically, in terms of cell number, GLT and Hu had the largest number of adipocytes, 3,724 and 5,335 respectively (Fig. 5D). In terms of gene expression, we performed differential analysis and screened out a total of 1,414 DEGs (Table S10), of which 797 were upregulated and 617 were downregulated (Fig. 5E).

The pathways significantly enriched in these 1,414 genes can be roughly divided into five categories, including cytoskeleton and cell structure-related pathways (adherens junction and regulation of cell shape), metabolism and energy-related pathways (fatty acid metabolism, glycerolipid metabolism and insulin signaling pathway), signal transduction-related pathways (PI3K-Akt signaling pathway, AMPK signaling pathway and HIF-1), extracellular matrix and cell adhesion-related pathways (extracellular matrix organization and ECM-receptor interaction), and immune-related pathways (TNF signaling pathway, cellular senescence and cancer) (Fig. S1C) (Table S11).

We crossed the DEGs obtained by the three omics and found that there were three genes (SESN1, RPRD1A and RASGEF1B) identified by all three omics, and 213 genes identified by both omics (Fig. 5F).

Fig. 5

Single-cell transcriptome landscape of tail adipose tissue in fat-tailed and thin-tailed sheep. (A) UMAP analysis of 19,471 single cells from GLT and Hu tail adipose tissue. In the UMAP map, 16 cell type clusters are marked with different colors. (B) Single-cell atlas of GLT and Hu tail adipose tissue identified a total of 10 cell types based on the expression of marker gene signatures. (C) From left to right: cell type information; number of cells of each type; number of unique molecular identifiers (UMIs) in each cell type; number of genes detected in each cell type and violin plots of marker gene expression in each cell type. (D) Comparison of the number of cells in each cell type in GLT and Hu tail adipose tissue. (E) DEGs in tail adipocytes of GLT and Hu. (F) Venn diagram showing the intersection of DEGs between fat-tailed and thin-tailed sheep at the genome, transcriptome, and single-cell levels

Expression levels of hub genes

We crossed the three omics methods in pairs and found that the genome and transcriptome (including differential analysis and WGCNA) jointly identified 7 candidate genes (such as ENPP1, MALL, MED12, MPRIP, PXN, SBNO2, SOCS2), RASGEF1B was excluded from WGCNA, while RPRD1A and SESN1 were recognized by all methods (Fig. 6A). Subsequently, in order to explore the transcription patterns of these hub genes (identified by three omics: SESN1, RPRD1A, RASGEF1B), we identified their expression levels in transcriptomes and single cells. The results showed that the expression levels of these three genes in fat tails were significantly lower than those in thin tails (Fig. 6B), and in GLT adipocytes were significantly lower than those in Hu (Fig. 6C).

In addition, the differential expression trends of SESN1 in CD8 + T cells, HSC CD34+, NK cells and Pre B cell CD34- were the same as those in adipocytes; the differential expression trends of RPRD1A in HSC CD34 + and Pre B cell CD34- were opposite to those in adipocytes; the expression trends of RASGEF1B in HSC CD34 + and NK cells were the same as those in adipocytes; the expression trends of SESN1 and RASGEF1B in HSC CD34+/Pre B cell CD34- were the same, and opposite to those of RPRD1A; the expression patterns of RPRD1A and RASGEF1B were opposite in HSC CD34+, and the same in Pre B cell CD34- (Fig. 6C).

Functional enrichment analysis

Subsequently, in order to clarify the biological processes and functions of hub genes that regulate fat deposition in sheep tails, we performed GO and KEGG enrichment analysis using 216 hub DEGs identified by at least two omics. GO results showed that these genes were significantly enriched in “extracellular matrix organization”, “protein phosphorylation”, “phosphorylation”, “barbed-end actin filament capping”, “positive regulation of endodermal cell differentiation”, “negative regulation of myofibroblast differentiation”, “cytoplasmic translation”, “collagen-containing extracellular matrix”, “extracellular region”, “mRNA binding”, “actin filament binding”, “protease binding”, “calcium ion binding”, “phosphodiesterase I activity”, “protein binding” and “peptidase activator activity”. KEGG results showed that these genes were significantly enriched in “VEGF signaling pathway”, “Insulin resistance”, “Focal adhesion”, “Cellular senescence” and “AMPK signaling pathway”. These pathways are significantly associated with fat deposition and metabolism (Fig. 6D).

Subsequently, we investigated the functions of the hub genes identified by the three omics in sheep and humans through multiple databases. In sheep, SESN1 regulates the “p53 signaling pathway” and “longevity regulation pathway” (Fig. S3A-B), has the function of initiating “L-leucine binding” and “oxidoreductase activity, acting on peroxide as acceptor”, and is also involved in “cellular response to L-leucine”, “cellular response to leucine starvation” and “negative regulation of TORC1 signaling”. RPRD1A regulates “Transcription Machinery” (https://www.kegg.jp/brite/oas03021+101109318), has the function of enabling “RNA polymerase II C-terminal domain binding”, and is involved in “mRNA 3’-end processing” (https://www.ncbi.nlm.nih.gov/gene/101109318). RASGEF1B enables “guanyl-nucleotide exchange factor activity” and participates in “Ras protein signal transduction” and “positive regulation of GTPase activity” (https://www.ncbi.nlm.nih.gov/gene/101114193). In humans, SESN1 is involved in the “p53 transcriptional gene network” (Fig. S3C), RPRD1B, as a homolog of RPRD1A, is involved in the “regulation of TCR signaling and T cell activation” (Fig. S3D), and RASGEF1A, as a homolog of RASGEF1B, is involved in the “breast cancer integrated pathway” (Fig. S3E).

Fig. 6

Expression and functional enrichment of hub genes regulating sheep fat tail formation. (A) Upset plot shows the intersection of DEGs in multi-omics analysis. (B) The expression of genes that regulate differential fat deposition in sheep tail at the transcriptome level, P < 0.05 indicates significant difference. (C) The expression of genes that regulate differential fat deposition in sheep tail at the single-cell level, * indicates P < 0.05, ** indicates P < 0.01, *** indicates P < 0.001. (D) Functional enrichment of genes with significant differences in two or more omics

Pseudo-time trajectories of sheep tail adipocytes

In order to avoid the influence of other types of cells on adipocytes, we re-clustered the adipocytes and identified four cell clusters (Fig. S2B). In order to more intuitively check the location of each subtype, we also displayed the location of these cells in the original map (Fig. 7A). Subsequently, we performed pseudo-time trajectory analysis on adipocytes. The results showed that as time progressed, the order of trajectory differentiation was cluster 1, cluster 3, cluster 2, and cluster 0 (Fig. 7B-D). The hub genes that regulate the difference in tail fat deposition showed different expression conditions as the differentiation state changed. The heatmap drawn allowed us to observe the common changes of different gene modules in pseudo-time-dependent genes at different times, and these genes followed similar dynamic trends (Fig. 7E-G, Fig. S4A). Combining the cell differentiation trajectory and the marker genes between the clusters (Table S12), it was determined that cluster 1 was a adipose stem cells (CD44), cluster 3 was a pre-adipocyte (PPARγ), and clusters 0 and 2 were different types of mature adipocyte subtypes. The bubble chart shows the marker genes of each subtype (Fig. S2C).

Communication network of sheep tail adipocytes

To explore whether the ligand-receptor interactions and signaling pathways between various types of cells differ in sheep with different tail types, this study used the CellChat tool to conduct an in-depth analysis of the single cell data of GLT and Hu. The results showed that the intercellular interaction strengths of GLT and Hu tail lipids were different (Fig. 7H). A total of 2533 significant ligand-receptor pairs were detected in GLT, involving 58 signaling pathways; a total of 2851 significant ligand-receptor pairs were detected in Hu, involving 60 signaling pathways (Table S13).

Fig. 7

Adipocytes exhibit temporal heterogeneity in their differentiation trajectories. (A) 9,059 adipocytes were subjected to UMAP analysis, and 4 clusters were marked with different colors. (B) Pseudo-time processes in adipocytes. (C-D) Adipocytes differentiation trajectory in a pseudo-time process. (E) Expression of genes regulating sheep tail fat deposition in pseudo-time trajectory. (F) Heatmap showing dynamic expression changes of genes in cell clusters. (G) Pseudo-time trajectories of adipocytes classified by state and cell clusters. (H) Communication network between GLT and Hu adipocytes

Potential regulatory network for tail fat deposition in sheep

There are 54 pathways that exist simultaneously in the communication between GLT and Hu, and 4 pathways that exist exclusively in GLT: CDH signaling pathways, MIF signaling pathways, NT signaling pathways and RELN signaling pathways; the pathways that exist exclusively in Hu are ADGRG signaling pathways, FN1 signaling pathways, GRN signaling pathways, HGF signaling pathways, Netrin signaling pathways and THBS signaling pathways. When adipocytes were used as ligand, LAMININ signaling pathways, FGF signaling pathways, ADGRA signaling pathways, and NECTIN signaling pathways were strongly expressed in both sheep, COLLAGEN signaling pathways, CD99-related signaling pathways, TENASCIN signaling pathways, MPZ signaling pathways, HSPG signaling pathways, ANGPTL signaling pathways, CXCL signaling pathways, and THY1 signaling pathways were strongly expressed only in GLT, and SEMA3 signaling pathways and BMP signaling pathways were strongly expressed only in Hu. When adipocytes were used as receptors, PTN signaling pathways, PERIOSTIN signaling pathways, THY1 signaling pathways, and VWF signaling pathways were strongly expressed in both sheep, BMP signaling pathways, CD99-related signaling pathways, MPZ signaling pathways, ADGRA signaling pathways, and NT signaling pathways were strongly expressed only in GLT, and FGF signaling pathways was strongly expressed only in Hu (Fig. S4B).

In the ADGRA pathway, Adipocytes in Hu as a ligand has a stronger effect than GLT (Fig. S5). In GLT, the signals sent by adipocytes as ligand are mainly received by adipocytes and monocytes, and the signals received by adipocytes as receptors are mainly sent by adipocytes and dendritic cells. In Hu, the signals sent by adipocytes as ligand are mainly received by B cells and monocytes, and the signals received by adipocytes as receptors are mainly sent by T cells and dendritic cells.

Discussion

Adipose tissue plays multiple roles in mammals. As an energy storage organ, adipose tissue provides energy buffering by storing lipids when food is abundant, and releases fatty acids for systemic use in states of starvation or high energy demand [20]. In addition, adipose tissue can secrete a variety of hormones (leptin and adiponectin, etc.) and cytokines, which are involved in regulating systemic metabolism, immune response and appetite control [21]. Therefore, adipose tissue is not only a single energy depot, but also a key endocrine organ for maintaining body homeostasis. In humans, body fat distribution is of great clinical importance in relation to a variety of metabolic diseases, with fat stored in the trunk being more pathogenic than fat stored in other compartments [22], and accumulation of adipose tissue in the upper body (abdominal region) being more pathogenic than accumulation of adipose tissue in the lower body (gluteofemoral region) [23]. In sheep, fat accumulation shows great differences in fat tails and thin tails. This differentiation directly affects the sheep’s body shape, production performance and immune function, making it an ideal model for studying lipid metabolism and gene regulation.

Adipose tissue expansion occurs through hyperplasia and hypertrophy [24]. In this study, genome-wide selection signals between fat-tailed and thin-tailed sheep were screened by three methods. The selective sweep test detected a group of genes related to adipocyte hyperplasia, including candidate genes that have been identified as involved in the fat-tail phenotype of sheep in previous studies (such as BMP2 [25], PDGFD [26], GLIS1 [27], and VEGFA [28]). These genes have also been detected in previous studies [1], but in this study, these genes were only detected by two selection signal analysis methods. This is due to two reasons that led to slightly different research results. The first reason is due to the difference in sample size and the population used. This study used 663 large-tailed sheep and 248 thin-tailed sheep, which has a larger sample size and breed number than previous studies (221 large-tailed sheep and 304 thin-tailed sheep [1]), which improves the representativeness and statistical power of the research results. The second reason is that previous studies used whole genome resequencing data [1], while this study used chip data. Chip data is more suitable for large-scale screening and comparison of common genes in multiple populations and has higher universality. In addition, the selective scanning test also screened out some candidate genes involved in the terminal differentiation of adipocytes, such as IGFBP-3 (interfering with PPARγ) [29], INSIG2 (involved in human adipocyte metabolism and body weight regulation) [30], JAK2 (JAK2/STAT3 pathway regulates C/EBPβ transcription) [31], ALG3 (involved in the biosynthesis of N-glucose precursors) [32] and Nrip1 (the loss of Nrip1 can reduce cell proliferation, prevent cell apoptosis, and inhibit adipogenesis) [33]. These genes and pathways are involved in the terminal adipogenic differentiation of adipocyte hyperplasia, hypertrophy and early adipogenesis, resulting in great differences in the volume of fat and the number of adipocytes in sheep tails.

In this study, 359 genes related to fat tail traits were screened by WGCNA. All of these genes were DEGs, including candidate genes that have been identified as involved in the fat tail phenotype of sheep in previous studies (such as VEGFA, SOCS2) [1]. In addition, some genes involved in adipogenesis and lipid droplet dynamics were screened, such as MAP2K1 (involved in cell proliferation, differentiation and metabolism) [34], NOS3 (regulating vascular function, adipocyte loss can enhance high-fat diet-induced hypertension) [35]. In the study of transcriptome analysis of Kazakh sheep (fat tail) and Tibetan sheep (thin tail), it was found that most of the DEGs were enriched in pathways such as fat digestion and absorption, glycine, serine and threonine metabolism [36]. In the comparison of Iranian fat-tailed sheep (Lori-Bakhtiari) and thin-tailed sheep (Zel) breeds, in addition to pathways affecting lipid metabolism, they were also enriched in MAPK signaling pathways, Wnt signaling pathways and ECM receptor interactions, which may contribute to the fat deposition in the sheep tail [37], which is consistent with our results. Among the seven genes screened simultaneously by the genome and transcriptome, SOCS2 [1], SBNO2 [38] and ENPP1 [39] have been identified in previous studies as being associated with fat deposition in sheep, while MALL (interacts with lipid raft structures in the cell membrane and participates in cell signaling) [40], MED12 (affects RNA polymerase activity) [41], MPRIP (involved in cytoskeleton regulation) [42], and PXN (cell adhesion and signal transduction) [43] affect lipid metabolism in different ways.

The single cell atlas shows the gene expression, differentiation status, and interactions of all cell types in the sheep tail fat tissue. We divided the cells in the sheep tail fat into 10 cell types, including adipocytes and 9 immune cells, represented by GLT and Hu, which are similar to the results in the single cell atlas of human and mouse adipose tissue [16]. Differential analysis of its adipocytes obtained a total of 1,414 DEGs, which directly led to pathways related to fat metabolism and lipid storage. Among them, fatty acid metabolism is directly related to the synthesis and decomposition of fat, affecting the energy metabolism process of adipocytes; the PPAR signaling pathway regulates lipid metabolism and the differentiation of adipocytes, affecting fat production and decomposition [44]. The insulin signaling pathway and insulin resistance regulate glucose and lipid metabolism and have an important impact on fat storage and degradation [45].

UMAP graphs usually reflect the differences in gene expression and functional responsibilities of cell types. In this study, some cell types are relatively far away from their precursor cells (such as B cells and Bcell_CD34- cells). This is because the transcriptional expression pattern of cells has changed after differentiation. For precursor cells or different cell subtypes, they may tend to be more like another cell type (or assume functions similar to another cell type) after differentiation [46]. Of course, this differentiation cannot be manipulated by omics or algorithms, it is the result of the inevitable selection of cells according to the environment. Different types of cells have assumed different biological functions. These reasons lead to the position of each cell in the UMAP graph. There are great differences in cell differentiation and function between pre B cells and B cells. B cells are mature and functional lymphocytes responsible for the production of antibodies and immune responses in the body. Pre B cells do not yet have the characteristics of mature B cells [47]. At the same time, both Pre-B cells and HSC are cell types in the hematopoietic system. HSC is the progenitor cell of all blood cells (including B cells) and the starting point of the hematopoietic system. Pre-B cells are the intermediate stage of HSC differentiation and belong to the early stage of B cell differentiation. Moreover, there is partial overlap in the surface markers and functions of Pre-B cells and HSC [48]. Similarly, T cells have various subtypes (such as helper T cells and regulatory T cells). In the local microenvironment of adipose tissue, the expression pattern of T cells may be affected (under conditions of obesity or metabolic diseases, adipose tissue will attract a large number of immune cells), and the expression of inflammatory or metabolic regulatory genes is also similar [49]. CD8 + T cells are often referred to as cytotoxic T cells (Cytotoxic T Lymphocytes), which are mainly responsible for killing abnormal cells such as virus-infected cells and tumor cells. Their functions are shared with NK cells (such as cytotoxicity and the ability to kill target cells) [50]. NK cells are between CD8 + T cells and T cells, reflecting their intermediate characteristics in function and gene expression profile, which determine their position in UMAP. Monocytes are myeloid precursor cells and are immature immune cells. In the tissue microenvironment, monocytes undergo migration, differentiation, and functional remodeling, and eventually form macrophages to play a role in tissue repair, pathogen phagocytosis, and immune regulation. There is a direct immune regulatory network between macrophages, NK cells, and CD8 + T cells. Chemokines secreted by macrophages (such as CCL5 and CXCL9) attract NK cells and CD8 + T cells [51]. Due to the differences in function and transcriptional patterns, these cells are in a special position in UMAP.

We performed cross-analysis on the DEGs obtained from the three omics and obtained three hub genes (SESN1, RPRD1A and RASGEF1B). We detected the expression of these hub genes at the transcriptome and single cell levels and found that their expression trends were consistent. In addition, these genes still have significant differences in certain immune cells, and their expression trends are the same or opposite to those in adipocytes, which indicates that adipocytes and immune cells in sheep tail fat work together to regulate fat accumulation and metabolism [52]. SESN1, RPRD1A, and RASGEF1B play multiple roles in adipose tissue. Specifically, SESN1 is an oxidative stress-sensing protein that activates antioxidant enzymes (such as superoxide dismutase SOD), reduce the level of reactive oxygen species, and thus protect adipocytes from damage [53]; At the same time, it also regulates PPARγ and other factors to promote the transformation of precursor adipocytes into mature adipocytes [54]. SESN1 can participate in regulating the body’s energy metabolism. When energy is insufficient, SESN1 promotes the oxidation of fatty acids to increase energy supply, while helping to maintain the energy balance of adipose tissue by enhancing the activity of the AMPK signaling pathway [55]. RPRD1A is a cyclin-dependent kinase 2B inhibitor-related protein that plays an important role in the splicing of pre-mRNA. It regulates the metabolism and energy balance of fat cells by regulating the processing of mRNA and interacting with multiple signaling pathways (such as insulin signaling pathway and the AMPK pathway) [56]. As a promoter of RAS (small GTPase), RASGEF1B (RasGEF domain family member 1B) can promote adipocyte differentiation through the PI3K/AKT and MAPK pathways. Changes in the activation of the RAS signaling pathway affect the energy metabolism of adipocytes (such as the enhancement of RAS signals promotes the synthesis and storage of fatty acids) [57]. Since the RAS signaling pathway plays a key role in the inflammatory response, the activity of RASGEF1B may affect the expression of inflammatory cytokines in adipose tissue, thereby promoting the occurrence of chronic low-grade inflammation caused by obesity.

The molecular mechanism of fat deposition differences involves the precise regulation of multiple signaling pathways and transcription factors [58]. The results of the comprehensive hub gene enrichment analysis show that these genes play a multi-level regulatory role in the deposition and metabolism of adipose tissue. VEGF and AMPK signaling pathways can effectively regulate energy balance and adipocyte differentiation, while extracellular matrix remodeling and focal adhesion may regulate fat deposition by affecting cell-cell interactions and cell morphology. Their enrichment in pathways such as insulin resistance and cell senescence suggests their potential role in regulating adipocyte proliferation, metabolic health, and adipose tissue function. In addition, multiple entries closely related to cell morphology and signal transduction indicate that these genes may provide a favorable cellular environment for fat deposition by regulating the microstructure and mechanical properties of cells [59].

Pseudo-time differentiation trajectories and dynamic expression patterns of hub genes reveal the differentiation process of different adipocyte subtypes in tail fat deposition. Pseudo-time differentiation trajectories and dynamic expression patterns of hub genes reveal the differentiation process of different adipocyte subtypes in tail fat deposition. During the differentiation of progenitor cells into mature adipocytes, the expression of different gene modules is gradually activated or inhibited over time, forming a continuum of differentiation within adipocytes. In addition, the differentiation pathways of these cell clusters suggest that different types of mature adipocytes may have unique functions in the regulation of fat deposition. The differentiation endpoints of cluster 0 and cluster 2 cells showed a division of labor in metabolic activity and energy storage, corresponding to the different characteristics of tail fat. In the communication network of adipocytes, pathways such as LAMININ signaling pathways, FGF signaling pathways, and ADGRA signaling pathways were strongly expressed in both tail types, indicating that these signaling molecules may be common regulators across tail types during tail fat formation [60]. In humans, the functional states of macrophages in obesity and normal are significantly different. The M2 anti-inflammatory state helps promote lipid storage and tissue repair, while the M1 pro-inflammatory state inhibits fat accumulation and promotes lipolysis. In sheep, the ratio of inflammatory M1 macrophages/anti-inflammatory M2 macrophages in the thin-tailed group is higher than that in the fat-tailed group [1]. In addition, other immune cells (such as T cells) also show different distribution and activation states, affecting the metabolic activity of adipocytes and lipid homeostasis. This is also the focus of our subsequent research.

Combining previous studies, we integrated the three main regulatory ways of adipocyte differentiation. (1) Insulin and IGF-1 activate the PI3K/Akt pathway through their receptors, further stimulating the expression of Necidin and GATA2 and enhancing the process of adipocyte differentiation [61]. Necidin and GATA2 jointly promote the activation of adipogenic genes, and ultimately promote the maturation and functionalization of adipocytes through the action of PPARγ [62, 63]. (2) The Wnt/β-catenin pathway is first activated through the Frizzled receptor, inhibiting the induction of adipogenic genes [64], and then activates the PI3K/Akt signaling pathway through IGF-1 and insulin receptors, thereby indirectly activating PPARγ and promoting adipocyte differentiation and maturation [61]. (3) TGF-β activates Smad3 through its receptor, and Smad3 enters the nucleus to inhibit the expression of adipogenic genes and limit adipocyte differentiation [65]. However, despite the inhibitory effect of TGF-β signals, PPARγ, as a core transcription factor for adipocyte differentiation, can still promote the maturation of adipocytes under appropriate conditions by regulating the expression of adipogenic genes [66]. In addition, SENS1 affects the energy status of cells by regulating the activation of AMPK. When cellular ATP levels decrease, AMPK is activated, initiating a series of metabolic reactions to restore energy balance [67]. Activation of AMPK promotes adipocyte proliferation and hypertrophy, and regulates adipocyte energy storage and metabolic processes [68]. RPRD1A and RASGEF1B also regulate adipocyte activity by affecting pathways such as lipid metabolism [69, 70] (Fig. 8).

Fig. 8

Potential regulatory network for tail fat deposition in sheep

Of course, this study also has some limitations. Single cell RNA sequencing data are derived from a large number of individual cells, representing repeated measurements at the technical level rather than independent biological replicates. Given this characteristic of scRNA-seq data, this study did not include biological replicates in the traditional sense. Therefore, a larger sheep sample cohort is needed for further confirmation. Secondly, although we conducted a joint analysis of large-scale genomic and transcriptomic data and used rigorous statistical methods to screen key genes at the single-cell level to enhance the reliability of the results, functional experiments are still needed in the future to further reveal the potential regulatory role of these hub genes.

Conclusions

In this study, we characterized the genomic characteristics of sheep with different tail types, detected the transcriptomic differences in tail adipose tissue between fat-tailed and thin-tailed sheep, and established a single-cell atlas of sheep tail adipose tissue. We found that the differential mechanism of sheep tail fat deposition not only involves adipocyte differentiation and proliferation, but is also closely related to cell-specific communication networks, including interactions with immune cells and tissue remodeling, to drive the typing of tail fat. The results of this study provide insights into further revealing the molecular mechanisms of differences in tail fat deposition in sheep with different tail types, and provide a deep explanation for the development and functional regulation of adipocytes. In addition, potential molecular markers for regulating sheep tail type provide a basis for the formulation of new breeding strategies and genetic improvement for specific needs.

Methods

Genomic data Preparation and preprocessing

We collected Ovine Infinium HD SNP BeadChip (600 K) data of 911 sheep (10 breeds, 663 fat-tailed sheep and 248 thin-tailed sheep) from the isheep database (https://ngdc.cncb.ac.cn/isheep/) [71]. Details of the data are provided in Table S1.

Subsequently, we used PLINK v1.9 to perform quality control on the SNP dataset [72], and SNPs that met any of the following conditions were removed: SNP missing rate higher than 10% (--geno 0.1), minimum allele frequency (MAF) lower than 5% (--maf 0.05). After quality control, 536,156 SNPs were retained for subsequent analysis.

Selection signal analysis

To identify potential selection signatures between fat-tailed and thin-tailed populations, we used three methods to scan the genome for selective sweeps. First, we used VCFtools v0.1.13 [73] to calculate the genome-wide pairwise F_ST values and − log₂(π_fat tail/π_thin tail) values using a sliding window method (50 kb sliding windows with 25 kb steps). Second, we performed a cross-population composite likelihood ratio test (XP-CLR) using the XP-CLR v1.1.2 [74] program with the same sliding window and step size settings. We then intersected the genes annotated by the three methods.

Transcriptome data and preprocessing

We collected transcriptome data of 60 samples (13 fat-tailed sheep and 47 thin-tailed sheep) from the NCBI database (https://www.ncbi.nlm.nih.gov/). Detailed information of the samples is given in Table S2.

Differential expression analysis

We used the R package “limma [75]” to compare the expression profile data of fat-tailed samples and thin-tailed samples to identify the DEGs between the two groups. The threshold was|log2FoldChange| >1, and P-adjust < 0.05 was considered significant, resulting in a total of 1311 DEGs.

Weighted gene co-expression network analysis

We performed weighted gene co-expression network analysis (WGCNA) on 27,054 genes in sheep samples using the R package “WGCNA [76]”. To ensure that the constructed co-expression network was close to a scale-free distribution, we selected 9 as soft power when performing WGCNA on the samples, obtained 3 modules, and then calculated their relationship with the samples. Finally, we determined to select genes from the blue module, and the 359 genes in the blue module.

Collection of single-cell adipose tissue and Preparation of cell suspension

In this study, a 3-month-old weaned male lamb from Guangling Large-tailed sheep and Hu sheep was selected, and samples of tail adipose tissue were collected. Animal surgical procedures were performed in accordance with the guidelines of the College of Animal Science and Veterinary Medicine, Shanxi Agricultural University (Taigu, China). All trial protocols were reviewed and approved by the institution (ethics committee approval reference number: SXAU-EAW-2022 S.UV.010009). After obtaining the adipose tissue, the tail adipose tissue of Guangling big-tailed sheep and Hu sheep was digested using type II collagenase, and then resuspended in culture medium to make a single-cell suspension that can be used for the single-cell sequencing process.

Formation of single-cell gel Bead-in emulsions and construction of library

Perform cell counting and viability tests on the single cell suspension to ensure that the cell viability is higher than 80%, and adjust the cell concentration to 1000 cells/µL. Subsequently, the prepared cell suspension is used to encapsulate beads with cell barcodes and cells in droplets using a microfluidic chip. The cells are lysed in the droplets, and the mRNA in the cells is connected to the cell barcode on the beads to form single cell gel bead-in-emulsions (GEMs). The reverse transcription reaction is carried out in the droplets, and then the emulsion is broken to construct the cDNA library. The cell barcode on the library sequence is used to distinguish which cell the target sequence comes from, and the sample index on the sequence is used to distinguish which sample the target sequence comes from. (Fig. S6A shows the single-cell RNA-Seq workflow of the GemCode platform.)

Cells and reaction reagents are placed in one channel on the microfluidic chip, and beads are placed in another channel, forming GEMs together. Reverse transcription is performed independently in each GEM, and then the labeled cDNAs are mixed and amplified for library construction. Fig. S6B shows a schematic diagram of the 10X library construction process. The single-cell sequencing was completed by Shenzhen BGI Co., Ltd. in this study.

Single cell sequencing data preprocessing and quality control

We used Cell Ranger software (https://support.10xgenomics.com/single-cell-gene-expression/downloads/latest) to align single-cell data with the ARS-UI_Ramb_v3.0 (GCF_016772045.2) reference genome (https://www.ncbi.nlm.nih.gov/datasets/genome/GCF_016772045.2/), followed by library splitting, cell splitting, and output of expression quantification matrix. Subsequently, we used the “Seurat [77]” package for quality control and subsequent analysis. This study excluded genes expressed by less than 10 cells in the sample, cells expressing less than 200 genes, and cells whose UMI counts of mitochondrial genes accounted for more than 20% of the total UMI counts.

After quality control, 19,471 cells were retained (10,191 cells in GLT sheep and 9,280 cells in Hu sheep). The function “NormalizeData” was used to correct the number of gene reads for each cell, and the corrected values were logarithmized. Then, the function “FindVariableFeatures” was used to screen out the top 2,000 highly variable genes (HVGs), and the function “ScaleData” was used to normalize the gene expression. On this basis, this study further used HVGs for principal component analysis to ensure that most of the variation information was retained while reducing the data dimension.

Analysis pipeline of sheep tail fat single-cell sequencing data

Remove batch effects

In order to eliminate the batch effects between different data sets, this study used the Runharmony function in the Harmony software to iteratively fine-tune similar cells in different batches [78], thereby correcting the batch effects and maintaining true biological variations, integrating cells in different data sets into the same space. This study used lambda = 1 to prevent over-correction of the data, and used the first twenty principal components for subsequent UMAP analysis.

Data dimensionality reduction and cell clustering

We used the UMAP [79] method to reduce the dimension and visualize the cells after removing the batch effect, and performed unsupervised clustering on the cells through the FindNeighbors and FindClusters functions, setting the resolution to 0.5 and dividing the cells into 16 clusters.

Marker gene identification and cell type annotation

In order to accurately determine the biological identity of each cell cluster, this study used the “FindAllMarkers” function to perform differential analysis on each cell cluster to identify the marker genes of each cluster. Based on the use of the “SingleR [80]” package and the CellMarker database [81], combined with existing literature for reference [82, 83, 84], the type of each cell was annotated to improve the accuracy of cell type identification. Subsequently, we extracted 9059 adipocytes, used the “FindVariableFeatures” function to screen out the top 2,000 HVGs, used “RunUMAP” to re-cluster the adipocytes, and then used the “FindAllMarkers” function to perform differential analysis on each cell cluster to identify adipocyte subtypes. In addition, we used the dotplot function to draw bubble plot of marker genes for each cell type and adipocyte subtype to enhance the reliability of the research results.

Differential expression analysis of adipocytes

In order to explore the differences between the tail fat of GLT and Hu at the single-cell level, we extracted the fat cells of the two species separately for differential analysis, and P-adjust < 0.05 was considered to be significantly different.

Pseudo-time trajectory analysis and cell communication analysis

In order to reveal the differentiation trajectory of sheep tail adipocytes, this study used the Monocle tool for pseudo-time series analysis [85]. First, we extracted the expression matrix information, gene information, and cell type information of adipocytes, and estimated the size factor through the function “estimateSizeFactors” to standardize the difference in mRNA between cells. The function “estimateSizeFactors” was used for subsequent differential analysis, and low-expression genes were filtered to reduce noise. Subsequently, the DEGs of each subtype were used as the basis for sorting, and the reverse graph embedding algorithm was used to reduce the dimension of the data. Then, the “orderCells” function was used to construct the pseudo-time trajectory and arrange the cells according to pseudo-time.

In this study, the CellChat tool [86] was used to analyze intercellular communication between different cell types in sheep tail adipose tissue to construct a probabilistic network of intercellular signaling and to investigate the ligand-receptor interactions between them.

Gene expression analysis

We crossed the DEGs selected from the three omics and then analyzed the expression of the hub genes. Specifically, in transcriptomics, we obtained hub gene expression data and tail type grouping data from sheep samples, performed significant difference analysis on them using R software, and drew violin plots. In single cell transcriptomics, we used the data after quality control and normalization process mentioned above to extract hub gene expression information and cell type information, and used R software to perform significant difference analysis and draw violin plots. It is worth mentioning that in order to ensure the scientificity, reliability and accuracy of the experimental results, we eliminated cells with 0 expression of each gene before performing significant difference analysis on this gene. Gene loss or non-expression will increase the cell base number and redundant information, leading to errors or deviations in the statistical process. In order for a gene to function, it must comply with the central dogma of replication, transcription, and translation. This process cannot be completed when the expression level is 0. Therefore, this study excluded them when conducting the significance analysis and did not consider these cells. In the analysis of significance of differences, P < 0.05 was considered statistically significant.

Functional enrichment analysis

We performed functional enrichment analysis on the differentially expressed genes screened in the three omics, and then we also performed enrichment analysis on the genes co-existing in the three omics. Specifically, we used DAVID (https://david.ncifcrf.gov/) online software to perform gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis on DEGs of different ranges, where a P value less than 0.05 was considered significant. In addition, in order to determine the three genes simultaneously identified by the three omics, we used the KEGG database (https://www.kegg.jp/), the Wikipathway database (https://www.wikipathways.org/), and the NCBI database (https://www.ncbi.nlm.nih.gov/) to perform functional annotation on the three genes simultaneously identified by the three omics.

Statistical analysis

All statistical analyses and visualizations were performed using R v4.2.0 software. Statistical analysis of multiple groups of data was performed using analysis of variance [38], and comparisons between two groups of data were performed using T-test (with parameters) and Wilcoxon rank-sum test (without parameters). For all statistical analyses, P < 0.05 was considered statistically significant.