RNA-Seq Analysis of Peripheral Whole Blood from Dairy Bulls with High and Low Antibody-Mediated Immune Responses—A Preliminary Study

Simple Summary Genetic selection for immune response features may be a useful strategy for enhancing animal health. In this study, we identified the DEGs for dairy bulls with high- and low-AMIR. Results revealed that several DEGs were substantially associated with the regulation of locomotion, tissue development, immune response, and detoxification. In addition, the results of the KEGG pathway analysis showed that most DEGs were enriched in pathways related to disease, inflammation, and immune response, providing insights for genomic selection. These findings help to better comprehend the immune response mechanism in dairy bulls. Abstract Enhancing the immune response through breeding is regarded as an effective strategy for improving animal health, as dairy cattle identified as high immune responders are reported to have a decreased prevalence of economically significant diseases. The identification of differentially expressed genes (DEGs) associated with immune responses might be an effective tool for breeding healthy dairy cattle. In this study, antibody-mediated immune responses (AMIRs) were induced by the immunization of hen egg white lysozyme (HEWL) in six Chinese Holstein dairy bulls divided into high- and low-AMIR groups based on their HEWL antibody level. Then, RNA-seq was applied to explore the transcriptome of peripheral whole blood between the two comparison groups. As a result, several major upregulated and downregulated genes were identified and attributed to the regulation of locomotion, tissue development, immune response, and detoxification. In addition, the result of the KEGG pathway analysis revealed that most DEGs were enriched in pathways related to disease, inflammation, and immune response, including antigen processing and presentation, Staphylococcus aureus infection, intestinal immune network for IgA production, cytokine–cytokine receptor interaction, and complement and coagulation cascades. Moreover, six genes (BOLA-DQA5, C5, CXCL2, HBA, LTF, and COL1A1) were validated using RT-qPCR, which may provide information for genomic selection in breeding programs. These results broaden the knowledge of the immune response mechanism in dairy bulls, which has strong implications for breeding cattle with an enhanced AMIR.

DEGs associated with the adaptive immune response in ca le; (2) study the key signaling pathways of the immune response and the changes of various regulatory molecules, thus elucidating the immunomodulatory mechanisms related to disease resistance in dairy cattle and providing new insights for improving animal health and welfare from an immunological perspective.

Materials and Methods
All procedures involving experimental animals were approved by the Animal Welfare Commi ee of Ningxia University (NXUC20210801). All efforts were made to reduce the suffering and discomfort of experimental animals.

Experimental Animals and Immunization Protocol
We randomly selected six healthy dairy bulls aged between 4.9-6.7 years and weighing 1000-1100 kg for this study. The six bulls were Holstein breed with similar genetic merit. They were of various genetic backgrounds, according to the pedigree information. Additionally, principal component analysis (PCA) was carried out to analyze the genetic relationships among these animals based on genotypes, showing no significant clustering according to Supplementary Materials File S1: Figure S1. All experimental animals received no veterinary treatment during the course of the experiment. Immunization protocols were conducted from September 2021 to October 2021, and AMIR indicators were evaluated by modified ELISA based on a previous study [6]. On d0 and d14, the bulls received an intramuscular injection of 0.5 mg of hen egg white lysozyme (HEWL; Sigma L4919, Shanghai, China), 0.5 mg heat-killed Candida albicans (HKCA-#tlrl-hkca, Invi-voGen, San Diego, USA), and 0.5 mg of Quil-A adjuvant (#vacquil, InvivoGen, Toulouse, France) dissolved in 1 mL of phosphate-buffered saline (PBS, pH 7.4). On both sides of the neck, an intramuscular injection of 1 mL was given using a single-use needle. The workflow of this study is summarized in Figure 1.

HEWL Serum Antibody Detection
Blood was collected via jugular venipuncture on d0, d7, d14, d21, and d28 to evaluate serum antibodies to the HEWL. Fetal calf serum served as the negative control, and d21 sera collected after two immunizations served as the positive control. Flat-bottomed 96-well polystyrene plates were coated with 0.05 mg/mL HEWL (Sigma L4919) in CoatingBuffer (0.05 M sodium bicarbonate, pH 9.6). We added 100 µL of diluted coating antibody to each well and incubated it at 4 • C overnight. The coating protein was removed, and the plates were washed five times using 200 µL of ELISA wash solution (50 mM Tris, 0.14 M NaCl, 0.05% Tween 20, pH 8.0). The plates were then blocked with blocking solution (50 mM Tris, 0.14 M NaCl, 0.05% Tween 20, pH 8.0) for 1 h at 37 • C and washed five times again. Sera were diluted in sample/conjugate diluent (50 mM Tris, 0.14 M NaCl, 0.05% Tween 20, pH 8.0) at 1/50 and 1/200 and incubated for 1 h at 37 • C. The plates were washed five times, and the secondary antibody was dissolved in sample/conjugate diluent and incubated at room temperature for 0.5 h followed by five washes; we used HRP-conjugated sheep anti-bovine IgG1 detection antibody (Bethyl Laboratories Inc., Montgomery, AL, USA, Bovine IgG1 ELISA Quantitative Set, Cat. No. E10-116). A total of 100 µL of TMB Substrate Solution was added to each well. The reaction was stopped by adding 100 µL of Stop Solution to each well after the plate had been developed for 15 min in the dark at room temperature. Antibody concentration was measured with an automatic ELISA microplate reader set to an absorbance of 450 nm and 630 nm and expressed in optical density units (OD). The sample to positive (S/P) ratio [6,19] and AMIR 14d (the ratio of optical density values at d14 of the immunization protocol to those at d0) were used to determine the AMIR based on the phenotypic deviation of serum antibody responses, and the formula is shown below. Accordingly, experimental bulls were classified into high-and low-AMIR groups.

RNA-Seq and Transcriptome Quantification
Given that the test animals produced a very small amount of antibodies on d7, transcriptome analysis was not performed on blood samples at this time point. From the blood samples collected on d0, d14, d21, and d28 of the immunization protocol by PAXgene Blood RNA tubes, total RNA was stabilized and extracted as directed by the manufacturer. RNA integrity was determined on a Bioanalyzer 2100 (Agilent Technologies, CA, USA) by a Nano 6000 Assay Kit with an acceptable RIN value above 8.0. Total RNA was used for the library construction following the Illumina NEB Next Ultra™ RNA protocol. The Agilent Bioanalyzer 2100 system was used to evaluate the library's quality. According to the manufacturer's instructions, the TruSeq PE Cluster Kit v3-cBot-HS (Illumina, San Diego, CA, USA) was used to cluster the index-coded sample data on a cBot Cluster Generation System. The library preparations were sequenced on an Illumina Novaseq platform after cluster generation, producing 150 bp paired-end reads.

Analysis of DEGs
The identification of DEGs between high-and low-AMIR groups was based on the expression level of each transcript. An exact test based on quantile-adjusted conditional maximum likelihood (qCML) was completed for DEG screening using the edgeR R package [23]. The threshold, fold change ≥2 and <0.05 for the alpha of the false discovery rate (FDR), was set according to a previous transcriptome study [24,25]. Furthermore, the similarity level of each sample was determined using hierarchical clustering of DEGs.
To generalize the functions and pathways of DEGs, the sequences were functionally annotated and classified by comparing them to the gene ontology (GO) and the Kyoto Encyclopedia of Genes and Genomes (KEGG) databases. By using the "clusterProfiler" package [26], which corrects for gene length bias, KEGG pathways and GO terms were enriched. DEGs with a p-value less than 0.05 were found to be highly enriched in GO terms and KEGG pathways.
The STRING (http://string-db.org/, accessed on 20 March 2023) database was used to investigate a PPI network among the DEGs to better comprehend and foresee the biological activity of the detected DEGs based on GO and KEGG pathway enrichment analysis. The cutoff threshold was set at a confidence score > 0.7 and a p < 0.05.
The RNA-Seq data were confirmed using RT-qPCR on six randomly selected genes. Primer 5.0 was utilized to design the primers, and General Biosystems Co., Ltd. (Anhui, China) was employed to synthesize the primers. Total RNA was reversed to cDNA by the TaKaRa Reverse Transcription Kit (TaKaRa, Dalian, China) according to the instructions provided by the manufacturer. Quantitative PCR was carried out on a 7500 Fast Real-Time PCR System (Thermo Fisher Scientific, Waltham, MA, USA) with the TB Green ® Premix Ex TaqTM Kit (TaKaRa). GAPDH was selected as an internal control. The2 −∆∆CT method was used to examine all data from three biological replicates for each sample [27].

Gene Expression Pattern Profiling
Fifty-six genes with significant differences between the high-and low-AMIR groups were selected for time series analysis. These genes are related to immune response, including the BOLA family, the CXCL chemokine family, the hemoglobin family, the interleukin family, the tumor necrosis factor family, etc. The Mfuzz R package's fuzzy c-means algorithm was used to profile DEGs based on their expression patterns. For each time point (d0, d14, d21, and d28) of the immunization protocol, the average FPKM value for each gene was utilized as the input. Following standardization, each gene was assigned to a distinct cluster based on its membership value.

Experimental Animals and Grouping
The obtained HEWL antibody titers on d0, d7, d14, d21, and d28 of immunization are shown in Figure 2A. The antibody response showed a gradually increasing trend within 21 days of immunization. The initial antibody titer was close to zero in the peripheral blood, indicating that the selection of experimental animals was reasonable. After d7 immunization, only a slight increase in the antibody response was induced to the HEWL immunogen. Subsequently, more antibodies were produced from d14. The antibody titer was highest at d28 and was significantly higher than that on d0 and d7 (p < 0.05).
The six animals were separated into high-and low-AMIR groups based on the S/P ratio and AMIR14d. The two evaluation criteria obtained similar results. For the high-AMIR group, there were three identified animals, H1, H2, and H3, whereas two animals, L1 and L2, were part of the low-AMIR group. The sixth individual was discarded because the serum antibody detection was not successful.
We also identified the trend of changes in the titer of antibody against HEWL in the two groups ( Figure 2B,C). There was an upward trend in both groups. In the high-AMIR group, except for the significantly higher antibody titer on d28 compared to d0 and d7, there were no significant differences among samples at any other time points, and the Animals 2023, 13, 2208 6 of 16 antibody titer on d14 reached a relatively high level ( Figure 2B). In the low-AMIR group, the antibody titers on d21 and d28 were remarkably higher than any other time points, and the antibody titer at d14 remained at a low level ( Figure 2C). The AMIR14d values and S/P of each sample are shown in Figure 3. The two indicators, AMIR 14d and S/P, for animals in the high-AMIR group were significantly higher than those of the low-AMIR group. there were no significant differences among samples at any other time points, and the antibody titer on d14 reached a relatively high level ( Figure 2B). In the low-AMIR group, the antibody titers on d21 and d28 were remarkably higher than any other time points, and the antibody titer at d14 remained at a low level ( Figure 2C). The AMIR14d values and S/P of each sample are shown in Figure 3. The two indicators, AMIR14d and S/P, for animals in the high-AMIR group were significantly higher than those of the low-AMIR group.  of the immunization protocol to those at d0; (B) sample-to-positive (S/P) ratio of HEWL antibody level from serum IgG1 at d14 of the immunization protocol. Red pillars represent high-AMIR animals, and blue pillars represent low-AMIR animals.

Differentially Expressed Genes
For each sample, paired-end transcriptome sequencing produced at least 6.0 Gbp of raw data. According to the trimmed results, the Q30 (%) of clean data for all samples collected on d14 was greater than 93.87%, and the GC contents ranged from 52.40 to 53.22% (Table 1). Approximately 89.97-93.12% of clean reads were uniquely aligned to the reference genome. A detailed description of the RNA-seq data for blood samples on d0, d21, and d28 is provided in Supplementary Materials File S1: Table S1.   there were no significant differences among samples at any other time points, and the antibody titer on d14 reached a relatively high level ( Figure 2B). In the low-AMIR group, the antibody titers on d21 and d28 were remarkably higher than any other time points, and the antibody titer at d14 remained at a low level ( Figure 2C). The AMIR14d values and S/P of each sample are shown in Figure 3. The two indicators, AMIR14d and S/P, for animals in the high-AMIR group were significantly higher than those of the low-AMIR group.  of the immunization protocol to those at d0; (B) sample-to-positive (S/P) ratio of HEWL antibody level from serum IgG1 at d14 of the immunization protocol. Red pillars represent high-AMIR animals, and blue pillars represent low-AMIR animals.

Differentially Expressed Genes
For each sample, paired-end transcriptome sequencing produced at least 6.0 Gbp of raw data. According to the trimmed results, the Q30 (%) of clean data for all samples collected on d14 was greater than 93.87%, and the GC contents ranged from 52.40 to 53.22% (Table 1). Approximately 89.97-93.12% of clean reads were uniquely aligned to the reference genome. A detailed description of the RNA-seq data for blood samples on d0, d21, and d28 is provided in Supplementary Materials File S1: Table S1.

Differentially Expressed Genes
For each sample, paired-end transcriptome sequencing produced at least 6.0 Gbp of raw data. According to the trimmed results, the Q30 (%) of clean data for all samples collected on d14 was greater than 93.87%, and the GC contents ranged from 52.40 to 53.22% (Table 1). Approximately 89.97-93.12% of clean reads were uniquely aligned to the reference genome. A detailed description of the RNA-seq data for blood samples on d0, d21, and d28 is provided in Supplementary Materials File S1: Table S1. PCA and cluster analysis were performed on the gene expression levels of blood samples on d0, d14, d21, and d28 of the immunization protocol. The clustering results of blood samples on d0 were scattered and had no obvious regularity (Supplementary Materials File S1: Figure S2). The high-and low-AMIR groups differed significantly in RNA-seq counts for blood samples collected on d14 of the immunization protocol based on PCA and clustering structure (Figure 4). According to the cluster results of blood samples on d21 and d28 of immunization (Supplementary Materials File S1: Figures S3  and S4), samples could not be clustered regularly. Thus, blood samples from the high-and low-AMIR groups on d14 were used for DEG analysis. PCA and cluster analysis were performed on the gene expression levels of blood samples on d0, d14, d21, and d28 of the immunization protocol. The clustering results of blood samples on d0 were sca ered and had no obvious regularity (Supplementary Materials File S1: Figure S2). The high-and low-AMIR groups differed significantly in RNA-seq counts for blood samples collected on d14 of the immunization protocol based on PCA and clustering structure (Figure 4). According to the cluster results of blood samples on d21 and d28 of immunization (Supplementary Materials File S1: Figures S3 and S4), samples could not be clustered regularly. Thus, blood samples from the high-and low-AMIR groups on d14 were used for DEG analysis. DEGs were identified using the qCML method. Based on the threshold level (fold change ≥ 2 and FDR < 0.05), 1109 genes were identified as differentially expressed ( Figure  5). Compared with the high-AMIR group, 637 genes were upregulated, while 472 genes were downregulated in the low-AMIR group.  DEGs were identified using the qCML method. Based on the threshold level (fold change ≥ 2 and FDR < 0.05), 1109 genes were identified as differentially expressed ( Figure 5). Compared with the high-AMIR group, 637 genes were upregulated, while 472 genes were downregulated in the low-AMIR group. PCA and cluster analysis were performed on the gene expression levels of blood samples on d0, d14, d21, and d28 of the immunization protocol. The clustering results of blood samples on d0 were sca ered and had no obvious regularity (Supplementary Materials File S1: Figure S2). The high-and low-AMIR groups differed significantly in RNA-seq counts for blood samples collected on d14 of the immunization protocol based on PCA and clustering structure (Figure 4). According to the cluster results of blood samples on d21 and d28 of immunization (Supplementary Materials File S1: Figures S3 and S4), samples could not be clustered regularly. Thus, blood samples from the high-and low-AMIR groups on d14 were used for DEG analysis. DEGs were identified using the qCML method. Based on the threshold level (fold change ≥ 2 and FDR < 0.05), 1109 genes were identified as differentially expressed ( Figure  5). Compared with the high-AMIR group, 637 genes were upregulated, while 472 genes were downregulated in the low-AMIR group.   HBE2, and HBA). Of all GO-MF terms, four GO terms were related to receptor activity, namely receptor ligand activity, signaling receptor activator activity, receptor regulator activity, and signaling receptor binding. The key genes enriched in GO-MF terms were TGFBP2, IL34, STC1, INHA, CXCL2, CCK, and GRO1 (Supplementary Materials File S2). The main purpose of our research is to explore biological functions, so we focus on biological processes, and the top 20 significant GO-BP terms are shown in Figure 6.  HBE2, and HBA). Of all GO-MF terms, four GO terms were related to receptor activity, namely receptor ligand activity, signaling receptor activator activity, receptor regulator activity, and signaling receptor binding. The key genes enriched in GO-MF terms were TGFBP2, IL34, STC1, INHA, CXCL2, CCK, and GRO1 (Supplementary Materials File S2). The main purpose of our research is to explore biological functions, so we focus on biological processes, and the top 20 significant GO-BP terms are shown in Figure 6. The KEGG enrichment analysis of DEGs identified 42 highly enriched pathways, mostly related to disease, inflammation, and immune response. Sixteen of all significantly enriched pathways (p < 0.05) were associated with various diseases. Graft-versus-host disease, ECM-receptor interaction, antigen processing and presentation, Staphylococcus aureus infection, intestinal immune network for IgA production, viral myocarditis, cytokine-cytokine receptor interaction, autoimmune thyroid disease, complement and coagulation cascades, allograft rejection, and TGF-beta signaling pathways were found to be enriched in canonical pathways related to immune response and inflammation ( Figure 7, Table 2  The KEGG enrichment analysis of DEGs identified 42 highly enriched pathways, mostly related to disease, inflammation, and immune response. Sixteen of all significantly enriched pathways (p < 0.05) were associated with various diseases. Graft-versus-host disease, ECM-receptor interaction, antigen processing and presentation, Staphylococcus aureus infection, intestinal immune network for IgA production, viral myocarditis, cytokine-cytokine receptor interaction, autoimmune thyroid disease, complement and coagulation cascades, allograft rejection, and TGF-beta signaling pathways were found to be enriched in canonical pathways related to immune response and inflammation ( Figure 7, Table 2      STRING analysis was performed to investigate the DEG potential interaction network (Figure 8). Among the upregulated genes, DCN, COL1A1, COL1A2, COL3A1, HBA, F2, and BOLA-DQA5 were situated in the network's core and connected to many other DEGs; for the downregulated genes, the key points included FBN1 and ITGA9, which linked to more genes. Additionally, not all DEGs displayed a link with others because their functions were either unclear or unrelated to one another, which were not included in the analysis.
Animals 2023, 13, x FOR PEER REVIEW 10 of 16 STRING analysis was performed to investigate the DEG potential interaction network ( Figure 8). Among the upregulated genes, DCN, COL1A1, COL1A2, COL3A1, HBA, F2, and BOLA-DQA5 were situated in the network's core and connected to many other DEGs; for the downregulated genes, the key points included FBN1 and ITGA9, which linked to more genes. Additionally, not all DEGs displayed a link with others because their functions were either unclear or unrelated to one another, which were not included in the analysis. The interaction was carried out with a confidence score of 0.7. Nodes represent genes, and lines between nodes refer to edges indicating various types of interactions, which are denoted by different colors and defined by the legends in the figure; stand-alone nodes without edges were removed.
RT-qPCR was performed to validate the gene expression identified by RNA-seq. The six genes selected were from various gene families. Thus, the relationship between each gene family and the immune response was further confirmed. Additionally, the results demonstrated that the expression levels of the six genes generally agreed with the RNA-Seq results based on read counts (Figure 9). The interaction was carried out with a confidence score of 0.7. Nodes represent genes, and lines between nodes refer to edges indicating various types of interactions, which are denoted by different colors and defined by the legends in the figure; stand-alone nodes without edges were removed.
RT-qPCR was performed to validate the gene expression identified by RNA-seq. The six genes selected were from various gene families. Thus, the relationship between each gene family and the immune response was further confirmed. Additionally, the results demonstrated that the expression levels of the six genes generally agreed with the RNA-Seq results based on read counts (Figure 9).
RT-qPCR was performed to validate the gene expression identified by RNA-seq. The six genes selected were from various gene families. Thus, the relationship between each gene family and the immune response was further confirmed. Additionally, the results demonstrated that the expression levels of the six genes generally agreed with the RNA-Seq results based on read counts (Figure 9).

Expression Pattern of Selected DEGs
Soft clustering was performed on the expression of the selected DEGs at different time points for the high-AMIR individuals. Since there were only a few genes, the number of clusters was manually set in this study. Finally, all genes were classified into five clusters. As shown in Figure 10, the expression of genes in Cluster 1 (IL5RA, RARRES2, F2, SPRY4, and ITGB4) increased after the immunization, maintaining a high level between d14 and d21, and then decreased slowly. The genes in Cluster 2 (IL34 and CCK) reached peak expression at d14 and then declined. In Cluster 3 (LTF, PGLYRP1, ANGPT2, C5, BOLA-DQA5, BOLA-DQB, HBA, HBE2, HBE4, HBM, HBG, and CXCL2), genes were downregulated from d0 to d14 of the immunization protocol, subsequently remaining relatively stable. The genes in Cluster 4 (SERPING1 and IL36A) showed irregular changes before d21, but were then upregulated. The expression pattern of genes in Cluster 5 (IL1RL2, COL1A1, and DCN) reached a peak at d21 after immunization. The gene list from each cluster and the FPKM values of each gene at different time points are presented in Supplementary Materials File S4. Of the six genes that were selected for RT-qPCR, BOLA-DQA5, C5, CXCL2, HBA, and LTF represented Cluster 3, while COL1A1 represented Cluster 5.

Expression Pa ern of Selected DEGs
Soft clustering was performed on the expression of the selected DEGs at different time points for the high-AMIR individuals. Since there were only a few genes, the number of clusters was manually set in this study. Finally, all genes were classified into five clusters. As shown in Figure 10, the expression of genes in Cluster 1 (IL5RA, RARRES2, F2, SPRY4, and ITGB4) increased after the immunization, maintaining a high level between d14 and d21, and then decreased slowly. The genes in Cluster 2 (IL34 and CCK) reached peak expression at d14 and then declined. In Cluster 3 (LTF, PGLYRP1, ANGPT2, C5, BOLA-DQA5, BOLA-DQB, HBA, HBE2, HBE4, HBM, HBG, and CXCL2), genes were downregulated from d0 to d14 of the immunization protocol, subsequently remaining relatively stable. The genes in Cluster 4 (SERPING1 and IL36A) showed irregular changes before d21, but were then upregulated. The expression pa ern of genes in Cluster 5 (IL1RL2, COL1A1, and DCN)

Experimental Design and Samples Information
Research has suggested that genetic improvement is dominated by the sire-daughter path compared to the dam-daughter path [28]. Another study also supports the argument that the sire has a more significant effect on the genetic improvement of dairy herds than the dam [29], and high immune response sires reduce the disease incidence in large commercial dairy populations in North America [30]. Likewise, the Immunity+™ sire line,

Experimental Design and Samples Information
Research has suggested that genetic improvement is dominated by the sire-daughter path compared to the dam-daughter path [28]. Another study also supports the argument that the sire has a more significant effect on the genetic improvement of dairy herds than the dam [29], and high immune response sires reduce the disease incidence in large commercial dairy populations in North America [30]. Likewise, the Immunity+™ sire line, which was identified with superior immune responsiveness by the Semex Alliance utilizing HIR Technology, may be expected to produce daughters with superior immunity or disease resistance [10]. Therefore, Holstein bulls as mature sires were selected as research subjects in the present study to contribute knowledge to the selection process of high-resistance breeding bulls.
Although there were only a few experimental animals, the results are encouraging and could be used to support future research based on a larger population. Because dairy cattle are bred by artificial insemination, there are a relatively small number of dairy bulls. Some bulls, especially Holstein bulls, may be closely related, which increased the difficulty of sample selection. In order to increase the reliability of the test results, six Holstein bulls were initially used to perform the experiment. Unfortunately, one bull sample was discarded due to uncertainty. However, serum antibody levels in the three high-AMIR animals were significantly higher than those in the low-AMIR group (Figure 3), indicating that the experimental animals were representative. Additionally, the selection of mature sires that have produced offspring may facilitate subsequent validation in the cow population. In addition, we referred to a published paper that also used Holstein dairy cattle as experimental animals, but only five samples were analyzed [31].
Given the advantage that the RNA-Seq technique can detect a higher proportion of DEGs than gene chips, it has been applied to map key genes for susceptibility or resistance to certain diseases, including ketosis and mastitis [12,13]. Instead of specific disease resistance, the current study, for the first time, reported RNA-Seq results for overall disease resistance and the adaptive immune response in dairy bulls. The ELISA results revealed that all animals showed an almost similar antibody titer, which was close to zero at d0. After the same treatment protocol, each individual produced an antibody reaction, and there was a significant difference in S/P ratio or AMIR14d based on antibody titer at d14 between the two groups. Therefore, the DEGs found in the transcriptome study were mostly attributed to the immune response of the animals to the HEWL immunogen. The method suggested here could detect differential gene expression profiles that may be associated with general disease resistance.

Differential Genes Screening and Functional Enrichment
The study of GO functional enrichment indicated upregulated DEGs related to the regulation of cell motility, tissue development, humoral immune response, and detoxification. LTF, PGLYRP1, C5, SERPING1, and RARRES2 were significantly enriched in the items related to immune response, namely the humoral immune response and antimicrobial humoral response. Lactotransferrin (LTF) is a gene in the transferrin family whose protein product is found in neutrophil secondary granules. It can regulate macrophage activity and stimulate lymphocyte synthesis by binding macrophages and lymphocytes, promoting lymphocyte proliferation by inducing specific immune responses. It has been proposed as a potential biomarker of mastitis susceptibility/resistance [32], subclinical endometritis (SCE) [33], and Mycobacterium avium subspecies-a paratuberculosis infection in dairy cattle [34]. Bovine PGLYRP1, found in naive neutrophil dense/large granules and neutrophil phagolysosomes [35], was the first member of PGRP identified for its bactericidal activity in vitro. Polymorphism in the gene was shown to be associated with mastitis resistance, and the H3H3 haplotype was regarded as a valuable genetic indicator of combination haplotypes for developing mastitis resistance in the Chinese Holstein [36]. Both C5 and SERPING1 are involved in the regulation of the complement system. During complement activation, C5 breaks down to produce C5a, which is a major allergic toxin involved in neutrophil chemotaxis and the release of pro-inflammatory cytokines [37]. SERPING1 encodes the C1 inhibitor, which is important in inhibiting complement component 1 (C1) and may be associated with the conventional complement activation pathway [38]. Chemerin (RARRES2) is an endogenous leukocyte chemoattractant and is frequently downregulated in several tumor types when compared to normal tissue [39]. This corresponds to the downregulated RARRES2 gene expression in low responders found in the present study.
DEGs associated with detoxification were mainly members of the hemoglobin (Hb) family, including HBE4, HBG, HBM, HBE2, and HBA. Long regarded as inert oxygen carriers, red blood cells are emerging as potential modulators of the innate immune response [40]. Hb and heme reportedly have immunomodulatory effects and could produce antimicrobial reactive oxygen species (ROS) to defend against invading microorganisms [40]. A recent study also indicated that Hb had immunomodulatory, signal transduction, and bactericidal effects [41]. Moreover, there is growing evidence that Hb is not only an oxygen transporter but also contributes substantially to the innate immune system. Bovine Hb is an animal protein that has been reported as a source of bioactive peptides. Nedjar-arroume et al. [42] extracted and discovered four antimicrobial peptides from bovine Hb that displayed antibacterial activity against Micrococcus luteus A270, Listeria innocua, Escherichia coli, and Salmonella enteritidis. In the present study, Hb family genes were considerably enriched in GO terms related to detoxification, further confirming the important role of Hb in immune regulation.
Among the KEGG pathways enriched in this study, nine pathways were consistent with the results obtained by the GWAS study conducted to assess overall adaptive immune response in cattle [43]; these included Graft-versus-host disease, antigen processing and presentation, systemic lupus erythematosus, intestinal immune network for IgA production, viral myocarditis, autoimmune thyroid disease, allograft rejection, asthma, and type I diabetes mellitus. In one study, the relationship between immune response and mastitis was examined, and it was found that high-AMIR cows had a significantly decreased incidence of clinical mastitis compared to average and low-AMIR cows [8]. Meanwhile, S. aureus has been considered one of the most important mastitis pathogens [44,45]. In the present study, S. aureus infection was identified as a KEGG pathway significantly related to the AMIR, in agreement with the research above.
The cytokine-cytokine receptor interaction is critical to health during immune and inflammatory responses to diseases [46]. In the present study, many DEGs were found to be enriched in the cytokine-cytokine receptor interaction pathway, including the interleukin family and CXC chemokine family genes. In addition, IL-34 expression has been shown to be dysregulated in a variety of immune-inflammatory diseases, infections, and metabolic and neurological disorders [47]. It has also been noted that there is a correlation between changes in IL-34 content and disease parameters in a variety of pathological conditions, including rheumatoid arthritis, systemic lupus erythematosus, heart failure, and viral infections. Genome-wide association analysis on cows showed that IL-34 was significantly associated with both clinical mastitis and SCS [48]. Also, the IL-34 level has been viewed as a helpful biomarker for predicting the progression of disease. CXCL2 is an inflammatory factor that directly participates in biological processes such as immune regulation, inflammatory responses, and angiogenesis by binding to the G protein-coupled receptor CXCR. Sharifi et al. [49] showed that CXCL2 could be used to distinguish mastitis from healthy samples with an accuracy of 83.85%.
Bovine lymphocyte antigen (BOLA) genes were also significantly enriched in multiple pathways, according to the KEGG pathway enrichment analysis. Known as the major histocompatibility complex (MHC) of cattle, BOLA genes are located on chromosome 23 in cattle [43] and have been intensively investigated by researchers because of their antigen presentation function. In addition, polymorphism in BOLA genes was found to be closely associated with bovine disease resistance and immunological traits [50]. According to the structure and function of the encoding products, BOLA genes can be classified into three categories: class I, class II, and class III. Among the DEGs identified significantly by KEGG analysis in the present study, BOLA-DQA5 and BOLA-DQB are class II genes. In cattle, polymorphism in class IIa genes affects the magnitude and epitope specificity of antigenspecific T-cell responses to FMD virus peptides. As indicated in Figure 5, the expression levels of BOLA-DQA5 and BOLA-DQB genes were extremely significantly different between high-and low-AMIR groups, confirming the importance of this gene in immune regulation. This is in accordance with previous studies.
The incorporation of DEGs from this study into animal genetic improvement programs such as selective breeding and genetic manipulation might improve cattle immunity.

Conclusions
In the present study, several critical DEGs were found to be associated with the high and low general AMIRs in dairy bulls, of which, six important genes (BOLA-DQA5, C5, CXCL2, HBA, LTF, and COL1A1) were validated by RT-qPCR. The DEGs, GO terms, and KEGG pathways identified in this study could help in further exploration of the underlying mechanisms of immune responses in dairy cattle. This work addressed a knowledge gap in the transcriptome differential analysis in the peripheral blood transcriptome of high-and low-AMIR bulls, and the identification of key genes could provide useful information for the breeding of enhanced immune response sires.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/ani13132208/s1. Figure S1: Principal component analysis of genotyped animals showing the diverse genetic background of the bulls; Table S1: An overview of RNA-Seq data for each sample on d0, d21, and d28; Figure S2: PCA, cluster dendrogram, and correlation matrix were performed for blood samples collected at d0 of the immunization protocol. Red points in the dot plot represent high-AMIR samples. Blue points in the dot plot represent low-AMIR samples; Figure S3: PCA, cluster dendrogram, and correlation matrix were performed for blood samples collected at d21 of the immunization protocol. Red points in the dot plot represent high-AMIR samples. Blue points in the dot plot represent low-AMIR samples; Figure S4: PCA, cluster dendrogram, and correlation matrix were performed for blood samples collected at d28 of the immunization protocol.

Institutional Review Board Statement:
The animal study protocol was approved by the Animal Welfare Committee of Ningxia University (NXUC20210801).