Differential mRNA Expression Profiling Reveals the Role of MiR-375 in Inflammation of Bovine Mammary Epithelial Cells

Simple Summary Bovine mammary epithelial cells (bMECs) are often used as cell models for mammary gland research. They are the most important cells for mammary gland function and the first line of defense for pathogen identification. MicroRNAs (miRNAs) are important regulatory factors involved in many physiological and pathological processes. Here, we examined a transcriptome profile of bovine mammary epithelial cell lines transfected with miR-375 inhibitor or negative control (NC) inhibitor, and further reveal the potential role of miR-375 in bMECs by differentially expressed mRNA analysis. We found that miR-375 potentially promotes inflammation in the mammary gland through the MAPK signaling pathway. Abstract MicroRNAs (miRNAs) are a class of small non-coding RNAs that regulate post-transcriptional gene expression and several biological processes. Bovine mammary epithelial cells (bMECs) mediate critical immune responses in the mammary gland and the occurrence of mastitis. Current research focuses on miRNA regulation of bMECs, but the miR-375 regulatory mechanism in bMECs is unclear. This study explored the role of miR-375 by profiling the transcriptome of miR-375-silenced bMECs using RNA-seq and identifying differentially expressed mRNAs (DIE-mRNAs). There were 63 DIE-mRNAs, including 48 down-regulated and 15 up-regulated mRNAs between miR-375-silenced bMECs and the controls. The Kyoto encyclopedia of genes and genomes (KEGG) and Gene Ontology (GO) functional analysis showed that the DIE-mRNAs enriched nuclear receptor subfamily 4 group A member 1 (NR4A1) and protein tyrosine phosphatase non-receptor type 5 (PTPN5) anti-inflammatory genes of the mitogen-activated protein kinase (MAPK) signaling pathway. However, they showed an opposite trend to the expression of miR-375 silencing, suggesting that miR-375 promotes bMEC inflammation through the MAPK signaling pathway. The findings of this study provide a new reference for understanding the regulation of bMEC inflammation and cow mastitis.


Introduction
Mastitis causes significant economic losses to the agricultural sector by reducing milk production and quality [1]. Many factors affect mastitis, including the genetics of the animal, which influences the susceptibility or resistance of the animal to the disease [2]. The cow mammary gland primarily comprises the bovine mammary epithelial cells (bMECs) involved in the synthesis and secretion of milk [3]. Furthermore, it is the first barrier preventing pathogens from invading cow mammary glands by secreting several immune regulatory factors [4]. Therefore, bMECs are important in the milk production and mammary gland immunity of dairy cows. However, many genes regulate the growth activity and function of bMECs [5][6][7].

Sample Collection and RNA Extraction
The MAC-T cells were collected 48 h after transfection. Total RNA was extracted from the transfected MAC-T using the TriZol kit (Takara Biomedical Technology Co., Ltd., Beijing, China) following the manufacturer's protocol. The concentration and quality of RNA were determined using agarose gel electrophoresis and a Multi-Mode Reader (BioTek, Winooski, VT, USA) (Table S1). High-quality RNA samples were used to construct cDNA libraries as described below.

Constructing cDNA Libraries and High-Throughput Sequencing
The sequencing library was constructed using the ABclonal mRNA-seq Lib Prep Kit (ABclonal, Wuhan, China), following the manufacturer's protocol. The mRNA was enriched by magnetic beads with Oligo (dT). The mRNA was randomly interrupted by the fragmentation and converted to double-stranded cDNA. Then, the double-stranded cDNA was purified using AMPure XP beads (Beckman Coulter, Pasadena, CA, USA). Further, the purified double-stranded cDNA was end repaired, A-tail enriched, and connected to sequencing joints. Then, AMPure XP Beads were used for fragment size selection. The selected fragments were PCR enriched and used to obtain the final cDNA library. Agilent Bioanalyzer 4150 (Agilent Technologies, Palo Alto, CA, USA) evaluated the library before qPCR validation. Then, the high-throughput sequencing was performed using the Illumina Novaseq 6000 (Illumina, San Diego, CA, USA) platform in Shanghai Applied Protein Technology (Shanghai, China).

Data Processing and Quality Control
After high-throughput sequencing, the raw data were stored in the FASTQ file format. We conducted an overall assessment of the quality of the raw reads by Perl, including the sequencing error rate, ATGC content, raw data composition, and comparative analysis with the reference genome to ensure the reliability of the sequencing results. Clean reads were obtained by removing the connector sequences and filtering out low-quality reads (the number of bases with a base mass value ≤25 accounted for over 60% of the total reads). The proportion of N (undetermined bases) was <5%. The obtained clean data were mapped to the bovine reference genome UMD 3.1 (http://oct2018.archive.ensembl. org/Bos_taurus/Info/Index (accessed on 25 March 2021)) using HiSAT2 software (http: //daehwankimlab.github.io/hisat2/ (accessed on 25 March 2021)). Finally, clean reads were used for subsequent analysis.

Analysis of Differentially Expressed mRNAs
The expression levels of each gene in each sample, expressed as Fragments Per Kilobase of transcript sequence per Million base pairs sequenced (FPKM), were calculated using FeatureCounts (http://subread.sourceforge.net/ (accessed on 26 March 2021)). Pearson's correlation analysis was performed between each sample. DESeq2 (http://bioconductor. org/packages/release/bioc/html/DESeq2.html (accessed on 26 March 2021)) software was used to analyze the differentially expressed mRNA in the miR-375 inhibitor group and the inhibitor NC group. The significant differentially expressed mRNAs were identified following the p-value < 0.05 and |log 2 FoldChange | > 1 threshold.

qPCR
We inquired the sequence information of miR-375, designed a neck ring structure primer specific for reverse transcription of miR-375 (5 -GTCGTATCCAGTGCAGGGTCCGAGGT ATTCGCACTGGATACGACTCACGCGA-3 ) and conducted reverse transcription reaction of miR-375. In addition, eight differentially expressed mRNAs (DIE-mRNAs) (3 up-regulated and 5 down-regulated mRNAs) were randomly selected for validating the RNA-seq results. Briefly, the 1 µg total RNA in bMECs was reverse transcribed into cDNA using the PrimeScriptTM RT Reagent Kit with gDNA Eraser (Takara Biomedical Technology Co., Ltd., Beijing, China) following the manufacturer's instructions. The 2×M5 HiPer SYBR Premix EsTaq (withTli RNaseH) fluorescence quantitative detection kit (Mei5 Biotechnology Co., Ltd., Beijing, China) detected the expression of DIE-mRNAs in transfected MAC-T cells using a CFX-96 Touch Real-Time PCR instrument (BioRad, Hercules, CA, USA). The qPCR reaction mix included 10.0 µL of 2×M5 HiPer SYBR Premix EsTaq (with Tli RNaseH), 0.8 µL of 10 nmol/µL forward and reverse primers, 100 ng of cDNA template, and adding sterilized deionized water to 20.0 µL. The reaction conditions included initial denaturation at 95 • C for 30 s, subsequent 40 cycles of denaturation at 95 • C for 5 s, and annealing at 60 • C for 30 s. All the experiments were repeated three times. The qPCR primers are listed in Table 1. All qPCR experimentation and analysis were performed following the minimum information for publication of quantitative real-time PCR experiments (MIQE) guidelines [27]. The amplification efficiency of the primers are 95-105%. The data were analyzed using SPSS software version 25.0 using GAPDH and RPS18 [28] as internal controls. The relative expression levels of the genes were calculated using the 2 −∆∆Ct method [29] and expressed as the mean ± standard deviation (X ± SD). The SPSS 25.0 software was used to test the significance of differential expression between groups through independent-samples t-tests, and p < 0.05 was considered statistically significant.

Functional Annotation of DIE-mRNAs
Gene cluster analysis was performed to display the gene expression pattern in different samples visually. Gene Ontology (GO) analysis for the molecular function (MF), biological process (BP), and cellular component (CC) of the DIE-mRNAs was performed using the clusterProfiler package in the R software (version: 4.0.3). Kyoto encyclopedia of genes and genomes (KEGG) analysis was also performed to identify pathways and the biological function of genes regulated by miR-375. qPCR findings (p < 0.01) (Figure 1b). The results showed that the silencing effect of miR-375 in MAC-T cells was good, and subsequent experiments could be carried out.
genes and genomes (KEGG) analysis was also performed to identify pathways and the biological function of genes regulated by miR-375.

The Silencing Efficiency of MiR-375 in MAC-T
Inverted fluorescence microscope revealed that miR-375 in MAC-T cells was efficiently silenced after transfection with the corresponding inhibitor (Figure 1a), consistent with qPCR findings (p < 0.01) (Figure 1b). The results showed that the silencing effect of miR-375 in MAC-T cells was good, and subsequent experiments could be carried out.

Sequence Data Analysis
We sequenced cDNA libraries using the Illumina high-throughput sequencing platform to study the miR-375 inhibitory effect on bMECs. The sequence data from the experimental and NC groups are listed in Table 2. The Q30 was over 92.5%, indicating a high base quality. Moreover, high-quality reads were obtained after removing the joint and low-quality reads ( Table 2). In general, these findings demonstrated the accuracy and reliability of the sequence data. Table 2. Quality assessment and comparison of sequencing data with reference genomes.

Sequence Data Analysis
We sequenced cDNA libraries using the Illumina high-throughput sequencing platform to study the miR-375 inhibitory effect on bMECs. The sequence data from the experimental and NC groups are listed in Table 2. The Q30 was over 92.5%, indicating a high base quality. Moreover, high-quality reads were obtained after removing the joint and low-quality reads ( Table 2). In general, these findings demonstrated the accuracy and reliability of the sequence data.

Overall Distribution of mRNA Expression
The expression levels of the overall genes in each sample are presented in Table S2. The results of Pearson's correlation analysis between samples revealed that the gene expression patterns were highly correlated among the samples (Figure 2a), further validating the reliability of the experimental design. Further analysis showed that the distribution of FPKM in each block diagram was similar, indicating that the overall gene expression abundance in each sample was similar (Figure 2b).

The Differently Expressed mRNAs in MAC-T after MiR-375 Inhibition
Inhibition of miR-375 dysregulated the expression of 63 mRNAs. A total of 48 mRNAs were significantly down-regulated, whereas the remaining 15 were significantly overexpressed (Table S3). We found that inhibition of miR-375 promotes expression of nuclear Animals 2022, 12, 1431 6 of 12 receptor subfamily 4 group A member 1 (NR4A1) and protein tyrosine phosphatase nonreceptor type 5 (PTPN5). Cluster analysis showed that the high and low expressed genes in the samples were clustered together, indicating that DIE-mRNAs regulate critical processes in bMECs (Figure 3).

Overall Distribution of mRNA Expression
The expression levels of the overall genes in each sample are presented in Table S2. The results of Pearson's correlation analysis between samples revealed that the gene expression patterns were highly correlated among the samples (Figure 2a), further validating the reliability of the experimental design. Further analysis showed that the distribution of FPKM in each block diagram was similar, indicating that the overall gene expression abundance in each sample was similar (Figure 2b).

The Differently Expressed mRNAs in MAC-T after MiR-375 Inhibition
Inhibition of miR-375 dysregulated the expression of 63 mRNAs. A total of 48 mRNAs were significantly down-regulated, whereas the remaining 15 were significantly overexpressed (Table S3). We found that inhibition of miR-375 promotes expression of nuclear receptor subfamily 4 group A member 1 (NR4A1) and protein tyrosine phospha-

Validation of the DIE-mRNAs Using qPCR
The expression pattern of eight DIE-mRNAs (three up-regulated and five downregulated genes) based on RNA-seq was validated using qPCR analysis (Figure 4). The qPCR results were consistent with the RNA-seq results, confirming the accuracy and reliability of the RNA-seq results.

GO Enrichment and KEGG Analyses of the DIE-mRNAs
The 63 DIE-mRNAs were annotated to molecular function (MF), biological process (BP), and cellular components (CC) through GO enrichment analysis (Table S4). These DIE-mRNAs enriched nuclear membrane and nuclear envelope under CC. For MF, the DIE-mRNAs enriched DNA-binding transcription activator activity, DNA-binding transcription factor activity, and growth factor receptor binding. However, the enriched BP included regulation of epithelial cell proliferation, intracellular receptor signaling, and extracellular structure organization (Figure 5a). Meanwhile, the KEGG pathway enrichment analysis revealed that the DIE-mRNAs mainly enriched the mitogen-activated protein kinase (MAPK) signaling pathway, retinol metabolism, and cortisol synthesis and secretion (Table S5 and Figure 5b). tase non-receptor type 5 (PTPN5). Cluster analysis showed that the high and low expressed genes in the samples were clustered together, indicating that DIE-mRNAs regulate critical processes in bMECs (Figure 3).

Validation of the DIE-mRNAs Using qPCR
The expression pattern of eight DIE-mRNAs (three up-regulated and five down-regulated genes) based on RNA-seq was validated using qPCR analysis (Figure 4). The qPCR results were consistent with the RNA-seq results, confirming the accuracy and reliability of the RNA-seq results.

GO Enrichment and KEGG Analyses of the DIE-mRNAs
The 63 DIE-mRNAs were annotated to molecular function (MF), biological process (BP), and cellular components (CC) through GO enrichment analysis (Table S4). These DIE-mRNAs enriched nuclear membrane and nuclear envelope under CC. For MF, the DIE-mRNAs enriched DNA-binding transcription activator activity, DNA-binding transcription factor activity, and growth factor receptor binding. However, the enriched BP included regulation of epithelial cell proliferation, intracellular receptor signaling, and extracellular structure organization (Figure 5a). Meanwhile, the KEGG pathway enrichment analysis revealed that the DIE-mRNAs mainly enriched the mitogen-activated protein kinase (MAPK) signaling pathway, retinol metabolism, and cortisol synthesis and secretion

Discussion
This study found that silencing miR-375 down-regulated and up-regulated the expression of 48 and 15 mRNAs, respectively, in MAC-T cells. Regulation of miRNA is a complex process that affects gene expression in the entire cell. MiRNA is directly regulated for some genes, and indirectly regulated for a large number of genes [9]. Therefore, this may be the reason why the down-regulated genes were more abundant than the up-regulated ones after silencing of miR-375. GO revealed that the 63 DIE-mRNAs regulated several Animals 2022, 12, 1431 9 of 12 biological processes, including regulation of epithelial cell proliferation, the intracellular receptor signaling pathway, and extracellular structure organization. KEGG enrichment analysis further revealed that the DIE-mRNAs regulate the MAPK signaling pathway, retinol metabolism, and cortisol synthesis and secretion. The MAPK signaling pathway regulates the growth and differentiation of cells, adaptation to environmental stress, inflammation, and other important cellular physiological and pathological processes [30][31][32]. In addition, the MAPK signaling pathway participates in the occurrence of mastitis in dairy cows [33]. Other signaling pathways were not closely related to the function of bMECs.
KEGG analysis further revealed that inhibiting miR-375 up-regulated the expression of PTPN5 and NR4A1 genes, all regulated via the MAPK signaling pathway. NR4A1 (also called Nur77) is a member of the orphan nuclear receptor family 4A (NR4A), which regulates inflammation and immunity [34,35]. The anti-inflammatory gene NR4A1 is rapidly expressed in the early stages of inflammatory upon entry of stimuli such as lipopolysaccharide (LPS), or secretion of cytokines such as interleukin-1β (IL-1β), and tumor necrosis factor α (TNF-α) [36]. It indicated that NR4A1 might be an important mediator in the early inflammation. NF-κB plays a critical role in regulating inflammation [37,38]. In vascular endothelial cells, NR4A1 up-regulates IκBα expression but inhibits the activation of NF-κB by binding to the IκBα promoter [39]. NR4A1 modulates the expression of NF-κB by directly interacting and blocking the binding of p65 to its κB, inhibiting the secretion of pro-inflammatory cytokines [40]. Modulating NR4A1 expression induces NF-κB dependent activation of macrophages [41]. The interaction between NR4A1 and NF-κB/p65 in microglia alleviates brain injury caused by cerebral ischemia, thus inhibiting neurogenic inflammation [42]. NR4A1 inhibited LPS-induced inflammation in acute liver injury by directly binding to TRAF6 [43]. Therefore, NR4A1 may be a potential target for regulating and preventing inflammation in bMECs. Meanwhile, miR-375 modulates mastitis in dairy cows by regulating the expression of NR4A1.
Reversible phosphorylation of tyrosine residues plays a key role in many signaling pathways [44], catalyzed by protein tyrosine phosphatases (PTPs) [45]. In humans, the PTP genes have been linked to several diseases and, thus, potential therapeutic targets in such complications [46]. PTP genes are divided into four families. The non-receptor protein tyrosine phosphatase (PTPN) belongs to the class I family [47]. The PTPN gene family members regulate numerous physiological processes and participate in the development and pathogenesis of numerous diseases [47][48][49]. PTPN5 (also known as striatum-enriched protein tyrosine phosphatase (STEP)) mainly participates in regulating neuronal signal transduction, and abnormal expression of this protein impairs motor control and cognitive function [50,51]. PTPN5 binds and reduces the affinity to MAPK substrates, negatively regulating the activity and cell localization of MAPK family members. These events block the kinase nuclear translocation of some cellular functions, such as inflammation [52]. In vivo studies revealed that PTPN5 inhibits the growth of breast tumors by blocking the epidermal growth factor (EGF)-induced MAPK signaling pathway [53]. Additionally, high PTPN5 activity decreases with aging [54], whereas PTPN5 deficiency induces neuronal inflammation and exacerbates ischemic brain injury [55]. Inhibition of PTPN5 reverses cognitive deficit impairment in mouse models with Alzheimer's disease [56]. Our RNA-seq results revealed that PTPN5 modulates inflammation by inhibiting the MAPK signaling pathway in the mammary gland of dairy cows with mastitis. Therefore, silencing miR-375 alleviates mastitis in cows by promoting PTPN5 expression and while inhibiting the MAPK signaling pathway.
It is worth noting that IL-6 is a key inflammatory cytokine. The infected mammary glands can promote the secretion of IL-6 and initiate the inflammatory response and body immunity by activating various signaling pathways [57]. It has been reported that the expression of IL-6 is positively correlated with the severity of mastitis in dairy cows [57,58]. We found that inhibition of miR-375 down-regulated IL-6 expression. Therefore, the downregulation of IL-6 in the miR-375 inhibition group suggested that inhibition of miR-375 might alleviate the inflammatory response.

Conclusions
In summary, miR-375 silencing dysregulated the expression of 63 mRNAs in bMECs. Additionally, miR-375 silencing increased the expression of NR4A1 and PTPN5 genes, all anti-inflammatory genes, via the MAPK signaling pathway. Given silencing of miR-375 significantly up-regulates NR4A1 and PTPN5 gene expression, miR-375 potentially promotes inflammation in the mammary gland through the MAPK signaling pathway. The findings of this study provide a new perspective on treating mastitis in cows.