Construction of the ceRNA Regulatory Network Associated with Milk Fat Metabolism

TL;DR

This study identifies circular RNAs and their regulatory networks involved in milk fat metabolism in dairy cows, offering new ways to improve milk quality through molecular breeding.

Contribution

The study constructs and validates a ceRNA regulatory network linked to milk fat metabolism in dairy cows.

Findings

Differential expression of circ_0009058, circ_0004021, circ_0011934, and circ_0008056 was confirmed in bovine mammary epithelial cells.

The ceRNA networks circ_0004021/bta-miR-541/PTPN6 and bta-miR-10175-3p/RHOA reduce lipid accumulation, while bta-miR-2309 increases it.

The study provides a theoretical foundation for improving milk quality through molecular breeding.

Abstract

The fat content and composition of milk directly influence its flavor and nutritional value. This study analyzed circular RNAs (circRNAs) in bovine mammary epithelial cells with high and low milk fat percentages, identifing significant differential expressions of circ_0009058, circ_0004021, circ_0011934, and circ_0008056. Experiments confirmed the presence of these circRNAs in the mammary tissues of dairy cow and suggested their potential regulatory roles in milk fat metabolism via the competing endogenous RNA (ceRNA) mechanism. A ceRNA regulatory network was further constructed, and the existence of the circ_0004021/bta-miR-541/PTPN6, circ_0008056/bta-miR-2309/ERBB3, and bta-miR-10175-3p/RHOA networks was experimentally validated. Among these, bta-miR-541 and bta-miR-10175-3p were found to reduce lipid accumulation in mammary epithelial cells, while bta-miR-2309 exhibited the opposite effect. These findings provide new insights for improving milk quality through molecular breeding.

Milk fat composition and content are crucial indicators that influence the taste, flavor and nutritional value of milk, and they are regulated by various non-coding RNAs. To enhance milk quality, it is essential to investigate the regulatory mechanisms underlying milk fat metabolism. This study analyzed significantly up-regulated circ_0009058 and circ_0004021 and significantly down-regulated circ_0011934 and circ_0008056 in bovine mammary epithelial cells (BMECs) with high and low milk fat percentages (MFP). RNase R digestion assays and Sanger sequencing were conducted to confirm their presence and high expression in the mammary tissue of dairy cows. Subcellular localization indicated that these RNA may exhibit regulatory functions through the competing endogenous RNA (ceRNA) network. Consequently, a ceRNA regulatory network was constructed, identifying six hub target genes in the ceRNA network (VAV1, PTPN6, PIK3R1, RHOA, ERBB3, PIK3CG) using the CytoHubba and MCODE plugins in Cytoscape (version 3.9). Quantitative real-time PCR (RT-qPCR) and dual luciferase reporter gene assays validated the existence of the interactive regulatory networks circ_0004021/bta-miR-541/PTPN6, circ_0008056/bta-miR-2309/ERBB3, and bta-miR-10175-3p/RHOA in BMECs. Functional validation of miRNAs within the ceRNA network demonstrated that bta-miR-541 and bta-miR-10175-3p reduced triglyceride and cholesterol levels while inhibiting lipid droplet secretion in BMECs. Conversely, the function of bta-miR-2309 was found to be opposite to that of bta-miR-541 and bta-miR-10175-3p. This study identified several critical candidate ceRNA networks involved in fatty acid metabolism in dairy cows through data analysis and a series of experiments. These functional and mechanistic studies provide a theoretical foundation for improving milk quality.

Figures

Interaction regulation network of candidate DE-circRNAs and enrichment pathways of target genes in the network. (A,C,E,G) Interacting regulatory networks for circ_0009058, circ_0004021, circ_0011934 and circ_0008056, respectively. Red hexagonal nodes represent circRNAs, blue square nodes represent miRNAs, and purple circular nodes represent mRNAs. Solid lines indicate predicted interactions based on bioinformatic analyses. (B,D,F,H) KEGG pathway enriched for target genes from the interactions network of circ_0009058, circ_0004021, circ_0011934, and circ_0008056, respectively.

The crucial sub-networks and hub genes obtained by the cytoHubba and MCODE algorithms in Cytoscape. (A,B,F) Three crucial sub-networks of the target gene screened from the PPI network, pink color represents the hub target gene. (C–E,G) Candidate key ceRNA networks constituted by hub target genes, and schematic illustration of circRNA-miRNA-mRNA binding sites.

Full text

1. Introduction

Milk is widely recognized for its rich nutritional value, appealing to consumers who consistently seek health and wellness benefits. The three primary nutrients in milk are milk fat, milk protein, and lactose. Milk fat comprises approximately 400 fatty acids, with the content and composition of these fatty acids being the main factors influencing the taste and nutritional value of milk [1]. In milk, 99% of the milk fat exists in the form of lipid droplets, with triglycerides (TAG) constituting about 98% of these droplets [2,3]. Milk fat serves as an efficient source of essential fatty acids for the body and acts as a carrier for fat-soluble vitamins, including vitamins A, D, and E [4,5]. The metabolism of milk lipids involves a complex regulatory process that encompasses de novo synthesis, uptake, transport of fatty acids, TAG synthesis, and secretion of lipid droplets [6]. Significant progress has been made in understanding the regulation of milk fat metabolism in the bovine mammary gland [7]. However, the advent of molecular techniques has revealed a complexity that surpasses previous understandings [8].

Circular RNAs (circRNAs) represent a recently identified class of endogenous RNAs that are formed through the back-splicing of pre-mRNAs into closed-loop structures. These molecules exhibit resistance to RNase R degradation, which confers a high level of stability [9]. CircRNAs are evolutionarily conserved and demonstrate spatiotemporally restricted, tissue-specific expression patterns [10]. They are localized in both the nucleus and cytoplasm, where they play roles in regulating transcription, splicing, and post-transcriptional processes [11,12]. The primary mechanisms of circRNA biogenesis include lariat-driven cyclization, intron pairing-driven cyclization, and RNA-binding protein (RBP)-mediated cyclization [13,14]. Functionally, circRNAs can act as microRNA sponges [15], influence the stability of other RNAs [16], and bind to RBPs, thereby altering their activity or localization [17,18]. Additionally, some circRNAs contain internal ribosome entry site (IRES) elements and can be translated into peptides or proteins under specific conditions [19].

Recent studies have demonstrated that circRNAs play significant regulatory roles in mammary epithelial cells. For instance, circRNA-006258 regulates triacylglycerol (TAG) synthesis by adsorbing miR-574-5p in sheep mammary epithelial cells [20]. Additionally, circ11103 and circ01592 positively influence milk lipid metabolism in bovine mammary epithelial cells (BMECs) [21,22]. Furthermore, circRNA-001091 has been shown to interact with several microRNAs in sheep mammary tissues, with miR-29 specifically regulating the secretion of milk proteins, TAG, and lactose in BMECs [23,24]. Wang et al. [25] identified multiple circRNAs in the mammary glands of dairy cows that regulate the lipid metabolism gene CD36 through miRNA sponge effects. Collectively, these findings indicate that circRNAs play crucial regulatory roles in mammary tissue and mammary epithelial cells.

To investigate the role of circular RNAs (circRNAs) in milk fat metabolism, we conducted RNA sequencing to profile circRNA expression in bovine mammary epithelial cells (BMECs) isolated from dairy cows exhibiting high and low milk fat percentages (MFP). The most significantly differentially expressed circRNAs were selected for further analysis. By integrating bioinformatics predictions with experimental validation, we systematically identified and functionally characterized key circRNA-associated competitive endogenous RNA (ceRNA) networks. This study aims to elucidate experimentally verified regulatory axes involved in lipid metabolism within the bovine mammary gland, thereby providing mechanistic insights into the synthesis of milk fat.

2. Materials and Methods

2.1. Sample Collection

We selected 245 primiparous Holstein cows with identical feeding and management backgrounds (Table S1). From this group, four cows with chronically high MFP and four cows with chronically low MFP were chosen based on a comprehensive year-long dairy herd improvement (DHI) dataset from the farm. To ensure that these eight cows had somatic cell counts (SCC) below 100,000/mL and exhibited significant variations in MFP, milk samples were collected from them in the morning, noon, and evening, proportionate to their yield, and mixed in a ratio of 4:3:3 for DHI analysis (Table S2). Additionally, fresh milk samples were aseptically collected from each cow and transferred into 50 mL centrifuge tubes. These tubes were promptly sealed, sterilized, and stored at 37 °C for transport back to the laboratory for the isolation of mammary epithelial cells. The isolation, culture, and characterization of bovine mammary epithelial cells (BMECs) have been completed by our research group in previous studies [26].

2.2. CircRNA Library Construction and Sequencing

The circRNA sequencing method utilized in this study mirrors the approach detailed in the existing literature [26]. In summary, the procedure involves the extraction of total RNA, the removal of ribosomal RNA (rRNA), and the digestion of linear RNA using RNase R. Subsequently, strand-specific libraries are constructed and subjected to Illumina PE150 sequencing. After quality control and filtering, the raw sequencing data are aligned to the ARS-UCD1.2 bovine reference genome using Bowtie2 (version 2.5.3).

2.3. CircRNAs Data Analysis

The analysis of circRNA data encompassed circRNA identification, differential expression analysis, and functional enrichment analysis. CircRNA was detected and identified using the find_circ and CIRI2 tools [27,28]. The differential expression analysis of transcript count matrices from high and low MFP in BMECs was conducted using the DESeq2 package in the R programming language [29]. The DAVID database (https://davidbioinformatics.nih.gov/, accessed on 5 January 2026) was utilized to identify significant enrichment based on the KEGG pathway, with a significance threshold set at p < 0.05.

2.4. RNase R Digestion Assay and Sanger Sequencing

Total RNA was digested with RNase R (5 U/µg, 37 °C for 30 min), purified, reverse transcribed, and subjected to RT-qPCR to compare the abundance of circRNA and linear RNA. Candidate circRNAs underwent 40 cycles of PCR amplification. The products were recovered from gels and subjected to Sanger sequencing using an ABI3730 gene analyzer (Applied Biosystems), following the methodology previously described [30]. Peak plots were obtained and aligned using BioEdit (v.7.2) software to identify the head and tail splice sites of the circRNA. The primer sequences used in this study are listed in Table S3.

2.5. RT-qPCR

The tissues and cells utilized in the present experiments were pre-preserved by our research group. Total RNA was extracted from tissues and BMECs using Trizol reagent (Takara, Dalian, China). Complementary DNA (cDNA) was synthesized through reverse transcription using the PrimeScript RT kit (Takara, Dalian, China) in accordance with the manufacturer’s instructions. For circRNAs, reverse transcription was conducted using the random primers provided in the kit. For miRNAs, the stem-loop method was employed to design specific reverse transcription and quantitative polymerase chain reaction (qPCR) primers. mRNA primers were designed utilizing Primer Premier 5.0 (Premier Biosoft International, Palo Alto, CA, USA). U6 was used as an internal reference for miRNAs, while GAPDH served as an internal reference for circRNAs and mRNAs. Detection was carried out on a CFX96 assay system (Bio-Rad) using the SYBR@ Green Premix Pro Taq HS qPCR kit (Accurate Biotechnology (Hunan) Co., Ltd., Changsha, China). All qPCR reactions were performed with three technical replicates for each of the three independent biological replicates (n = 3). Detailed primer information for qPCR is provided in Table S4.

2.6. CeRNA Network Construction and Hub Gene Analysis

The binding sites for circRNA-miRNA and miRNA-mRNA were predicted using Targetscan (version 7.2) and miRanda (version 3.3a) software. The ceRNA pairs were then screened based on Context+ scores and free energy. Subsequently, the target genes were imported into STRING 11.0 to construct protein–protein interactions (PPIs). Hub genes and core subnetworks were identified using Cytoscape’s CytoHubba (version 3.9), which employs four algorithms: Degree, EPC, MCC, and MNC, along with the MCODE plugin. These methods are effective for identifying hub genes within the PPI network [31].

2.7. Cell Cultures and Transfections

BMECs were cultured in Dulbecco’s Modified Eagle Medium (DMEM/F-12) (BI, Beijing, China) supplemented with 20% fetal bovine serum (FBS) (CellMax, Beijing, China) and 1% penicillin/streptomycin (Solarbio, Beijing, China) at 37 °C in a 5% CO_2_ atmosphere. To induce lactation, the medium was further supplemented with 5 μg/mL insulin (MedChem Express, Monmouth Junction, NJ, USA), 5 μg/mL hydrocortisone (MedChem Express, Monmouth Junction, NJ, USA), and 30 ng/mL prolactin (ProSpec, Jerusalem, Israel). HEK-293T cells were cultured in high-glucose DMEM (Gibco, New York, NY, USA) containing 10% FBS and 1% penicillin/streptomycin under the same conditions as BMECs. When cell confluence reached 70–80%, miRNA transfection was conducted following the manufacturer’s instructions for Lipofectamine 3000 (Invitrogen, Carlsbad, CA, USA) [32].

2.8. Dual-Luciferase Reporter Assay

The miRNA mimics and their corresponding negative controls were procured from Tsingke Biotechnology Co., Ltd. (Beijing, China). Wild-type (WT) and mutant (MUT) sequences of circRNA and mRNA were cloned into the 3′ end of the firefly luciferase reporter gene (luc2) within the pmirGLO vector. Prior to conducting the Dual-Luciferase reporter gene assay, HEK-293T cells were maintained in 24-well plates, and the reporter vectors along with miRNA mimics (or their negative controls) were co-transfected using Lipofectamine 3000 when the cell confluency reached 70–80%. Each group consisted of three independent biological replicates (n = 3), with three technical replicates performed for each measurement. Following 48 h of transfection, the luciferase activity of both firefly and Renilla was measured using the Dual-Luciferase^®^ Reporter Assay System (Promega, Madison, WI, USA) in accordance with the manufacturer’s instructions [32].

2.9. TAG and Cholesterol Content Determination

The TAG and total cholesterol (TC) levels in transfected BMECs were assessed using the Cellular TAG and TC Enzymatic Assay Kit (Beijing Pulilai Gene Technology Co., Ltd., Beijing, China), following the manufacturer’s guidelines. Each experimental group, comprising either the mimic or negative control (NC), consisted of three independent biological replicates (n = 3). In brief, the cells were washed three times with phosphate-buffered saline (PBS) and subsequently collected by adding cell lysate, allowing them to stand for 10 min. The cells were then centrifuged at 2000× g for 5 min at room temperature, and the resulting supernatant was utilized for enzymatic assays.

2.10. Oil Red O and BODIPY Staining

Lipid droplets in BMECs were stained using the Oil Red O Staining Kit (Solarbio, Beijing, China). Staining was conducted on three independent biological replicates per group (n = 3), with representative images displayed for each replicate. Initially, the cells were washed with phosphate-buffered saline (PBS) and fixed using Oil Red O fixative for 30 min. Following this, the cells were washed with distilled water, and 60% isopropanol was added for 20 to 30 s. Subsequently, the cells were stained with freshly prepared Oil Red O working solution for 30 min and then thoroughly washed with distilled water. Nuclei were stained with hematoxylin for 1 to 2 min, after which the cells were washed again with distilled water. Finally, the cells were incubated with Oil Red O buffer and photographed.

Treated BMECs were stained with BODIPY dye (Invitrogen, CA, USA) following the manufacturer’s instructions. The BMECs were fixed in 4% paraformaldehyde for 30 min and subsequently washed three times with PBS. The cells were then incubated with BODIPY dye for 30 min. After incubation, the cells were washed with PBS three times, and the nuclei were stained with DAPI (G-CLONE, Beijing, China) for 10 min. Finally, anti-fluorescence attenuation sealer tablets (GCLONE, Beijing, China) were applied, and images were captured using a fluorescence microscope.

2.11. Statistical Analysis

Data analysis and visualization were conducted using Origin 8.0 software. The results are presented as mean ± standard error of the mean. An unpaired two-tailed Student’s t-test was employed to evaluate differences between the two groups. Statistical significance was established at a p-value threshold of <0.05, with *, **, and *** indicating p < 0.05, p < 0.01, and p < 0.001, respectively.

3. Results

3.1. Screening and Identification of Differentially Expressed circRNAs

A circRNA library of BMECs from cows exhibiting high and low MFP was constructed, followed by sequencing and circRNA identification. A total of 309 differentially expressed circRNAs (DE-circRNAs) were identified between the high and low MFP groups (Table S5). Among these, the two most significantly upregulated DE-circRNAs were circ_0009058 (log2FoldChange = 7.5216, p-value = 0.0263) and circ_0004021 (log2FoldChange = 6.6562, p-value = 2.11 × 10^−5^), and the two most significantly down-regulated DE-circRNAs were circ_0011934 (log2FoldChange = −6.5955, p-value = 0.0077) and circ_0008056 (log2FoldChange = −5.9305, p-value = 0.0002). Notably, circ_0009058 was clipped from exon4, exon5 and exon8 of CPEB4 (Figure 1A), while circ_0004021 was back-spliced from exon9 of PABPC1 (Figure 1B). In addition, circ_0011934 was back-spliced from exon20-24 of SORCS1 (Figure 1C) and circ_0008056 is back-spliced. from exon1, exon2, exon4, exon5, and exon6 of ST3GAL6 (Figure 1D). The circular structure of these four circRNAs was confirmed through RNase R digestion experiments and Sanger sequencing. The four candidate circRNAs (circ_0009058, circ_0004021, circ_0011934, circ_0008056) exhibited no significant reduction in expression following RNase R digestion when compared to the untreated control. In contrast, the corresponding linear RNAs demonstrated a decrease in abundance exceeding 90% after digestion. Subsequently, Sanger sequencing was performed on the amplification products of these four circRNAs, and the sequencing results were analyzed using BioEdit software to compare with reference sequences, thereby confirming the head-to-tail linkage structure of the circRNAs. Consequently, their circRNA structures were further validated.

3.2. Subcellular Localization and Tissue Expression of Candidate DE-circRNAs

The potential functions of circ_0009058, circ_0004021, circ_0011934, and circ_0008056 in the mammary gland were investigated. RT-qPCR was conducted to assess the expression levels of these circRNAs in various tissues of dairy cows. Notably, circ_0009058, circ_0004021, circ_0011934, and circ_0008056 exhibited relatively high expression levels in the mammary tissues of dairy cows (Figure 2A–D). Subsequently, RT-qPCR assays were employed to determine the specific localization of these circRNAs, revealing that circ_0004021, circ_0011934, and circ_0008056 were predominantly localized in the cytoplasm, whereas circ_0009058 was expressed in both the nucleus and cytoplasm (Figure 2E). Cytoplasmic circRNAs have the potential to function as molecular sponges for miRNAs by competitively binding to them. This interaction modulates the inhibitory effects of miRNAs on their target genes, thereby influencing the expression of associated genes. Such mechanisms are essential to the regulatory processes involved in mammary gland development and lactation. Consequently, these four circRNAs may play significant regulatory roles in the development and lactation of the mammary glands in dairy cows through the ceRNA network.

3.3. CeRNA Network Construction for Candidate DE-circRNAs

Considering the potential role of candidate DE-circRNAs as ceRNAs, we utilized TargetScan (version 7.2) and miRanda (version 3.3a) software to predict their target miRNAs and analyze the binding sites, thereby screening for high-scoring targeted miRNAs associated with the candidate DE-circRNAs (Table S6). Following this, we predicted the target genes of the high-scoring miRNAs using the same methodology (Table S7). We constructed the ceRNA regulatory networks for the four DE-circRNAs (Figure 3A,C,E,G). Subsequently, we conducted KEGG pathway enrichment analysis for the target genes within the ceRNA networks. The findings indicated that the target genes within the circ_0009058 regulatory network were predominantly enriched in pathways related to lipid metabolism, including AMPK signaling, metabolic pathways, oxytocin signaling, inflammatory mediator regulation of TRP channels, and insulin secretion (see Figure 3B). Conversely, target genes in the circ_0004021 regulatory network were primarily enriched in pathways associated with ether lipid metabolism, glycerophospholipid metabolism, cholesterol metabolism, fat digestion and absorption, and the PI3K-AKT signaling pathway (Figure 3D). The target genes within the circ_0011934 regulatory network were enriched in the phospholipase D signaling pathway, insulin resistance and secretion, sphingolipid signaling pathway, prolactin signaling pathway, PI3K-AKT signaling pathway, regulation of lipolysis in adipocytes, and the mTOR signaling pathway (Figure 3F). Lastly, target genes in the circ_0008056 regulatory network were enriched in the phospholipase D signaling pathway, glucagon signaling pathway, sphingolipid metabolism, AMPK signaling pathway, and glycerolipid metabolism pathway (Figure 3H). These results further suggest that candidate DE-circRNAs may play a significant regulatory role in milk fat metabolism through the ceRNA network.

3.4. Screening for Candidate Key ceRNA Networks

To explore potential key ceRNA networks, we adopted bioinformatics methods to obtain hub target genes, as hub target genes typically regulate other genes within related pathways. Three sub-networks were screened from the PPI network of target genes using the MCODE plugin in Cytoscape (Figure 4A,B,F). Subsequently, the CytoHubba plugin in Cytoscape was utilized to obtain the hub genes. Ultimately, six candidate hub genes from three sub-networks were identified, including VAV1, PTPN6, PIK3R1, RHOA, ERBB3, and PIK3CG. These hub target genes are proposed to constitute putative ceRNA networks, which we hypothesize represent key candidate regulatory modules in bovine milk fat metabolism. Our in silico analyses suggest the following potential regulatory axes: circ_0004021 is inferred to regulate PTPN6 and VAV1 genes through bta-miR-541 and bta-miR-2330-5p, respectively (Figure 4C); circ_0011934 is putatively linked to PIK3R1 regulation through bta-miR-2449 (Figure 4D); circ_0009058 is predicted to target RHOA through bta-miR-10175-3p (Figure 4E); and bta-miR-2309, potentially sponged by circ_0008056, is computationally associated with the regulation of ERBB3 and PIK3CG (Figure 4G). The base sequence of the 3’UTR seed region was found to be complementarily paired by predicting their targeted binding sites each other. Subsequently, its authenticity will be verified through experiments.

Interaction regulation network of candidate DE-circRNAs and enrichment pathways of target genes in the network. (A,C,E,G) Interacting regulatory networks for circ_0009058, circ_0004021, circ_0011934 and circ_0008056, respectively. Red hexagonal nodes represent circRNAs, blue square nodes represent miRNAs, and purple circular nodes represent mRNAs. Solid lines indicate predicted interactions based on bioinformatic analyses. (B,D,F,H) KEGG pathway enriched for target genes from the interactions network of circ_0009058, circ_0004021, circ_0011934, and circ_0008056, respectively.

The crucial sub-networks and hub genes obtained by the cytoHubba and MCODE algorithms in Cytoscape. (A,B,F) Three crucial sub-networks of the target gene screened from the PPI network, pink color represents the hub target gene. (C–E,G) Candidate key ceRNA networks constituted by hub target genes, and schematic illustration of circRNA-miRNA-mRNA binding sites.

3.5. RT-qPCR Identifies Interplay Regulatory Relationships in Candidate ceRNA Networks

Based on the ceRNA theory, to further determine the interactive regulatory relationships among the molecules in the candidate key ceRNA networks, The mimics of bta-miR-541, bta-miR-2330-5p, bta-miR-10175-3p, bta-miR-2309, and bta-miR-2449 were transfected into the BMECs, respectively, and the transfection efficiency of these miRNAs was detected by RT-qPCR after 48 h. The results showed that the expression levels of miRNAs were significantly higher (p < 0.01) in the mimic group compared to the NC group (Figure 5A–E). The effectiveness of transfection and expression of miRNAs in cells was demonstrated. The effect of miRNAs in the candidate key ceRNA network on the expression levels of circRNAs and their target genes was examined by RT-qPCR. The results revealed that bta-miR-541 significantly decreased the expression levels of circ_0004021 and its target gene PTPN6 (p < 0.01) (Figure 5F,G). Notably, bta-miR-2330-5p negatively regulated circ_0004021 but did not inhibit the expression of the target gene VAV1. In addition, bta-miR-10175-3p significantly inhibited the expression levels of circ_0009058 and target gene RHOA (Figure 5J,K). Bta-miR-2309 negatively regulated target-bound circ_0008056 and ERBB3 (Figure 5L,M). However, the expression of circ_0011934 was not detected in bta-miR-2449 mimic-treated cells and its target gene PIK3R1 was not inhibited (Figure 5N). Collectively, these experimental results revealed that the interaction networks circ_0004021/bta-miR-541/PTPN6, circ_0009058/bta-miR-10175-3p/RHOA, and circ_0008056/bta-miR-2309/ERBB3 were key ceRNAs in milk fat regulation in dairy cows.

3.6. Dual Luciferase Reporter Assay Validates Interactive Regulatory Relationships of Key Cerna Networks

Dual luciferase reporter gene assay revealed that overexpression of bta-miR-541 significantly attenuated the relative luciferase activity of circ_0004021 and PTPN6 WT groups. However, the relative luciferase activities of MUT circ_0004021 and PTPN6 were not significantly altered (Figure 6A,B). Overexpression of bta-miR-2309 significantly decreased the relative luciferase activity of WT group circ_0008056 and ERBB3, which was effectively abrogated by mutating the 3’UTR sequences of circ_0008056 and ERBB3 (Figure 6C,D). Dual-luciferase activity assay was performed to detect circ_0009058/bta-miR-10175-3p/RHOA interactions, indicating that bta-miR-10175-3p mimic had no effect on luciferase activity when co-transfected with circ_0009058-WT (Figure 6E). Moreover, luciferase activity was significantly decreased when co-transfected with RHOA-WT (Figure 6F). These results confirm the existence of the circ_0004021/bta-miR-541/PTPN6, circ_0008056/bta-miR-2309/ERBB3, and bta-miR-10175-3p/RHOA networks in BMECs.

3.7. Impact on Lactation of miRNAs in the Key ceRNA Network

Furthermore, the regulatory functions of bta-miR-541, bta-miR-2309, and bta-miR-10175-3p in milk lipid metabolism in BMECs were explored. Notably, bta-miR-541 and bta-miR-10175-3p significantly decreased TAG and TC contents in BMECs (Figure 7A), whereas bta-miR-2309 promoted TAG and TC synthesis (Figure 7B). Milk fat exists in the form of lipid droplets, and the effects of miRNAs for lipid droplet secretion were observed using oil red O and BODIPY staining. The results revealed that bta-miR-541 and bta-miR-10175-3p inhibited the secretion of lipid droplets, while, bta-miR-2309 promoted the secretion of lipid droplets (Figure 7C,D). These results fully demonstrate that miRNAs in the circ_0004021/bta-miR-541/PTPN6, circ_0008056/bta-miR-2309/ERBB3, and bta-miR-10175-3p/RHOA interworking networks play potential regulatory roles in dairy cow milk fat metabolism.

3.8. Regulation for the Fatty Acid Synthesis Pathway of miRNA in the Key ceRNA Network

The synthesis pathway of fatty acid (FA) involves ab initio synthesis of FAs, uptake of FAs from the blood, intracellular transport and desaturation of FAs, synthesis of TAG, and secretion of final lipid droplets. Hence, the effects of bta-miR-541, bta-miR-2309, and bta-miR-10175-3p on the FA synthesis pathway were investigated. The expression levels of the relevant genes in the FA synthesis pathway were analyzed using RT-qPCR. The results indicated that bta-miR-541 and bta-miR-10175-3p significantly decreased FA concentration by significantly inhibiting the expression levels of the FA de novo synthesis marker gene ACACA (Figure 8B), the FA uptake marker gene ACSL1 (Figure 8C), the FA transport and desaturation marker genes FABP3 and SCD (Figure 8A), the TAG synthesis marker genes DGAT1 and LPIN1 (Figure 8D), and expression levels of the lipid droplet secretion marker genes XDH and PLIN2. All the marker genes related to the FA synthesis pathway were significantly up-regulated by bta-miR-2309 (Figure 8A–E). Therefore, bta-miR-2309 ultimately increased FA secretion by promoting the entire pathway of FA synthesis.

4. Discussion

With the development of high-throughput sequencing technology, researchers have found that circRNAs are ubiquitous in a variety of species and have important physiological regulatory functions [33]. Several previous studies have reported that circRNAs play critical regulatory roles in the mammary glands of dairy cows [21,34,35]. These studies have fully demonstrated the existence of a large number of regulatory circRNAs in the mammary glands of dairy cows, which are worthy of further exploration and research.

In this study, we selected the two circRNAs with the most significant upregulation (circ_0009058 and circ_0004021) and the two with the most significant downregulation (circ_0011934 and circ_0008056) as the primary candidate genes influencing milk fat metabolism in dairy cows. This selection was based not only on their extreme differential expression but also on their high expression in mammary tissue, cytoplasmic localization, and the association of their parental genes with lipid metabolism pathways.

Growing evidence suggests that circRNAs can regulate the expression of their parental genes. [11,36]. Additionally, circRNAs can also act as sponges for miRNAs targeting their parental genes, thereby maintaining parental gene expression [15,27,37].Thus, in this study, we predicted the potential roles of circRNAs by the function of the parent genes of circRNAs. Specifically, blockade of the expression of CPEB4, the parent gene of circ_0009058, prevented diet-induced weight gain and reduced adipose tissue enlargement, lipid accumulation, and the pro-inflammatory capacity of macrophages [38]. This suggests that circ_0009058 may influence body weight and fat metabolism by regulating CPEB4 expression. Meanwhile, the parental gene PABPC1 of circ_0004021 is involved in the regulation of lipid droplet secretion and intracellular lipid transport [39,40]. This suggests that circ_0004021 may play an important role in lipid metabolism. In addition, the parental gene SORCS1 of circ_0011934 was involved in the processes of energy homeostasis, insulin secretion, and fat deposition in the organism [41]. For circ_0008056, it is not clear whether its parent gene ST3GAL6 directly regulates lipid metabolism. However, ST3GAL6 mediates the regulation of miR-26a through the AKT/mTOR signaling pathway [42]. Moreover, miR-26a is considered a critical regulator of lipid metabolism [43]. In summary, the present study reveals the potential mechanism of these circRNAs in the regulation of lipid metabolism through the function of their parent genes.

Subcellular localization and tissue expression revealed that circ_0009058, circ_0004021, circ_0011934, and circ_0008056 are predominantly present in the cytoplasm and highly expressed in mammary tissues, suggesting that they can competitively bind to miRNAs to regulate mammary gland development and lactation. Therefore, their ceRNA regulatory network was constructed. Considering the numerous target genes predicted in the ceRNA network, CytoHubba and MCODE plug-ins in Cytoscape were utilized to identify six hub target genes (VAV1, PTPN6, PIK3R1, RHOA, ERBB3, PIK3CG) from the target gene’s PPI network. RT-qPCR and dual luciferase reporter gene assays were performed to examine the interactive regulatory roles of molecules in candidate ceRNAs composed of hub target genes. Finally, circ_0004021/bta-miR-541/PTPN6, circ_0008056/bta-miR-2309/ERBB3, and bta-miR-10175-3p/RHOA were screened as the key regulatory networks for exploring the mechanism of milk fat regulation in dairy cows.

In the circ_0004021/bta-miR-541/PTPN6 regulatory network, PTPN6, also known as SHP-1, is considered a major regulator of inflammation and glucose homeostasis [44,45]. PTPN6 was found to be localized to lipid rafts and interacted with EGFR [46,47]. Furthermore, the activation of the EGFR signaling pathway can lead to up-regulation of lipid metabolism [48]. This suggests that PTPN6 may indirectly regulate lipid metabolism by affecting EGFR signaling. PTPN6 also interacts with the lipid metabolism gene CD36, and the transcriptional regulator of HDL cholesterol metabolism TCF1 [49,50,51], further suggesting its critical role in lipid metabolism. In addition, bta-miR-541 in this network effectively inhibited TAG, TC synthesis, and lipid droplet secretion as a sponge that upregulates circ_0004021. Therefore, the signaling axis circ_0004021/bta-miR-541/PTPN6 may have potential regulatory roles in lipid metabolism in BMECs.

In the present study, we also identified the regulatory axis bta-miR-10175-3p/RHOA, which has a potential role in the regulation of milk fat metabolism in dairy cows. RHOA, a small GTPase of the RHO family, plays a key role in lipid metabolism. It was found that activation of RHOA significantly increased EGFR expression [52], whereas decreased RHOA activity upregulated PPARγ mRNA expression and induced lipogenic differentiation [53]. In addition, RHOA regulated the activity and expression of the lipid transporter protein ABCA1 [54], and its downregulation led to a decrease in the expression of HMGCR and SREBF2 mRNA, key regulators of cellular cholesterol metabolism, as well as an increase in intracellular cholesteryl ester content and apolipoprotein B (APOB) concentration [55]. Notably, bta-miR-10175-3p significantly inhibited milk fat synthesis. This finding suggests that bta-miR-10175-3p/RHOA regulatory axis has an important regulatory role in milk fat metabolism in dairy cows.

In the circ_0008056/bta-miR-2309/ERBB3 regulatory network, ERBB3, as a member of the EGFR family, interacts with the PIK3 regulatory subunit and transmits signals to the mTOR pathway to regulate insulin secretion and intracellular glucose metabolism [56]. In addition, ERBB3 is a potent activator of the PI3K/AKT pathway, forming a complex with ERBB2 to activate the downstream PI3K/AKT/mTOR pathway [57,58]. While the PI3K/AKT pathway plays an important role in the initiation of lactation, its downstream mTOR signaling plays a key role in the regulation of protein and lipid synthesis and cell growth [59,60]. In this network, bta-miR-2309 acts as a competitive binding factor for the down-regulation of circ_0008056 and was able to significantly up-regulates TAG and TC content as well as lipid droplet secretion in BMECs. In summary, this study concluded that the circ_0008056/bta-miR-2309/ERBB3 regulatory network may be involved in regulating the PI3K/AKT/mTOR pathway in the mammary gland, which affects lactation and lipid metabolism in dairy cows. The ceRNA networks we identified provide novel insight into understanding the molecular mechanisms of milk fat metabolism in dairy cows.

While this study provides compelling evidence for the involvement of specific ceRNA networks in milk fat metabolism, certain limitations should be acknowledged. Firstly, the RNA-seq analysis for screening differentially expressed circRNAs was conducted using a relatively small cohort (n = 4 per group), which, although sufficient for an exploratory study, may affect the generalizability of the circRNA profile. Future studies with larger and independent populations are warranted to confirm the robustness of these biomarkers. Secondly, the regulatory roles of the identified circRNAs were inferred through network analysis and miRNA manipulation. Direct functional validation via circRNA knockdown or overexpression in BMECs would further solidify their causal roles within the proposed ceRNA axes. Addressing these points in future work will deepen our understanding of the circRNA-mediated regulatory landscape in bovine lactation.

5. Conclusions

This study identified circ_0009058, circ_0004021, circ_0011934, and circ_0008056 in bovine mammary epithelial cells as regulators of milk fat metabolism. It established regulatory networks involving circ_0004021/bta-miR-541/PTPN6, circ_0008056/bta-miR-2309/ERBB3, and bta-miR-10175-3p/RHOA, which are integral to the regulation of milk fat metabolism (Figure 9). These findings offer novel insights and provide a theoretical foundation for the regulation of fatty acid synthesis and metabolism in milk.