MicroRNAs modulate immunological and inflammatory responses in Holstein cattle naturally infected with Mycobacterium avium subsp. paratuberculosis

MicroRNAs (miRNAs) regulate the post-transcriptional expression of genes by binding to their target mRNAs. In this study, whole miRNA sequencing was used to compare the expression of miRNAs in ileocecal valve (ICV) and peripheral blood (PB) samples of cows with focal or diffuse paratuberculosis (PTB)-associated lesions in gut tissues versus (vs) control cows without lesions. Among the eight miRNAs differentially expressed in the PB samples from cows with diffuse lesions vs controls, three (miR-19a, miR-144, miR32) were also down-regulated in cows with diffuse vs focal lesions. In the ICV samples, we identified a total of 4, 5, and 18 miRNAs differentially expressed in cows with focal lesions vs controls, diffuse lesions vs controls, and diffuse vs focal lesions, respectively. The differential expression of five microRNAs (miR-19a, miR-144, miR-2425-3p, miR-139, miR-101) was confirmed by RT-qPCR. Next, mRNA target prediction was performed for each differentially expressed miRNA. A functional analysis using the predicted gene targets revealed a significant enrichment of the RNA polymerase and MAPK signaling pathways in the comparison of cows with focal vs no lesions and with diffuse vs focal lesions, respectively. The identified miRNAs could be used for the development of novel diagnostic and therapeutical tools for PTB control.


Ethical statement
All methods were carried out in accordance with relevant guidelines and regulations.The study is reported in accordance with ARRIVE guidelines (https:// arriv eguid elines.org).The Animal Ethics Committee of the Servicio Regional de Investigation y Desarrollo Agroalimentario (SERIDA) approved the procedures on the animals included in this study.All procedures were authorized by the Regional Consejería de Agroganadería y Recursos Autóctonos of the Principality of Asturias (approval code PROAE 29/2015 and PROAE 66/2019) and were carried out in accordance with the European Guidelines for the Care and Use of Animals for Research Purposes (2012/63/EU).PB, gut tissues, and fecal samples were collected by trained personnel and in accordance with good veterinary practices.

Animals and PTB infectious status
The fourteen Holstein cows animals included in this study belonged to a single farm in Asturias (Spain).The infectious status of these animals was determined by histopathological analysis of gut tissues, ELISA for the detection of anti-MAP antibodies, and fecal and gut tissues PCR and bacteriological culture as previously described 21,36 .According to their location and extension, inflammatory cell type, and mycobacterial load, PTBassociated histopathological lesions were classified in focal and diffuse as previously described 14 .The average age of the animals without lesions and with focal and diffuse lesions was 5.45, 5.09, and 4.38 years old, respectively.

RNA extraction, RNA-Seq library preparation, sequencing, and differential expression analysis
Total RNA isolation, RNA-Seq library preparation, and sequencing were previously performed 21 .The raw reads were filtered by length (> 75 bp) and by the percentage of ambiguous bases (< 10%) using FastQC v0.11.9 37 and had PCR primers and Illumina adapters trimmed with Trim Galore 38 .The fastq reads (NCBI Gene Expression Omnibus database; Acc.Number GSE137395) were mapped against the most recent Bos taurus reference genome sequence (ARS-UCD1.2) using the Spliced Transcripts Alignment to a Reference aligner (STAR 2.5.3a) 39.The alignments files (.bam) from STAR were used to generate a table of counts for each gene using the FeatureCounts function from Rsubread 2.12.0 40 .Gene counts were then normalized with the mean-of-ratios method of DESeq2 1.38.0 41 .In the current study, DESeq2 was also used for the differential gene expression analysis.An mRNA was differentially expressed if its false discovery rate (FDR)-adjusted p-value was lower than 0.05.sRNA extraction, sRNA-Seq library preparation, and sequencing.
MiRNAs were isolated from PB samples using the PAXGene miRNA isolation kit following the instructions of the manufacturer (Qiagen).Total RNA from gut tissue samples was extracted with the miRNeasy Mini Kit according to the manufacturer's protocol (Qiagen).For each sample, 150 ng of total RNA was used to prepare a sRNA-Seq library using the NEB-Next Small RNA Sample Preparation following the manufacturer´s protocol (New England Biolabs, Ipswich, MA, US).Each of the generating libraries was run through a mini-PROTEAN TBE Precast gel 5% (BioRAD, Madrid, Spain), and fragments from 140 to 325 bp of the library were excised from the gel and purified.A unique pool containing all the purified libraries was prepared and sequenced in 101 bp single-read in a NovaSeq 6000 sequencer (Illumina) at the Genomic Unit of the Scientific Park of Madrid, Spain.

miRNAs differential expression analysis
Quality control of the raw reads was performed using FastQC v0.11.9 37 .Trimmed reads were analyzed using sRNAtoolbox, which consists of independent tools for sRNA analysis 42,43 .sRNAbench, a tool from sRNAtoolbox, was used to filter and align the trimmed reads to the Bos taurus genome (ARS-UCD1.2).Reads with less than 15 bp and with a mean phred score below 20 were filtered out of the analysis, and the maximum number of allowed mismatches in the mapping process was 1. miRbase 22.1 database was used to annotate the miRNAs, and other ncRNAs were annotated using RNAcentral 20.0 and Ensembl release 104 databases.For the normalization and differential expression analysis of the miRNAs, the alignments obtained from sRNAbench were fed directly to sRNAde, another web-based tool from sRNAtoolbox 42 .sRNAde uses different algorithms (DESeq2 2.12 and EdgeR 1.38.0) to test for differential expression.A miRNA was differentially expressed if its FDR-adjusted p-value Figure 1.Bioinformatic pipeline.sRNA-Seq and RNA-Seq from PB and ICV samples of infected and control cows was performed, and the resulting reads were trimmed and aligned to the Bos taurus reference genome.Then, miRNA and mRNA differential expression was performed, and the results were integrated to identify a negative correlation between each differentially expressed miRNA and its predicted mRNA target(s).

Reverse transcription quantitative PCR (RT-qPCR) for RNA-Seq validation
Changes in the expression of five of the differentially expressed miRNAs were validated by RT-qPCR.For specific and highly sensitive RT-qPCR detection of a mature miRNA, 10 ng of total RNA was reverse transcribed into cDNA using the miRCURY LNA RT kit (Qiagen) according to the manufacturer´s instructions.For the RT reactions, a total volume of 10 µl (2 × MiRCURY RT Enzyme mix, reaction buffer, target RNA, and RNasefree water) was incubated for 60 min at 42 °C to perform the RT and 5 min at 95 °C to inactivate the reverse transcriptase.RT control reactions without the enzyme mix were included.The resultant cDNAs (1:40 dilution) were amplified on a LightCycler ® 480 thermal cycler using the miRCURY SYBR ® Green PCR Kit (Qiagen) in the presence of miRNA-specific and locked nucleic acid (LNA)-enhanced forward and reverse primers (Qiagen).PCR conditions were: 95 °C for 2 min (heat activation), followed by 40 cycles of 95 °C for 5 s (denaturation) and 56 °C for 30 s (annealing and extension), and a final melting curve.Appropriate controls (no template) were included.RT-qPCR experiments were performed by duplicate and in each experiment three replicates were run for each combination of primers and samples.
To obtain Ct values, the 2nd derivative max method of the LightCycler ® 480 Software was used.Using the results of the samples from cows without lesions or with focal lesions as controls, fold changes in expression were calculated using the 2 -ΔΔCt method.Normalization was performed using the U6 snRNA and 5S rRNA endogenous reference genes.Two comparisons were made: samples of cows with diffuse lesions vs controls and with diffuse vs focal lesions.The results were expressed as fold changes and were standardized by log 2 transformation to be comparable to the sRNA-Seq differential expression results.To test if the differences between the groups were statistically significant, an unpaired t-test was run using R 45 .The Pearson correlation coefficient between the sRNA-Seq and RT-qPCR quantitative results was calculated using R. Differences and correlations were considered statistically significant if the p-value was ≤ 0.05.

miRNA-gene target prediction
The putative target mRNAs for each miRNA were predicted using the tool mRNAconsTarget from sRNAtoolbox with the miRanda 46 and PITA 47 algorithms with default parameters.In parallel, the prediction database Tar-getScan 48,49 was used to predict mRNA targets with a score ≥ 0.2.To reduce potential false positive targets, only targets predicted by the three algorithms were considered for further analysis.Next, the correlations between the miRNAs and predicted mRNA expression levels were tested with the Spearman and Pearson correlation tests using R 45 .Using both correlation tests, mRNAs with a negative correlation (R < 0 and p-value ≤ 0.05) with respect to their related miRNA were considered as putatively affected by the miRNA expression.

Functional analysis
Using the predicted gene targets for each comparison, an enrichment analysis of gene ontology (GO) processes and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways was performed using ClusterProfiler 4.0 50 and String 12.0 51 .An FDR-adjusted p-value lower than 0.05 was considered significant.The function of the candidate genes was searched in GeneCards by searching their gene symbol.

sRNA-seq data
The animals included in the study were tested by histopathological analysis of gut tissues, acid-fast stain (Zielh-Neelsen), specific antibody response against MAP measured by ELISA, fecal and tissue qPCR, and bacteriological culture as previously described 21,36 .The infectious status of the 14 animals included in the study are presented in Table 1.The four control animals had no histopathological lesions in gut tissues and had a negative result for ELISA, fecal and tissue bacteriological culture, and qPCR.All the infected animals had focal or diffuse lesions in gut tissues.The five animals with focal lesions in gut tissues had a positive tissue qPCR, none of them had heavy bacterial load (> 50 cfu/g) in feces and tissues, and none of them was culled due to PTB-associated clinical signs.Only one of the animals with focal lesions was fecal PCR positive.All the animals with diffuse lesions had positive ZN and ELISA results, positive fecal and tissue qPCR results, and positive bacteriological culture from gut tissues.Four of the five animals with diffuse lesions had heavy bacterial load in gut tissues.
sRNA-Seq libraries were prepared from total RNA extracted from PB and ICV samples of the 14 animals which were classified into three groups according to the results of the histopathological analysis: animals with focal and diffuse lesions and animals without lesions in gut tissues and regional lymph nodes (controls).An average of 40,954,530 and 36,762,469 raw reads per library were obtained after sequencing the PB and ICV libraries, respectively, and an average of 35.8 million reads per library (92.4%) passed the quality control process.The alignment of the filtered reads resulted in 34.2 million mapped reads (94.6%) per library.Mapped reads were annotated using miRBase v22 to identify known miRNAs in PB and ICV samples, and 31.92 and 15.51 million reads mapped to known miRNAs in the miRBase v22, respectively.The remaining reads that were not identified as miRNAs were lncRNAs, transfer RNAs (tRNAs), small nuclear RNAs (snRNAs), small nucleolar RNAs (snoRNAs), ribosomal rRNA (rRNA), and unassigned RNAs (Fig. 2).

miRNA differential expression analysis
A total of 420 and 358 mature miRNAs identified in the PB and ICV samples were analysed for differential gene expression with the EdgeR and DESeq2 algorithms.In general, there were more miRNAs with decreased than increased expression levels in all the comparisons (Tables 2 and 3).As expected, the number of differentially www.nature.com/scientificreports/expressed miRNAs was larger in the ICV than in the PB samples.None of the miRNAs was differentially expressed in both the PB and ICV samples and, therefore, our results showed PB and ICV-specific miRNA expression.A list of the differentially expressed miRNAs identified with the EdgeR and/or DESeq2 algorithms in the PB samples of the infected animals is presented in Table 2.In the comparison of focal vs controls, we did not identify any dysregulated miRNA.In the PB samples from animals with diffuse lesions vs controls, the consensus of the two algorithms revealed five miRNAs; two upregulated (bta-miR-2425-3p and bta-miR-139) and three downregulated (bta-miR-27a-5p, bta-miR-144, and bta-miR-181a).In this comparison, DESeq2 was able to identify an additional bta-miR-32 (fold = −2.62)and the Edge2 identified the bta-miR-27a-5p (fold = −3.68)and bta-miR-181a (fold = −1.86).The bta-miR-2425-3p was the most upregulated miRNA in our study with a  www.nature.com/scientificreports/fold of 5.76.Among the eight miRNAs differentially expressed in PB samples from cows with diffuse lesions vs controls, three downregulated miRNAs (bta-miR-19a, bta-miR-144, bta-miR-32) were also found downregulated in PB samples of cows with diffuse vs focal lesions.The list of the differentially expressed miRNAs identified in the ICV samples with the EdgeR and DESeq2 algorithms is presented in Table 3.The bta-miR-2478 (fold = 3.23) was the most highly upregulated miRNA in the tested ICV samples.In the ICV samples, we identified a total of 4, 5, and 18 miRNAs differentially expressed in cows with focal lesions vs controls, diffuse lesions vs controls, and diffuse vs focal lesions, respectively.In the comparison of focal lesions vs controls, DESeq2 was able to identify two upregulated (bta-miR-2478 and bta-miR-150) and two downregulated (bta-miR-23a and bta-miR-23b-3p) miRNAs.In the ICV samples of cows with diffuse lesions, both DESeq2 and EdgeR identified two downregulated miRNAs (bta-miR-135b and bta-miR-433) when compared with control cows.In this comparison, DESeq2 was able to identify two additional downregulated miRNAs (bta-miR-146a and bta-miR-99a-5p) and EdgeR identified the bta-miR-215 (fold = 2.60).Some of the miRNAs that were differentially expressed in this comparison (bta-miR-433, bta-miR-146a, bta-miR-99a-5p) were also dysregulated in the comparison of diffuse vs focal lesions.Two miRNAs (bta-miR-23a and bta-miR-23b-3p) were downregulated in the comparisons; focal lesions vs controls and diffuse vs focal lesions.

Validation of the differentially expressed miRNAs in PB samples by RT-qPCR
The expression changes observed in five of the differentially expressed microRNAs (miR-19a, miR-144, miR-2425-3p, miR-139, miR-101) were validated by RT-qPCR using total RNA extracted from PB samples of the animals included in the study (Fig. 3).The selected miRNAs were all identified with both DESeq2 and EdgeR algorithms.The bta-miR-2425-3p and bta-miR-139 were selected because they were the only miRNAs significantly upregulated in PB samples.The bta-miR-19a and bta-miR-144 were selected because they were found downregulated in all the tested comparisons.The results of the qRT-PCR are presented in Fig. 3 and were expressed as fold changes and standardized by log 2 transformation to be comparable to the miRNA-Seq differential expression results analyzed with DESeq2.Using multiple reference genes for normalization minimizes technical variation 52 .Therefore, the normalization of the RT-qPCR results was performed using the expression of the U6 snRNA and 5S rRNA as the endogenous reference genes.As seen in Fig. 3 similar results were obtained when the U6 snRNA or 5S rRNA were used as endogenous genes for data normalization.The RT-qPCR results showed a statistically significant upregulation of the bta-miR-2425-3p in the comparison of cows with diffuse lesions vs controls (p-value ≤ 0.05).In the comparison of diffuse vs focal lesions, a statistically significant downregulation of the bta-miR-19a was observed (p-value ≤ 0.05).A similar trend (upregulation or downregulation) was observed in the results obtained with both RNA-Seq and RT-qPCR methods.Indeed, the Pearson correlation coefficient calculated for the miRNA-Seq and RT-qPCR data was equal to 0.978 and 0.959 when the U6 snRNA or 5S rRNA were used as endogenous genes; p-values = 0.003 and 0.0006, respectively.

miRNA-gene target prediction
Potential target mRNAs of the differentially expressed miRNAs were predicted using three different algorithms (miRanda, PiTA, TargetScan).All the consensus mRNA targets were then investigated for negative Spearman www.nature.com/scientificreports/and Pearson correlations with corresponding miRNA using mRNA sequencing data from the same animals.The number of differentially expressed miRNAs and their predicted mRNA targets for each comparison is presented in Table 4.In the PB samples, 35 and 32 mRNA targets were identified in the comparisons: diffuse lesions vs controls and diffuse vs focal lesions, respectively.In ICV samples, 33, 2, and 198 mRNA targets were identified in the comparisons focal lesions vs controls, diffuse lesions vs controls, and diffuse vs focal lesions, respectively.The heatmaps representing the changes in the expression of the differentially expressed miRNAs in the different comparisons along with the changes in the expression of their target mRNAs are presenting in Figs. 4 and  5.As seen in Fig. 4A and B, 32 mRNA targets identified in the two comparisons of PB samples (diffuse lesions vs controls and diffuse vs focal lesions) were upregulated and potentially controlled by only two downregulated miRNAs, bta-miR-144 and bta-miR-19a.In the ICV samples, 33 mRNA targets would be controlled by two upregulated (bta-miR-150 and bta-miR-2478) and two downregulated (bta-miR-23a and bta-miR-23b-3p) miRNAs in the comparison focal lesions vs controls (Fig. 4C).In the comparison of cows with diffuse lesion vs controls, the upregulation of bta-miR-215 (fold = 2.60) would result in the downregulation of the mRNAs encoding the G Protein-Coupled Receptor 22 (GPR22) and the BMX Non-Receptor Tyrosine Kinase (BMX).Finally, 198 mRNAs were found dysregulated and would be controlled by 16 DE miRNAs in the comparison of cows with diffuse vs focal lesions (Fig. 5).

Functional enrichment analysis
Functional analysis of all the predicted gene targets was performed by enrichment analysis of GO processes and KEGG pathways (Table 5).Using the predicted gene targets identified in the PB samples, no GOs or pathways could be identified.Using ClusterProfiler, two KEGG pathways, the RNA polymerase (bta:03020) and the mitogen-activated protein kinase (MAPK) signaling pathway (bta:04010) were found enriched in the ICV samples of cows with focal lesions vs controls and with diffuse vs focal lesions, respectively.In the comparison of focal lesions vs controls, the RNA Polymerase I Subunit E (POLR1E) was positively associated with the bta-miR-23a and the RNA Polymerase III Subunit G (POLR3G) was negatively associated with the bta-miR-2478, the most highly upregulated miRNA in the tested ICV samples (fold = 3.23).In the comparison of diffuse vs www.nature.com/scientificreports/focal lesions, downregulation of bta-miR-214, bta-miR-23a, bta-miR-99a-5p, bta-let-7c, bta-miR-204, bta-let-7e caused upregulation of 11 genes belonging to the MAPK signaling pathway (bta04010).In this comparison, the regulation of immune system processes (GO:0002682) and the nucleolus cellular component (GO:0002682) were significantly enriched with 29 and 22 enriched genes in each pathway, respectively.In addition, several KEGG pathways related to the induction of cytokines and inflammatory immune response (WP997) were enriched including the MAPK signaling pathway, pathways in cancer, Kaposi sarcoma-associated herpesvirus infection, Hepatitis C, Measles, Epstein-Barr virus infection, pancreatic cancer, P53 signaling pathway, lysosome, PI3K-Akt, and JAK-STAT signaling pathways.

Discussion
Previous studies analyzed the differential expression of miRNAs in Holstein cattle experimentally infected with MAP 33,34 .A comparison of the serum miRNome of MAP-challenged IFNɣ responders to unchallenged controls six months after infection did not identify significant differences in miRNA expression 33 .Similarly, Shaughnessy et al. did not detect differentially expressed miRNAs at either the early 6 months or late (43, 46, and 49 months) intervals across seropositive and seronegative animals 34 .Using naturally infected cattle previous studies compared  www.nature.com/scientificreports/ the miRNA profiles of cattle positive and negative for MAP antibodies by ELISA 32 or compared miRNA profiles according to the stage of MAP infection defined by fecal shedding, ELISA, and clinical signs 53 .In cattle, however, signs of infection such as the appearance of focal lesions in gut tissues can be detected before fecal shedding and ELISA, and therefore, examining the transcriptome of animals according to the histopathological lesions facilitates the identification of animals in the subclinical stage of the infection when the amount of MAP and antibodies can be undetected with current pre-mortem diagnostic methods.In the present study, the absence or presence of PTB-associated lesions (focal or diffuse) was used to define the severity of MAP infection and to compare the miRNA profiles of animals in early and more advanced stages of the infection vs uninfected cattle without lesions in gut tissues.To our knowledge only one study compared miRNA expression in the serum of cattle classified according to histopathological data, serology, and MAP bacterial culture into four disease groups: control, mild, moderate, and severe 54 .Instead of using RNA-Seq, Gupta et al. used NanoString technology (Seattle, WA, US) to detect bovine miRNAs using probe sets specific for human miRNAs.Careful should be taken in generalizing human results to cattle as the miRNome and targets of miRNAs can be organism specific.Indeed, none of the miRNAs detected by Gupta et al., was detected in the current study.
Most previous studies analyzed circulating miRNA profiles in serum samples.In our study, we examined whether the PB miRome recapitulated the miRome of the ICV, the primary site of MAP colonization.None of the identified miRNAs was differentially expressed in both the PB and ICV samples and, therefore, our results showed PB and ICV-specific miRNA expression.Some miRNAs identified in or study matched with the ones identified in previous studies describing miRNA responses to mycobacterial infections including the bta-miR-27a-5p, bta-miR-32, bta-miR-23a, bta-miR-135b, bta-miR146a, bta-miR146b, bta-miR-433, bta-let-7e, and btalet-7c 32,35,55,56 .Using the predicted gene targets of the miRNAs identified in the PB samples, no GOs or pathways could be identified.However, by analyzing the predicted gene targets of the miRNAs identified in the ICV samples, we were able to distinguish pathways associated with different stages of MAP infection.To our knowledge, our study is the first to analyze miRNAs profiles in tissue samples of MAP-infected cattle.
Our analysis showed that specific miRNAs were dysregulated in cows with different PTB-associated lesions.Specifically, in PB samples, we identified eight miRNAs in the comparison of cows diffuse lesions vs controls and three in the comparison of cows with diffuse vs focal lesions.Importantly, the miRNA-Seq analysis's results and the RT-qPCR results correlated.Among the eight miRNAs differentially expressed in PB samples from cows with diffuse lesions vs controls, three (bta-miR-19a, bta-miR-144, bta-miR-32) were also found downregulated in PB samples of cows with diffuse vs focal lesions.In agreement with these findings, the bta-miR-32 was also found downregulated in a previous study where the miRNA profiles of Holstein cattle positive and negative for MAP antibodies were compared 32 .
In the ICV samples, we identified a total of 4, 5, and 18 miRNAs differentially expressed in cows with focal lesions vs controls, diffuse lesions vs controls, and diffuse vs focal lesions, respectively.This suggests that the miRNA expression in ICV changes more as PTB-associated lesions become more severe.In the comparison of cows with focal lesions vs controls, DESeq2 was able to identify two upregulated (bta-miR-2478 and bta-miR-150) and two downregulated (bta-miR-23a and bta-iR-23b-3p) miRNAs.Some of the differentially expressed miRNAs in the comparison diffuse lesions vs controls (bta-miR-433, bta-miR-146a, bta-miR-99a-5p) were also dysregulated in the comparison diffuse vs focal lesions.The bta-miR-146a with a conserved human ortholog acts as an inhibitor of the pro-inflammatory immune response 57 .In our study, upregulation of miR-146a would negatively regulate the TNF Receptor Associated Factor 6 (TRAF6) and control the pro-inflammatory immune response 58 .Two miRNAs (bta-miR-23a and bta-miR-23b-3p) were downregulated in both comparisons, focal lesions vs controls and diffuse vs focal lesions.We also integrated miRNA and mRNA expression data to ascertain whether the differences in miRNA profiles would be reflected in the mRNA data set.The 32 common mRNA targets identified in the two comparisons of PB samples were upregulated in all the comparisons and would be controlled by only two downregulated miRNAs, bta-miR-144 and bta-miR-19a.In ICV samples, 33, 2, and 198 mRNA potential targets were identified in the comparisons focal lesions vs controls, diffuse lesions vs controls, and diffuse vs focal lesions, respectively.Functional analysis using the predicted gene targets allowed the identification of two enriched KEGG pathways, the RNA polymerase (bta:03020) and the mitogen-activated protein kinase (MAPK) signaling pathway (bta:04010), in the ICV samples of cows with focal lesions vs controls and with diffuse vs focal lesions, respectively.
A more thorough examination of the gene targets in each enriched pathway allowed us to obtain a better understating of the regulation of the innate and inflammatory responses by miRNAs (Fig. 6).In the comparison of cows with focal lesions vs controls (Fig. 6A), the RNA Polymerase I Subunit E (POLR1E) would be positively regulated by bta-miR-23a, and the RNA Polymerase III Subunit G (POLR3G) negatively regulated by the bta-miR-2478, the most highly upregulated miRNA in the tested ICV samples (fold = 3.23).POLR3G is part of the RNA polymerase III.POLR3G is involved either in the recruitment and stabilization of the subcomplex within RNA polymerase III, or in stimulating catalytic functions of other subunits during initiation.The RNA polymerase III acts as a nuclear and cytosolic non-self dsDNA sensor, is involved in the positive regulation of innate immune responses and Type I IFN and Nuclear factor (NF-κβ) activation through the RIG1 pathway and plays a key role in sensing and limiting infection by intracellular bacteria and DNA viruses.The POLR3G downregulation might be responsible for the blockade of the innate immune responses generally observed upon MAP infection.This blockage might allow MAP to persist and establish infection within macrophages during the long-term subclinical phase of the infection.The downregulation of POLR3G might be a way for MAP to prolong its survival within infected macrophages by decreasing DNA-driven innate immune responses.
In the comparison of cows with diffuse lesions vs controls (Fig. 6B), the upregulation of bta-miR-215 (fold = 2.60) would result in the downregulation of the mRNAs encoding the G Protein-Coupled Receptor 22 (GPR22) and the BMX Non-Receptor Tyrosine Kinase (BMX).BMX is expressed in hematopoietic cells of the myeloid lineage, granulocytes and monocytes, and has a significant role in cytokine signaling and inflammation 59 .More specifically, BMX is required for the phosphorylation and activation of the Myeloid Differentiation Primary Response protein (MyD88), Mal, and transforming grow factor-β activated kinase (TAK1), which leads to the activation of the NF-кβ and mitogen-activated protein kinase (MAPK) signaling pathways.Impaired production of BMX might therefore result in a control of type I IFN and pro-inflammatory cytokines production.Interestingly, a recent study showed that the downregulation of miR-27a in MAP-infected macrophages significantly inhibited pro-inflammatory responses in MAP infected macrophages by targeting the expression of TAK1 binding proteins 2 and 3 (TAB2/3), important components of the MAPK signaling pathway 60 .In the comparison of cows with diffuse lesions vs controls, we observed that the upregulation of bta-miR-146a was negatively associated with the expression of TRAF6 which in normal conditions would ensure a check-in of the pro-inflammatory signalling of NF-κβ and MAPK.Mir146a expression is induced both in macrophages and in mice after mycobacterial infections, and further supress the iNOs expression and NO generation via NF-кβ and MAPK signalling 61,62 .
In the comparison of cows with diffuse vs focal lesions (Fig. 6C), we observed that the downregulation of bta-miR-214, bta-miR-23a, bta-miR-99a-5p, bta-let-7c, bta-miR-204, bta-let-7e was associated with the upregulation of 12 genes belonging to the MAPK signaling pathway.The MAPK signaling pathway plays a significant role in the induction of pro-inflammatory responses that occur in MAP-infected cattle in advanced stages of the infection 63 .Among the 12 potential gene targets controlled by the differentially expressed miRNAs identified in the comparison diffuse vs focal lesions we found the Tumor Necrosis Factor Receptor Type 1-associated DEATH Domain Protein (TRADD), Rac Family Small GTPase 2 (RAC2), and Colony Stimulating Factor 1 (CSF1).TRADD mediates programmed cell death, RAC2 is involved in the generation of reactive oxygen species, and CSF1 promotes the release of pro-inflammatory chemokines, and thereby plays an important role in the stimulation of inflammatory processes.In our study, TRADD, RAC2 and CSF1 expression regulation was associated with the bta-mir-214.Besides the MAPK signaling pathway, the software String revealed that the target genes identified in the comparison diffuse vs focal lesions were associated with the regulation of the immune system process (GO:0002682) and with the induction of cytokines and inflammatory immune responses (WP997) in several diseases including cancer (bta05200), Kaposi sarcoma-associated herpesvirus infection (bta05167), Hepatitis C (bta05167), Measles (bta05162), Epstein-Barr virus infection (bta05169), pancreatic cancer (bta05212).Signaling pathways associated with inflammation such as the P53 signaling pathway (bta04115), lysosome (bta04142), PI3K-Akt (bta04151), and JAK-STAT (bta04630) were also found significantly enriched.

Conclusions
Through modulation of the expression of miRNAs, MAP infection modifies host gene expression and modulates the innate and inflammatory responses in distinct stages of the infection.In early stages of the infection, the bta-miR-247 would downregulate POLR3G, an activator of the RNA polymerase III that functions as a nonself dsDNA sensor and is implicated in the activation of DNA-driven innate immune responses.The POLR3G downregulation might be responsible, at least in part, for the blockade of the innate immune responses that allows MAP persistence within infected macrophages during the long-term subclinical stage of the infection.As the infection progresses, MAP grows within granulomas but the upregulation of the bta-miR-215 and bta-miRNA-146a would control the activation of an excessive pro-inflammatory response by targeting BMX and TRAF6, respectively.In clinical stages of the infection, however, the dysregulation of bta-miR-214, bta-miR-23a, bta-miR-99a-5p, bta-let-7c, bta-miR-204, bta-let-7e would upregulate several genes of the MAPK signaling pathway causing the activation of a strong pro-inflammatory response.Aside from the utility of the miRNAs identified in the current study as biomarkers, modifying or targeting the expression of specific miRNAs might serve as a novel strategy for PTB control. https://doi.org/10.1038/s41598-023-50251-9

Figure 2 .
Figure 2.Categories of the identified ncRNAs in peripheral blood and ileocecal valve samples.The reads that were not identified as miRNAs were long-non coding RNAs (lncRNAs), transfer RNAs (tRNAs), small nuclear RNAs (snRNAs), small nucleolar RNAs (snoRNAs), ribosomal rRNA (rRNA), and other and unassigned RNAs.

Figure 3 .
Figure 3. Validation of the differentially expressed miRNAs in PB samples by RT-qPCR.Validation of the differences in expression of five miRNAs (bta-miR-101-3p, bta-miR-19a-3p, bta-miR-139, bta-miR-144, bta-miR-2425-3p) was performed by RT-qPCR.Changes in the expression levels (expressed as the log 2 of the fold change) of the selected differentially expressed miRNAs in the comparisons diffuse vs controls and with diffuse vs focal lesions was calculated using RNA-Seq (in yellow) and RT-qPCR.RT-qPCR results using the U6 snRNA or 5S rRNA as endogenous genes are presented in grey and purple, respectively.Significant changes in the miRNA expression (p-value ≤ 0.05) for each comparison are indicated with an asterisk.

Figure 4 .
Figure 4. Negative correlations between miRNAs and predicted target mRNAs expression levels.Heatmaps representing changes in the expression of the differentially expressed miRNAs in peripheral blood samples in the comparisons: diffuse lesions vs controls (A) and diffuse vs focal lesions (B) along with the changes in expression of their predicted target mRNAs.(C) Heatmap representing changes in the expression of the differentially expressed miRNAs in ICV samples in the comparison focal vs controls along with the changes in expression of their predicted target mRNAs.

Figure 5 .
Figure 5. Negative correlations between miRNAs and predicted target mRNAs expression levels.Heatmap representing changes in the expression of the differentially expressed miRNAs in ICV samples in the comparison diffuse vs focal lesions along with the changes in expression of their predicted target mRNAs.

Figure 6 .
Figure 6.Potential role of the identified miRNAs in the regulation of innate and inflammatory responses in response to MAP infection.(A) In the comparison of cows with focal lesions vs controls, the POLR3G would be negatively regulated by the bta-miR-2478.POLR3G is a member of the RNA polymerase III that functions as a non-self dsDNA sensor.The POLR3G downregulation might be responsible for the blockade of the innate immune responses that allow MAP persistence within infected macrophages.(B) In the comparison of cows with diffuse vs controls, the upregulation of bta-miR-215 would reduce the expression of BMX and the activation of MAPK and NF-кβ.Upregulation of bta-miR-146a would negatively regulate TRAF6 and the pro-inflammatory response during the subclinical stage of MAP infection.(C) In the comparison of cows with diffuse vs focal lesions, however, we observed that the dysregulation of bta-miR-214, bta-miR-23a, bta-miR-99a-5p, bta-let-7c, bta-miR-204, bta-let-7e would upregulate the expression of several genes of the MAPK signaling pathway causing the activation of a strong pro-inflammatory response commonly observed in the clinical stages of the infection.BMX Non-Receptor Tyrosine Kinase (BMX), transforming grow factor-β activated kinase (TAK1), mitogen-activated protein kinase (MAPK), TNF Receptor Associated Factor 6 (TRAF6), interleukin 1 receptor-associated kinase (IRAK), Toll-like receptors (TLR), cyclic GMP-ANP Synthase (cGAS), stimulator of interferon response CGAMP interactor 1 (STING), RNA Polymerase III Subunit G (POLR3G), Interferon regulatory factor 3 (IRF3), nuclear factor (NF-кβ), RNA sensor RIG1, interferon (IFN), tumor necrosis factor (TNF), Interleukin 1 Receptor Associated Kinase 4 (IRAK1/4).Created with BioRender.com.

Table 1 .
Results of the histopathological analysis, ZN stain, ELISA, PCR, and bacteriological culture from the 14 animals included in the study.Bacterial load was classified as low (< 10 CFU), medium (between 10 and 50 CFU) or heavy (> 50 CFU).N Ziehl-Neelsen, CFU colony forming units.

Table 2 .
List of the differentially expressed miRNAs identified with the EdgeR and DESeq2 algorithms in the peripheral blood samples of animals with diffuse lesions vs controls and animals with diffuse vs focal lesions.

Table 3 .
List of the differentially expressed miRNAs identified with the EdgeR and DESeq2 algorithms in ileocecal valve samples of infected animals with focal or diffuse lesions vs controls and with diffuse vs focal lesions.

Table 4 .
Number (No.) of differentially expressed (DE) miRNAs in peripheral blood (PB) and ileocecal valve (ICV) for each comparison and their predicted mRNA targets identified with PiTA, miRanda, and TargetScan algorithms and with a negative Spearman and Pearson correlation.

Table 5 .
Gene ontologies (GOs) and canonical pathways associated with the predicted gene targets for the differentially expressed miRNAs in each comparison.