Comparative Characterization of Immune Response in Sheep With Caseous Lymphadenitis Through Analysis of the Whole Blood Transcriptome

1. Introduction

Caseous lymphadenitis (CL) is a contagious chronic bacterial disease of global relevance that affects small ruminants, horses, cattle, pigs, wildlife, and even humans. It is also known as pseudotuberculosis or "cheesy gland." The disease is caused by a gram-positive facultative intracellular parasite called Corynebacterium pseudotuberculosis (Cp). It belongs to the vitally important CMNR group encompassing Corynebacterium, Mycobacterium, Nocardia, and Rhodococcus [1,2,3]. It is often referred to as "iceberg disease" because it can cause a gradual decline in physical condition, thin ewe syndrome, fertility problems, behavioral changes, decreased milk and wool production, and poor carcass quality [4,5,6,7]. Unfortunately, the disease is challenging to eradicate due to difficulties in detecting subclinically affected animals, poor response to antibiotic treatment, and the ability of C.pseudotuberculosis to persist in different environments while evading the host's immune system [8,9].

Sheep infected with Cp experience a slow progression of encapsulated pyogranulomas, which eventually lead to the formation of larger abscesses in superficial or visceral lymph nodes and organs, such as the lungs, liver, and kidneys [10,11]. The main factors contributing to the virulence of Cp are its cell wall lipids and exotoxin phospholipase D, which promote bacterial persistence and spread within the host and result in macrophage death and chronic inflammation [12,13,14]. Researchers are currently conducting genomic and proteomic studies to identify Cp's antigenic potential and pathogenicity islands involved in evading the host's immune system. These studies aim to produce more effective vaccines for preventing infection and limiting emerging antibiotic resistance to reduce CL outbreaks on farms [5,15,16,17]. However, the cost of commercial vaccines remains high, and they have yet to provide full protection against infection or effective control of CL.

The reasons why some sheep and goats are more susceptible to Cp than others, whether due to physiological or genetic factors, are still not fully understood. However, studies suggest that a cellular immune response, specifically a Th1 response, plays a significant role in controlling the infection, as the microorganism is facultatively intracellular. The formation of pyogranulomas, which are characteristic of C. pseudotuberculosis infection, involves both humoral and cell-mediated immunity, a complex process that requires the activation of macrophages and T-lymphocytes [18,19,20]. Recent advances in high-throughput genomic technologies have enabled the first comprehensive profiling of the goat spleen transcriptome and proteome in response to experimental infection with Cp [21,22,23]. Currently, genome-wise approaches and the generation of big data are considered promising selection-associated techniques that can provide animal farming industries with the ability to cope with the issues caused by diseases through breeding for health traits [24].

Cases of CL have increased in the Czech Republic over the last decade, affecting an average of 20% of animals in approximately 30 tested sheep and goat flocks. To better understand the disease and host defence, we studied global gene expression in sheep infected with Cp. The study's main goals were the comprehensive characterization of the ovine whole blood transcriptome in the later course of the infection in naturally infected ewes. Moreover, we performed a functional comparison of the differentially expressed genes between healthy and affected animals based on their differing serological statuses. Thus, this research provides new insights into the host immune response against Cp in sheep.

2. Materials and Methods

2.1. Animals and sampling procedure

The study included nine five-year-old females of the East Frisian breed from two different farms located in Czechia. The animals had been maintained under comparable housing conditions and nutritional regiments and were regularly dewormed. A commercial ELISA serological test (Elitest CLA, HYPHEN Biomed, Neuville-sur-Oise, France) was carried out to detect any IgG antibodies against phospholipase D, which could indicate the presence of disease. The confirmation test for Cp presence was performed by classic bacterial cultivation. The first experimental group consisted of the three serologically negative ewes exposed to Cp infection, which did not show CL clinical signs (group exposed - negative; EN). The second experimental group from the same flock consisted of the three diseased, positive ewes (group exposed – positive; EP) with visible abscesses. Finally, the third control group included three healthy females from the other noninfected farm (control group – CN). Peripheral blood samples (3 ml) were drawn from the jugular vein into Tempus™ Blood RNA Tubes (Thermo Fisher Scientific, Waltham, MA, USA), immediately cooled and then stored at -80 °C until RNA isolation.

2.2. RNA isolation, library preparation and sequencing of the transcriptome

Total RNA was extracted from the cryopreserved stabilized whole blood of ewes using a Tempus™ Spin RNA Isolation Kit (Thermo Fisher Scientific) according to the manufacturer's instructions. RNA quantity and quality were assessed by a NanoDrop 1000 spectrophotometer and Qubit 4 fluorometer (Thermo Fisher Scientific). The RNA integrity was verified using an Agilent Bioanalyzer 2100 (Agilent Technologies, Santa Clara, CA, USA). All samples fulfilled the criteria of a 260/280 ratio between 2.0 and 2.2 and RNA integrity numbers greater than 7; the average yield was 114.7 ng/µl of RNA. The mRNA-seq library was constructed per the manufacturer's instructions using the NEBNext®Ultra™ II Directional RNA Library Preparation Kit (New England Biolabs, Ipswich, MA, USA). Samples were sequenced on an Illumina NextSeq 550 system following the vendor’s protocol using Mid Output Kit v2.5 and PhiX Sequencing Control Kit v3 (Illumina Inc., San Diego, CA, USA). On average, 32 million 150-bp paired-end reads were obtained per sample.

2.3. Transcriptome data processing, differential expression analysis and Gene Ontology classification

The raw reads in FASTQ format were cleaned from low-quality reads and adapter sequences using Trim Galore software. The remaining reads were aligned to the NCBI Ovine transcriptome reference (GCF_002742125.1_Oar_rambouillet_v1.0_rna) using Bowtie2 (v2.3.4.3). Unmapped reads were extracted using SAMtools (v0.1.19-44428cd). Furthermore, the obtained contigs were compared to the Ovis aries transcriptome and genome reference and to the human genome, viral, bacterial and rRNA database (SILVA db) using blast (v2.5.0+). The cleaned contigs were annotated using Trinotate (v3.1.1). Mapping and quality statistics were performed using QualiMap (v.2.2.2-dev) in the Multisample BAM QC mode [25]. The count of mapped reads was extracted and saved into a matrix using the featureCounts (v 1.6.5) program [26].

Transcriptomic data were divided into three groups of sheep (EN, EP, and CN) based on their serological status. We used two different approaches, DESeq2 (v 1.26.0) and Edge R (v 3.28.0), to perform differential expression analysis. DESeq2 corrects for library size and RNA composition by using a method called median of ratios, while Edge R uses a normalization method called trimmed mean of M-values (TMM). We compared the transcript expression between all animal groups (CN x EN, CN x EP, and EN x EP) and used the Benjamini‒Hochberg (BH) method for p value adjustment in DESeq2 or the false discovery rate (FDR) in Edge R. Transcripts with a log2 fold change greater than 1.5 or less than -1.5 and a padj < 0.05 or FDR < 0.05 were considered significantly differentially expressed (DE).

The sets of DEGs were interpreted using the Gene Ontology (GO) system with the online software tool WEB-based GEne SeT AnaLysis Toolkit, WebGestalt, which allowed for classification of enriched cellular components, biological processes, and pathways implemented in the KEGG (Kyoto Encyclopedia Genes and Genomes), Reactome and Wiki Pathway databases [27]. Due to a noticeable lack of functional data for Ovis aries, human or cattle genomes need to be used as references. The BH correction for a padj (FDR) value of 0.05 and the top 10 significance levels were set as the thresholds for significantly enriched GO terms and pathways. Redundancy reduction was realized as the Weighted Set Cover approach, which finds top gene sets while maximizing gene coverage.

2.4. Reverse-transcriptase real-time qPCR validation

To verify the accuracy of the high-throughput sequencing data, four DEGs with higher opposite expression, MFA5 (microfibril-associated protein 5), MEDAG (mesenteric estrogen-dependent adipogenesis), HMCN1 (hemicentin 1), and LY6G6C (lymphocyte antigen 6 family member G6C), were chosen for analysis via real-time quantitative PCR to confirm their transcriptional levels. Reverse transcription was carried out using the High-Capacity cDNA Reverse Transcription Kit according to the manufacturer's instructions (Thermo Fisher Scientific). The relative mRNA levels were analysed using 2-step RT qPCR with specific predesigned TaqMan Gene Ovine Expression Assays (Thermo Fisher Scientific) consisting of primers and fluorescently labelled probes (Supplementary Table S1). The short TaqMan probes were designed to span exon‒exon junctions to avoid contamination of genomic DNA. As an optimal reference, B2M (beta 2 microglobulin), POLR2A (RNA polymerase II subunit A), and PPIB (peptidylprolyl isomerase B) genes were selected (geNorm average stability value M = 0.25). qPCR measurements were performed in triplicate using a LightCycler 480 II Instrument (Roche Diagnostics, Basel, Switzerland, EU). The reaction mixture was prepared in an 8 µl volume and consisted of the TaqMan Gene Expression Assay, 480LightCycler Probes Master Mix (Roche Diagnostic), and 2 µl of 5 x diluted cDNA. Experimental data analysis was performed using qBase+ Premium software (Biogazelle, Ghent, Belgium). The relative gene expression was estimated using the generalized qBase model of delta-delta-Ct, which takes into account gene-specific PCR efficiency correction and multiple reference gene normalization [28].

3. Results

3.1. Sequencing quality overview

We used high-throughput RNA deep sequencing to analyse the global expression profiles of sheep with different CL serological statuses in the affected flock compared to the control, healthy animals. The basic sequencing statistics of the reads are shown in Table 1. After removing adapters and sequences of low quality, an average of 31.94 mil reads per library was finally obtained. Alignment to the ovine reference yielded mean values per library of 26.5 million reads in positive ewes, 24.2 mil reads in negative ewes and 19.9 mil reads in uninfected controls, representing 63.64 to 78.29% of the cleaned reads. The rate of uniquely mapped reads was lower, reaching a maximum mean of 57% in the control group, but the sample coverage indicated that every base of reference was, on average, represented more than 21 times in the reads. The reads mapping to multiple regions were excluded from downstream DE analysis. The GC content was determined to be between 41% and 51%, meeting the necessary analytical criteria. We obtained 27,430 transcripts and 12,230 genes, covering 54.8% and 43.3% of known transcripts and genes, respectively, in the ovine transcriptome and genome reference. The quality scores were excellent and sufficient for the analysis of transcriptional levels and the determination of differentially expressed transcripts and genes between the groups.

3.2. Characterization of the ovine whole blood transcriptome

We analysed the transcriptome of healthy and diseased ewes from ordinary flocks to understand how sheep respond to late-stage natural Cp infection. Differential gene expression was determined using DESeq2 and Edge R software, which gave similar results, indicating their validity. Compared with the controls, we identified 253 differentially expressed transcripts (DETs) in the EP, of which 118 were downregulated and 135 were upregulated, and 385 DETs in the EN, of which 226 were downregulated and 159 were upregulated. Data on the differentially expressed transcripts are provided in Supplementary Table S2. The identified transcripts were assigned to the 191 annotated DEGs in the EP (102 upregulated and 89 downregulated) and 256 in the EN (106 upregulated and 150 downregulated). The number of differentially regulated transcripts observed in the EN group was more than 13% above the level in the EP group. Regarding the log2 fold change, the CL-positive tested sheep had higher average values (FC = 2.75) and CL-negative sheep had lower (FC = 1.99) values. Under the given experimental conditions, we also analysed ten differentially expressed transcripts of five genes between the EN and EP groups.

The numbers of uniquely expressed and overlapping common genes are depicted in a Venn diagram (Figure 1). A total of 86 DEGs were found to be shared between the EN and EP groups, of which 52 were downregulated and 34 were upregulated. Alternatively, we observed 168 unique DEGs in EN animals (70 upregulated vs. 98 downregulated) and 102 unique DEGs in the EP group (68 upregulated vs. 34 downregulated). The volcano plot (Figure 2) shows expression profiling of the transcripts related to ewe infection with C. pseudotuberculosis. The Venn diagram and Volcano plots illustrate that EN animals had more unique DEGs, and the share of downregulated genes was greater than that in the EP group. In contrast, in the EP group, the number of upregulated genes was slightly higher.

3.3. Functional significance of unique and common DEGs

The unique and common DEGs observed among the top 10, ranked by their log2 fold-change, are presented in Table 2 and Table 3, respectively. In addition, the functional importance of the DEGs and their biological role were screened using the NCBI Gene Resource and UniProt Knowledgebase considering bovine or human orthologues.

3.3.1. Unique genes with the highest DE

The top unique upregulated transcripts in the blood of the diseased ewes included the most abundant killer cell lectin-like receptor subfamily B member 1 (KLRB1, CD161, CLEC5B), which is involved in the native immune response (IR). Other highly expressed genes encode the enzyme kelch-like family member 18 (KLHL18), which regulates proteosomal protein ubiquitination and protein glucosyl galactosyl hydroxylysine glucosidase (PGGHG). This top group further consists of kidney-specific chloride voltage gated channel (CLCNKA) and immunoglobulin superfamily DCC subclass member 3 (IGDCC3). Vascular cell adhesion protein 1 (VCAM-1) participates in various processes of adaptive immunity, and phospholipid phosphatase related 5 (PLPPR5) regulates neurite biosynthesis.

The top unique downregulated genes in diseased ewes were represented by extracellular matrix (ECM) proteins and organizators, namely, asporin (ASPN), which showed the lowest expression (FC = -6.7), secreted protein acidic and cysteine rich (SPARC), and secreted phosphoprotein 1 (SPP1). Decreased expression was further measured in twist family bHLH transcription Factor 1 (TWIST1), which regulates tissue development and fat metabolism, and solute carrier family 7 member 11 (SLC7A11), which enables the transport of small amino acids.

In the exposed group, the top unique genes with elevated expression were more functionally diversified and showed considerably less expression change (FC of 2.41 to 3.45) than in diseased animals. The most abundant transcripts were identified as a ubiquitous transcription factor, zinc finger protein 227 (LOC101109384), a multidomain ankyrin repeat domain 7 protein (ANKRD7) whose function is not yet fully recognized, and an eye component crystallin gamma S (CRYGS). In addition, immune-relevant highly expressed transcripts were identified as ubiquitin ligase tripartite motif containing 45 (TRIM 45), which may act as a transcriptional regulator in the mitogen-activated protein kinase signalling pathway (MAPK) and a pathogen pattern recognition receptor of T lymphocytes - antigen WC1.1.

Conversely, the top downregulated genes showed higher FC in the negative sheep than in the positive sheep, ranging from 3.8 to 7.44. The unique least expressed genes were large structural proteoglycan 4 (PRG4), an activator of Wnt signalling secreted frizzled related protein 1 (SFRP1), and a heme protein myeloperoxidase (MPO). Among the low-abundance transcripts, we further observed two ECM components containing collagen: tubulointerstitial nephritis antigen-like 1 (TINAGL1) and collagen type VI alpha 1 chain (COL6A1).

3.3.2. Common genes with the highest DE

The common and highly differentially expressed genes are presented in Table 3. The highest expression with the same direction was observed for these immune-relevant DEGs: MHC class I polypeptide-related sequence B (MICB), lymphocyte antigen 6 family member G6C (LY6G6C), related to MHC class III, and hemicentin 1 (HMCN1). On the other hand, the microfibril-associated protein 5 (MFAP5), proteoglycan syndecan 2 (SDC2), and mesenteric estrogen-dependent adipogenesis (MEDAG) transcripts were identified as being the least numerous in the tested blood samples of the CL-affected animals, and had the lowest gene expression values in the study in terms of FC. Only two common genes with the opposite direction of expression between diseased and exposed CL groups were mapped and predicted as MHC class I receptor - leukocyte immunoglobulin-like receptor subfamily A member 6 (LILRA6), involved in adaptive IR, and immuno-associated protein GTPase IMAP family member 7 (GIMAP7).

3.4. Enrichment of gene ontology (GO) terms and pathways in the investigated ovine groups

3.4.1. Cellular localization of differentially expressed genes

To explain the molecular function of identified DEGs that may be related to the ovine CL, the particular gene clusters connected with different levels of expression within the EP and EN ewes were interpreted using the WebGestalt tool. The results of the cellular component enrichment are depicted in Figure 3. In the positive group, the gene list contains 191 DEG IDs, of which 184 IDs were unambiguously mapped to unique Entrez genes and 134 gene IDs (70.2%) could be annotated to the enriched cellular components (CC). Similarly, in the negative group, the list consists of 256 gene IDs, of which 246 IDs were mapped and 178 of them (69.5%) were annotated. In terms of the affected cellular components, the DE genes in diseased ewes were localized mainly in the cell surface (GO:0009986; FDR 0.003) and endocytic vesicle lumen (GO:0071682; FDR 0.003), followed by the ECM (GO:0031012; FDR 0.005), extrinsic component of plasma membrane (GO:0019897; FDR 0.006) and endoplasmic reticulum (GO:0005783; FDR 0.05). In exposed ewes, we recorded significant differential gene expression in collagen-containing ECM (GO:0062023; FDR 0.0006), ECM (GO:0031012; FDR 0.004) and supramolecular fibre (GO:0099512; FDR 0.02).

3.4.2. Influenced biological processes and role of DEGs in the host immune response

The examined ewes did not show significant enrichment of the classified biological processes with ongoing Cp infection. The top-ranked categories, according to their FDR value, are presented in Figure 4. In the case of the EP group, the GO Slim summary was based on 160 DEGs (85.1% of all DEGs of the EP group) that could be successfully mapped to the genome reference and annotated to the categories. In the EN group, 215 DEGs (85.3%) could be classified. After redundancy reduction, the most enriched in the positive ewes were immune response (GO:0006955; FDR 0.2) and regulation of immune system process (GO:0002682; FDR 0.3). Among the potentially most influenced processes, we further identified the organonitrogen compound catabolic process (GO:1901565; FDR 0.6), response to cytokine (GO:0034097; FDR 0.6), and ECM organization (GO:0030198; FDR 0.6). In negative ewes, we mainly observed enrichment of the cell cycle (FDR 0.07; GO:0007049), immune response (FDR 0.5), regulation of cell death (GO:0010941; FDR 0.7) and positive regulation of transport (GO:0051050; FDR 0.7), particularly norepinephrine transport (GO:0015874; FDR 0.4).

The ewes' immune response may include up to 51 immunomodulating DEGs in the EP and EN groups; 20 genes are shared between the groups. Using GO analysis, we further refined their most likely biological function during the interaction with Cp in the infected flock. The involvement of genes in specific biological processes, ranked according to their enrichment ratio, is detailed in Table 4.

In positive, diseased animals, we observed the influence of biological processes that lead, on the one hand, to activation and positive regulation of the immune cells and molecules and, on the other hand, to the deregulation of some effector components of the sheep defence system. In particular, regulation of antigen receptor-mediated SP and B-cell receptor SP, adaptive IR based on somatic recombination of immunoglobulin receptors, response to other organism, and some processes of innate IR involved mainly more transcribed DEGs. Humoral IR, complement activation and homeostasis of calcium ions showed lowered transcriptional activity. The momentary cytokine-mediated inflammatory response measured in the blood transcriptome of the diseased sheep might be lower, but cytokine regulation itself was well balanced. Cytokine signalling accounted for approximately half of all DEGs assigned to the immune response in diseased sheep.

In the exposed sheep, we observed almost the same number of IR genes that were differentially regulated as in the positive animals. Unlike the positive animals, the majority of immune DEGs showed limited transcriptional activity in the corresponding BP. Reduced gene expression was recorded for part of innate IR and regulation of protein metabolic process. Upregulated DEGs preferentially act in the positive regulation of natural killer cell-mediated cytotoxicity. Furthermore, they positively participate in I-kappa B kinase/NF-kappa B signalling and gamma-delta T-cell activation. Apoptosis and stress response were both positively and negatively regulated in exposed animals. Unlike positive, diseased animals, higher transcription of the genes driving cellular calcium ion homeostasis was observed.

3.4.3. Differential gene expression in biological pathways influenced by C.pseudotuberculosis

To understand the main potentially affected pathways, the Kyoto Encyclopedia of Genes and Genomes (KEGG), Reactome [29] and WikiPathways [30] databases were applied under the same conditions as the GO analysis. Analysis was subsequently performed on all 447 predicted DE genes across the comparisons among the affected and healthy control sheep. As the table demonstrates, functional enrichment of the pathways following the results obtained in the BP classification of both experimental groups corresponded with the revealed CC. The share of successfully annotated DEGs in the three tested databases was rather low, from 39.9% to 56.4% in the EP group and 41.3% to 52.4% in the EN group. In this respect, Reactome provided the best coverage of the terms through all comparisons. An overview of the main pathways with altered gene expression that met at least the p < 0.005 criteria is highlighted in Table 5. We investigated additional pathways that CL affects in the positive, diseased ewe that are often connected with the immune system. The obtained results show significant DE of the “Binding and Uptake of Ligands by Scavenger Receptors” pathway. Highly different transcript abundances were observed for genes maintaining cell and tissue structure and dynamic cooperative signalling in the “ECM organization” and “ECM proteoglycans” pathways. Transcriptional regulation of the “Immunoregulatory interactions between a Lymphoid and non-Lymphoid cell” pathway seems to be substantially changed in the EP as well as in the EN group. Other important, possibly highly altered pathways in the EP group include “Cytokine signalling” and “Immune System”.

In the exposed ewes, significant DE was observed for p53 signal transduction events (FDR 0.04) and correlated TP53 regulation of cell death gene pathways. TP53 (tumour suppressor protein 53) controls the transcription of genes involved in the cell cycle and regulates the expression of a number of genes involved in the intrinsic apoptosis pathway triggered by cellular stress. Consistent with previous DE results in the cellular component category, differential regulation of gene expression in the “Cell cycle” and “Syndecan interactions” was identified.

3.5. Validation of RNA-Seq data by RT‒qPCR

To confirm the transcriptome data, four target DEGs with the opposite expression (MFA5, MEDAG, HMCN1, and LY6G6C) and three reference genes with very good geNorm expression stability (B2M, POLR2A and PPIB) were analysed by RT‒qPCR. The same RNA samples used for RNA-Seq were used for qRT‒PCR. We evaluated the normalized relative expression levels measured in the downregulated MEDAG and MFAP5 genes and upregulated HMCN1 and LY6G6C genes. The comparative analysis revealed the consistent trend of the results obtained by both methodological approaches and confirms the validity of DE measured in the blood of experimental and control sheep (Figure 5).

4. Discussion

We first describe the potential influence of gene expression on difficult-to-treat caseous lymphadenitis through the analysis of the blood transcriptome in sheep. In contrast to previous animal model studies conducted under controlled experimental conditions in mice, goats, or cell cultures [18,31,32,33,34], we measured gene expression in sheep raised under field conditions. This approach yielded unique insights into the host‒pathogen relationship but also has certain limits. Despite the comprehensive data, transcriptome analysis can only offer limited information about metabolic and health status during sampling. In addition, soluble immune molecules often act briefly, locally, at low concentrations and can therefore be challenging to analyse in the blood. It is important to note that whole blood is a complex mixture of various cells and may not necessarily reflect molecular events in organs [35]. Nevertheless, the alteration of gene expression in peripheral blood may be specific with disease progression and used for pathophysiology identification [36]. Understanding how the host and pathogen interact in the initial phases of a disease is crucial to comprehending the infection's mechanism. However, transcriptomic analysis of the later phases may guide changes occurring in the preclinical and clinical phases of the disease [37].

In this study, we highlight the value of RNA-seq in monitoring genes that could modify the immune response during infection. Our findings show that the blood transcriptome differences between the two experimental groups were less pronounced. Several authors have observed that the full onset of CL can be preceded by an incubation period in which antibodies are not detectable, and some animals may not show them at any phase after the infection [38,39,40]. There is also evidence of a slight difference in the blood transcriptome between infected Mycobacterium avium subsp. paratuberculosis (MAP) and exposed animals belonging to paratuberculosis-positive herds [37,41]. The phylogenetically related genera Mycobacterium and Corynebacterium probably share some pathological features due to their intracellular nature [4,7].

A survey of DEGs observed in the experimental groups discussed below is presented in Table 2, Table 3 and Table 4. Compared to the previous transcriptome measure of early exposure to Cp in goats, [21] the number of differentially expressed genes was substantially lower in our study. Additionally, the changes in metabolic processes and pathways in the later phase of infection were less pronounced. The host immune mechanisms appeared to be relatively balanced in controlling the adverse effects of infection, especially in diseased ewes. The number of differentially expressed genes in infected animals often decreases with time [42].

A member of the Ig superfamily, VCAM-1 (vascular cell adhesion molecule), was identified as the main differentially expressed unique gene between diseased and exposed ewes. It is a cell surface glycoprotein expressed by the cytokine-activated endothelium that helps regulate inflammation-associated vascular adhesion and the transendothelial migration of macrophages and T cells. The expression of VCAM-1 is activated by proinflammatory cytokines and, together with integrins, it plays a central role in leukocyte recruitment during inflammation [43]. Its high concentration in the blood of diseased animals may indicate active regulation of the inflammatory response to ongoing Cp infection.

The most abundant and unique DE gene investigated in the EP group compared to healthy controls was KLRB1 (killer cell lectin-like receptor subfamily B member 1, CD161), which encodes a receptor expressed on natural killer (NK) cells and peripheral blood memory T cells. It defines a functionally distinct subset of proinflammatory NK cells and marks cells that have retained the ability to respond to innate cytokines during their differentiation. It is interesting to note that upregulation of KLRB1 was also found in whole blood of MAP-positive cattle [44]. The most downregulated gene in diseased animals detected using RNA-seq was a common gene for microfibril-associated protein 5, MFAP5. In cattle, it is involved in definitive haemopoiesis, and its deficiency in mice is connected with decreased levels of neutrophils in circulation [45].

The highly abundant transcript among exposed animals was hemicentin, HMCN1. This gene encodes a conserved large ECM protein, a member of the Ig superfamily [46]. It is involved in cell adhesion, cytoskeleton organization and response to bacteria. On the other hand, mesenteric estrogen-dependent adipogenesis, MEDAG, was examined as the most downregulated transcript. This gene is believed to promote adipocyte differentiation, lipid accumulation, and glucose uptake in mature adipocytes, but its role in immunity is yet to be determined [47]. It is also noteworthy that the transcriptional activity of the PRG4 gene was deficient in the EN group. It has been proposed that PRG4 is regulated by the inflammatory response, typically in a negative fashion, and once inflammation has subsided, PRG4 returns to homeostatic levels to maintain tissue health [48].

The infection primarily interfered with the expression of many sensing surface receptors. Various other differentially expressed C-lectin receptors of NK cells and macrophages were identified, such as CLECL1 and CLEC2D and CLEC4E (Mincle). Recently, it was discovered that Mincle (macrophage inducible Ca2+-dependent lectin) acted as a receptor for corynomycolic acids of the corynebacterial cell wall. It probably cooperates with other C-lectin receptors in Cp recognition, activation of immune cells and modulation of the proinflammatory response to prevent possible tissue damage [49,50,51]. In addition to the C-lectin receptors, RNA-Seq analysis identified several G protein-coupled receptors and six different transcripts of LIR receptors. Ig-like receptors LIR, LILRA5 and LILRA6 of NK cells showed DE in both experimental groups, which was more pronounced in exposed ewes. They drive leukocyte activity and play an essential role in regulating the immune response. The significance of LIR receptors in the innate immune response has been recently broadened to include the ability to develop antigen-specific immune memory [52].

The category of peptide presenting antigens to NK receptors was almost highly upregulated in affected ewes. It consists of two classical BOLA class I histocompatibility antigens, OVAR and OLA-I, and MHC-III lymphocyte antigen 6 family member G6C, LY6G6C, and nonclassical MHC proteins and cooperating molecules, MHC class I polypeptide-related sequence B, MICB, UL16-binding protein 2-like, ULBP2 and butyrophilin subfamily 2, member A2, BTN2A2. Many studies have confirmed the pivotal role of MHC molecules in antigen presentation required to elicit an immune response against invading pathogens. The expression of MHC genes is altered during mycobacterial infection in cattle [42,53,54]. In transcriptome analysis of the Cp-infected spleen of dairy goats, the MHC class I protein complex was identified as one of the most enriched terms [21]. The genes for CD1A and CD1E molecules structurally related to MHC, which present lipids of microbial origin to T-cell receptors, were differentially expressed only in the transcriptome of the diseased ewes. The CD1 system is considered a versatile player in the immune response, sitting at the crossroads of innate and adaptive immunity, which may be involved in numerous infectious and inflammatory responses [55].

Our study found that both experimental groups had a decreased abundance of GZMA, GZMB, and GZMH transcripts in their blood. These transcripts are part of the granzyme serine protease family, which is specific to cytolytic T lymphocytes and NK cells. Proteases directly contribute to immune defence by inducing cell cytolysis and pathogen clearance in innate immunity [56]. GZMA is a potential biomarker of human Mycobacterium tuberculosis (MT) infection and disease, as its serum levels are significantly lower in patients with active MT disease [57]. Additionally, GZMA may be captured and internalized within mycobacteria-infected monocytes in influenced organs to inhibit bacterial growth and prevent infection [58]. Based on our findings, we hypothesize that granzyme levels may fluctuate in the peripheral blood of CL-affected animals and that its lower transcription can be attributed to a particular disease stage. The family of GIMAP enzymes (specifically GIMAP1, GIMAP7, and GIMAP8) was expressed in either the same or different directions in both diseased and exposed sheep. GTPases of immunity-associated proteins are regulators of lymphocyte survival and homeostasis and are linked to inflammatory and autoimmune diseases [59,60]. Recent studies suggest that GIMAP proteins may interact with each other and be involved in the movement of the cellular cargo along the cytoskeletal network [61].

Several authors have emphasized the significance of the cellular immune response and proinflammatory cytokines in limiting bacterial growth and dissemination inside the host in the initial phase of Cp infection [18,62]. Through transcriptome analysis of ovine blood samples, we confirmed the involvement of T lymphocyte-mediated immunity. Gamma-delta T-cell activation predominates in exposed ewes, along with the expression of their specific WC1.1 antigen. γδ T lymphocytes are thought to play crucial roles in immunosurveillance and host defence, particularly against mycobacteria [63]. Glycoprotein WC1.1 is highly represented in ovine peripheral blood and acts as both a γδ T-cell activating coreceptor and pattern recognition receptor [64]. In goat cell culture, only WC1(+) γδ T cells were able to produce proinflammatory cytokines [65]. We also observed increased regulation of innate proinflammatory NF-kappa B signaling in EN group. In contrast, a detailed analysis of the goat transcriptome revealed downregulation of the NF-κB pathway in Cp-infected individuals [21]. However, bovine MAP-infected macrophages produce larger amounts of NF-κB during the initial phase of infection [66]. Many DEGs (approximately 30) in exposed animals are involved in programmed cell death and response to cellular stress as part of normal physiological processes and native immunity. The regulation of their expression tends to be balanced, probably to maintain homeostasis during infection and to prevent excessive cell damage. It has been reported that apoptosis is an essential defence mechanism for host resistance to pathogen invasion. However, Cp infection likely contributes to elevated programmed cell death by apoptosis [21]. Accordingly, death of host cells colonized by corynebacteria has been observed to occur more frequently during the early stage of the pathogenic process [67].

In our study of positive sheep with caseous lymphadenitis, we observed an increase in the regulation of the adaptive immune response through signal transduction and cytokine signalling. The pathogenesis of ovine caseous lymphadenitis was shown to be associated with the production of cytokines at the pyogranuloma level. However, the local cytokine patterns associated with different courses of infection and the time following infection were not distinguished [2,34]. We found both positively and negatively regulated host responses to cytokines but partly lowered cytokine and chemokine production and its regulation. During this phase of the disease, lower expression of chemokine CXCL10, chemokine receptor CXCR6, complement protein C4BPA, complement receptor C3AR1, and anti-inflammatory interleukin IL13 suggests a limited inflammatory response, complement activation, chemotactic activity and lymphocyte migration towards the inflammation site. Likewise, MAP has the ability to downregulate the host complement, hence increasing the chance of survival of bacteria when outside the host cells [53]. Recently, C-X-C motif chemokine ligand 10 has been proposed as a possible biomarker of bovine tuberculosis with the potential to differentiate between latent and active disease. Due to the granulomatous nature of mycobacteria and also related corynebacteria, chemokine recruitment of leukocytes may be a host response to contain invading bacteria, and the restriction of this process by bacteria can subvert the host immune response and establish a latent infection [54]. We suppose that the mechanism behind this phenomenon may be the ability of phospholipase D to reduce leukocyte migration towards sheep serum [68].

Although the transcriptional regulation of the humoral IR appeared to decrease, B-cell activation and B-cell receptor signalling orchestrated by the highly expressed receptor TNFRSF13C (TNF receptor superfamily member 13C) were upregulated in diseased ewes. This receptor is considered the principal receptor required for B-cell-activating factor (BAFF)-mediated mature B-cell survival [69]. The initial immune response to Cp has been shown to include a strong humoral component, but little is known about the kinetics of protective antibodies [62,70]. After experimental infection in goats, the humoral IR showed individual intensity variations and regularly declined over time [20]. The positive regulation of the adaptive immune response based on somatic recombination of immune receptors built from Ig superfamily domains indicates the potential of the host to develop both clonal-specific antibodies against Cp and protective immunological memory during disease progression. We suggest that this can contribute to more effectively eliminating bacteria that can escape from granulomas and spread through the lymphatic system.

We identified significant enrichment of several pathways that may be altered in relation to ongoing infection. In diseased animals, the only statistically significant pathway of Binding and Uptake of Ligands by Scavenger Receptors was recognized. Scavenger receptors bind free extracellular ligands as the initial step in their clearance from the body. They represent an important part of the innate immune defence by acting as bacterial pattern-recognition receptors. The receptors can modulate the inflammatory response as coreceptors of Toll-like receptors, and some intracellular pathogens utilize them to enter host cells [71]. Pathways of ECM proteoglycans and ECM organization containing significant amounts of DEGs were also highly differentially regulated in our research, especially in a negative way. ECM, as a major component of the cellular microenvironment, influences cell adhesion and migration, and its composition is highly heterogeneous and dynamic, being constantly remodelled [72]. Therefore, we suppose the current infection can degrade the essential immune and adhesive functions of ECM and its large protein structures. In accordance with our results, transcriptional profiling of Holstein cows yielded similar findings on ECM modifications in response to natural MAP infection [73].

We also observed substantial expression differences in the immunoregulatory interactions between lymphoid and nonlymphoid cell pathways, which play a critical role in modifying the response of lymphoid cells to pathogenic organisms. Lymphoid cells are able to regulate their location and movement in accordance to their state of activation, and home in on tissues expressing the appropriate complementary ligands [74]. The pathway involves mainly identified NK immune receptors (LILRs and KIRs) and adhesion molecules (i.e., VCAM-1) necessary for antigen presentation by MHC and cooperating molecules (i.e., MICB) and maintaining immunological synapses [75]. It is transcriptionally differentially regulated in both experimental groups facing infection compared to healthy animals.

The p53 pathway was found to be the only one significantly affected in exposed ewes. This finding aligns with the significant transcriptional regulation of genes related to cell death and cell cycle. The p53 transcription factor, known as the tumour suppressor protein, is the key gatekeeper in the cellular response to stress signals. It can help exposed animals respond adequately to stress and prevent DNA damage that is potentially caused by infectious pressure. It also allows cells to recover from damage and survive, preventing premature cell death.

In conclusion, we want to draw attention to some of the limitations or shortcomings of this study that need to be considered. First, to improve the accuracy of the results, it may be necessary to work with a larger sample size during CL infection and exposure. Additionally, due to the incomplete functional annotation of sheep genes, it was necessary to rely on bovine or human homologues, which prevented full interpretation of some discovered genes. Last, the differences in gene expression levels we observed between the investigated ewe groups could be affected by other factors, such as the animal's physiological status, management, or environmental differences.

5. Conclusions

We performed the first comprehensive analysis of the peripheral blood transcriptome in sheep naturally infected with Corynebacterium pseudotuberculosis from an ordinary flock. Bacterial infection interfered with host gene expression, which was reflected preferentially in the transcriptional regulation of genes that modulate the immune response. Two investigated groups (exposed - negative group and diseased - positive group) of adult East Frisian ewes displayed similar transcriptome patterns during longer contact with pathogens but differences in gene expression compared to healthy animals. Hundreds of sheep DEGs associated with the pathogenesis of caseous lymphadenitis, major cellular components, biological processes, and pathways potentially impacted by the disease were identified. Although it is difficult to unambiguously ascribe these differential gene products as contributors to pathogenesis of caseous lymphadenitis, this study has identified genes not previously associated with CL exposure in sheep. The most differentially expressed genes were involved in lymphocyte function, MHC antigen presentation, and extracellular matrix organization. Several major biological processes appear to be triggered in response to CL, including increased activation of the adaptive immune response, altered cytokine production and signalling, and suppression of the inflammatory response. In addition, exposed ewes tended to have notably influenced innate immune processes and programmed cell death. This study suggests that the pathway of antigen binding scavenger receptors might play a substantial role in the protective immune response against Cp. It is not precisely known whether the pathogen, the host, or both drive the changes in gene expression described in Cp-infected sheep. It would be interesting to hypothesize that the differences in expression may provide information about biomarkers that distinguish the preclinical stage of CL; however, further research will be necessary to explore the roles of particular genes in the disease course and the involved mechanisms.