Screening and Joint Analysis of Key lncRNAs for Milk Fat Metabolism in Dairy Cows

Background: Long noncoding RNAs (lncRNAs) play an important regulatory role in various biological processes as a key regulatory factor. However, there are largely unknown for the function and expression prole of lncRNAs in milk fat synthesis of dairy cows. Results: In this study, RNA sequencing (RNA-seq) was used to research the whole genome expression of lncRNAs and mRNA transcripts in bovine mammary epithelial cells (BMECs) of dairy cows with high and low milk fat percentage (MFP), and joint analysis was carried out. We identied a total of 47 differentially expressed genes (DEGs) and 38 differentially expressed lncRNAs (DELs, Padj < 0.05), 11 candidate DEGs that may regulate milk fat metabolism were screened by enrichment analysis. Downregulated differential gene ENPP2 and upregulated differential gene BCAT1 are more likely to participate in the milk fat metabolism, and its function needs further experiments verication. The enrichment analysis of target genes predicted by DELs identied 7 cis (co-localization) and 10 trans (co-expression) candidate target genes related to milk lipid metabolism, corresponding to a total of 18 DELs. Among them, the targeting relationship between long intervening/intergenic noncoding RNA (lincRNA) TCONS_00082721 and FABP4 gene that predicts milk fat metabolism by co-localization and co-expression is worthy of attention. Based on the expression information of DELs, differential microRNAs (miRNAs), and lipid metabolism-related target genes, 156 competing endogenous RNAs (ceRNAs) interaction regulation networks related to milk fat metabolism were constructed. The regulatory network centered on miR-145 will be the focus of subsequent experimental research. The ceRNAs regulatory network related to TCONS_00082721 and TCONS_00172817 are more likely to be involved in milk fat synthesis. Conclusions: These results will provide new ways to understand the complex biology of dairy cow milk fat synthesis and provide valuable information for the breed improvement of Chinese Holstein cattle.


Background
Milk quality is affected by many factors, such as population genetic structure, reproductive performance, feeding and management level, and is closely related to the main raw material composition of milk. Milk fat is an important component in butter and yogurt, and its content and composition are the main reference elements in milk quality evaluation. Nowadays, milk fat content is not only one of the important indicators of core competitiveness of dairy industry, but also the main target feature of dairy cow breeding [1]. To a certain extent, it also plays an important role in nutrition and metabolism during human growth and development [2]. Therefore, the exploration of the theory and ways of milk fat formation and regulation, and the improvement of milk fat content have become the focus of lactation biology research in the world.
In the past, many scholars have extensively studied on the complex regulatory mechanisms of breast development and elucidated the main ways of milk fat synthesis (including ab initio synthesis and absorption of FA in the blood) [3]. Quantitative real-time reverse transcription PCR (qRT-PCR) can better discover the transcriptional regulation steps of milk fat synthesis. High-throughput sequencing has identi ed a series of effective genes and regulatory factors in milk fat metabolism [4]. As we all know, the study of gene regulation mainly depends on protein-coding genes. However, more and more evidence suggests that non-coding RNAs (ncRNAs), which are mistaken for "transcriptional noise", are also important factors in regulating complex developmental processes in organisms [5]. Long noncoding RNAs (lncRNAs) are newly identi ed long-stranded non-coding RNAs in the breast, with a length of more than 200 nucleotides and more tissue-speci c than coding genes [6]. The mode of action of lncRNAs on genes is more complex than that of microRNAs (miRNAs). At present, transcriptional interference is a more mature regulatory model [7]. Functional studies in mammals have found that lncRNAs regulate gene expression at epigenetic, transcriptional and post-transcriptional levels, and also play key regulatory roles in a variety of biological processes [8]. In regulating milk fat metabolism in goats, Yu et al. [9] found that compared with early lactation, the expression trend of lncRNAs TCONS-00144434, TCONS-00148514 and TCONS-00055666 during mature lactation was the same as that of FASN, LPL, GPAM and MSMO1 genes with the function of fatty acid biosynthesis and cholesterol storage (The full names of all genes in this paper are shown in Supplementary Table 1). However, the speci c regulatory mechanism remains to be further studied. Chao et al. [10] found that the target genes of TCONS_00162862 in dairy cow were involved in lipid transporter protein activity, ligase activity, and fatty acid transport, and mainly enriched in MAPK, PI3K-Akt, and insulin signaling pathways. Zheng et al. [11] performed transcriptome sequencing on four breast tissues of Holstein cows during peak lactation and late lactation, and found that the target genes of differentially expressed LNC-XLOC_274111 were FABP3 and FABP4, which would encode fatty acid-binding proteins in the bovine mammary gland. The target genes of XLOC_000752 were signi cantly enriched in PPAR, AMPK, and insulin signaling pathways, as well as glycerol metabolism pathways.
So far, research on lncRNAs in milk fat metabolism of ruminants are relatively lacking. The functions and expression pro les of lncRNAs in the bovine mammary gland are still unclear, and there are even fewer studies to screen differentially expressed lncRNAs (DELs) for differences in high and low milk fat percentage (MFP). In addition, the dynamic pattern of the interaction of ceRNA in milk fat synthesis remains unknown. To better understand the metabolic process of milk fat synthesis during lactation, Holstein cows with extreme differences in MFP were screened, fresh milk was collected and bovine mammary epithelial cells (BMECs) were isolated. The differentially expressed genes (DEGs) and DELs related to milk fat synthesis in the same lactation stage were screened by transcriptomics technology, and the prediction of lncRNAs targeting miRNAs to construct a ceRNA regulatory network of genes, miRNAs and lncRNAs related to milk fat metabolism, It is an important step in understanding the complex biological processes involved in milk production,and provides valuable information for the improvement of Chinese Holstein cattle breeds.

Results
Isolation, culture, and identi cation of BMECs from milk As shown in Fig.1, the fresh milk (per 40 mL) of dairy cows with high MFP and low MFP was centrifuged at 1000 r/min for 10 min. It is found that the milk fat layer area of high MFP dairy cows are higher than that of low MFP dairy cows, and there are also different milk fat layer thickness (a, b). The isolated BMECs is a unique "pebble" epithelial cells shape, and c-f is the ow chart of the isolation BMECs, genetics, growth, and cell secretory characteristics are normal [12]. The cell growth curve conforms to the "S" rule, indicating that the isolated cells have grown well [13]. Keratin 18 belongs to the intermediate brin family, is participated in the formation of the cytoskeleton, and is the marker of epithelial cells [13]. In this study, immuno uorescence showed that keratin 18 was positive (g, h), and the expression of keratin 8 (Epithelial cell-speci c protein) was detected in both high and low groups (Fig. 2a), indicating that the isolated cells had speci c epithelial cell characteristics [12]. Chromosome analysis showed that the cells were in a normal diploid con guration and contained 60 chromosomes, which was consistent with the chromosome number of dairy cows (i). BMECs TAG content and expression level of adipogenic genes in high and low milk fat groups Shen et al. [14] showed that the isolated BMECs maintained cell-speci c genetic, structural, and biological functions in the 2 nd and 5 th passages. To con rm whether the BMECs isolated from the high-milk fat group and low-milk fat group have the differences in milk fat synthesis at the cellular level, we detected the changes in the TAG content and adipogenic gene expression of F2 generation BMECs. The results found that there was a great difference in TAG content, and the expression levels of fat-related genes SCD, PPAR γ, and FABP3 were higher in high-milk fat group than those in low-milk fat group, which further con rmed the difference between high and low MFP groups at the molecular level. It has laid a reliable material foundation for subsequent transcriptome sequencing and other experiments (Fig. 2).
CPC/PFAM/CNCI software predicted that there were 25 702 novel_mRNA with coding potential and 7602 candidate novel_lncRNAs without coding potential. Among all the candidate novel_lncRNAs, the lncRNAs located in the intergenic region account for 51.8%, the antisense contains 24.75%, the sense overlapping contains 23.50%, and no sense intronic. After obtaining the FPKM of all transcripts, the distribution of transcripts expression levels of different samples was shown by block diagram (Fig. 3a), indicating that the BMECs expression level in cows with high MFP was slightly higher than that in cows with low MFP. The correlation of gene expression level between samples is an important index for testing the experimental reliability and sample selection. In this study, the correlation coe cients of the intra-group and inter-group samples were calculated based on the FPKM values of all genes in each sample and drawn into a heat map (Fig. 3b). As can be seen from the gure, intra-group correlation coe cients are all above 0.96, graeter than the ideal sampling and experimental conditions (R 2 = 0.92), indicating a high degree of similarity in expression patterns among samples and a difference between groups. Therefore, we think that the transcriptome sequencing results are reliable and can be used for subsequent analysis.

Differential expression and clustering analysis of DEGs and DELs
Among all the differential expressed mRNAs screened, a total of 15 573 mRNAs were found in the highmilk fat group and the low-milk fat group, 717 in the high-milk fat group, and 614 in the low-milk fat group, However, the number of lncRNAs shared between the high-milk fat group and the low-milk fat group was relatively small (4 666), with 372 speci c lncRNAs in the high-milk fat group and 365 lncRNAs in the low-milk fat group. A total of 47 signi cant DEGs and 38 DELs were screened from BMECs with high and low MFP. The FPKM of DEGs and DELs were transformed by log 10 (FPKM+1) and then clustered by heat map (Fig. 4). DEGs or DELs with similar expression patterns cluster together and differ signi cantly in color, indicating that these genes share a common function or are involved in a common metabolic pathways. Compared with the low-milk fat group, 24/27 DEGs/DELs expression in the highmilk fat group was signi cantly upregulated, and 23/11 DEGs/DELs expression was signi cantly downregulated.
DEGs enrichment analysis and screening of DEGs related to lipid metabolism Based on the functional enrichment analysis of 47 DEGs, a total of 173 signi cantly enriched GO items were identi ed, as shown in Fig. 5a, showing the rst 20 signi cantly enriched GO items (the same below). Among the signi cantly enriched GO items, the BP related to lipid metabolism is the regulation of phospholipid translocation, positive regulation of phospholipid translocation, and glycerophospholipid catabolic process. The MF item is triglyceride lipase activity, lipase activity, aminophospholipid transporter activity, phosphodiesterase I activity, sterol esterase activity, and carboxylic ester hydrolase activity. KEGG enrichment analysis showed that 6 signaling pathways were signi cantly enriched (Fig.   5b), and the pathways related to lipid metabolism were Wnt, Rap1 and 2-Oxocarboxylic acid metabolism. 11 DEGs that may regulate milk fat metabolism were screened out ( Table 1). The enrichment of GO/KEGG is shown in Supplementary Table 2. PDGFD, BCAT1, and APOL3 genes expression was upregulated, while ATP8A2, PTPRR, KCNMA1, ZFYVE28, ENPP2, DKK1, CES4A, and CTSH genes expression was downregulated. ENPP2 gene is signi cantly enriched to BP, MF, and KEGG pathway related to lipid metabolism, and it is inferred that it is more likely to regulate milk fat metabolism.

Target genes prediction and functional analysis of DELs
The target genes (cis) of DELs predicted by co-localization were enriched to 258 signi cant GO items (P < 0.05), as shown in Fig. 6a. In the signi cantly enriched GO items, the BP is related to lipid metabolism that regulates the secretion of arachidonic acid, and MF is lipid phosphatase activity. KEGG enrichment analysis showed that 7 signal pathways were signi cantly enriched (Fig. 6b), among which Fc gamma Rmediated phagocytosis was related to lipid metabolism. Through enrichment analysis, 7 candidate target genes that may regulate milk fat metabolism were screened, namely VARS2, ITGA6, ATP4A, PPAP2C, FGF1, AMPH, and SYK. The enrichment of GO/KEGG is shown in Supplementary Table 3.
The co-expression predicted DELs target genes (trans) were signi cantly enriched to 509 GO items (P < 0.05), as shown in Fig. 7a. Among them, BP related to lipid metabolism includes lipid metabolic process, long-chain fatty-acyl-CoA biosynthetic process, fatty acid extension, fatty acid transport and negative regulation of triglyceride biosynthetic process, and so on. MF items includes fatty acid synthase activity, fatty acid transporter activity, long-chain fatty acid-binding, and lysophospholipid transporter activity, and so on. The target genes were signi cantly enriched to 17 signal pathways (P < 0.05) (Fig. 7b), among which the pathways related to lipid metabolism include PPAR, cAMP, Fatty acid metabolism and Fc epsilon RI. Under the condition that the correlation coe cient between DELs and target gene expression was greater than 0.97, 10 target genes that might regulate milk fat metabolism were screened out, namely PLA2G4E, INPP4B, FABP3, ACADSB, FABP4, OXSM, FABP7, GPX1, CYP27A1, and ALOX12 genes. The GO/KEGG enrichment is shown in Supplementary Table 4. Among the 10 lipid metabolism-related target genes, the expression of PLA2G4E, FABP4, FABP3, and ACADSB genes was upregulated (P < 0.05), while the expression of FABP7 and CYP27A1 genes was downregulated (P < 0.05). The upregulation of PLA2G4E, FABP4 and OXSM genes is signi cantly enriched in lipid metabolism-related BP, MF and KEGG. Therefore, these three genes and their corresponding DELs need more attention in the follow-up functional mechanism research. The 17 lipid metabolism-related target genes were predicted and screened by co-localization and coexpression. All DELs that regulate these target genes are shown in Table 2. A total of 18 DELs that may be related to milk fat metabolism were screened and their corresponding relationships with target genes are shown in Supplementary Tables 5 and 6. Among the 18 DELs that may be related to milk fat metabolism, 12 were upregulated and 6 were downregulated. Bold font is the DELs predicted by both colocalization and co-expression that may regulate milk fat metabolism. Italic is the DELs corresponding to the target genes (PLA2G4E, FABP4, and OXSM), which is signi cantly enriched in BP, MF, and KEGG. These are also the lncRNAs that need more attention in further experiments.
The interaction between 18 DELs and lipid metabolism-related target genes is shown in Fig. 8. A lncRNA may interact with multiple lipid metabolism-related genes, and a gene may also be regulated by multiple lncRNAs. Among the target genes related to lipid metabolism, the expression trend of signi cantly upregulated or downregulated genes is consistent with the targeted regulation of DELs, which is not only in line with the general mode of action of lncRNAs on genes but also the focus of subsequent experiments.

Construction of lncRNA_miRNA_mRNA interaction regulation network
The differential miRNAs used in this section were derived from the small RNA sequencing results of the research team. LncRNAs can act as precursor molecules of miRNAs. When looking for lncRNAs that interact with miRNAs, it is necessary to lter out lncRNAs that may be precursors of miRNAs. Based on the homology between lncRNAs and miRNAs precursors, we use Blastn software to search for such lncRNAs, and the results are listed in Supplementary Table 7. Then, the lncRNAs corresponding to differential miRNAs were predicted by miRanda software, and miRNAs regulated by the DELs are screened. The target genes of differential miRNAs were predicted and cross-differentiated by using two software, miRanda, and RNAhybrid. Based on ceRNA theory, We searched the lncRNAs_target gene pairs with the same miRNAs binding sites, construction of lncRNA_miRNA _mRNA regulatory relationship with lncRNAs as a decoy, miRNAs as core and mRNAs as the target, and calculation of Pearson correlation coe cients for the expression of lncRNA_miRNA and miRNA_mRNA (Fig. 9).
As shown in the regulatory network diagram, all lipid metabolism-related target genes SYK and PPAP2C (cis) are also the target genes of lipid metabolism DELs, constituting the regulatory relationships of TCONS_00054231_bta-miR-455-3p_SYK, TCONS_00184475_bta-miR-455-3p_SYK, and TCONS_00027906_bta-miR-331-3p_PPAP2C. The upregulated genes DIS3L2, SRRM2, NDRG1 and the downregulated genes ANO1 and MTM1 are the DEGs screened in this study, and 18 ceRNAs regulatory networks contain these DEGs. Also, lipid metabolism-related DELs are involved in a total of 76 ceRNAs regulatory networks, including 10 lncRNAs, 8 miRNAs, and 32 genes. There are 80 ceRNAs regulatory networks involved in lipid metabolism-unrelated DELs, consisting of 8 lncRNAs, 6 miRNAs, and 24 genes, are enriched for target genes GO/KEGG as shown in Supplementary Table 8. Downregulated bta-miR-2387 is a differential miRNAs with the largest number of potential lipid metabolism-related target genes (13) and only interacts with lipid metabolism-related differential lncRNA TCONS_00172817. Non-lipid metabolism-related differences lncRNA TCONS_00004679, TCONS_00091027, and TCONS_00083670 may all play a regulatory role in bta-miR-2387.
It has been shown that the downregulation of bta-miR-145 among many differential miRNAs can regulate the lipid droplet formation and TAG accumulation [15]. Therefore, the regulatory relationship centered on bta-miR-145 is one of the focuses of our subsequent experiments. Notably, TCONS_00027906, TCONS_00054231, TCONS_00082721, TCONS_00045874, and TCONS_00083670 act together on bta-miR-145, which in turn regulates the expression of FASN, which was a potential regulatory network for probing the mechanism of milk fat regulation of the potential regulatory network. The downregulated expression of bta-miR-421, bta-miR-331-3p and bta-miR-197 was a lipid metabolism DELs that speci cally targets lipid metabolism-related differential miRNAs, while bta-miR-222 was a non-lipid metabolism differential lncRNA speci cally targets lipid metabolism-related differential miRNAs.
Interestingly, among the numerous lipid metabolism-related differential miRNAs, only bta-miR-370 is an upregulated differential miRNA, while the rest are downregulated.

Protein interaction network analysis
Protein interaction network analysis of DEGs (11) screened the possible regulation of milk lipid metabolism, target genes predicted by the co-localization and co-expression of DELs (17), and lipid metabolism-related genes corresponding to differential miRNAs targeted by DELs (34), the results are shown in Fig 10. There are 49 nodes and 113 interaction link sets in the protein interaction network diagram, with an enriched P-value of 1.0e-16 and a con dence level of 0.4 (moderate). The different colored bonds represent the source of evidence for protein interactions, and the more bonds involved in biological functions, the more likely they are. The PLTP PPP2R5A SOCS5 APOC3 VARS2 ITGA6 PDGFD BCAT1 ATP8A2 KCNMA1 ZFYVE28 CES4A and CTSH genes are not represented in the gure because they do not interact with each other or with other genes. The interaction among HADHA ACAT1 ACAA1 HADHB ACADS and ACADSB genes, FABP3, FABP4, and FABP7 genes, FASN, OXSM, and MCAT genes are strongly supported by the data with very high intensity as can be seen from the gure. Secondly, SYK and PLCG1, TNF and JUN, PPAP2C and SPHK2, TBL1X and RXRA may also have shared functions between genes. Notably, MAPK3 genes has a strong interaction with several lipid metabolism-related genes (JUN, SMAD2, GH1, SPHK2, FASN, PTPRR, PLA2G4E, PLA2G4F, PLA2G4D, FGF1, and PLCG1), respectively, and may jointly play an important role in regulating milk lipid metabolism. The remaining genes also interacted with each other. Although the data support is relatively small, they still have some reference value.
The expression levels of DELs and DEGs were veri ed by qRT-PCR We randomly selected 8 DELs and 8 DEGs, analyzed their relative expression levels by qRT-PCR in high and low groups (three repeats), and converted the difference multiplea by log 2 FoldChange. As shown in Fig. 11, the results of the qRT-PCR were consistent with those of RNA-seq, which con rmed the accuracy of transcriptome sequencing.

Discussion
Candidate DEGs and DELs in milk fat metabolism It has been con rmed that many genes involve in the lactation initiation, maintenance and growth and development of the mammary gland through direct or indirect regulation [16][17][18]. In this study, a total of 11 candidate DEGs related to lipid metabolism were screened and enriched in PI3K-Akt, MAPK, and Wnt lipid metabolism-related pathways [19,20]. The downregulated ENPP2 gene is signi cantly enriched in BP, MF, and lipid metabolism-related pathways, indicating that it has a greater possibility of regulating milk fat metabolism. Moreover, studies have shown that the lack of ENPP2 has a signi cant protective effect on hepatic steatosis [21]. Dong et al. [22] found that the expression of PDGFD in adipose tissue of obese individuals was higher than that of lean individuals. In this study, PDGFD was also highly expressed in BMECs with high MFP meanwhile, the downregulated gene KCNMA1 also showed a downregulated trend in white adipose tissue and hypertrophied adipocytes of mice on the high-fat diet [23], which con rmed the reliability of sequencing, suggesting that PDGFD and KCNMA1 genes may play an important role in milk fat metabolism. BCAT1 is a branched-chain aminotransferase that has been shown to mainly affect MFP, milk protein percentage, and milk yield of dairy cows [24]. CES4A has also been shown to be mainly involved in lipid metabolism [25]. However, how BCAT1 and CES4A regulate milk fat synthesis has not yet been reported This will be the focus of our future research. It is worth noting that the latest literature reportes that the differentially downregulated DKK1 related to lipid metabolism seems to play a signi cant role in promoting obesity in mice [26], which is inconsistent with the results of this study. It may be related to different species or tissues in the study. Other lipid metabolism-related DEGs APOL1, ATP8A2, PTPRR, ZFYVE28, and CTSH have been reported to associate with kidney disease [27], nerve disease [28], ovarian cancer [29], glomerular ltration barrier [30] and atherosclerosis [31].
LncRNAs can regulate gene expression within 100 kb (cis), and can also change the expression of distant mRNAs (trans) [32]. The functional enrichment analysis of 38 target genes of DELs showed that 17 target genes (10 co-expression and 7 co-localization) were enriched in lipid metabolism-related pathways.
Among them, FABP3, FABP4, and FABP7 are trans-target genes. Studies have con rmed that these genes encode highly abundant fatty acid-binding proteins in bovine mammary glands and transport long-chain fatty acids to the endoplasmic reticulum to synthesize TAG [33]. Protein interaction network analysis showed that the interactions between FABP3, FABP4, and FABP7 were also strongly supported by data, which further indicated that they may have similar functions and participate in common BP.
Based on bioinformatics analysis, a total of 18 candidate DELs that may affect milk fat synthesis in dairy cows were screened out. A lncRNA may target multiple genes related to lipid metabolism, and a gene may be regulated by multiple lncRNAs [11]. These results indicated that lncRNAs were diverse in regulating the function of target genes. In addition, some of the lipid metabolism-related target genes are also signi cantly upregulated or downregulated. The target genes expression trends are consistent with the DELs that may be targeted and regulated, which is consistent with the general mode of action of lncRNAs on genes [34], and further re ects the accuracy of the sequencing results. It is worth noting that FABP4 is a gene encoding a highly rich fatty acid-binding proteins in the mammary gland, which may have a targeting relationship with TCONS _00082721 and TCONS _00172817 [22]. Among all lncRNAs related to lipid metabolism, TCONS_00082721, TCONS_00119434, TCONS_00191498, TCONS_00007612, TCONS_00118412, and TCONS_00163391 are long intervening/intergenic noncoding RNAs (lincRNAs).
Because lincRNAs do not overlap with exon sequences in organisms and show a higher degree of tissue speci city [35] and disease speci city [36]. They are involved in chromatin modi cation, epigenetic regulation, genomic imprinting, transcriptional control, and mRNA processing before and after translation [37]. Therefore, lincRNAs related to milk fat metabolism will be the focus of our subsequent studies. The lncRNAs that might regulate milk fat metabolism are predicted by co-localization and co-expression to be TCONS_00104103, TCONS_00162586, TCONS_00172817, TCONS_00082721, and TCONS_00143117.
Interestingly, the lincRNA TCONS_00082721 with FABP4 as the candidate target gene is also a DELs that regulates milk fat metabolism predicted by the co-localization and co-expression. To clarify this hypothesis, it needs to be further veri ed in functional experiments. These ndings will be of great signi cance for further research on the mechanism of milk fat synthesis in dairy cows.

Protein interaction and LncRNA_miRNA_mRNA regulatory network analysis
As an important regulatory factor, miRNAs can change the expression of target mRNA and involve widely in various biological processes. In recent years, more and more studies have shown that lncRNAs have the function of sponge adsorption of miRNAs and can shield the inhibition or degradation of target genes by competitive binding of miRNAs. And rapidly upregulate its expression [38,39]. In this study, we rst predicted the targeted differential miRNAs of DELs and then predicted the target genes of differential miRNAs. Thirty-four target genes related to lipid metabolism of differential miRNAs with DELs targeting relationship were screened out by GO/KEGG enrichment analysis.
Protein interaction network analysis showed that the target genes HADHA, HADHB, ACAA1, ACADS, ACADSB, and ACAT1 had a strong interaction relationship and might have common functions. HADHB is involved in the β-oxidation of fatty acids and catalyzes the second to fourth steps of the β-oxidation of fatty acid. The other ve genes are acyltransferases or dehydrogenases, which play important regulatory roles in the β-oxidation of fatty acids [40]. Also, the interaction between FASN, OXSM, and MCAT genes is strongly supported by the data. FASN is an enzyme involved in fatty acid (FA) synthesis and plays a signi cant role in the formation of new adipogenesis in mammals. It has been found that the variation of FASN is related to a variety of dairy traits, including milk fat content, total milk solid content and peak milk yield, etc [41] Shi et al. [42] found that FASN variation can be used as a genetic marker to improve milk fat content and milk solid level of yaks. The co-overexpression of SCAP + SREBP1 will also lead to the increase in the abundance of FASN gene mRNA and the formation of lipid droplets [43]. In short, the current study on the FASN gene in the milk fat synthesis has been relatively mature at present. OXSM is a 3-oxyacyl -[acyl carrier protein] synthase. Xiong et al. [44] found that OXSM has inhibitory effects on fat catabolism, fat anabolism, and fatty acid oxidation, and promotes glycerol transport and polyunsaturated fatty acid synthesis. MCAT is also an important gene involved in lipid metabolism [45]. Therefore, there is an inevitable interaction between FASN, OXSM, and MCAT genes, and they all play a potential role together in the regulation of lipid metabolism. MAPK3 plays a major role in autophagy, and there is evidence that MAPK3 is involved in lipid metabolism [46]. We also found that the MAPK3 gene has a strong interaction with multiple lipid metabolism-related genes (JUN, SMAD2, GH1, SPHK2, FASN, PTPRR, PLA2G4E, PLA2G4F, PLA2G4D, FGF1, and PLCG1). It is speculated that this gene may play a crucial role in the lipid metabolism process, and the detailed mechanism needs further study.
Through functional analysis, a total of 156 milk fat metabolism-related ceRNAs were constructed in this study. Among the miRNAs targeted by DELs, miR-145 has been shown to regulate lipid drop-formation and TAG accumulation in goat mammary epithelial cells [15]. The regulatory network centered on miR-145 will be the focus of subsequent experimental research. Since TCONS_00082721 is a differential lincRNA predicted by co-localization and co-expression that may regulate milk fat metabolism, the competitive regulatory relationship between TCONS_00082721_miR-145_FASN and TCONS_00082721_miR-145_TNF will be further studied and cell experiments will be carried out to verify the interaction and function between them. The target gene predicted by TCONS_00172817 is FABP4 [22], which encodes a highly abundant fatty acid-binding protein in the mammary glands of dairy cows. These lncRNAs may regulate milk fat metabolism predicted by the co-localization and co-expression. Therefore, the ceRNAs regulatory network related to TCONS_00172817 is likely to be involved in the process of milk fat synthesis and needs attention. For non-lipid metabolism DELs,TCONS_00045874 and TCONS_00083670 may play a role in miR-145, regulating the expression of lipid metabolism-related genes FASN [43] and TNF [47]. Among all the DELs related to lipid metabolism targeted differential miRNAs, bta-miR-2387 has the most abundant predicted target genes, which may be an important gene in the process of milk fat synthesis in dairy cows. Therefore, the ceRNAs regulatory network with bta-miR-2387 as the core will have more choices in functional veri cation of its interaction and function, and the possibility of success will also be greater. To sum up, the ceRNAs network analysis was performed on the transcriptome data (mRNA, miRNAs, lncRNAs) obtained from BMECs with high and low MFP. It provides a new way to understand the function of genes in the biological process of milk fat metabolism, and it also provides a new perspective for studying the lactation process of dairy cows.

Conclusion
In this study, we isolated BMECs of Holstein cows with extreme differences MFP from fresh milk and measured lncRNAs and mRNA transcriptome pro les. Through functional enrichment analysis, 11 DEGs and 18 DELs related to milk fat metabolism were screened out, and 156 lncRNA_miRNA_mRNA interaction regulation networks related to milk fat metabolism were constructed. The results of this study provide a new way to understand the gene function during milk fat synthesis and provide many new insights into the structure of molecular pathways during lactation.

Ethics statement
In this study, BMECs were isolated and cultured from milk samples of dairy cows. 200 mL milk samples were collected from each dairy cow, and aseptic operation was strictly carried out during the sampling process, so there was no harm to dairy cows.
Experimental cow screening and milk sample collection 245 Holstein cows were selected from Ningxia Nongkeng Heilanshan Maosheng dairy farm with the same breeding background, similar lactation (150-220 days) and age. Early, middle and late milk samples of each cow were collected for dairy herd improvement (DHI) determination. Holstein cows with high MFP (n = 4) and low MFP (n = 4) with somatic cell count (SCC) less than 100000 / mL were selected (Table 3). Fresh milk samples were aseptically collected into 50 mL centrifuge tubes, and four tubes per cow were collected for a total of 200 mL. Tighten the cap of the bottle, seal and sterilize it, put it in a thermos ask containing 37 ℃ sterile water, and return it to the laboratory for later use.

Isolation and identi cation of primary BMECs
The composition of the complete culture medium and culture method were the same as those of Duan et al. [48]. When the cells grew to 80% (F 2 generation), total cellular RNA (TaKaRa; D9108A) was extracted.
To exclude genomic DNA contamination, DNAase (TaKaRa; 00079459) was used for reverse transcriptional preprocessing (TaKaRa; BK3702). Primers were designed using Primer Premier 5.0 and PCR ampli cation was performed for genes speci cally expressed in BMECs related to milk fat synthesis and keratin 8 (TaKaRa; DR100A). PCR reaction procedure: pre-denaturation at 94°C for 4 min, denaturation at 94 °C for 30 s, annealing (Supplementary Table 10  The basic identi cation of BMECs is the same as that of Shen et al. [14]. The triglyceride (TAG) content of BMECs was determined by using the triglyceride enzymatic assay kit (E1025) of the high-fat samples. Then the protein concentration of the sample was measured by using the BCA protein content detection kit, and the TAG content was adjusted for the protein concentration per mg.
Library construction and transcriptome sequencing The Trizo method was used to extract total RNA from BMECs (F2) of high and low MFP cows (4 heads), and the quality of RNA samples was strictly controlled. The integrity of RNA and DNA contamination were analyzed by agarose gel electrophoresis. The concentration and purity of RNA (OD260/280) were detected by Nanodrop, and the integrity of RNA was accurately detected by Agilent 2100 bioanalyzer. The 260/280 ratio of all samples was approximately 2.0, and the RNA integrity index (RIN) was ≥8.0. After determining the RNA purity and quality, ribosomal RNA was removed to construct strand-speci c library.
After passing the library inspection, Illumina PE150 sequencing was performed. After the original data is obtained, the reads with adapter, N (undetermined base information) ratio greater than 0.002, and lowquality bases with a read length of more than 50% are removed. After sequencing error rates (Q20 and Q30) and GC content distribution checks, clean reads for subsequent analysis were obtained. tophat2 (http://tophat.cbcb.umd.edu), Hisat2 (http://ccb.jhu.edu/software/hisat2), and STAR (http://code.google.com/p/rna-star) software were used to compare and analyze the RNA sequencing(RNA-seq) data.

DEGs, DELs screening, and functional enrichment analysis
The level of gene expression was quanti ed using fragments per kilobase of exon per millions of reads (FPKM) value. Expression difference signi cance analysis was performed using edgeR software. The corrected P-value (Padj) was used to determine the signi cance level, and Padj < 0.05 was used as the standard of differential signi cance, while |log2FoldChange|>1.5 was another important criterion for screening DELs and DEGs. Co-localization of lncRNAs with protein-coding genes was used to nd genes within 100 kb upstream and downstream of lncRNAs, and co-expression was used to predict the target genes of lncRNAs. The gene ontology (GO) database was used to enrich cellular composition (CC), biological process (BP), and molecular function (MF) of DEGs and target genes of DELs. Kyoto encyclopedia of genes and genomes (KEGG) was used to identify the major biochemical metabolic pathways and signal transduction pathways involved in genes. The protein interaction network between DEGs and DELs target genes was analyzed by STRING database (https://string-db.org/) and further visualized using Cytoscape (http://www.cytos-cape.org/). MiRNA target prediction miRanda database was used with the default parameters to identify conserved miRNA target sites in the 3'UTR of lncRNAs. MiRNAs target gene prediction was based on the intersection of miRanda (Total score ≥ 140, Total energy ≤ -15 Kcal/mol) and RNAhybrid (mfe ≤ -20 Kcal/mol, P < 0.05).
qRT-PCR validation of sequencing data To con rm the sequencing results, eight DEGs and DELs were randomly selected for qRT-PCR validation, respectively. First-strand cDNA synthesis was performed using the PrimeScript RT kit (Takara, Dalian, China). qRT-PCR (three replicates) was performed by SYBR Premix Ex Taq™II (TaKaRa) on the Bio-Rad CFX96 Touch™ Real-Time PCR Detection System (Bio-Rad, Hercules, CA, USA). Ampli cation procedure: 95°C for 30 s, 95°C for 5 s, annealing for 30 s, 40 cycles. qRT-PCR primers were designed using Primer Premier 5.0 with primers spanning at least one intron, and the sequence and annealing temperature of each primer are shown in Supplementary Table 10

Consent for publication
Not applicable.
Availability of data and materials statement All data generated or analyzed in this study are included in this article [and its supplementary information le], and the datasets have been submitted to the SRA database with the accession number PRJNA730595. Access to the data of permanent link to https://www.ncbi.nlm.nih.gov/sra/PRJNA730595.

Competing interests
The authors declare that we have no competing interests.
Financial support statement  Tables   Table 1 Candidate DEGs related to milk fat metabolism. Log2FoldChange is logarithm of fold change with a base of 2, Padj is the corrected signi cance test probability value, gene location is the speci c location of the gene on the chromosome. Table 2 DELs related to lipid metabolism predicted by co-expression and co-localization. Italic is the DELs corresponding to the target genes (PLA2G4E, FABP4, and OXSM), which is signi cantly enriched in BP, MF, and KEGG. Bold font is the DELs predicted by both co-localization and co-expression that may regulate milk fat metabolism. transcript location is the speci c location of the DELs on the chromosome. Padj is the corrected signi cance test probability value.