Cattle domestication has led to the development of various breeds with unique genomic architecture [1], which is phenotypically expressed not only through physical appearance and color [2] variations but also through valuable traits such as milk productivity, and especially milk composition [3]. The ratio of milk components can change because of certain factors, such as diet and season of the year, but these changes differ among breeds [4]. The stage of lactation is another criterion that determines the metabolite composition of milk. Milk is produced by epithelial cells of the mammary gland, which under the influence of both exogenous and endogenous factors undergo numerous physiological changes throughout lactation. Ultimately, these changes determine the intensity of milk component synthesis and, hence, milk composition [5]. Functional genomics studies based on gene network analysis revealed the complexity of molecular adaptation of the mammary gland to lactation. This is assumed to be determined by changes in the transcriptome of mammary epithelial cells during lactogenesis and galactopoiesis [6], whereas the expression level of key candidate genes may vary in animals of different breeds [7]. The emergence of new technologies (DNA microarrays, genome and transcriptome sequencing, etc.) provides a better understanding of mechanisms controlling the regulation of milk production. The use of genomic information has increased selection efficiency in cattle, although the variability of inherited traits is not fully understood [8]. Understanding gene regulation and interactions at the post-transcriptional level can be a useful tool for identifying novel epigenetic markers of productive traits. The complex interaction of genetic and environmental factors determines epigenomic changes, which significantly affect trait expression level [9]. MiRNAs are essential epigenetic components that participate in fundamental processes, including proliferation, embryonic development, tissue differentiation, and apoptosis, and influence lipogenesis, hematopoiesis, and immunity. As a part of the non-coding RNA class, they can regulate up to 60% of gene expression at the post-transcriptional level by binding to complementary RNA molecules, which leads to translation repression or mRNA degradation and, thus, to changes in cellular protein level in different tissue cells [10]. MiRNAs are secreted by body cells and found in all body fluids as part of stable protein or lipid complexes [11]. Being expressed by mammary epithelial cells, they are involved in intracellular communication and signaling pathways at the cellular level, thereby determining the functioning of the mammary gland itself and ultimately the nutrient ratio in milk [12]. Milk as a non-invasive source of miRNAs is an excellent target for studying the mammary transcriptome. Milk fat-derived miRNAs were reported to accurately map the miRNAome of breast tissue [13]. In a recent study [14], the stability of human breast milk miRNAs was experimentally confirmed by treating milk samples with RNase, a low-pH solution, and a triple freeze-thaw cycle (−20°C). Several studies showed that the expression level of miRNA in cow milk depends on housing and environmental conditions, feed ration [15,16], age and breed [17], and the physiological condition of the udder, including mastitis [18]. According to other authors, miRNA expression differs not only in lactating and dry cows but also in groups of cows with high and low milk fat and protein contents [19].

It is well known that Holstein cattle have a high potential for milk productivity, whereas high milk fat content is a distinctive feature of Ayrshire cattle. Modern populations of Ayrshire and Holstein cattle are distinguished by their unique genomic architecture, formed as a result of long-term breeding and artificial selection [20,21]. In this regard, it is relevant to study in detail the contribution of some key milk miRNAs as epigenetic regulators of lactopoiesis and galactopoiesis. MiRNAs including miR-106b (BTA 25; MI0009724), miR-191 (BTA 22; MI0005034), and miR-30d (BTA 14; MI0004747), are presumably involved in the processes regulating the synthesis of protein-fat components of milk, as shown in other studies [22,23]. Overexpression of miR-106b in the mammary gland tissues of Holstein cows resulted in downregulation of the CDKN1A gene and alteration of protein synthesis pathways [24]. Although miR-191 is one of the most prevalent miRNAs in bovine mammary gland tissues, its function in lipid and protein synthesis remains unclear. However, its human homologue regulates transcription factors, chromatin remodelers, and cell cycle genes involved in proliferation, apoptosis, differentiation, and migration [25]. Furthermore, miR-191 was identified as a diagnostic marker for breast cancer in women [26]. MiR-30d is a universal miRNA that interacts with many target genes and performs various biological roles. For instance, miR-30d indirectly participates in blood glucose level regulation by activating insulin transcription [27].

It is essential to study the expression patterns of miR-106b, miR-191, and miR-30d in cow milk throughout lactation. To accomplish this aim, we conducted a comprehensive analysis of the aforementioned miRNA expression in the milk of Holstein and Ayrshire cows at different stages of lactation, considering milk composition.

Effective breeding of dairy cattle requires considering the peculiarities of milk composition in different breeds. Milk composition, according to Schwendel et al [31], depends on a variety of paratypical and genetic factors, so the determination of paratypical factors is of great practical importance for industrial production. In addition, breed and housing system (including feed ration) constitute the minimum set of factors to be considered when comparing milk samples. Therefore, the statistical analysis takes into account a complex factor that includes both breed and housing system. Breed can be a decisive factor that largely accounts for differences between housing systems [32]. Although we obtained differences between cow groups in protein and casein content, the “farm” factor did not significantly affect the content of these components in milk. This result, apparently, can be explained by the fact that both farms provided conditions sufficient for the same manifestation of the studied traits, as well as by the fact that cows were at the same stage of lactation and received a balanced mono fodder in accordance with their physiological status. As our data show (Table 3), the content of protein-fat constituents of Holstein cows’ milk changed significantly during 10 months of lactation, and only the level of PUFA and TFA remained relatively stable. In Ayrshire cows, only the values of milk protein, casein, C14:0, C18:0, and TFA content significantly changed (at least at p<0.05). In general, both breeds showed an increase in the content of all investigated milk constituents until late lactation (month 7 to 8), which is probably due to physiological changes caused by late lactation and decreased milk yield. ANOVA analysis revealed that the factor “month of lactation” had no significant effect on C18:1, LCFA, PUFA, and SCFA content in milk of both breeds during the first lactation (Table 2). The combined effect of the two factors “farm+month of lactation” was significant for all analyzed milk constituents in both breeds.

The milk composition characteristics of Holstein and Ayrshire cattle may be defined by differences in their genomic architecture [21]. This factor determines the presence of different sets of candidate genes, which are phenotypically expressed through milk traits [33]. Higher fat, protein and casein content in milk, compared to Holstein cows, is a specific breed trait of Ayrshire cattle (p<0.01 to 0.001) [34]. Similar interbreed differences were observed in our study, where milk from Ayrshire cows contained on average 1.82% more fat (p<0.001), 0.06% more protein and 0.08% more casein (p<0.05). Significant differences in milk fat content appear to cause differences in fatty acid content as well (p<0.001). PCC analysis for both breeds showed clustering of fatty acids with high correlation coefficients (p<0.05 to 0.001). For the Ayrshire breed, PCC analysis demonstrated clustering of protein and casein traits relative to fat and fatty acids, with no significant correlations between these traits. For Holstein cows, a unidirectional correlation between protein, casein, and fat components of milk was established. The identified differences in the ratio of protein-fat components in these breeds presumably determine the technological properties of milk [34].

The results obtained in this research, based on the real-time stem-loop RT-qPCR data, revealed different expression patterns of miR-106b, miR-191, and miR-30d in milk samples during lactation in Ayrshire and Holstein cows. Against the background of considerable fluctuations in protein-fat components (p<0.05 to 0.001), except for PUFA, the expression levels of miR-106b, miR-191, and miR-30d in milk samples of Holstein cows were reduced in the first two months, whereas in Ayrshire cows they were reduced by the fourth month. In contrast, the milk composition of the latter was relatively stable throughout lactation, except for protein, casein, C14:0, C18:0, and TFA (at least p<0.05).

Correlation analysis demonstrated that both breeds exhibited positive relationships between the expression levels of miR-106b, miR-191, and miR-30d with each other, while displaying negative correlations with several fatty acids, which is consistent with the generally accepted statement that miRNA is a negative post-translational regulator. Moreover, among the components characterizing milk fat composition, only C18:1, C18:0, and LCFA were negatively correlated with each of the three miRNAs in milk from Holstein cows. The Ayrshire breed also exhibited negative correlations between milk fat components and miRNAs, including miR-106b and miR-30d with C18:0, and miR-191 with TFA, which was absent in milk samples from Holstein cows. The Ayrshire breed specificity was the positive correlation of miR-106b and miR-30d with milk protein components. According to the literature, miRNAs do not always act as negative post-translational regulators. Previously, while studying the functional role of 11 miRNAs in the processes of milk fat synthesis regulation in epithelial cells of the goat mammary gland, it was found that increased expression of three miRNAs (miR-23a, miR-103, and miR-200a) led to activation of transcription and translation of their target genes [35].

The study of miRNAs is important not only for understanding galactopoiesis, but also for human health. Fatty acids (FA) associated with the expression of the studied miRNAs are important for healthy human diet. Stearic acid (C18:0) exhibits a neutral or beneficial effect on blood cholesterol levels. It can also be desaturated to oleic acid that, in turn, is associated with lower risk of cardiovascular disease, diabetes, and obesity, as it can improve the lipid profile, insulin sensitivity, and inflammatory markers [36].

SCFA can influence the gut-brain axis, which is the bidirectional communication between the gastrointestinal tract and the central nervous system, by affecting the production and release of neurotransmitters, hormones, and cytokines. SCFA can also modulate the energy balance, glucose homeostasis, lipid metabolism, and immune function by activating specific receptors, such as G-protein coupled receptors (GPR41, GPR43, and GPR109A) and free fatty acid receptors (FFAR2 and FFAR3), or by inhibiting histone deacetylases [37].

TFA have been associated with increased risk of cardiovascular disease, diabetes, obesity, and cancer, as they can adversely affect the lipid profile, insulin resistance, inflammation, and oxidative stress [38]. However, the effects of TFA may vary depending on their type and dose. Industrial TFA (iTFA), on the other hand, have been shown to have more detrimental effects on health, compared to ruminant TFA (rTFA), as they can increase the levels of LDL cholesterol, triglycerides, and lipoprotein(a), decrease the levels of HDL cholesterol, and impair the endothelial function [39].

SFA can also influence the gut-brain axis, energy balance, glucose homeostasis, lipid metabolism, and the immune function, depending on their chain length and concentration [40].

Cow’s milk protein is a complex mixture of bioactive peptides and proteins, which can modulate various physiological systems in humans, such as the immune, cardiovascular, gastrointestinal, and nervous systems [41]. Casein, which constitutes about 80% of the total protein content of cow’s milk, influences the physical properties of milk and dairy products, such as viscosity, stability, and texture, which are relevant for processing and storage of dairy products [42].

Therefore, predicting the component composition of cow’s milk is crucial for enhancing the quality of human nutrition. ROC curves data suggest that miRNAs can be used as a prognostic marker in the early prediction of cow’s milk composition. Thus, in our study, miR-106b had high predictive values for C18:1, SFA, and SCFA (AUC = 1) in Holstein cows, and for protein, casein, and C18:0 (AUC = 1) in Ayrshire cows. In contrast, miR-30d had high predictive values (AUC = 1) for protein, casein and C18:0 only in Ayrshire cows. In the presence of some reliable correlations in both cow groups, miR-191 showed low predictive value, and therefore cannot be recommended as a predictive marker of protein-fat composition of milk.

The miRNAs’ biological features, specifically their ability to target dozens of genes, suggest their involvement in many signaling pathways. Understanding the mechanisms of genetic and epigenetic regulation of processes affecting the synthesis of milk components may improve milk productivity in cows. More than 6,000 genes regulate milk synthesis processes, and expression of these genes is observed in both mammary gland tissues and other types of tissues [43].

In our study, 1,531 target genes involved in nearly 200 reactome pathways were identified for miR-106b, miR-191, and miR-30d (Table 5). The target gene network plot identified a cluster of 5 genes that were targets for all miRNAs analyzed in the study, i.e. TAOK1, FSD1L, PPARGC1B, PPP1R16B, and CLMN.

The TAOK1 gene, serine/threonine-protein kinase TAO1, is a MAP3K protein kinase that, through regulation of the mitogen-activated protein kinase pathway, modulates a significant number of cellular processes. TAOK1 was previously identified as one of the functional candidate genes for immunoglobulin G (IgG) and IgM immunoglobulin content in the colostrum and serum of Holstein cows [44], confirming its involvement in the synthesis of milk components.

The FSD1L gene, fibronectin type III and SPRY domain-containing protein 1, encodes type 2 cystatins, members of the cystatin family of intracellular and extracellular protease inhibitors. The function of this gene is poorly understood; however, in a study by Tahir et al [45], a single nucleotide polymorphism (SNP) located near the genomic region comprising the FSD1L gene (BTA8) was associated with heifer fertility traits, which are known to be negatively correlated with cow milk production traits.

The PPARGC1B gene, peroxisome proliferator-activated receptor gamma coactivator 1-beta isoform X2, is known to be involved in the Insulin resistance biological pathway (bta04931), as well as to regulate glucose homeostasis and mitochondrial biogenesis. Its participation in the regulation of milk fat synthesis was previously confirmed [46]. Earlier in a GWAS study, the PPARGC1B gene was identified as significant for milk production traits in Sahiwal-Tharparkar cows [47]. In another study, addition of high oleic sunflower seeds to the cows’ ration was associated with an increase in the proportion of C18:0 and C18:1 in milk and an increase in PPARGC1B transcriptional activity in the mammary gland [48].

The PPP1R16B gene, protein phosphatase 1 regulatory inhibitor subunit 16B, is involved in biological processes such as endothelial barrier establishment and positive regulation of blood vessel endothelial cell proliferation involved in sprouting angiogenesis. According to Pszczola et al [49], this gene was associated with cows’ methane production level, which, according to the authors, is caused by the presumed impact of this gene on digestion and nutrient assimilation processes. At the same time, when taking phenotypic traits into account, the researchers also made adjustments to the fat and protein content of milk. Earlier studies have linked fatty acid levels in milk to methane production in cows; in particular C18:1 level had a negative correlation with CH₄ [50].

The CLMN gene, calmin isoform X1, has been identified as one of the genes sensitive to the vitamin A metabolite all-trans retinoic acid (atRA), which regulates the development of the nervous system during embryonic development [51]. Although the role of CLMN gene in regulating the synthesis of protein-fat constituents of milk is not clear, there is evidence supporting its involvement in lipid metabolism. GWAS analysis revealed an association of SNP in the CLMN gene with changes in total cholesterol and plasma LDL-cholesterol levels in patients taking statins [52]. In cattle studies, the CLMN gene haplotype was associated with meat marbling score in Korean beef cattle Hanwoo [53].

Some target genes detected by reactome pathways analysis (Table 5) were associated with certain milk constituents. Polymorphic variants of GPAM gene were associated with milk fat, protein and dry matter content in Holstein cows’ milk [54]. A method combining GWAS with a gene-centric approach revealed that the PYCR1 gene is included in genomic regions associated with fatty acid profile and milk fat content in cow milk [55].

Our data provide new insights into the relationship between miRNA expression and milk composition and suggest that miR-106b, miR-191, and miR-30d play a potentially significant role in galactopoiesis in both Holstein and Ayrshire cattle. The differences observed in the expression patterns of the analyzed miRNAs may be caused by specific genomic architectures of the studied breeds. This is probably caused by the fact that in a number of candidate genes the phenotypic effect varies both by lactation number and dairy cattle breed.

Since different breeds may have different genetic and epigenetic backgrounds that affect miRNA expression and milk composition, it would be advisable to compare the results with other dairy cattle breeds (Jersey, Danish Red) and combined productivity breeds (Brown Swiss, Simmental). Identification of the target genes and pathways of these miRNAs in mammary tissue and milk would also be useful, as well as confirmation of their action by functional analyses such as knockdown or overexpression experiments. This will help to identify the molecular mechanisms of lactation, as well as the epigenetic mechanisms controlling milk composition and quality in the studied breeds, which may benefit the dairy industry and consumers.

Analysis of miRNA expression levels can also be used to predict milk composition and quality or to select cows with desired milk characteristics. In addition, external factors like diet and environment, can alter miRNA expression levels to enhance milk composition and quality. These applications will require further validation and optimization in future studies.

Our results showed that milk from Holstein and Ayrshire breeds had different protein-fat component ratios. The expression levels of miR-106b, miR-191, and miR-30d in milk samples differed depending on the lactation period and breed. The correlation analysis did not reveal a correlation between changes in miRNA expression level and milk fat content. However, a negative relationship with several milk fatty acids in both breeds was observed. A positive correlation between the expression levels of miR-106b and miR-30d and the total protein and casein content was found in the milk of Ayrshire cows. This finding suggests that different epigenetic mechanisms regulating the synthesis of milk components in mammary gland tissues might exist in different breeds. The calculated AUC values in ROC curves analysis showed that changes in the expression of miR-106b and miR-30d can be a prognostic marker for the milk composition in terms of C18:0, protein and casein content in Ayrshire cows. In Holstein cows, the expression level of miR-106b can also be a prognostic factor, but for fatty acid composition of milk (C18:1, SFA, and SCFA).

The analysis of the relative expression levels of key miRNAs in dairy cow milk contributes to a better understanding of epigenetic factors affecting milk composition and expands our knowledge of lactogenesis and galactopoiesis. The proposed real-time RT-qPCR method using stem-loop primers and the designed TaqMan probe is a universal tool for milk miRNA analysis, which other researchers can use for further characterization of these miRNAs.