Introduction
Milk is a biologically active secretion containing essential nutrients and immune-related components that support neonatal growth and dairy productivity. Beyond its nutritional importance, milk quality is influenced by physiological factors that reflect mammary gland function. Among these, somatic cell count (SCC) is a key indicator of mammary health and subclinical inflammation, making it a central parameter in dairy research (Yakan et al., 2021; Karagözlü et al., 2024).
In cows, somatic cells naturally present in milk are largely composed of mammary epithelial cells and leukocytes (Vitenberga-Verza et al., 2022). An increase in the number of these cells in milk is associated with inflammation and infections, and it is known that high SCC is related to a loss in milk quality (Yakan et al., 2021). SCC is used as a marker for determining mammary health and milk quality, but the threshold values for evaluating SCC in cow's milk differ between developed and developing countries. In Sweden, an SCC above 200,000 cells/mL in milk is considered a sign of disease (Frössling et al., 2017), while in several EU countries and Australia this threshold is 400,000 cells/mL, and in Turkiye, it is 500,000 cells/mL (Hisira et al., 2023; Darbaz et al., 2023).
Although husbandry practices such as breed, lactation period, and milking frequency cause fluctuations in SCC, under normal physiological conditions, SCC is low but increases due to a rise in leukocytes. The status of mammary health can be assessed through direct measurement of SCC and the California Mastitis Test (CMT). However, in subclinical mastitis, SCC measurement and CMT application can limit the accuracy of diagnosis (Abed et al., 2021). Although microbiological analysis is considered the gold standard for diagnosing subclinical mastitis, limitations associated with this method have led to an increased search for molecular biomarkers in recent research, with microRNAs (miRNAs) coming to the forefront (Özkan et al., 2024).
MicroRNAs were first identified in 1993 in studies conducted on C. elegans. These RNA molecules, 18-24 nucleotides in length, are highly conserved across species (Lee et al., 1993). These molecules exhibit high stability in body fluids, play roles in intercellular communication, and can be actively released into the extracellular environment. Their ability to be differentially expressed in abnormal conditions such as disease and infection, or in changing physiological states, makes them potential biomarkers for evaluating various physiological and pathological conditions (O’Brien et al., 2018; Do et al., 2021). miRNAs found in milk are resistant to acidic environments, RNase digestion, temperature, and multiple freeze/thaw cycles, and their stability suggests their potential as diagnostic biomarkers for diagnosing subclinical mammary infections and assessing milk quality (Cai et al., 2018; Kok et al., 2018).
According to miRBase (2025), over 1,000 mature miRNAs have been identified in cattle. Studies have shown that miR-27a-3p is involved in pathways such as inflammation, apoptosis, and cell proliferation, and it regulates lipid metabolism by affecting genes like Peroxisome Proliferator-Activated Receptor Gamma (PPARG) in bovine mammary epithelial cells. Additionally, genes involved in energy metabolism, such as Carnitine Palmitoyltransferase 1B (CPT1B) and Retinoid X Receptor Alpha (RXRα), are also targets of miR-27a-3p (Lin et al., 2013; Luo et al., 2020; Wang et al., 2021). miR-29B-2, which plays significant regulatory roles in cellular processes such as apoptosis, cell cycle, cell aging, cell differentiation, and cell proliferation, has been reported to have potential as a biomarker in mastitis in cattle (Srikok et al., 2020; Horita et al., 2021). In another study, it was found that miR-142-5p levels in milk decreased when SCC was below 200,000 cells/mL (Stefanon et al., 2023). miR-142-5p, is more abundant in milk compared to other miRNAs and is involved in fatty acid metabolism, has also been shown to regulate inflammatory pathways in cattle (Chartoumpekis et al., 2012; Cai et al., 2018; Tzelos et al., 2022). Moreover, miR-148a plays roles in the inflammatory response during infection, adipogenesis, and the regulation of numerous genes during lactation (Okumura et al., 2021; Li et al., 2022).
In this study, the expression patterns of miR-27a-3p, miR-29B-2, miR-142-5p, and miR-148a were investigated in Holstein cows with low (SCC <200,000 cells/mL) and high SCC (SCC >200,000 cells/mL). The relationships between these target miRNAs and SCC, as well as basic milk compositional and quality parameters (protein, lactose, fat, fat-free dry matter (FFDM), electrical conductivity, and freezing point), were also examined. Additionally, identification target genes, protein-protein interactions (PPI) for the relevant miRNAs and pathway analysis were investigated using bioinformatics tools.
Materials and methods
This study utilized animal material consisting of multiparous Holstein cows (parity 2-4) in mid-lactation (120-180 days in milk), housed in a private enterprise in Kayseri province (Türkiye) with over 200 milking cows. Prior to sample collection, the udder lobes were cleaned with 70 % ethyl alcohol and wiped with disposable paper towels. After discarding the first stream of milk, samples were collected under aseptic conditions and subjected to the CMT (Özkan et al., 2024). Milk samples collected before the morning milking were divided into two groups based on SCC: those with SCC below 200,000 cells/mL (n=20) and those with SCC above 200,000 cells/mL (n=20). The group with SCC below 200,000 cells/mL was defined as the low SCC group (LSCC), while the group with SCC above 200,000 cells/mL was defined as the high SCC group (HSCC). From each animal, 100 mL of milk sample (2x50 mL falcon tubes) was collected under sterile conditions. One of the falcon tubes containing the milk samples was used for determining SCC, pH, and basic milk quality and compositional parameters (protein, lactose, fat, fat-free dry matter, electrical conductivity, and freezing point), while the other was used for molecular analyses. All cows were kept under identical management and nutritional conditions, receiving the same standardized total mixed ration (TMR) formulated according to NRC recommendations (NRC, 2001).
Determination of SCC, pH, and compositional parameters of milk samples
All milk samples were analyzed immediately after collection under cold-chain conditions. Prior to measurements, each device was calibrated according to the manufacturer's instructions. The somatic cell counter (Lactoscan SCC 6010, Bulgaria) was used with its self-calibration, which was routinely verified through the device’s internal control system. The pH meter (Hanna HI83141, USA) was calibrated using manufacturer-supplied buffer solutions before each set of measurements. Milk composition parameters were measured using the Milkotester Master Classic M2 P1 (Bulgaria), which was calibrated following the daily and weekly cleaning and calibration routines recommended by the manufacturer (Özkan et al., 2024).
RNA Isolation, polyadenylation, and cDNA synthesis
Samples transported to the laboratory under cold chain conditions were centrifuged at 3000 xg for 10 minutes at +4 ℃, followed by incubation at -20 ℃ for 10 minutes. The resulting cream layer was aseptically removed using a spatula. RNA isolation was performed using 250 µL of the skim milk portion obtained after centrifugation, following a modified Trizol method (Özkan et al., 2024). The obtained RNA pellets were resuspended in 15 µL of nuclease-free water and assessed for purity and concentration using a nucleic acid quantifier (Merinton-SMA 1000, China).
Following isolation, poly(A) tails were added to the miRNAs present in the total RNA using a Poly(A) Polymerase kit (ABM, Canada, Cat. no: E017). For this purpose, 500 ng of RNA was mixed with 5 μL of 5x Poly(A) Polymerase Yeast Reaction Buffer, 2.5 µL of MnCl2 (25 mM), 1.25 μL of ATP (10 mM), 1 μL of Poly(A) Polymerase Yeast (1 U/μL), and 0.3 µM RNA, and the volume was adjusted to 25 µL with nuclease-free water on a cold block. After preparation, the mixture was briefly centrifuged and incubated at 37 ℃ for 15 minutes. The poly(A) reaction was terminated by incubation at 65 °C for 20 minutes.
cDNA synthesis was performed using the OneScript® Plus cDNA Synthesis Kit (ABM, Canada, Cat. no: G236). For this, 10 µL of poly(A)-tailed samples were mixed with 4 µL of 5x RT Buffer, 1 µL of dNTP, 1 µL of oligo(dT) primer, and 1 µL of RTase enzyme. The volume was then adjusted to 20 µL with nuclease-free water. cDNA synthesis was carried out at 50 °C for 15 minutes, followed by incubation at 85 °C for 5 minutes to terminate the reaction. The resulting cDNA was diluted 10-fold and stored at -20 ℃ until qPCR analysis.
qPCR application
Amplification of target miRNAs was performed using qPCR (Rotor Gene Q MDx 5plex HRM, Qiagen, USA) and a SYBR Green I dye-containing kit (Power SYBR Green PCR Master, Thermo Fisher Scientific, USA, Cat No: 4368702). Each sample was run in duplicate. The reaction consisted of an initial denaturation at 95 °C for 10 minutes, followed by 50 cycles of 95 ℃ for 15 seconds, 57.0-63.2 ℃ for 60 seconds, and 72 ℃ for 30 seconds (Table 1). Additionally, melting curve analysis was performed to verify the specificity of the amplified regions. The temperature increased from 62 °C to 99 °C in 1 °C increments, with fluorescence recorded at each step to determine the melting temperature (Tm) and confirm product specificity. Because of U6 snRNA has commonly been preferred as a reference in recent milk and milk related miRNA studies, U6 was used as the internal control for normalization (Gao et al., 2024; Ma et al., 2024; Ünal et al., 2025). Primer design was based on the miRbase database
(www.mirbase.org).
Table 1. Primer sequences of the studied miRNAs and snoU6
AT: Annealing temperature
Identification of miRNAs target genes
The target genes of the identified homologous miRNAs were analyzed using the miRNet (https://www.mirnet.ca/) and miRTarBase databases (Chang and Xia, 2022). The targeted genes were visualized using the miRNet database. The identified miRNA target genes were further analyzed and visualized using a Venn diagram generated through the Venny tool (v.2.1, https://bioinfogp.cnb.csic.es/tools/venny/index.html) available on the BioinfoGP platform (Oliveros, 2007).
Protein-protein interaction network analysis
Gene interactions were visualized using the STRING database (v.12.0, https://string-db.org/), including proteins targeted by at least two miRNAs. Among these, only proteins with an interaction score ≥ 0.150 were considered, while proteins with more than five interactions were excluded from visualization (Lan et al., 2022). The analysis was conducted by considering multiple parameters, including experimental data, database annotations, co-expression, neighborhood relationships, gene fusions, and co-occurrence (Szklarczyk et al., 2023). The data obtained were visualized using the Cytoscape software (v.3.10.3), and proteins that did not exhibit interactions were excluded from further analysis. As part of the PPI network analysis, the MCODE algorithm was employed to identify the most highly interconnected gene clusters (Bader and Hogue, 2003). The MCODE analysis was performed using default parameters (Degree Cutoff: 2, Node Score Cutoff: 0.2, K-Core:2, Max. Depth:100).
Functional and pathway enrichment analysis
Functional and pathway analyses were conducted to annotate the predicted target genes from the PPI network analysis. This analysis was performed using genes obtained from PPI interactions through the STRING database. In this context, proteins exhibiting interactions were assessed through the Gene Ontology (GO) database and categorized into biological processes, molecular function, and cellular component for functional enrichment analysis. Pathway analysis was performed using the Kyoto Encyclopedia of Genes and Genomes (KEGG), WikiPathways, and Reactome databases. These analyses are exploratory and rely on database-predicted target genes rather than experimentally generated gene expression data. Accordingly, the reported FDR values represent database-based overrepresentation scores of the predicted targets (Benjamini and Hochberg, 1995). The obtained results were visualized using the R programming language in RStudio (v.2024.12.0), employing the ggplot2 (v.3.5.1) and dplyr libraries (v.1.1.4) for graphical representation.
Statistical analysis
Prior to the study, a literature review was conducted to determine the minimum required sample size (de Souza et al., 2025). Using a type I error probability (α) of 0.05 and a power (1-β) of 0.80, and assuming an effect size (d) of 0.91 for the differences in milk quality parameters between groups and the relationships between milk quality parameters and miRNA levels, it was calculated that a minimum of 40 milk samples from cows would be appropriate for inclusion in the study. Power analysis was performed using PASS 11 and G*Power Version 3.1.9.2 statistical software. miRNA expression results were determined using the 2-ΔΔCt method and expressed as fold changes (Livak and Schimittgen 2001).
Before conducting significance tests on the obtained variables, normality was assessed using the Shapiro-Wilk test, and homogeneity of variances was evaluated using Levene's test. In the study, the differences in milk quality parameters between groups were analyzed using the Student's t-test for variables meeting the assumptions, and the Mann-Whitney U test for variables not meeting the assumptions. To control for potential type I error inflation arising from multiple comparisons, p-values obtained from group-wise tests of milk composition traits were additionally adjusted using the Benjamini–Hochberg False Discovery Rate (FDR) method. Adjusted p-values did not differ meaningfully from raw p-values and did not change the interpretation of statistical significance.
The strength and direction of the relationship between milk quality parameters and miRNA levels were determined using Spearman's correlation coefficient. All statistical analyses were performed using the IBM SPSS 23.0 statistical package program, and results were evaluated considering a significance level of p<0.05.
Results and discussion
Milk quality parameters and SCC-related findings
In the HSCC, all samples tested positive for the CMT, with a mean CMT score of 1.91±0.17. The overall SCC of the collected samples was approximately 320,000 cells/ml (318.05±44.34 × 10³). In the LSCC, the SCC was approximately 60.000 cells/mL (60.16±56.07×10³), whereas in the HSCC, it was around 540.000 cells/mL (540.77±196.77×10³) (Figure 1).
LSCC (n=20), HSCC (n=20), and LSCC+HSCC (Overall, n=40). Values were presented as Mean±SD
Figure 1. Distribution of somatic cell counts in groups
Leukocytes present in milk increase in response to any inflammation in the mammary tissue, leading to an elevation in the SCC. SCC is a crucial parameter associated with mammary gland health and milk quality. Essential amino acid content, caseins, whey protein fractions, and milk proteins are among the key compositional parameters contributing to milk quality (Gai et al., 2021). Although mastitis has been reported to adversely affect milk protein content and composition, the protein content of HSCC and LSCC milk samples exhibited similar patterns (Tiantong et al., 2023).
The compositional parameters, including protein, fat, lactose, and FFDM ratios, as well as pH and freezing point, were found to be similar across the groups. However, the electrical conductivity was significantly higher in the HSCC group (adj. p<0.001) (Figure 2).
LSCC: low somatic cell count group (n=20), HSCC: high somatic cell count group (n=20), Values were presented as Mean ± SEM. Fat was compared with Mann-Whitney U test. Other parameters were compared with Student-t test. ***: adj. p<0.001
Figure 2. Compositional parameters in groups
While milk fat is a primary determinant in assessing milk quality and pricing, the fat content in both groups was found to be similar but lower than that in healthy cow milk. Milk fat content is influenced by genetic and environmental factors (Marumo et al., 2022). The unexpectedly low milk fat content in milk samples may be attributed to factors such as diet composition and management practices. Indeed, the ratio of roughage to concentrate, diet composition, particle size, feeding frequency, and lactation stage are known to affect milk fat content (Ponnampalam et al., 2024). The observed negative correlation between milk fat and FFDM, protein, lactose was considered an expected finding (Yakan et al., 2021).
Although some studies have suggested that an increase in SCC may be associated with a decrease in milk lactose content, the present study found that lactose levels in the HSCC were comparable to those in the LSCC (Alessio et al., 2021; Tomanić et al., 2024). The similarity in FFDM content between groups may be attributed to the comparable proportions of fat, protein, and lactose in milk.
A near-neutral pH is indicative of healthy and high-quality milk, whereas a lower pH suggests increased milk acidity (Poghossian et al., 2019). The average pH of healthy and high-quality cow milk is known to range between 6.60 and 6.80 (Tadesse et al., 2023). Consistent with the literature, the findings of the present study fell within this reference range, with no significant differences detected between groups. The positive correlation observed between pH and FFDM and lactose content, as well as the expected negative correlation with freezing point, were deemed consistent findings. Another important parameter, electrical conductivity, was found to be higher in the HSCC. In addition to compositional parameters, minerals and ions in milk contribute to its electrical conductivity at varying levels (Yakan et al., 2021). Several studies have reported an association between mammary gland health and milk electrical conductivity (Yakan et al., 2021; Neculai-Valeanu & Ariton, 2022). The positive correlation between SCC and electrical conductivity observed in this study indirectly supports this relationship.
Soluble substances present in milk lower its freezing point (Ilie et al., 2010). Milk fat, which exists in an emulsified state, does not significantly affect the freezing point of milk. Similarly, the protein composition of milk does not influence its freezing point (Pesce et al., 2016). However, lactose and minerals, which impact the osmotic pressure of milk, are known to affect its freezing point (Zagorska & Ciprovica, 2013; Pesce et al., 2016). The freezing point of raw bovine milk, one of its most stable physical properties, typically ranges between -0.522°C and -0.540°C, and the findings of the present study were in this reference range. Moreover, no significant differences were observed between groups.
Expressions of miRNAs, correlations results and database analysis
The A260/A280 ratios of the isolated RNA samples were determined to be >1.80, and RNA concentrations were approximately 150 ng/µL (152.63 ± 17.12). In the HSCC group, miR-27a-3p levels were found to be upregulated by approximately 3-fold compared to the LSCC group (p<0.05). Similarly, miR-29B-2 levels in the HSCC group were upregulated by more than 3-fold (p<0.05). The expression levels of miR-142-5p were found to be similar between the groups, whereas miR-148a levels were upregulated by more than 6-fold in the HSCC group compared to the LSCC group (p<0.05).
HSCC: high somatic cell count group (n = 20). Values are shown as Mean ± SEM. *p<0.05.
Figure 3. Expression patterns of miRNAs in groups. LSCC: low somatic cell count group (n = 20)
A significant positive correlation was identified between the expression levels of miR-27a-3p in milk, assessed according to SCC, and miR-29B-2 (r = 0.583; p<0.001), miR-148a (r = 0.629; p<0.001), as well as electrical conductivity (r = 0.316; p<0.05). Additionally, miR-29B-2 was positively correlated with miR-148a (r = 0.605; p<0.001), SCC (r = 0.366; p<0.05), fat content (r = 0.319; p<0.05), and electrical conductivity (r = 0.390; p<0.05). Moreover, miR-142-5p activity in milk was negatively correlated with milk protein content (r = -0.311; p<0.05), while miR-148a exhibited a positive correlation with SCC (r = 0.441; p<0.01). Significant positive and negative correlations were also observed among key compositional parameters, with the correlation findings presented in the correlogram in Figure 4.
SCC: somatic cell count, FFDM: fat-free dry matter, FP: freezing point, EC: electrical conductivity
Figure 4. Correlation between miRNAs and other parameters studied
miR-27a-3p is known to be involved in molecular pathways such as inflammation and apoptosis (Luo et al., 2020). Although this miRNA was positively correlated with electrical conductivity, no direct correlation with SCC was identified. It has been reported that the activity of this miRNA varies in bovine mammary epithelial cells under different environmental conditions, such as heat stress (Wang et al., 2023). Furthermore, miR-27a-3p has been suggested to increase particularly in mastitis cases caused by Gram-negative bacteria and could serve as a potential biomarker for pathogen-induced subclinical mastitis (Wang et al., 2021; Özmen et al., 2025). The present findings align with previous literature, indicating that CMT and SCC alone may not be sufficient for assessing mammary gland health (Wang et al., 2021; Sangwa et al., 2025). Although some studies suggest that miR-27a-3p is associated with compositional parameters such as milk protein and fat, and thus may play a crucial role in milk quality, no significant correlation was observed between miR-27a-3p and protein, fat, or lactose in this study (Wang et al., 2021; Wang et al., 2023). Consequently, further research is needed to investigate this miRNA by considering genetic and environmental factors in cows.
miR-29B-2 has been reported to be expressed in milk and exhibits increased activity in inflammatory conditions such as mastitis (Srikok et al., 2020). A study suggested that miR-29B-2 could be utilized as a biomarker for detecting mastitis in cattle (Srikok et al., 2020). In this study, miR-29B-2 was found to be more than 3-fold upregulated in the HSCC and positively correlated with miR-27a-3p, which is also known to be upregulated in pathogen-specific mastitis cases. Additionally, the present study identified a positive correlation between miR-29B-2 and miR-148a, SCC, fat, and electrical conductivity. These findings suggest that miR-29B-2 may be associated with subclinical mastitis characterized by high SCC and milk quality, indicating its potential as a biomarker for evaluating these parameters.
A study demonstrated that miR-142-5p affects the expression levels of genes involved in pathways such as cell proliferation and apoptosis in bovine mammary epithelial cells and may exhibit pro-inflammatory activity in mastitis by targeting genes like BAG5 (B-cell Activation Gene 5) (Lu et al., 2021). In dairy cattle, miR-142-5p expression has been reported to decrease significantly when SCC is below 200,000 cells/mL (Stefanon et al., 2023). However, in the present study, miR-142-5p exhibited similar expression patterns in groups. Given the SCC levels in the LSCC and HSCC, it is suggested that this miRNA should be evaluated over a broader SCC range. Additionally, miR-142-5p, which has been implicated in energy metabolism through pathways such as adipogenesis, was found to be negatively correlated with milk protein (Chartoumpekis et al., 2012; Mobuchon et al., 2017). Based on both the present and previous studies, it is suggested that miR-142-5p may be SCC-dependent in its activity, and SCC levels could be a critical factor influencing this relationship (Tzelos et al., 2022).
miR-148a has been reported to exhibit anti-inflammatory activity and to undergo expression changes in response to infections with Escherichia coli and Staphylococcus aureus in mammary tissue (Hasankhani et al., 2023). However, another study reported that miR-148a expression remains unchanged in subclinical mastitis caused by Staphylococcus spp., whereas it is upregulated in Streptococcus spp.-induced subclinical mastitis (Özkan et al., 2024). Srikok et al. (2020) also noted that miR-148a can be downregulated in mastitic animals. Consistent with these findings, Petracci et al. (2023) reported that miR-148a expression may decrease in response to elevated somatic cell content in milk. In contrast, the present study demonstrated a positive correlation between miR-148a and SCC, suggesting that increasing somatic cell content may enhance miR-148a expression. This pattern aligns with reports indicating that miR-148a may participate in regulatory immune responses within mammary epithelial cells, potentially through pathways associated with FOXP3-mediated modulation of immune tolerance (Wang et al., 2024). Since miR-148a is the most abundantly expressed miRNA in bovine milk (Cai et al., 2018), the heterogeneous findings reported across studies suggest that shifts in milk composition and inflammatory status may differentially modulate its expression. Therefore, the positive association observed between miR-148a and SCC in the present study highlights the need for additional mechanistic studies to clarify the biological basis of this relationship.
The target genes of miRNAs associated with milk somatic cell count and milk quality were visualized using the miRNet database (Figure 5A). A total of 1,020 target genes were identified for miR-142, miR-148a-3p, miR-29b-3p, and miR-27a-3p. Among these, the PTEN (Phosphatase and Tensin Homolog) gene was found to be a common target of all four miRNAs, while SMAD3 (SMAD Family Member 3), TGFB2 (Transforming Growth Factor Beta 2), OTUD4 (OTU Deubiquitinase 4), and RPS6KA5 (Ribosomal Protein S6 Kinase A5) were regulated by three of these miRNAs. Additionally, 82 genes were identified as being co-regulated by two miRNAs in this study. The distribution of these target genes was visualized using a Venn diagram, which was presented in Figure 5B. A comprehensive list of genes targeted by miRNAs in this study is provided in Appendix S1.
Figure 5. Visualization of miRNA Target Genes. (A) miRNet-based network visualization of target genes for miR‐142, miR‐148a-3p, miR‐29b-3p, and miR‐27a-3p, showing that PTEN is targeted by all four miRNAs, and SMAD3, TGFB2, OTUD4, and RPS6KA5 by three. (B) Venn diagram displaying the overlap among target genes, with 82 genes co-regulated by at least two miRNAs
The target genes of the miRNAs investigated were retrieved from the STRING database and visualized using Cytoscape. The analysis revealed that the 84 targeted proteins exhibited 402 edges (Figure 6A). A detailed list of all interacting proteins, along with their corresponding interaction scores, is included in Appendix S1. The initial MCODE analysis identified a cluster of 14 genes with a score of 10.759, involving 70 edges. Among these, genes such as GATA3 (GATA Binding Protein 3), STAT3 (Signal Transducer and Activator of Transcription 3), and SMAD2 (SMAD Family Member 2), which play roles in host response regulation, were identified. Additionally, CREBBP (CREB Binding Protein), a key protein involved in lipid metabolism, was found to exhibit significant interactions (Figure 6B). In a subsequent MCODE analysis, a distinct gene cluster was identified, including BCL2 (BCL2 Apoptosis Regulator), MCL1 (MCL1 Apoptosis Regulator, BCL2 Family Member), NFE2L2 (NFE2 Like BZIP Transcription Factor 2), TP53 (Tumor Protein P53), and the SMAD gene family, which are involved in inflammation, apoptosis, and oxidative stress metabolism, exhibiting significant interactions with a score of 4.545 (Figure 6C). Furthermore, genes such as DNMT1 (DNA Methyltransferase 1), KDM6B (Lysine Demethylase 6B), and HMGB1 (High Mobility Group Box 1), associated with epigenetic modifications, were also found to exhibit significant interactions.
Figure 6. (A) Interaction network of the targeted proteins, (B) Genes representing the highest interaction cluster identified by MCODE analysis with a score of 10.759, (C) Genes representing a highly interactive cluster identified by MCODE analysis with a score of 4.545
Functional enrichment analyses were performed to explore the roles of the interacting genes within the Gene Ontology (GO) framework, specifically in the categories of biological processes, molecular functions, and cellular components (Figure 7).
Figure 7. Gene ontology (GO) analysis results presenting database-based categorization of predicted target genes into biological processes (GO: BP), molecular functions (GO:MF), and cellular components (GO:CC)
The interactions of the identified genes within biological pathways were assessed using the KEGG, Reactome, and WikiPathways databases and visualized in Figure 8. The analysis indicated database-enriched associations with pathways related to immune system regulation and infectious diseases, particularly those related to somatic cell transfer into milk, such as TLR4 and IL-17 signaling pathways. All terms with an FDR < 0.05 represent database-derived overrepresentation of the predicted target genes through the Gene Ontology, KEGG, Reactome, and WikiPathways resources, and the full list is provided in Appendix S1.
Figure 8. Pathway annotations of predicted target genes derived from KEGG, Reactome, and WikiPathways databases
Bioinformatic analyses suggested that the four miRNAs under investigation (miR‐27a, miR‐29b, miR‐142, and miR‐148a) may target the PTEN gene. In the study by Xu et al. (2022), it was demonstrated that in bovine mammary epithelial cells, PTEN suppresses autophagy triggered by microbial agents, thereby mitigating the inflammatory response. In this context and considering that somatic cell count is the most widely used biomarker in the diagnosis of subclinical mastitis (Sangwa et al., 2025), it is postulated that these miRNAs may regulate PTEN expression, consequently influencing somatic cell dynamics and playing a role in subclinical mastitis. Additionally, the genes TGFB2, SMAD3, RPS6KA5, and OTUD4, which are commonly targeted by three of the miRNAs, have been identified as potential candidates involved in regulating inflammatory responses in mammary epithelial cells. Although previous studies have suggested that the TGF‑β family, via SMAD signaling, may affect somatic cell proliferation (Panahipour et al., 2021), no study to date has directly examined the relationship between RPS6KA5 and OTUD4 expression and SCC in bovine milk. However, recent research has demonstrated that the OTUD4 and RPS6KA5 genes play significant roles in mammary health (Ci et al., 2024).
Protein-protein interaction analyses revealed that the high interaction density among proteins belonging to the HIF1A, TP53, NFE2L2, BCL2, and SMAD families underscores the critical mechanistic roles of fundamental cellular processes in regulating milk quality parameters. Previous studies conducted on bovine models have similarly reported the effects of these genes on milk quality, thereby corroborating our findings (Zhuang et al., 2023; Li et al., 2024). Furthermore, the association of MAPK1 with increased milk protein synthesis, along with the high interaction exhibited by CREBBP suggests that these signaling pathways significantly impact milk production and quality (Zhang et al., 2022).
Recent studies have reported that the expression of the GATA3 gene in mammary epithelial tissue is associated with critical functions in mammary and milk health (Bacha et al., 2024; Sandström et al., 2024). In this context, it is posited that GATA3, which has not been extensively studied in bovine models, may play a potential role in the regulation of milk quality. Moreover, functional enrichment analyses indicate that the targeted genes may exert a systemic effect on milk regulation via their involvement in various biological processes and metabolic pathways. Both in vivo and in vitro studies have demonstrated that the TLR4 and IL‑17 signaling pathways directly regulate the migration of somatic cells in mammary epithelial cells, while analyses of biological processes and cellular components have shown that the formation of transferase complexes modulates cell migration (Glynn et al., 2014; Vitenberga-Verza et al., 2022). Data obtained from the WikiPathways database further indicates that the targeted genes are directly associated with lipid and carbohydrate metabolism in milk through leptin and insulin signaling. Collectively, these findings offer new perspectives on elucidating the molecular basis of milk quality and underscore the need for further experimental studies.
Conclusion
miR-27a-3p, miR-29B-2, and miR-148a are significantly upregulated in Holstein cows with high SCC, suggesting their potential involvement in immune response and regulation of milk composition and quality. Correlations between these miRNAs, SCC, and EC highlight their relevance as biomarkers for mammary gland health in lactating cows. Target gene analysis in this study suggested PTEN as a potential key regulator, along with SMAD3, TGFB2, and STAT3, which are associated with inflammation and immune pathways. Additionally, genes linked to apoptosis, oxidative stress, and epigenetic modifications were highlighted by the bioinformatic analyses. These results provide database-derived functional annotations of the predicted targets and suggest biological pathways that could be prioritized for future experimental validation.