Skip to main content

Mining key circRNA-associated-ceRNA networks for milk fat metabolism in cows with varying milk fat percentages

Abstract

Background

Cow milk fat is an essential indicator for evaluating and measuring milk quality and cow performance. Growing research has identified the molecular functions of circular RNAs (circRNAs) necessary for mammary gland development and lactation in mammals.

Method

The present study analyzed circRNA expression profiling data in mammary epithelial cells (MECs) from cows with highly variable milk fat percentage (MFP) using differential expression analysis and weighted gene co-expression network analysis (WGCNA).

Results

A total of 309 differentially expressed circRNAs (DE-circRNAs) were identified in the high and low MFP groups. WGCNA analysis revealed that the pink module was significantly associated with MFP (r = − 0.85, P = 0.007). Parental genes of circRNAs in this module were enriched mainly in lipid metabolism-related signaling pathways, such as focal adhesion, ECM-receptor interaction, adherens junction and AMPK. Finally, six DE-circRNAs were screened from the pink module: circ_0010571, circ_0007797, circ_0002746, circ_0003052, circ_0004319, and circ_0012840. Among them, circ_0002746, circ_0003052, circ_0004319, and circ_0012840 had circular structures and were highly expressed in mammary tissues. Subcellular localization revealed that these four DE-circRNAs may play a regulatory role in the mammary glands of dairy cows, mainly as competitive endogenous RNAs (ceRNAs). Seven hub target genes (GNB1, GNG2, PLCB1, PLCG1, ATP6V0C, NDUFS4, and PIGH) were obtained by constructing the regulatory network of their ceRNAs and then analyzed by CytoHubba and MCODE plugins in Cytoscape. Functional enrichment analysis revealed that these genes are crucial and most probable ceRNA regulators in milk fat metabolism.

Conclusions

Our study identified several vital circRNAs and ceRNAs affecting milk fat synthesis, providing new research ideas and a theoretical basis for cow lactation, milk quality, and breed improvement.

Peer Review reports

Introduction

Milk fat an important nutrient and key evaluation indicator for milk. Conjugated linoleic acid—which is abundant in milk fat—is essential for cholesterol downregulation [1] and low-density lipoprotein levels [2] in humans and for the defense against atherosclerosis [3]. Milk fat contains essential minerals and fat-soluble vitamins for humans [4,5,6]. Cheese—a further processed milk product—holds a significant position in the global premium market for dairy products. However, milk fat has an important influence on cheese flavor, not only because the fatty acids used to synthesize cheese are flavor substances in their own right [7], but also as precursors to flavor substances, such as methyl ketones, secondary alcohols, lactones, and esters, in cheese [8]. Therefore, the safe and effective increase of milk fat content in milk is one of the necessary tools to strengthen the dairy industry’s core competitiveness in the global market.

The advancement of sequencing technology and bioinformatics algorithms has enabled the identification of a class of potentially functional circular RNAs (circRNAs) [9]. These circRNAs are linked together through the 5’ and 3’ ends of the parent gene to form a circular structure [10]. Studies have demonstrated that circRNAs have important molecular functions for mammary gland development and lactation in mammals. circRNA-006258 is closely related to mammary epithelial cell (MEC) growth and milk synthesis in goats [11]. circ_015343 reduces milk production and milk fat synthesis and inhibits MEC growth in sheep [12]. circ01592 and circ09863 increase the levels of triglycerides, cholesterol and unsaturated fatty acids content in MECs of dairy cows [1314]. circRNA8220 promotes the proliferation and synthesis of β-casein and triglycerides in the MECs of goat [15]. Considering the importance of circRNAs to MECs, we discovered more circRNAs affecting MEC growth and lactation using differential expression analysis and weighted gene co-expression network analysis (WGCNA). WGCNA is an effective systems biology method for analyzing RNA-seq data, including mRNA [1617], miRNA [1819], lncRNA [2021], and circRNA [2223]. WGCNA canlocate core genes faster using gene connectivity information. Weak effector genes can also be mined, elucidating the biological mechanisms underlying traits. The differential expression analysis and WGCNA methods complement each other and help in the rapid and comprehensive identification of the DE-circRNAs that regulate MFP.

Given the potential of circRNAs to indirectly regulate gene expression in MECs, it is necessary to identify and characterize circRNAs in MECs of cows with different milk fat percentage (MFP), since this circRNAs may be involved in epigenetic and genetic regulation of mammary function. In this study, RNA-seq was used to obtain the expression levels of circRNA in MECs of lactating cows with significantly differential MFP. DE-circRNAs mediating milk lipid metabolism were mined by using differential expression analysis and WGCNA, and we explore the effects of circRNA-mediated regulatory networks on mammary gland development and lactation in dairy cows. This study provides a new research idea and a theoretical basis for future studies on the mechanism of circRNA-regulated milk fat metabolism in dairy cows.

Results

Combined differential expression analysis and WGCNA screening for candidate DE-circRNAs

We constructed a circRNA library of MECs from the high- and low-MFP groups, then sequenced and identified circRNAs. A total of 309 DE-circRNAs were found at a significantly higher level in the high-MFP group than in the low-MFP group (Fig. 1A; Table S1). Subsequently, WGCNA was constructed using transcriptome sequencing data. A total of 18 co-expression modules were obtained after merging modules with a similarity more significant than 75% (Fig. 1B). The number of circRNAs contained in each module ranged from 20 (grey) to 198 (turquoise), among which there were eight modules with over 100 circRNAs, including the turquoise, blue, brown, and yellow modules (Fig. 1C). Module-trait correlation analysis demonstrated that multiple modules were associated with MFP (Fig. 1D), among which the pink module was significantly negatively correlated with MFP (r = − 0.85, P = 0.007) (Fig. 1D-E). The pink module contained 101 circRNAs (Table S2). Functional annotation of 101 circRNAs revealed that the significantly enriched GO terms included regulation of the triglyceride metabolic process, regulation of the fatty acid metabolic process, positive regulation of MEC proliferation, and the relation of insulin and phosphatidylinositol. The significantly enriched Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways included focal adhesion, ECM-receptor interaction, and adherens junction. These findings suggests that the pink module might contain critical circRNAs that regulate lipid metabolism (Fig. 1F). Subsequently, the gene significance (GS) and module membership (MM) values were calculated to screen the key circRNAs (Fig. 1G), resulting in 11 key circRNAs (|GS| ≥ 0.90 and|MM| ≥ 0.60). Among them, circ_0010571, circ_0007797, circ_0002746, circ_0003052, circ_0004319, and circ_ 0012840 belonged to DE-circRNAs (Table S3); therefore, these circRNAs were further validated as candidates.

Fig. 1
figure 1

Differential expression analysis and WGCNA of circRNAs in MECs of cows with different MFP. (A) DE-circRNAs in the HMF and LMF groups. (B) Clustered tree diagram of circRNAs in MECs of dairy cows, each color represents a module. (C) The number of circRNAs clustered in different modules. (D) Heatmap of correlation between modules and MFP (each module contains correlation coefficient and corresponding P-value). (E) Module significance values ​​for co-expression modules associated with MFP (module significance values ​​represent a summary of circRNA significance for all circRNAs in each module, with different column colors representing different modules). (F) Functional annotation of circRNAs in the pink module. (G) A scatterplot of GS for MFP vs. MM in the pink module respectively (a dot represents a circRNA in the pink module)

Validation of the circular structure of candidate circRNAs

To demonstrate the circular structure of the six candidate circRNAs, we used divergent primers and convergent primers to detect the resistance of circRNAs and linear RNAs to linear RNase R. The results showed that circ_0002746, circ_0004319, and circ_ 0012840 resisted linear RNase R digestion (Fig. 2A-D), but circ_0010571, circ_0007797 could not resist linear RNase R digestion (Fig. S1). Subsequently, the head-to-tail splice sites of circ_0002746, circ_0003052, circ_0004319, and circ_ 0012840 were also confirmed by through Sanger sequencing (Fig. 2E-H).

Fig. 2
figure 2

Circular structure identification of circRNAs. (A-D) We identified by RT-qPCR using divergent and convergent primers and found that circ_0002746, circ_0003052, circ_0004319 and circ_0012840, but not linear_0002746, linear_0003052, linear_0004319, and linear_0012840, were resistant to RNase R digestion. (E-H) Backsplice sites of circ_0002746, circ_0003052, circ_0004319 and circ_0012840 were confirmed by Sanger sequencing

Tissue expression and subcellular localization of candidate DE-circRNAs

To investigate the potential functions of circ_0002746, circ_0003052, circ_0004319 and circ_0012840 in the mammary glands of dairy cows, we used RT-qPCR to examine the expression levels of these circRNAs in various tissues of dairy cows. The results demonstrated that circ_0002746, circ_0003052, and circ_0004319 were highly expressed in mammary tissues compared with other tissues (Fig. 3A-C). In mammary tissue, the expression abundance of circ_ 0012840 was second only to that of the small intestine (Fig. 3D). To identify the specific locations where these circRNAs functioned, we isolated and examined the nucleus and cytoplasm of MECs using RT-qPCR. circ_0002746, circ_0003052, circ_0004319, and circ_0012840 were expressed at significantly higher levels in the cytoplasm than in the nucleus (Fig. 3E), suggesting that these four circRNAs may have potential regulatory functions on mammary gland development and lactation in dairy cows through the ceRNA network.

Fig. 3
figure 3

(A-D) The relative expression levels of circ_0002746, circ_0003052, circ_0004319 and circ_0012840 in different tissues of dairy cows, respectively. (E) The expression abundances of circ_0002746, circ_0003052, circ_0004319 and circ_0012840 in the cytoplasm and nucleus, with U6 and GAPDH as internal controls for nuclear and cytoplasm, respectively

Construction and screening of ceRNA networks for candidate circRNAs

We used the Targetscan (v7.2) [24] and miRanda (v3.3a) [25] software to predict circRNA/miRNA and miRNA/mRNA interactions and to screen the top 5 miRNAs that bind to each candidate circRNA (Table S4). We screened 372 target genes of the top 5 miRNAs bound by circRNAs according to the context+ ≤ − 0.20 and free energy ≤ − 20 of the miRNA/mRNA interaction relationship (Table S5). Although many target genes existed in the ceRNA network, we further screened the hub target genes. Hub genes are vital in biological processes, and their regulation often influences other genes within related pathways. We first identified four crucial sub-networks from the PPI network of target genes using the MCODE plugin in Cytoscape (Fig. 4A, C, E, H), and then, we performed gene ontologies (GO) and KEGG analysis to clarify the role of target genes in these subnetworks. Subnetwork 1’s main functions were focused on metabolic processes such as V-type ATPase for proton transport and ATP hydrolysis and synthesis, and the enriched pathways was mainly oxidative phosphorylation (Fig. 4B). Subnetwork 2 was enriched in G-protein coupled receptor signaling pathway, lipid catabolic process, phosphatidylinositol phospholipase C activity, and other GO terms. KEGG analysis was involved in lipid metabolism-related signaling pathways, such as calcium signaling pathway, insulin secretion, oxytocin signaling pathway, and phosphatidylinositol signaling system (Fig. 4D). Subnetwork 3’s functions were focused mainly on the glycosylphosphatidylinositol (GPI)-anchor biosynthetic process of biological process (BP). In cellular component (CC) and molecular function (MF), Subnetwork 3 was mainly involved in the synthesis of phosphatidylinositol-related proteins. The main enriched KEGG pathway was GPI-anchor biosynthesis (Fig. 4F). The target genes in subnetwork 4 performed their functions in the mitochondrial respiratory chain complex IV, and the main enriched pathway was also oxidative phosphorylation (Fig. 4G). We obtained seven hub genes from these four sub-networks: GNB1, GNG2, PLCB1, PLCG1, ATP6V0C, NDUFS4, and PIGH (Fig. 4A, C, E, H). These seven hub target genes compose key ceRNA network for exploring the mechanism of milk fat regulation in cows (Fig. 5).

Fig. 4
figure 4

The crucial sub-networks and hub genes obtained by the cytoHubba and MCODE algorithms in Cytoscape. (A) (C) (E) (H) Four crucial sub-networks of the target gene screened from the PPI network. (B) (D) (F) (G) Enrichment analysis results of sub-networks (A), (C), (E), and (H) respectively

Fig. 5
figure 5

The candidate key ceRNA network regulating milk fat metabolism

Discussion

WGCNA can combine gene expression with phenotypic data, making it more suitable for analyzing complex data. There are several applications for diverse omics data (e.g. transcriptomics, proteomics, and metabolomics) and various organisms (animal, plant, and microbial) [2627]. WGCNA provides a solution to the multiple testing problems by reducing the size of large networks into a small number of hub nodes, allowing comparison of external traits with a limited number of variables. And can cluster genes with the same function or pathway to form functional modules [28]. Recognizing the core modules helps to annotate the results of systems biology scale experiments, thus adding valuable biological information. Therefore, we used WGCNA to construct a co-expression network of circRNAs from high- and low-MFP Holstein cows and found that the pink module was significantly and negatively correlated with MFP. Focal adhesion was the most significantly enriched pathway for circRNAs in this module, which connects downstream to the mitogen-activated protein kinase (MAPK) signaling pathway, thereby affecting lipid metabolism [29]. Among other pathways, the extracellular matrix (ECM) can regulate the proliferation, apoptosis, and polarity of MECs [30,31,32]. The AMPK signaling pathway acts as an energy sensor that regulates metabolism in the body and cells, including lipid metabolism [33]. In epithelial cells, tight junctions are essential for cell adhesion and prevent the lateral diffusion of lipids and proteins. Cholesterol and long-chain fatty acids are abundant in its plasma membrane. The adherens junction increases cholesterol levels in the plasma membrane to facilitate tight junction formation [3435]. The enriched pathways suggest these circRNAs may have potential regulatory mechanisms for mammary gland development and milk fat metabolism. Finally, six candidate DE-circRNAs (circ_0010571, circ_0007797, circ_0002746, circ_0003052, circ_0004319, and circ_0012840) were screened from the pink module by combining differential expression analysis and WGCNA. The circular validation revealed that circ_0002746, circ_0003052, circ_0004319, and circ_0012840 belonged to circRNAs, which were the primary focus of this study.

Subcellular localization and tissue expression revealed that circ_0002746, circ_0003052, circ_0004319, and circ_0012840 were predominantly present in the cytoplasm and highly expressed in mammary tissue. These results provide strong evidence that the four DE-circRNAs to regulate mammary gland development and lactation by competitively binding miRNAs. In circRNA/miRNA and miRNA/mRNA interactions regulation, the target genes of circ_0002746, PPARD, ELOVL2, LSS, ACAA2 and PPARGC1A, were associated with the regulation of lipid metabolism [3637]. circ_0003052 as sponge of miR-2454-3p to regulate SLC27A6 with high expression in adipose tissue, and inhibiting SLC27A6 expression may significantly affect lipid metabolism pathways, including lipid biosynthesis, transport, and β-oxidation in mammary cells [3839]. Among the miRNAs that interacted with circ_0004319, miR-11,999 exhibited the highest number of target genes, among which ELOVL7, ACADS, and APC were lipid metabolism-related candidate genes [4041]. circ_0012840 and circ_0003052 interact with miR-7864 to regulate the expression of PLA2G2E, which promotes lipid accumulation in adipose tissue and liver [42]. circ_0012840 interacts with miR-214 and miR-761 to regulate VPS4A, an important regulator of endosomal cholesterol transport [43].

We conducted a PPI network analysis of target genes in the ceRNA regulatory network of four candidates DE-circRNAs and screened seven hub target genes from the PPI network. GNB1 (guanine nucleotide-binding protein (G protein), beta polypeptide 1) encodes the Gβ subunit of a heterotrimeric G protein complex that includes Gα and Gγ subunits. This complex function can transduce multiple intracellular signaling cascades [44]. This gene promotes lipolysis in adipose tissue and, when expressed, increases blood glycerol levels [45]. GNG2 (G protein subunit gamma 2) expression is positively correlated with adipocyte size [46], and its upregulation can directly activate PI3K IB, thereby activating the PI3K-Akt pathway [47]. This pathway is also regulated by the Gβγ subunits of the trimer G protein complex formed by GNB1 and GNG2 [4849]. Gβ1γ2 produces phosphoinositol by stimulating phospholipase Cβ, activating MAPK and Akt [50]. This evidence suggests that GNB1 and GNG2 are involved in lipid metabolic pathways. PLCB1 and PLCG1 belong to the phospholipase C (PLC) gene. PLC protein—a key enzyme for metabolizing inositol lipids—plays a key role in multiple transmembrane signal transduction pathways that regulate various cellular processes, including cell proliferation and mobility [51]. PLCB1 (phospholipase Cβ1) is involved in adipocyte differentiation [52]. PLCG1 (phospholipase Cγ1) has two major lobes: one contains the active site that modifies lipids, and the other sits on top of the active site to prevent lipids from reaching it [53]. PLCB1 and PLCG1 are important candidate genes for fat deposition [5455]. ATP6V0C and the genes in Subnetwork 2 are part of the vacuolar ATPase (V-ATPase), which is crucial in stimulating mitochondrial gluconeogenesis and insulin secretion in the body [56]. The tight binding between the lipid phosphatidylinositol 3,5-bisphosphate (PI(3,5)P2) and the membrane of V-ATPase can activate V-ATPase activity and proton pump [5758]. Cholesterol depletion significantly affects V-ATPase activity and the initial transfer [59]. ATP6V0C affects glucose metabolism through phosphorylation during glycolysis [60]. NDUFS4 is an auxiliary subunit of the mitochondrial membrane respiratory chain NADH dehydrogenase (complex I), and NDUFS4 functions in the later stages of complex I assembly [6162]. The NADH shuttle substantially maintains mitochondrial energy metabolism and glucose-induced insulin secretion [63]. The PIGH gene encodes an endoplasmic reticulum-associated protein involved in GPI anchor biosynthesis, which are glycolipids found in many blood cells that anchor proteins to the cell surface. The protein encoded by the PIGH gene is a subunit of GPI–N-acetylglucosamine transferase (GPI–GlcNAc transferase) that transfers GlcNAc to phosphatidylinositol (PI) lipids in endoplasmic reticulum cells [64]. In short, these hub target genes may be involved in energy metabolism, lipid metabolism, and mitochondrial function through ceRNA networks. The ceRNA network comprising the seven hub target genes as a key ceRNA network for exploring the mechanism of milk fat regulation in cows.

Conclusion

We screened four DE-circRNAs using differential expression analysis, WGCNA, and circular validation. Tissue expression and subcellular localization suggested that these DE-circRNAs may have potential regulatory functions on mammary gland development and lactation in dairy cows through ceRNA networks. Thus, we constructed the ceRNA regulatory network of candidate DE-circRNAs and screened out the key ceRNA networks regulating milk fat metabolism, which helped us further explore the regulatory mechanism of milk fat metabolism. This study also provides new clues for molecular breeding of dairy cows.

Methods

Selection of experimental animals and sample preparation

Based on year-round dairy herd improvement (DHI) measurements at Nongkeng Helanshan Maosheng dairy farm, we screened 245 mid–late lactation Holstein cows with similar average daily milk yield (35.21–37.21 kg) and consistent feeding and management backgrounds (Table S6). We then screened 4 long-term high-MFP and 4 long-term low-MFP cows from 245 cows and aseptically collected fresh milk samples from each cow. One portion of each sample was sent to the testing center for DHI determination, whereas the other portion was placed in sterile water at 37 °C and returned to the laboratory for MEC isolation. The isolation, culture, and identification of MECs were completed during the early stage of our research group [65]. Cow milk MECs had the characteristic “pebble” morphology of epithelial cells (Fig. S2A), and S-shaped growth curve, which was consistent with cell growth (Fig. S2B), and the expression of epithelial cell-specific keratin 18 was positive (Figs. S2C, D). Besides, there was a significant difference in the triglyceride content of cow milk MECs of the high- and low-MFP groups (Fig. S2E), and the expression levels of the lipogenic genes SCD, PPARγ, and FASN were higher in the high-MFP group than in the low-MFP group (Fig. S3). The aforementioned sample preparation details are described in a recent paper [66].

RNA-Seq library construction and sequencing

The sequencing in the present study belongs to the same batch as that in a recent study [66] and is, therefore, methodologically identical. Total RNA was extracted from cow milk MECs using the TRIzol method. RNA degradation and contamination were monitored on 1% agarose gels. RNA purity was checked using the NanoPhotometer® spectrophotometer (IMPLEN, CA, USA). RNA integrity was assessed using the RNA Nano 6000 Assay Kit of the Bioanalyzer 2100 System (Agilent Technologies, CA, USA). The 260/280 ratio of all samples ranged from 1.70 to 1.90, and the RNA Integrity Index (RIN) was ≥ 8.00. Sample RNA for circRNA sequencing was stripped of ribosomal RNA (Epicenter Ribozero™ rRNA Removal Kit, Epicentre, USA), and linear RNA was digested with RNase R (Epicentre, USA). Sequencing libraries were prepared according to the manufacturer’s instructions for the NEBNext® Ultra™ Directional RNA Library Prep Kit for Illumina® (NEB, USA). After passing the library inspection, Illumina PE150 sequencing was performed. After filtering the raw data, the obtained clean reads were aligned with the downloaded reference genome (https://bovinegenome.elsiklab.missouri.edu/downloads/ARS-UCD1.2) using the Bowtie2 software (v2.2.8).

Identification and differential expression analysis of circRNAs

Identification of circRNAs and differential expression analysis were performed according to our previously described methods [66]. The circRNA was detected and identified using find_circ [67] and CIRI2 [68]. Transcripts per million (TPM) were used to normalize known and novel circRNAs in each sample [69]; normalized expression levels = (readCount × 1,000,000)/libsize (libsize is the sum of circRNA read counts). Differential expression analysis of transcript count matrices of high and low MFP in cow milk MECs was performed using the R package “DESeq2” [70]. The resulting P-value was adjusted using Benjamini and Hochberg’s approach for controlling the false discovery rate. Genes identified using DESeq2 with an adjusted P-value < 0.05 were designated as differentially expressed.

Weighted gene co-expression network analysis (WGCNA)

WGCNA was used for network construction and identification of consensus modules. Weighted gene network construction requires the optimal selection of soft threshold power β to improve co-expression similarity and calculate the degree of adjacency. Therefore, the function “pickSoftThreshold” (based on the criterion of approximate scale-free topology) from the R package “WGCNA” [71] was used to pick out the optimal soft threshold power β. The function “blockwiseConsensusModules” was employ used to calculate the consensus topology overlap and produce consensus modules. Based on the WGCNA analysis parameters of Yang et al. [72], we set the following: power = soft threshold power β (when r = 0.80); modules containing 20 genes as a minimum number (minModuleSize = 20); the module detection sensitivity of 2 (deepSplit = 2); module merged cut height of 0.25 (mergeCutHeight = 0.25, i.e., merged into one module if the correlation coefficient of eigengenes within the module was greater than 0.75). To avoid rearrangement of eigengene within modules according to intramodular connectivity (KME), we set the following parameters: minKMEtoStay = 0, maxBlockSize = 10,000, and the remaining parameters followed the default values of the function. Subsequently, the co-expression network was built using the function “blockwiseModules”, with the following parameters: power = soft threshold power β; TOMType = “notations”; minModuleSize = 20; mergeCutHeigh t = 0.25; maxBlockSize = 20,000; pamRespectsDendro = FALSE; verbose = 3; the other parameters were set to default. This process produces co-expression modules that are significantly correlated with MFP. Finally, the GS and MM of the eigengenes in the consensus module were calculated using the function “corAndPvalue”, and the selection criteria for circRNAs were|MM| ≥ 0.90,|GS| ≥0.60.

Circular structure verification

To verify the circular structure of circRNA, we designed divergent primers and convergent primers for each circRNA. First, circRNAs and linear RNAs were tested for resistance to RNase R. We performed RNase R treatment on the total RNA. A portion of the total RNAs was added with 5U/µg RNase R and 2 µl of 10× RNase R Reaction Buffer, whereas the other portion was added with an equal amount of RNase-free water and 2 µl of 10× RNase R Reaction Buffer. After incubating at 37 °C for 30 min. the RNase R-treated RNA was purified using the RNeasy MinElute Cleanup Kit (QIAGEN, Germany). The RNA was then reverse transcribed to cDNA, and the expression levels of circRNA and linear RNA were detected using PCR and reverse transcription-quantitative polymerase chain reaction (RT-qPCR). Subsequently, the head-to-tail splice sites of the circRNAs were identified through Sanger sequencing. The circular structure was verified using our previous verification method [66]. Table S7 lists the primer sequences used in the present study.

Tissue expression and subcellular localization

The experimental methodology in this section is also consistent with our previous research methodology [66]. TRIzol reagent was used to extract the total RNA from tissues of the heart, liver, kidney, uterus, ovaries, small intestine, and mammary gland. RNA was isolated from the cell cytoplasm and nucleus using the Cytoplasmic and Nuclear RNA Purification kit (Norgen Biotek), and GAPDH and U6 were used as cytoplasmic and nuclear fractionation indicators, respectively. The first-strand cDNA was synthesized using the PrimeScript RT Reagent Kit with gDNA Eraser (Takara, Dalian, China). RT-qPCR was used to detect the expression of circRNAs in various tissues, the cytoplasm, and the nucleus.

Target relationship prediction

CircRNA/miRNA and miRNA/mRNA interaction pairs were predicted using miRanda (v3.3a) [24] and TargetScan (v7.2) [25] software. circRNA/miRNA and miRNA/mRNA interactions were analyzed based on TargetScan’s Context + and miRanda’s Free Energy criterion.

Analysis of hub genes in PPI networks

A PPI network analysis was performed following the methodology previously described by Yang et al. [72]. Protein network interactions were obtained using the Strings website (https://string-db.org/, v11.0), where a minimum interaction score of 0.90 deemed sufficient to obtain high-confidence protein network interactions. The MCODE plugin in Cytoscape was applied to identify critical subnetworks and seeds of nodes (or hub genes). The CytoHubba plugin in Cytoscape detects hub genes by four centrality methods—Degree, Edge Percolation Component (EPC), Maximum Cliff Centrality (MCC), and Maximum Neighborhood Component (MNC)—which are practical methods for identifying hub genes from PPI networks [73]. Subsequently, functional enrichment analysis was performed on the genes of key subnetworks. The function “enrichGO” was applied to the annotation of GO, including BP, MF, and CC. The function “enrichKEGG” was applied for the KEGG annotations to uncover relevant signaling pathways. All enrichment analysis results were visualized using the R package “ggplot2”.

Data availability

All data generated or analyzed in this study are included in this article (and its Supplementary Information file), and the datasets have been submitted to the SRA database with the accession number PRJNA730595. Data can be accessed at: https://www.ncbi.nlm.nih.gov/sra/PRJNA730595.

References

  1. Jun X, Mingyue Z, Lingjie LI, Hou X, Zeng W. Conjugated linoleic acid improves glucose and lipid metabolism in diabetic mice. J South Med Univ. 2019;39:740–6.

    Google Scholar 

  2. Kritchevsky D, Tepper SA, Wright S, Czarnecki SK. Influence of graded levels of conjugated linoleic acid (CLA) on experimental atherosclerosis in rabbits. Nutr Res. 2002;22:1275–9.

    Article  CAS  Google Scholar 

  3. McCarthy C, Lieggi NT, Barry D, Mooney D, de Gaetano M, James WG, McClelland S, Barry MC, Escoubet-Lozach L, Li AC, Glass CK, Fitzgerald DJ, Belton O. Macrophage PPAR gamma co-activator-1 alpha participates in repressing foam cell formation and atherosclerosis in response to conjugated linoleic acid. EMBO Mol Med. 2013;5:1443–57.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  4. Bach AC, Ingenbleek Y, Frey A. The usefulness of dietary medium-chain triglycerides in body weight control: fact or fancy? J Lipid Res. 1996;37(4):708–26.

    Article  CAS  PubMed  Google Scholar 

  5. Kasai M, Maki H, Nosaka N, Aoyama T, Ooyama K, Uto H, Okazaki M, Igarashi O, Kondo K. Effect of medium-chain triglycerides on the postprandial triglyceride concentration in healthy men. Biosci Biotechnol Biochem. 2003;67(1):46–53.

    Article  CAS  PubMed  Google Scholar 

  6. Haenlein G. Goat milk in human nutrition. Small Ruminant Res. 2004;51:155–63.

    Article  Google Scholar 

  7. Seth K, Bajwa U. Effect of acidulants on the recovery of milk constituents and quality of Mozzarella processed cheese. J Food Sci Technol. 2015;52(3):1561–9.

    Article  CAS  PubMed  Google Scholar 

  8. Pagthinathan M, Nafees MSM. Biochemistry of cheese ripening. AGRIEAST J Agricultural Sci. 2017;10:16.

    Article  Google Scholar 

  9. Kristensen LS, Andersen MS, Stagsted LVW, Ebbesen KK, Hansen TB, Kjems J. The biogenesis, biology and characterization of circular RNAs. Nat Rev Genet. 2019;20(11):675–91.

    Article  CAS  PubMed  Google Scholar 

  10. Salzman J, Chen RE, Olsen MN, Wang PL, Brown PO. Cell-type specific features of circular RNA expression. PLoS Genet. 2013;9(9):e1003777.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  11. Zhang M, Ma L, Liu Y, He Y, Li G, An X, Cao B. CircRNA-006258 sponge-adsorbs mir-574-5p to regulate cell growth and milk synthesis via EVI5L in goat mammary epithelial cells. Genes (Basel). 2020;11(7):718.

    Article  PubMed  Google Scholar 

  12. Wu X, Zhen H, Liu Y, Li L, Luo Y, Liu X, Li S, Hao Z, Li M, Hu L, Qiao L, Wang J. Tissue-specific expression of circ_015343 and its inhibitory effect on mammary epithelial cells in sheep. Front Vet Sci. 2022;9:919162.

    Article  PubMed  PubMed Central  Google Scholar 

  13. Chen Z, Cao X, Lu Q, Zhou J, Wang Y, Wu Y, Mao Y, Xu H, Yang Z. circ01592 regulates unsaturated fatty acid metabolism through adsorbing miR-218 in bovine mammary epithelial cells. Food Funct. 2021;12(23):12047–58.

    Article  CAS  PubMed  Google Scholar 

  14. Chen Z, Zhou J, Wang M, Liu J, Zhang L, Loor JJ, Liang Y, Wu H, Yang Z. Circ09863 regulates unsaturated fatty acid metabolism by adsorbing miR-27a-3p in bovine mammary epithelial cells. J Agric Food Chem. 2020;68(32):8589–601.

    Article  CAS  PubMed  Google Scholar 

  15. Zhu C, Jiang Y, Zhu J, He Y, Yin H, Duan Q, Zhang L, Cao B, An X. CircRNA8220 sponges miR-8516 to regulate cell viability and milk synthesis via Ras/MEK/ERK and PI3K/AKT/mTOR pathways in goat mammary epithelial cells. Anim (Basel). 2020;10(8):1347.

    Google Scholar 

  16. Sabino M, Carmelo VAO, Mazzoni G, Cappelli K, Capomaccio S, Ajmone-Marsan P, Verini-Supplizi A, Trabalza-Marinucci M, Kadarmideen HN. Gene co-expression networks in liver and muscle transcriptome reveal sex-specific gene expression in lambs fed with a mix of essential oils. BMC Genomics. 2018;19(1):236.

    Article  PubMed  PubMed Central  Google Scholar 

  17. Wang J, Sui J, Mao C, Li X, Chen X, Liang C, Wang X, Wang SH, Jia C. Identification of key pathways and genes related to the development of hair follicle cycle in cashmere goats. Genes (Basel). 2021;12(2):180.

    Article  CAS  PubMed  Google Scholar 

  18. de Oliveira PSN, Coutinho LL, Cesar ASM, Diniz W J D S, de Souza MM, Andrade BG, Koltes JE, Mourão GB, Zerlotini A, Reecy JM. Regitano L C A. co-expression networks reveal potential regulatory roles of mirnas in fatty acid composition of nelore cattle. Front Genet. 2019;10:651.

    Article  PubMed  PubMed Central  Google Scholar 

  19. Do DN, Dudemaine PL, Fomenky BE, Ibeagha-Awemu EM. Integration of miRNA weighted gene co-expression network and miRNA-mRNA co-expression analyses reveals potential regulatory functions of miRNAs in calf rumen development. Genomics. 2019;111(4):849–59.

    Article  CAS  PubMed  Google Scholar 

  20. Ling Y, Zheng Q, Sui M, Zhu L, Xu L, Zhang Y, Liu Y, Fang F, Chu M, Ma Y, Zhang X. Comprehensive analysis of LncRNA reveals the temporal-specific module of goat skeletal muscle development. Int J Mol Sci. 2019;20(16):3950.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  21. Wang J, Chai Z, Deng L, Wang J, Wang H, Tang Y, Zhong J, Ji Q. Detection and integrated analysis of lncRNA and mRNA relevant to plateau adaptation of Yak. Reprod Domest Anim. 2020;55(11):1461–9.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  22. Shen M, Li T, Chen F, Wu P, Wang Y, Chen L, Xie K, Wang J, Zhang G. Transcriptomic analysis of circRNAs and mRNAs reveals a complex regulatory network that participate in follicular development in chickens. Front Genet. 2020;11:503.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  23. Deng R, Cui X, Dong Y, Tang Y, Tao X, Wang S, Wang J, Chen L. Construction of circRNA-Based ceRNA network to reveal the role of circRNAs in the progression and prognosis of hepatocellular carcinoma. Front Genet. 2021;12:626764.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  24. Miranda KC, Huynh T, Tay Y, Ang YS, Tam WL, Thomson AM, Lim B, Rigoutsos I. A pattern-based method for the identification of MicroRNA binding sites and their corresponding heteroduplexes. Cell. 2006;126(6):1203–17.

    Article  CAS  PubMed  Google Scholar 

  25. Jan CH, Friedman RC, Ruby JG, Bartel DP. Formation, regulation and evolution of Caenorhabditis elegans 3’UTRs. Nature. 2011;469(7328):97–101.

    Article  CAS  PubMed  Google Scholar 

  26. Pei G, Chen L, Zhang W. WGCNA application to proteomic and metabolomic data analysis. Methods Enzymol. 2017;585:135–58.

    Article  CAS  PubMed  Google Scholar 

  27. Zhu M, Xie H, Wei X, Dossa K, Yu Y, Hui S, Tang G, Zeng X, Yu Y, Hu P, Wang J. WGCNA analysis of salt-responsive core transcriptome identifies novel hub genes in rice. Genes. 2019;10(9):719.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  28. Galán-Vásquez E, Perez-Rueda E. Identification of modules with similar gene regulation and metabolic functions based on co-expression data. Front Mol Biosci, 6:139.

  29. Yu X, Fang X, Gao M, Mi J, Zhang X, Xia L, Zhao Z, Albrecht E, Maak S, Yang R. Isolation and identification of bovine preadipocytes and screening of micrornas associated with adipogenesis. Animals. 2020;10(5):818.

    Article  PubMed  PubMed Central  Google Scholar 

  30. Streuli CH, Bailey N, Bissell MJ. Control of mammary epithelial differentiation: basement membrane induces tissue-specific gene expression in the absence of cell-cell interaction and morphological polarity. J Cell Biol. 1991;115(5):1383–95.

    Article  CAS  PubMed  Google Scholar 

  31. Strange R, Li F, Saurer S, Burkhardt A, Friis RR. Apoptotic cell death and tissue remodelling during mouse mammary gland involution. Development. 1992;115(1):49–58.

    Article  CAS  PubMed  Google Scholar 

  32. Pullan S, Wilson J, Metcalfe A, Edwards GM, Goberdhan N, Tilly J, Hickman JA, Dive C, Streuli CH. Requirement of basement membrane for the suppression of programmed cell death in mammary epithelium. J Cell Sci. 1996;109(3):631–42.

    Article  CAS  PubMed  Google Scholar 

  33. Kohjima M, Higuchi N, Kato M, Kotoh K, Yoshimoto T, Fujino T, Yada M, Yada R, Harada N, Enjoji M, Takayanagi R, Nakamuta M. SREBP-1c, regulated by the insulin and AMPK signaling pathways, plays a role in nonalcoholic fatty liver disease. Int J Mol Med. 2008;21(4):507–11.

    CAS  PubMed  Google Scholar 

  34. van Meer G, Simons K. The function of tight junctions in maintaining differences in lipid composition between the apical and the basolateral cell surface domains of MDCK cells. EMBO J. 1986;5(7):1455–64.

    Article  PubMed  PubMed Central  Google Scholar 

  35. Shigetomi K, Ono Y, Inai T, Ikenouchi J. Adherens junctions influence tight junction formation via changes in membrane lipid composition. J Cell Biol. 2018;217(7):2373–81.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  36. Cohain AT, Barrington WT, Jordan DM, Beckmann ND, Argmann CA, Houten SM, Charney AW, Ermel R, Sukhavasi K, Franzen O, Koplev S, Whatling C, Belbin GM, Yang J, Hao K, Kenny EE, Tu Z, Zhu J, Gan LM, Do R, Giannarelli C, Kovacic JC, Ruusalepp A, Lusis AJ, Bjorkegren JLM, Schadt EE. An integrative multiomic network model links lipid metabolism to glucose regulation in coronary artery disease. Nat Commun. 2021;12(1):547.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  37. Voillet V, San Cristobal M, Père MC, Billon Y, Canario L, Liaubet L, Lefaucheur L. Integrated Analysis of Proteomic and Transcriptomic Data highlights late fetal muscle maturation process. Mol Cell Proteom. 2018;17(4):672–93.

    Article  CAS  Google Scholar 

  38. Ruan D, Zhuang Z, Ding R, Qiu Y, Zhou S, Wu J, Xu C, Hong L, Huang S, Zheng E, Cai G, Wu Z, Yang J. Weighted single-step GWAS identified candidate genes Associated with Growth traits in a Duroc Pig Population. Genes (Basel). 2021;12(1):117.

    Article  CAS  PubMed  Google Scholar 

  39. Yen MC, Chou SK, Kan JY, Kuo PL, Hou MF, Hsu YL. New insight on Solute Carrier Family 27 Member 6 (SLC27A6) in Tumoral and non-tumoral breast cells. Int J Med Sci. 2019;16(3):366–75.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  40. Suto JI, Kojima M. Identification of quantitative trait loci that determine plasma total-cholesterol and triglyceride concentrations in DDD/Sgn and C57BL/6J inbred mice. Cholesterol. 2017;2017:3178204.

    Article  PubMed  PubMed Central  Google Scholar 

  41. Zhang P, He Q, Wang Y, Zhou G, Chen Y, Tang L, Zhang Y, Hong X, Mao Y, He Q, Yang X, Liu N, Ma J. Protein C receptor maintains cancer stem cell properties via activating lipid synthesis in nasopharyngeal carcinoma. Signal Transduct Target Ther. 2022;7(1):46.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  42. Sato H, Taketomi Y, Ushida A, Isogai Y, Kojima T, Hirabayashi T, Miki Y, Yamamoto K, Nishito Y, Kobayashi T, Ikeda K, Taguchi R, Hara S, Ida S, Miyamoto Y, Watanabe M, Baba H, Miyata K, Oike Y, Gelb MH, Murakami M. The adipocyte-inducible secreted phospholipases PLA2G5 and PLA2G2E play distinct roles in obesity. Cell Metab. 2014;20(1):119–32.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  43. Du X, Kazim AS, Dawes IW, Brown AJ, Yang H. The AAA ATPase VPS4/SKD1 regulates endosomal cholesterol trafficking independently of ESCRT-III. Traffic. 2013;14(1):107–19.

    Article  CAS  PubMed  Google Scholar 

  44. Ford CE, Skiba NP, Bae H, Daaka Y, Reuveny E, Shekter LR, Rosal R, Weng G, Yang CS, Iyengar R, Miller RJ, Jan LY, Lefkowitz RJ, Hamm HE. Molecular basis for interactions of G protein betagamma subunits with effectors. Science. 1998;280(5367):1271–4.

    Article  CAS  PubMed  Google Scholar 

  45. Hua C, Geng Y, Niu L, Chen Q, Cai L, Tao S, Ni Y, Zhao R. Stimulating lipolysis in subcutaneous adipose tissues by chronic dexamethasone administration in goats. Livest ence. 2018;214:62–7.

    Article  Google Scholar 

  46. Heinonen S, Saarinen L, Naukkarinen J, Rodríguez A, Frühbeck G, Hakkarainen A, Lundbom J, Lundbom N, Vuolteenaho K, Moilanen E, Arner P, Hautaniemi S, Suomalainen A, Kaprio J, Rissanen A, Pietiläinen KH. Adipocyte morphology and implications for metabolic derangements in acquired obesity. Int J Obes (Lond). 2014;38(11):1423–31.

    Article  CAS  PubMed  Google Scholar 

  47. Dong Z, Ba H, Zhang W, Coates D, Li C. iTRAQ-based quantitative proteomic analysis of the potentiated and dormant antler stem cells. Int J Mol Sci. 2016;17(11):1778.

    Article  PubMed  PubMed Central  Google Scholar 

  48. Foust DJ, Godin AG, Ustione A, Wiseman PW, Piston DW. Two-color spatial cumulant analysis detects heteromeric interactions between membrane proteins. Biophys J. 2019;117(9):1764–77.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  49. Stephens LR, Eguinoa A, Erdjument-Bromage H, Lui M, Cooke F, Coadwell J, Smrcka AS, Thelen M, Cadwallader K, Tempst P, Hawkins PT. The G beta gamma sensitivity of a PI3K is dependent upon a tightly associated adaptor. Cell. 1997;89(1):105–14.

    Article  CAS  PubMed  Google Scholar 

  50. Shi CS, Lee SB, Sinnarajah S, Dessauer CW, Rhee SG, Kehrl JH. Regulator of G-protein signaling 3 (RGS3) inhibits Gbeta1gamma 2-induced inositol phosphate production, mitogen-activated protein kinase activation, and akt activation. J Biol Chem. 2001;276(26):24293–300.

    Article  CAS  PubMed  Google Scholar 

  51. Rhee SG, Bae YS. Regulation of phosphoinositide-specific phospholipase C isozymes. J Biol Chem. 1997;272(24):15045–8.

    Article  CAS  PubMed  Google Scholar 

  52. Faenza I, Bavelloni A, Fiume R, Santi P, Martelli AM, Maria Billi A, Lo Vasco VR, Manzoli L, Cocco L. Expression of phospholipase C beta family isoenzymes in C2C12 myoblasts during terminal differentiation. J Cell Physiol. 2004;200(2):291–6.

    Article  CAS  PubMed  Google Scholar 

  53. Hajicek N, Keith NC, Siraliev-Perez E, Temple BR, Huang W, Zhang Q, Harden TK, Sondek J. Structural basis for the activation of PLC-γ isozymes by phosphorylation and cancer-associated mutations. Elife. 2019;8:e51700.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  54. Moreira GCM, Boschiero C, Cesar ASM, Reecy JM, Godoy TF, Trevisoli PA, Cantão ME, Ledur MC, Ibelli AMG, Peixoto JO, Moura A S A M T, Garrick D, Coutinho LL. A genome-wide association study reveals novel genomic regions and positional candidate genes for fat deposition in broiler chickens. BMC Genomics. 2018;19(1):374.

    Article  PubMed  PubMed Central  Google Scholar 

  55. Srivastava S, Srikanth K, Won S, Son JH, Park JE, Park W, Chai HH, Lim D. Haplotype-based genome-wide association study and identification of candidate genes associated with carcass traits in hanwoo cattle. Genes. 2020;11(5):551.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  56. Hirao J, Tojo A, Hatakeyama S, Satonaka H, Ishimitsu T. V-ATPase blockade reduces renal gluconeogenesis and improves insulin secretion in type 2 diabetic rats. Hypertens Res. 2020;43(10):1079–88.

    Article  CAS  PubMed  Google Scholar 

  57. Li SC, Diakov TT, Xu T, Tarsio M, Zhu W, Couoh-Cardel S, Weisman LS, Kane PM. The signaling lipid PI(3,5)P2 stabilizes V1-V(o) sector interactions and activates the V-ATPase. Mol Biol Cell. 2014;25(8):1251–62.

    Article  PubMed  PubMed Central  Google Scholar 

  58. Banerjee S, Clapp K, Tarsio M, Kane PM. Interaction of the late endo-lysosomal lipid PI(3,5)P2 with the Vph1 isoform of yeast V-ATPase increases its activity and cellular stress tolerance. J Biol Chem. 2019;294(23):9161–71.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  59. Costa GA, de Souza SB, da Silva Teixeira LR, Okorokov LA, Arnholdt ACV, Okorokova-Façanha AL. Façanha A R. Tumor cell cholesterol depletion and V-ATPase inhibition as an inhibitory mechanism to prevent cell migration and invasiveness in melanoma. Biochim Biophys Acta Gen Subj. 2018;1862(3):684–91.

    Article  CAS  PubMed  Google Scholar 

  60. Son SW, Chau GC, Kim ST, Um SH. Vacuolar H+-ATPase subunit V0C regulates aerobic glycolysis of esophageal cancer cells via PKM2 signaling. Cells. 2019;8(10):1137.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  61. Scacco S, Petruzzella V, Budde S, Vergari R, Tamborra R, Panelli D, van den Heuvel LP, Smeitink JA, Papa S. Pathological mutations of the human NDUFS4 gene of the 18-kDa (AQDQ) subunit of complex I affect the expression of the protein and the assembly and function of the complex. J Biol Chem. 2003;278(45):44161–7.

    Article  CAS  PubMed  Google Scholar 

  62. Pereira B, Videira A, Duarte M. Novel insights into the role of Neurospora Crassa NDUFAF2, an evolutionarily conserved mitochondrial complex I assembly factor. Mol Cell Biol. 2013;33(13):2623–34.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  63. Eto K, Tsubamoto Y, Terauchi Y, Sugiyama T, Kishimoto T, Takahashi N, Yamauchi N, Kubota N, Murayama S, Aizawa T, Akanuma Y, Aizawa S, Kasai H, Yazaki Y, Kadowaki T. Role of NADH shuttle system in glucose-induced activation of mitochondrial metabolism and insulin secretion. Science. 1999;283(5404):981–5.

    Article  CAS  PubMed  Google Scholar 

  64. Watanabe R, Inoue N, Westfall B, Taron CH, Orlean P, Takeda J, Kinoshita T. The first step of glycosylphosphatidylinositol biosynthesis is mediated by a complex of PIG-A, PIG-H, PIG-C and GPI1. EMBO J. 1998;17(4):877–85.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  65. Mu T, Hu H, Feng X, Ma Y, Wang Y, Liu J, Yu B, Wen W, Zhang J, Gu Y. Screening and conjoint analysis of key lncRNAs for milk fat metabolism in dairy cows. Front Genet. 2022;13:772115.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  66. Feng X, Cai Z, Gu Y, Mu T, Yu B, Ma R, Liu J, Wang C, Zhang J. Excavation and characterization of key circRNAs for milk fat percentage in Holstein cattle. J Anim Sci. 2023;101:skad157.

    Article  PubMed  Google Scholar 

  67. Memczak S, Jens M, Elefsinioti A, Torti F, Krueger J, Rybak A, Maier L, Mackowiak SD, Gregersen LH, Munschauer M, Loewer A, Ziebold U, Landthaler M, Kocks C, le Noble F, Rajewsky N. Circular RNAs are a large class of animal RNAs with regulatory potency. Nature. 2013;495(7441):333–8.

    Article  CAS  PubMed  Google Scholar 

  68. Gao Y, Zhang J, Zhao F. Circular RNA identification based on multiple seed matching. Brief Bioinform. 2018;19(5):803–10.

    Article  CAS  PubMed  Google Scholar 

  69. Zhou L, Chen J, Li Z, Li X, Hu X, Huang Y, Zhao X, Liang C, Wang Y, Sun L, Shi M, Xu X, Shen F, Chen M, Han Z, Peng Z, Zhai Q, Chen J, Zhang Z, Yang R, Ye J, Guan Z, Yang H, Gui Y, Wang J, Cai Z, Zhang X. Integrated profiling of microRNAs and mRNAs: microRNAs located on Xq27.3 associate with clear cell renal cell carcinoma. PLoS ONE. 2010;5(12):e15224.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  70. Love MI, Huber W, Anders S. Moderated estimation of Fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15(12):550.

    Article  PubMed  PubMed Central  Google Scholar 

  71. Zhang B, Horvath S. A general framework for weighted gene co-expression network analysis. Stat Appl Genet Mol Biol. 2005;4:17.

    Article  Google Scholar 

  72. Yang C, Ding Y, Dan X, Shi Y, Kang X. Multi-transcriptomics reveals RLMF axis-mediated signaling molecules associated with bovine feed efficiency. Front Vet Sci. 2023;10:1090517.

    Article  PubMed  PubMed Central  Google Scholar 

  73. Chin CH, Chen SH, Wu HH, Ho CW, Ko MT, Lin CY. cytoHubba: identifying hub objects and sub-networks from complex interactome. BMC Syst Biol. 2014;4(4):S11.

    Article  Google Scholar 

Download references

Acknowledgements

We thank all laboratory members for their comments and suggestions. We are grateful to Dr. Mingjie Chen from NewCore Biotech Co., Ltd. Shanghai, for helping us conduct bioinformatics analysis. We thank Helanshan Maosheng Dairy Farm for providing us with experimental animals and dairy herd improvement (DHI) measurements data.

Funding

This project is supported by the Key Research Project of the Ningxia Hui Autonomous Region (Grant No: 2022BBF02017) and the special breeding project of high-quality and high-yield dairy cows in the Ningxia Autonomous Region (Grant No: 2019NYYZ05).

Author information

Authors and Affiliations

Authors

Contributions

X.F. Data analysis and manuscript writing; L.T. Experimental verification; R.M. and L.M. Article grammar modification; J.Z. and Y.G. Revised the manuscript and provided reagents; T.M. and C.W. Isolation and culture of mammary epithelial cells; B.Y. and J.L. Conceptual analysis and article modification.

Corresponding author

Correspondence to Juan Zhang.

Ethics declarations

Ethics approval and consent to participate

All animal experiments involved in this study were carried out under the approval of the Laboratory Animal Welfare and Ethics Review Committee of Ningxia University, ethical approval number NXU-2023-068. All experimental procedures were conducted following the guidelines of the Laboratory Animal Welfare and Ethical Review Committee of Ningxia University. Consent was obtained from the dairy farmers for the use of animals in the study.

Consent for publication

Not applicable.

Competing interests

The authors declare no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Electronic supplementary material

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Feng, X., Tong, L., Ma, L. et al. Mining key circRNA-associated-ceRNA networks for milk fat metabolism in cows with varying milk fat percentages. BMC Genomics 25, 323 (2024). https://doi.org/10.1186/s12864-024-10252-y

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12864-024-10252-y

Keywords