RNA-seq analysis identifies key genes and signaling pathways involved in androgen promotion of sebaceous gland proliferation in Hetian sheep

Abstract

Sebaceous glands (SGs), essential elements of the skin barrier in sheep, are tightly regulated by the androgen signaling pathway. This study investigated how androgens influence SGs hyperplasia in sheep by combining morphological assessments with transcriptomic analysis. Histological examination revealed a significant increase (P < 0.05) in both the number and average size of SGs in the testosterone-treated groups across four concentration gradients compared to the control group. Transcriptome sequencing was performed on fifteen skin samples, revealing that, compared to the control group, 371 genes were upregulated and 115 were downregulated across the four treatment groups. The upregulated genes were predominantly enriched in lipid metabolism pathways, while the downregulated genes were mainly associated with keratin filament structures. Protein-protein interaction (PPI) network analysis identified ten hub genes involved in mitochondrial β-oxidation and lipid synthesis. These findings were further supported by gene set enrichment analysis (GSEA), which indicated significant activation of fatty acid metabolic pathways (P < 0.05). The results of this study suggest that androgens may be involved in SGs hyperplasia in sheep by regulating mitochondrial lipid metabolism-related pathways.

Keywords: GSEA; Hetian sheep; Masson staining; RNA-seq; Skin.

Fig. 1

Masson staining analysis showing the impact of treatment with different concentrations of testosterone on the skin structure of Hetian sheep. (a) The skin structure of CG, TG1, TG2, TG3 and TG4 groups. The magnification of the upper image is 100 times, and the magnification of the lower image is 400 times; (b) Comparative analysis of SG counts between subject groups; (c) Comparative analysis of SG area between subject groups; (d) Statistical analysis of collagen fiber area in each group. Arrowheads indicate SGs, while arrows indicate hair follicles. Scale bars: the upper image: 100 μm; the lower image: 50 μm. *significant difference (P < 0.05), **Highly significant difference (P < 0.01).

Fig. 2

Comparative analysis of gene expression patterns in the skin of Hetian sheep between the TG groups and CG. (a) Number of DEGs between each TG and CG; (b) Volcano plot showing changes in DEGs for each TG compared with CG.

Fig. 3

Analysis of DEGs, expression trends, and GO enrichment across the treatment groups. (a) Venn diagram of DEGs across the four treatment groups; (b) Heatmap showing the expression of the 486 identified genes; (c) Expression patterns of the 371 upregulated genes; (d) GO functional enrichment analysis of the 371 upregulated DEGs; (e) Expression patterns of the 115 down-regulated genes; (f) KEGG pathway enrichment analysis of the 115 down-regulated DEGs.

Fig. 4

Results of keratin filament gene expression. (a) List of keratin filament genes; (b) Heatmap of keratin filament expression levels across the experimental groups.

Fig. 5

Construction and analysis of the DEG PPI network. (a) Protein interaction network of the upregulated DEGs; (b)The color intensity (light to dark) corresponds to Degree values (from 1 to 25), with larger values indicated by darker colors. Higher Degree values also suggest greater biological relevance; (c) Protein interaction network from the A network of the down-regulated DEGs; (d) GO analysis of down-regulated DEGs; (e) KEGG pathway analysis of down-regulated DEGs; (f) GO analysis of upregulated DEGs; (g) KEGG pathway analysis of upregulated DEGs.

Fig. 6

Hub gene analysis of upregulated DEGs. (a) Network of the top 10 hub genes and their associated genes; (b) Sub-network of surrounding genes, color-coded according to Degree values which ranged from 1 to 16; (c) Ranking of the 10 hub genes based on Degree values (18 to 25); (d) Functional enrichment analysis of the 10 hub genes was conducted using the GO database. (e) Pathway enrichment analysis of the 10 hub genes was performed using KEGG.

Fig. 7

GO analysis and gene-pathway enrichment network analysis by GSEA for the TG2 group compared with the CG group. Significant up-regulation of (a) cellular lipid metabolic processes (P = 0.000152, NES = 1.421); (b) fatty acid metabolic processes (P = 0.000169, NES = 1.646); (c) lipid catabolic processes (P = 0.000169, NES = 1.833) and (d) lipid metabolic processes (P = 0.000148, NES = 1.306) in TG2; (e) Bar graphs were plotted using the 20 processes exhibiting the lowest FDR values; (f) Gene-pathway enrichment network analysis plots were constructed using the 10 processes exhibiting the lowest FDR values.

Fig. 8

Analysis of KEGG pathways and gene-pathway enrichment networks by GSEA for the TG2 group compared with the CG group. (a) Fatty acid metabolism was significantly upregulated in TG2 (P = 0.011810, NES = 1.576); (b) Hedgehog signaling pathway was significantly decreased in TG2 (P = 0.005917, NES=−1.590); (c) Peroxisome pathway was enhanced considerably in TG2 (P = 0.000376, NES = 1.717); (d) TGF-beta signaling pathway was significantly decreased in TG2 (P = 0.007490, NES=−1.508); (e) Bar graphs showing all the apparent pathways; (f) Gene-pathway enrichment network analysis plot.

Fig. 9

Verification of RNA-seq results by qRT-PCR. The left axis indicates the relative gene expression level (mean ± SEM), while the right axis shows RNA-seq Count values. Bar graphs and line charts show qPCR results and RNA-seq results, respectively.

Fig. 10

Immunofluorescence detection of ACSL1 in the SGs of treatment and control groups. Scale: 100 μm, magnification 100×.

Fig. 11

Immunofluorescence detection of CERS4 in the SGs of treatment and control groups. Scale: 100 μm, magnification: 100×.

Fig. 12

Fluorescence intensity analysis. (a) Relative fluorescence intensity of ACSL1 in the SGs of each group; (b) Relative fluorescence intensity of CERS4 in the SGs of each group. (mean ± SEM).