Identifying differentially expressed genes in goat mammary epithelial cells induced by overexpression of SOCS3 gene using RNA sequencing

The suppressor of cytokine signaling 3 (SOCS3) is a key signaling molecule that regulates milk synthesis in dairy livestock. However, the molecular mechanism by which SOCS3 regulates lipid synthesis in goat milk remains unclear. This study aimed to screen for key downstream genes associated with lipid synthesis regulated by SOCS3 in goat mammary epithelial cells (GMECs) using RNA sequencing (RNA-seq). Goat SOCS3 overexpression vector (PC-SOCS3) and negative control (PCDNA3.1) were transfected into GMECs. Total RNA from cells after SOCS3 overexpression was used for RNA-seq, followed by differentially expressed gene (DEG) analysis, functional enrichment analysis, and network prediction. SOCS3 overexpression significantly inhibited the synthesis of triacylglycerol, total cholesterol, non-esterified fatty acids, and accumulated lipid droplets. In total, 430 DEGs were identified, including 226 downregulated and 204 upregulated genes, following SOCS3 overexpression. Functional annotation revealed that the DEGs were mainly associated with lipid metabolism, cell proliferation, and apoptosis. We found that the lipid synthesis-related genes, STAT2 and FOXO6, were downregulated. In addition, the proliferation-related genes BCL2, MMP11, and MMP13 were upregulated, and the apoptosis-related gene CD40 was downregulated. In conclusion, six DEGs were identified as key regulators of milk lipid synthesis following SOCS3 overexpression in GMECs. Our results provide new candidate genes and insights into the molecular mechanisms involved in milk lipid synthesis regulated by SOCS3 in goats.


Introduction
Ruminant milk is an important source of dietary nutrition that is rich in fatty acids for the human body.Compared to cow milk, the fat globules in goat milk are smaller in diameter and are easily digested by the human body, especially infants and children (1).Goat milk is rich in unsaturated and short-and medium-chain fatty acids, which have preventive and auxiliary therapeutic effects on coronary heart disease and diabetes (2).However, the functions of the key genes and the mechanisms of fatty acid synthesis in goat milk remain unclear.Therefore, discovering important molecules associated with lipid metabolism will contribute to improving the gene regulation network of milk fat synthesis and provide an experimental basis and new molecular targets for the improvement of milk quality.
The suppressor of cytokine signaling 3 (SOCS3) is a negative regulator that blocks multiple cellular signaling processes that can be activated by various cytokines, growth factors, and hormones.During the dry lactation period, low expression of SOCS3 promotes mammary gland development and milk synthesis by enhancing prolactin signaling, suggesting that SOCS3 participates in the regulation of lactation in dairy cows (3).SOCS3 negatively regulates Janus kinase 2 (JAK2)/signal transducer and activator of transcription 5 (STAT5) signaling pathway to inhibit milk synthesis and cell proliferation in bovine mammary epithelial cells (BMECs) (4).JAK2/ STAT5, the mammalian target of rapamycin, and sterol regulatory element-binding protein (SREBP) signaling pathways regulate milk fat synthesis by inhibiting SOCS3 expression in BMEC lines (5).In particular, n-3 and n-6 fatty acids inhibit SOCS3 mRNA expression and downregulate nuclear factor-kappa B (NF-κB) and interleukin-6 levels in mice with high-fat diet-induced obesity (6).In addition, low STAT5 activity downregulates SOCS3 expression and inhibits proliferation and lactation in BMECs (7).
Previous studies have shown that SOCS3 is closely associated with milk production in dairy livestock.However, the underlying mechanism through which SOCS3 regulates milk lipid synthesis remains unclear.Therefore, in this study, we overexpressed SOCS3 in goat mammary epithelial cells (GMECs) and screened for downstream molecules associated with milk lipid synthesis using RNA sequencing (RNA-seq).Our results provide a reference for follow-up studies on SOCS3-mediated key candidate genes and regulatory mechanisms of milk lipid synthesis in dairy livestock.

Ethics statement
Wanlin white goats were obtained from Hefei Boda Animal Husbandry Science and Technology Development Co., Ltd.The animal study protocol was approved by the Animal Care and Use Committee of Anhui Agricultural University (SYXK 2016-007).

Construction of SOCS3 overexpression vector
The coding region sequence of goat SOCS3 (GenBank no.XM_018063683.1)was obtained from the National Center for Biotechnology Information database. 1The goat SOCS3 coding sequence was subcloned into a PCDN3.1 vector between the BamHI and EcoRI restriction enzyme sites to generate the PC-SOCS3 vector.The primer sequences (cleavage sites are underlined) used for PC-SOCS3 vector construction are as follows: forward, CCGGAATTCATGGTCACCCACAGCAAGT; reverse, CGCGGATCCCTAAAGCGGGGCATCGTACT.

Total RNA extraction and RNA-seq
After transfection for 48 h, total RNA from the SOCS3 overexpression (n = 3) and control groups (n = 3) were extracted using a Total RNA Extraction Kit (R1200, Solarbio, Beijing, China), and the samples were preserved on dry ice (−78.5°C).Library construction and RNA-seq were performed by Novogene Science and Technology (Beijing, China).RNA quality and integrity were assessed using an RNA Nano 6,000 Assay Kit on a Bioanalyzer 2,100 system (Agilent Technologies, CA, United States).The RNA concentration was >300 ng/μL, and the RNA integrity number value was >9 (Supplementary Table S2).Library fragments were purified using the AMPure XP system (Beckman Coulter, Beverly, MA, United States) to select 370-420 bp-long cDNA fragments.The high-quality RNA samples were sequenced using an Illumina NovaSeq 6,000 system, generating 150-bp double-end reads.The sequencing depth was approximately 6.62 G, and the GC content was approximately 51.93% (Supplementary Table S3).

RNA-seq data analysis
Raw reads obtained by RNA-seq were filtered to remove low-quality reads and obtain clean ones.Clean paired-end reads were mapped to Capra hircus reference genome ARS1.2 (GCF001704415.2) 2 using HISAT2 software.Gene expression was determined using the million mapped reads per kilobase value.R package (3.22.5) was used to calculate Pearson's correlation coefficients for the samples.DESeq2 software (1.20.0) was used to analyze the differential expression between the control and SOCS3 overexpression groups.The Benjamini-Hochberg algorithm was used to adjust the p-value (padj) to control for the false discovery rate.FoldChange >1.5 (| log2FoldChange | > 0.59) and adjusted p-value <0.05 were set as the significance threshold for differential expression.The Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses were performed using the ClusterProfiler R package (3.8.1), which corrected for gene length bias.Functional categories with an adjusted p-value <0.05 were considered statistically significant.Enrichment of GO terms and KEGG pathways with genes differentially expressed between treatments was analyzed based on the cumulative hypergeometric distribution.

RNA-seq data verification
Total RNA was extracted from the control and SOCS3 overexpression groups, and reverse transcription was performed using the PrimeScript RT Kit (RR047, Takara, Shiga, Japan).Quantitative real-time PCR (qPCR) was performed on an Applied Biosystems StepOne Plus system (Thermo Scientific, Waltham, MA, United States) using a TB Green kit (RR820, Takara, Shiga, Japan).Ubiquitously expressed transcript (UXT) and ribosomal protein S9 (RPS9) were used as reference genes (9).The primers used for qPCR are listed in Supplementary Table S1.Data were normalized to internal controls and analyzed using the relative quantification (2 -ΔΔCt ) method.Data are shown as mean ± standard error and are statistically analyzed by an unpaired t-test using SPSS 20.0 software (IBM, Chicago, IL, United States).Differences were considered statistically significant for p < 0.05 (*p < 0.05, **p < 0.01).

Milk lipid synthesis assay
The quantitative results from goat mammary tissue showed that the abundance of SOCS3 mRNA was highest in the early and peak 2 https://www.ncbi.nlm.nih.gov/lactation periods and lowest in the dry lactation period (Figure 1).SOCS3 gene expression decreased globally from the early to dry lactation periods.
Transfection efficiency was detected by qPCR, and the expression of SOCS3 in the overexpression group was 1,050-fold higher than that in the control group (Figure 2A).To understand the effect of SOCS3 on milk lipid synthesis, we quantified the levels of TAG, T-CHO, NEFA, and lipid droplets.We found that the levels of TAG, T-CHO, NEFA, and lipid droplets decreased in the SOCS3 overexpression group (Figures 2B-E).Therefore, SOCS3 inhibits milk lipid synthesis in GMECs.

Evaluation of RNA-seq data quality
The RNA-seq was used to explore the regulatory network of SOCS3 that affects milk lipid synthesis.When the six samples were grouped, the correlation coefficient (R 2 ) between the SOCS3 overexpression and control groups was approximately 0.99, indicating that the three biological samples in each group showed similarly high performance (Supplementary Figure S2; Figure 3).After filtering the raw reads with adapter sequences and low-quality reads, a total average of 44,091,933 clean reads were obtained with phred base quality scores of Q20 > 97% and Q30 > 94% (Supplementary Table S3).

Differentially expressed gene analysis
While verifying the difference in gene expression due to SOCS3 overexpression in GMECs, 430 DEGs were identified, including 204 upregulated and 226 downregulated genes.Volcano plots of DEGs between the PCDNA3.1 and PC-SOCS3 groups are shown in Figure 4.A heatmap of the hierarchical clustering analysis of DEGs between the two groups is shown in Supplementary Figure S3.Lipid synthesis is Suppressor of cytokine signaling 3 (SOCS3) gene expression in goat mammary gland of different lactation stages.Mammary gland samples were obtained from Wanlin white goats during early, peak, middle, late, and dry lactation stages.Data were normalized to the early lactation stage.Different lowercase letters represent significant differences in SOCS3 levels (p < 0.05).

Validation and dynamic expression of DEGs
To verify the DEGs involved in lipid synthesis, cell proliferation, and apoptosis and the RNA-seq results, we selected six DEGs, comprising STAT2, FOXO6, CD40, MMP11, MMP13, and BCL2, for validating the RNA-seq data by qPCR.The qPCR data corroborated the RNA-seq findings, indicating the reliability of the RNA-seq data (Figure 5).

Functional annotation of DEGs
The GO and KEGG enrichment analyses of DEGs were performed.The GO functional annotation was performed to identify the biological functions of the DEGs.A total of 430 DEGs were enriched in 404 GO terms, including 190 in biological process terms, 45 in cellular component terms, and 169 in molecular function terms (Figure 6).Pathway analysis of the DEGs associated with lipid synthesis was performed using the KEGG database.The upregulated genes were mainly enriched in chemical carcinogenesis-reactive oxygen species and prion disease, whereas the downregulated genes were enriched in Epstein-Barr virus infection, influenza A, nucleotide-binding oligomerization domain (NOD)-like receptor, and tumor necrosis factor (TNF) signaling pathways (Figure 7).Gene network analysis showed that multiple DEGs were enriched in the JAK/STAT and NF-κB signaling pathways, which were closely associated with lipid metabolism, cell proliferation, and apoptosis progress (Figure 8).

Discussion
Milk lipid synthesis in dairy livestock is regulated by various signaling molecules (10).As a negative regulator of signal transduction, SOCS3 affects milk lipid metabolism in ruminants (11).In this study, SOCS3 overexpression decreased milk lipid synthesis in GMECs.We then explored SOCS3 regulation of potential downstream signaling molecules involved in milk lipid synthesis by RNA-seq.In BMECs, SOCS3 binds to the phosphorylated tyrosine residues of JAK2 and inhibits milk lipid metabolism by blocking STAT5 activity (12).SOCS3 can inhibit the binding of leptin and its receptor during feedback, resulting in pathway blockage of the signal transduction of leptin, insulin, and growth hormone signals in mice (13).SOCS3 disrupts fatty acid synthesis by reducing SREBP activity, thus validating its involvement in milk lipid synthesis in mammary epithelial cells.In this study, the GO and KEGG analyses showed that several DEGs were enriched in lipid metabolism and cell survival following SOCS3 overexpression.STAT2 plays a critical role in immune responses to extracellular and intracellular stimuli and upregulates the expression and secretion of pro-inflammatory mediators (14).In particular, STAT2 promotes the transcription and expression of acetyl-CoA carboxylase, thus enhancing lipid synthesis in colorectal cancer cells (15).Fatty acid binding protein 4 upregulates macrophage inflammation-related interleukin-6 (IL-6) level and regulates lipid metabolism in atherosclerosis induced by activation of the JAK2/STAT2 pathway (16).Unphosphorylated STAT2 bridges the interferon-stimulated response and NF-κB elements in the IL-6 promoter and increases cancer cell survival by enhancing IL-6 expression (17).SOCS3 and SOCS1 lead to sustained cytokine activation of STAT2 and contribute to the progression of inflammation in mouse macrophages (18).SOCS3 is a negative regulator of the JAK/STAT pathway and is a cytokine-inducible SH2 domain-containing protein involved in cell proliferation, differentiation, apoptosis, and immune regulation (19).In this study, SOCS3 overexpression significantly inhibited lipid synthesis, which was speculated to be mediated by STAT2 downregulation.
FOXO6 is a transcription factor that mediates insulin signaling during glucose and lipid metabolism.FOXO6 constitutes a distinct route by which the liver orchestrates insulin-dependent regulation of gluconeogenesis and triglyceride-rich very low-density lipoprotein production.Increased hepatic FOXO6 activity activates glucose production and triglyceride secretion in response to fasting (20).Phosphorylated FOXO6 promotes resistance to oxidative stress and inhibition of pro-inflammatory mediators by facilitating the nuclear translocation of NF-κB induced by lipopolysaccharides in aged rats (21).FOXO6 activation induces hepatic lipid accumulation by increasing the expression and transcriptional activity of peroxisome proliferator-activated receptor gamma during endoplasmic reticulum stress in HepG2 cells (22).Furthermore, it triggers cellular triglyceride-mediated lipid accumulation in the livers of aged rats fed a high-fat diet, leading to hepatic steatosis and hyperglycemia.Conversely, the expression of genes involved in lipid synthesis is reduced in FOXO6-knockout mice, resulting in decreased hepatic lipid accumulation (23).FOXO6 promotes inflammation by activating cytokine IL-1β and induces lipid triglyceride accumulation in the mouse liver and human hepatocellular carcinomas (24).In our study, SOCS3 overexpression significantly inhibited the accumulation of triglycerides, which may be associated with FOXO6 downregulation.
Survival affects lactation and lipid metabolism in mammary epithelial cells.In this study, SOCS3 overexpression significantly upregulated the mRNA levels of BCL2, MMP11, and MMP13 and   downregulated the mRNA levels of CD40, which are involved in cell proliferation and apoptosis.As an inhibitor of apoptosis, the loss of BCL2 affects lipid metabolism, and BCL2 probably has a close relationship with lipid peroxidation and reactive oxygen species production (25).Furthermore, BCL2 is a downstream gene of the phosphatidylinositol-4,5-bisphosphate 3-kinase (PI3K)/protein kinase B (AKT) signaling pathway, which is involved in the cell cycle, metabolism, and apoptosis.In our study, KEGG enrichment analysis showed that the PI3K/AKT pathway was downregulated, which was consistent with the expression trend of BCL2.Activation of the PI3K/AKT pathway enhances cell proliferation in mice and GMECs (26).The PI3K/AKT pathway regulates milk lipid metabolism, cell proliferation, and fatty acid synthesis in GMECs (27).MMP11 overexpression in macrophages promotes the proliferation and migration of breast cancer cells and monocyte recruitment (28).MMP13 contributes to the proliferation, migration, and anchorage-independent clonogenicity of mouse mammary tumorigenic cell lines (29).The tumor necrosis factor Top 10 Gene Ontology (GO) terms of DEGs determined using the GO analysis.BP (biological processes), CC (cell components), and MF (molecular functions).receptor family member CD40 is critical for the immune system and induces tumor cell-specific apoptosis (30).However, how the regulation of BCL2, MMP11, MMP13, and CD40 mediated by SOCS3 overexpression affects lipid metabolism and the complex relationship between cell survival and lipid metabolism remains unclear.Interestingly, functional enrichment analysis showed that SOCS3 overexpression downregulated multiple apoptosis-related biological processes, such as the TNF-and NOD-like receptor signaling pathways.The TNF signaling pathway includes signal transducers and adaptor proteins.A transcriptome study showed that Fumonisin B1 induces apoptosis via the TNF signaling pathway in porcine kidney cells (31).Excessive inflammatory overload-associated activation of the TNF pathway ultimately induces hepatic lipid accumulation in high-fat diet-fed mice (32).NOD-like receptor X1 regulates multiple cellular processes, including antiviral immunity, apoptosis, reactive oxygen species generation, and mitochondrial metabolism (33).Mitochondrial NOD-like receptors control apoptosis in response to intrinsic and extrinsic cues in murine embryonic fibroblasts and human embryonic kidney 293 cells (34).However, whether the regulation of lipid metabolism by SOCS3 is associated with the TNF-and NOD-like receptor pathways needs to be further explored.

Conclusion
Our study identified 430 DEGs, including 226 downregulated and 204 upregulated genes, by overexpressing SOCS3 in GMECs.Among these, STAT2, FOXO6, BCL2, MMP11, MMP13, and CD40 may be key regulatory genes involved in milk lipid metabolism.These genes play crucial roles in lipid synthesis, cell proliferation, and apoptosis.These results provide a reference for follow-up studies on the functions of key candidate genes involved in milk lipid synthesis regulated by SOCS3 in goats.

FIGURE 2 SOCS3
FIGURE 2 SOCS3 overexpression inhibits lipid synthesis in goat mammary epithelial cells (GMECs).(A) Overexpression efficiency of SOCS3 gene.(B) Effect of SOCS3 overexpression on triacylglycerols.(C) Effect of SOCS3 overexpression on total cholesterol.(D) Effect of SOCS3 overexpression on nonesterified fatty acids.(E) Effect of SOCS3 overexpression on lipid droplets determined using Oil Red O staining.

FIGURE 3
FIGURE 3Heat map of the correlation coefficient between RNA sequencing (RNA-seq) samples.The numbers in the box represent the correlation coefficient of the two samples.

FIGURE 4
FIGURE 4Volcano plot of differentially expressed genes (DEGs) between the PCDNA and PC-SOCS3 groups.The red dots represent significant upward adjustments, the green dots represent significant downward revisions, and the gray dots indicate no difference.Filter criteria: Adjusted p-value <0.05 and FoldChange >1.5 (| log2FoldChange | > 0.59).Genes in the middle of the two dotted lines and below the horizontal dotted line are not DEGs.

FIGURE 5
FIGURE 5Validation of DEGs by quantitative real-time PCR (qPCR) in the SOCS3 overexpression group.

FIGURE 7 The
FIGURE 7The Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis of DEGs.(A) Upregulated gene enrichment pathways.(B) Downregulated gene enrichment pathways.Adjusted p-value <0.05 indicates the significance of pathway enrichment.