Control of Mycobacterium avium subsp. paratuberculosis load within infected bovine monocyte-derived macrophages is associated with host genetics

The genetic loci influencing individual resistance to Mycobacterium avium subsp. paratuberculosis (MAP) infection are still largely unknown. In the current study, we searched for genetic loci associated with resistance to MAP infection by evaluating the performance of monocyte-derived macrophages (MDMs) isolated from the peripheral blood of 75 healthy Holsteins cows and infected ex vivo with MAP. Bacterial load (log colony-forming units, log CFUs) within MDMs was quantified at 2 h and 7 days p. i. using a BACTEC MGIT 960 instrument. In addition, the expression levels of some genes with important roles in the innate immune response including epiregulin (EREG), complement component C3 (C3), galectin-9 (Gal9), and nitric oxide (NO-) were measured in the supernatant of the infected cells. DNA from peripheral blood samples of the animals included in the study was isolated and genotyped with the EuroG MD bead Chip (44,779 single nucleotide-polymorphisms, SNPs). Linear mixed models were used to calculate the heritability (h2 ) estimates for each indicator of MDM performance, MAP load within MDMs and EREG, C3, Gal9, and NO-expression. After performing a genome-wide association study, the only phenotypes that showed SNPs with a significant association were the bacterial load within MDMs at 2 h (h2 = 0. 87) and 7 days (h2 = 0.83) p.i. A total of 6 SNPs, 5 candidate genes, and one microRNA on the Bos taurus chromosomes BTA2, BTA17, BTA18, and BTA21 were associated with MAP load at 2 h p.i. Overlap was seen in two SNPs associated with the log CFUs at 2 h and 7 d p.i. The identified SNPs had negative regression coefficients, and were, therefore, associated with a low bacterial load within MDMs. Some of the identified SNPs were located within QTLs previously associated with longevity, reproductive, and udder health traits. Some of the identified candidate genes; Oxysterol Binding Protein Like 6, Cysteine and Serine Rich Nuclear Protein 3, and the Coiled-Coil Domain Containing 92 regulate cellular cholesterol trafficking and efflux, apoptosis, and interferon production, respectively. Taken together, our results define a heritable and distinct immunogenetic profile in MAP-infected macrophages designed to limit bacterial load early after infection.


Introduction
Bovine paratuberculosis (PTB) or Johne´s disease (JD) is a chronic enteritis of domestic and wild ruminants caused by Mycobacterium avium subsp. paratuberculosis (MAP). PTB is a major problem for animal health that compromises animal welfare and causes important economic losses to the dairy industry (1). In Europe and North America, PTB is endemic in dairy cattle, with herd prevalence estimates higher than 50% (2). PTB-associated losses include increased susceptibility to other diseases, increased somatic cell counts, increased incidence of clinical mastitis, reduced fertility, reduced milk production, costs of testing, and involuntary culling of cows (3). Animals are infected early in life through the fecal-oral route but clinical onset appears when animals are 18 months or older (4). Once ingested, MAP reaches the jejune and ileum and crosses the intestinal mucosa by binding to fibronectin b1 receptors present on M cells located in Peyer's patches (5,6). In the submucosa, MAP is phagocytized by sub-epithelial macrophages. Within infected macrophages, MAP has developed several survival mechanisms such as preventing the maturation and acidification of the phagosome and its fusion with the lysosome and preventing the presentation of antigens to the immune system (7). MAP has a zoonotic potential and it has been detected in samples of patients with Crohn´s disease (CD), ulcerative colitis, and idiopathic inflammatory bowel disease (IBD)-associated colorectal cancer (8)(9)(10)(11). Inactivated vaccines against MAP interfere with diagnostic tests for bovine tuberculosis and therefore are not allowed in many countries. One solution to control this disease is the identification of genetic loci determining variation in immunerelated traits that could be used to select PTB-resistant animals.
PTB is a multifactorial disease that arises as the result of the interaction of genetic, environmental, and microbial factors which affect the various disease outcomes. According to their extension in the intestine, cellular infiltrate, and amount of MAP; PTBassociated lesions were classified as focal, multifocal, and diffuse (diffuse paucibacillary or lymphoplasmacytic, diffuse intermediate, and diffuse multibacillary or histiocytic) (12, 13). Our research group has explored very recently the associations between PTBassociated pathology and host genetics (14). Overall, a total of 380 SNPs were found associated (FDR ≤ 0.05; P ≤ 5 10 -7 ) with antemortem (serum ELISA) and post-mortem (tissue PCR and culture and histopathology) PTB diagnostic definitions in a common set of Spanish Holstein cattle (N = 983) using whole-genome sequence data (WGS) data (14)(15)(16). However, we were unable to identify SNPs associated with resistance to MAP infection.
Reductionist approaches investigate a host biological subsystem such as a key cellular function whose performance is the criterion for the classification of the population (17). Recently, a strong effect of host genetics on the in vitro nitric oxide (NO -) production of monocyte-derived macrophages (MDMs) in response to two common bacterial pathogens of dairy cattle, Escherichia coli and Staphylococcus aureus, has been described (18). Macrophages are important host defense cells against MAP infection. Since MAPinfected macrophages' performance can be better analyzed under a controlled environment, the probability of identifying resistant animals is much higher using ex vivo macrophage models in comparison to holistic approaches based on field studies. In addition, it is generally hypothesized that mycobacteria growth in macrophages in vitro predicts the in vivo risk of disease or infection (19,20). Moreover, the direct mycobacteria growth inhibition assay (MGIA) has demonstrated cross-species potential (21). Recently, we have demonstrated that MAP-infected MDMs constitute a very useful model of infection to validate the functional impact of specific expression quantitative trait loci (cis-eQTLs) potentially implicated in susceptibility to MAP infection (22). Cis-eQTLs are genetic variants typically located in gene regulatory regions that alter gene expression in an allele-specific manner. Cis-eQTLs are functional links between genomic variants, gene expression, and ultimately phenotype. In the current study, we hypothesize that MAP-infected MDMs might help to identify resistant individuals that respond stronger than others to MAP infection and, as result, reduce the intracellular MAP burden significantly.
The measure of indicators of MDMs' performance in response to MAP infection is another strategy for the identification of resistant cattle. To date, however, no serum biomarkers have been associated with resistance to MAP infection. In our previous study (22), the integration of gene expression data (RNA-Seq) and genotypes (54,609 SNPs per animal) from a cohort of cows naturally exposed to MAP allowed the identification of 192 cis-eQTLs associated with the expression of 145 genes in peripheral blood samples. Four of these cis-eQTLs regulated the expression of several genes with important roles in the innate immune response; the complement component C3 (C3), epiregulin (EREG), galectin-9 (Gal9), and nitric oxide synthase (NOS1). In the current study, we tested whether the measurement of MAP load within infected MDMs and the expression in the supernatants of MAP-infected MDMs of several genes with important roles in the innate immune response could be used as correlates of MDMs performance. Our objective was to identify genetic loci associated with resistance to MAP infection using ex vivo MAP-infected MDMs from Spanish Holstein cows. For this purpose, we performed a genome-wide association study (GWAS) to identify SNPs and candidate genes associated with the MAP load within infected MDMs and with the expression of Gal9, ERG, C3, and NOin the supernatants of infected MDMs. Since there is evidence that some allelic variants may contribute to resistance to multiple pathogens, the identified SNPs and candidate genes were compared with QTLs and candidate genes for other bovine diseases and health and reproductive traits. In addition, the candidate genes identified in our study were compared with human candidate genes previously identified for CD, IBD, and colorectal cancer. The workflow of the study is presented in Figure 1.

Animals and blood samples
The reference population consisted of 75 Spanish Holstein cows from a single farm located in the Basque Country. In 2022, only five cows from this herd had an ELISA positive result and all the animals had a negative fecal PCR result. Only adult cows (2 years or older, 4.5 years mean age) were included in the study. All the animals included in this study were registered with the Spanish Federation of Holstein cattle (CONAFE; www.conafe.com.) and had negative fecal PCR results when blood samples were collected and in subsequent annual samplings. In addition, the animals used in this study were not diagnosed with any other pathogen. Blood samples were collected in groups of 16 per sampling day. The validation population consisted of 16 cows from the same farm.

Mycobacterial growth inhibition assays
MGIA were performed as previously described (22)(23)(24). Briefly, fifteen milliliters of peripheral blood were drawn from the tail vein of healthy Holstein cows into heparinized Vacutainer tubes (Becton, Dickinson and Company, Sparks, MD, USA) and diluted 1:2 in Hanks balanced salt solution (HBSS). Leucosep tubes were filled with 15 ml of Ficoll-Paque (1.084 g/cm 3 ) (GE Healthcare, Uppsala, Sweden) and centrifuged at 1,000 rpm for 30 seconds at room temperature. Subsequently, the diluted blood was overlaid on the top of the Ficoll-Paque and centrifuged at 800 g for 15 minutes at room temperature. The plasma layer was removed and the cell interphase containing peripheral blood mononuclear cells (PBMCs) was collected and transferred to a clean tube. PBMCs were washed twice in HBSS and centrifuged at 400 g for 10 minutes to remove platelets. Supernatants were aspirated and the purified PBMCs were resuspended in RPMI-1640 supplemented with 20 mM Lglutamine, 10% heat-inactivated bovine serum (Lonza, Spain), 100 U ml-1 penicillin G, and 100 mg ml-1 streptomycin sulfate (Lonza, Spain). PBMCs were cultured at a concentration of 1 x 10 6 cells/ml into 24-well plates and incubated at 37°C in a humidified 5% CO 2 incubator for 2 h. Non-adherent cells were removed by washing, and adherent cells were incubated in fresh medium for 7 days at 37°C to allow differentiation to MDMs. Differentiated MDMs were inoculated in triplicate with a single-cell suspension of MAP K10 strain at a multiplicity of infection (MOI) of 10:1 (bacteria:cells). After 2 h, the supernatant was removed, and the cells were washed twice with HBSS to remove extracellular bacteria. Infected MDMs were lysed at this time (2 h p. i.) or were cultured at 37°C for 7 days in fresh medium. At each time point, the supernatant was aspirated and infected MDMs were lysed by vigorous pipetting with 0.5 ml of 0.1% Triton X-100 (Sigma-Aldrich) in sterile water for 10 min.   Study design. The approach starts from the measures of MAP load within MAP infected-MDMs and the C3, EREG, Gal9, and NOexpression levels in the supernatants of MAP-infected-MDMs. These phenotypes are combined with the genotypes of the animals in a GWAS study for the calculation of variance components and heritability (h 2 ) estimates and the identification of associated SNPs and candidate genes.
MGIT 960. The tubes were incubated at 37 ± 2°C for up to 42 days in a BACTEC MGIT 960 instrument (Becton, Dickinson, and Company). The earliest instrumental indicator of positivity (i.e., time to detection [TTD]) for each tube was recorded. The predicted number of bacteria in each positive tube was calculated using standard curves that relate TTD (in days) to the estimated log CFUs (25). For the standard curves, a ten-fold dilution series of a MAP K10 strain cellular suspension (McFarland= 1) were prepared and the TTDs of 100 μl of each dilution in the BACTEC MGIT 960 system were detected. MAP load (log CFUs) were plotted versus TTDs and the mathematical equation was used to determine the estimated log CFUs for each sample. The log CFU ratios were calculated by dividing the estimated log CFUs at day 7 by that at 2 h p. i. MAP survival indexes were calculated according to  according to the manufacturer's instructions. The assay has an optimal range of 20-500 pmol nitrite, making it up to 50 times more sensitive than colorimetric methods utilizing the Griess reagent. Nitrates are analyzed after quantitative conversion to nitrites through enzymatic reduction; used in this manner, the assay provides effective quantitation of NO -. Briefly, the supernatants from the MAP-infected MDMs after 2 h of infection (40 μl) were brought to a volume of 50 μl by adding 10 μl of H 2 O and placed into 96-well black plates in triplicate. Subsequently, 100 μl of quantitation reagent were added and mixed by pipetting. The plate was incubated for 10 minutes at room temperature and 5 μl of quantitation developer were added to each well and gently mixed by pipetting. Right after quantitation developer addition, fluorescence was measured using a multimodal microplate reader Synergy ™ HTX (Biotek, Winooski, Vermont, US) at 365/450 nm (ex/em). A standard curve was prepared with the Measure-iT ™ nitrite quantitation reagent concentrate included in the kit (0, 2.75, 5.5, 11, 22, 33, 44, 55 μM). Assays were conducted in duplicate in a final volume of 110 μl by adding 100 μl of quantification agent, fluorescence was measured, and plotted versus picomoles of NO -. The equation was used to determine the NOconcentration for each sample.

Quantification of EREG, Gal9, and C3 protein levels
The EREG, Gal9, and C3 expression was assessed in the supernatants of MAP-infected MDMs collected at 7 days p. i. using quantitative sandwich ELISAs according to the manufacturer's instructions (MyBioSource, San Diego, US). The sensitivity of the EREG, Gal9, and C3 ELISA kits is 5 pg/ml (detection range, 31.2-1.000 pg/ml), 0.1 ng/ml (detection range, 0.625-20 ng/ml), and 2 μg/ ml (detection range, 15.6 mg/ml-500 mg/ml), respectively. Briefly, standards and samples (50 μl) were added in duplicate into a Microelisa Stripplate provided with each kit. One hundred microliters of horseradish peroxidase-conjugated antibody were added to each well. After incubation for 60 min at 37°C in the dark, the plate was washed four times with 350 μl of wash solution and incubated with 50 μl of 3, 3′, 5, 5′-Tetramethylbenzidine for 15 min at 37°C in the dark. After adding 50 μl of stop solution into each well, the OD values were measured in an ELISA reader at 450 nm (Thermo Scientific Multiskan, US). We average the duplicate readings for each standard and sample and subtract the average OD of the blank. A standard curve was generated by plotting the mean OD values of each standard on the vertical axis and the corresponding concentration on the horizontal axis. The levels of EREG, Gal9, and C3 in each sample were interpolated from the standard curve. The correlations between the intracellular MAP load within MDMs at 2 h and 7 d p. i., and the quantification of EREG, Gal9, and C3 protein levels were analyzed using the Spearman's rank correlation coefficient (r) implemented in R 4.1.2. considering a coefficient with a P-value less than or equal to 0.05 as significant.

Genotyping
Peripheral blood samples were collected in EDTA tubes and DNA was extracted using the QIAmp DNA Blood Mini Kit according to the manufacturer's instructions (Qiagen, Hilden, Germany). Purified genomic DNA was genotyped with the EuroG MD Bead Chip at the molecular genetic laboratory service of the Spanish Federation of Holstein Cattle (CONAFE). The InfiniumTM iScan software (Illumina, San Diego, CA) was used for allele assignation. All the SNPs passed a call rate > 0.80. After filtering out SNPs with minimum allele frequency (MAF)< 0.01, the final marker set included 44,779 SNPs.

Variance components and heritability (h 2 ) estimates
The intracellular MAP load at 2 h and 7 d p. i. and the expression of NO -, EREG, Gal9, and C3 were the quantitative phenotypes analyzed. The variance components and h 2 explained by all the SNPs were calculated using the genome-wide complex trait analysis (GCTA) software 1.93.2, according to the following formula where s 2 G Is the variance explained by all the SNPs and s 2 e is the residual variance (28).

Genome-wide association study
Genotypes and the quantitative phenotypes were analyzed using the mixed linear model association analysis of the GCTA 1.93.2 software. Briefly, the model is y = a + bx + g + e, where y is the phenotype, a is the mean term, b is the additive effect (fixed effect) of the candidate SNP to be tested for association, x is the SNP genotype indicator variable coded as 0, 1 or 2, g is the polygenic effect (random effect) assumed to be distributed as N~(0 s 2 G ),and e is the residual assumed to be distributed as N~(0 s 2 e ). Age was included as a covariate in the analysis. To account for multiple testing, a 5% genome-wide false discovery rate (FDR) was used (29). The inflation factor (l) and quantile-quantile plots were used to compare the observed distributions of -log (P-values) to the expected distribution under the no-association model. l values close to 1 suggest appropriate adjustment for potential substructure and l > 1.2 suggest population stratification. The regression coefficients (b-values) were calculated using the GCTA 1.93.2 software.

Identification of SNPs and candidate genes
The EuroG MD Bead chip is based on the University of Maryland UMD 3.1 Bos taurus reference genome; hence, the coordinates of the GWAS-identified SNPs were converted to the ARS-UCD1.2 genome by using liftOver software (https://genome.ucsc.edu). The location of the identified SNPs (e.g. upstream or downstream of a transcript, in the coding sequence, in non-coding RNA, in regulatory regions) in the ARS-UCD1.2 genome was determined using the Ensembl Variant Effect predictor (VEP). None of the identified SNPs were within 500,000 base pairs of each other or on linkage disequilibrium. The candidate genes located within 50,000 base pairs to each side of the SNPs were identified using Ensembl (https://www.ensembl.org). The function of all the identified genes was searched in GeneCards (http:// www.genecards.org) by searching their gene symbol. Since there is evidence that some allelic variants may contribute to resistance to multiple pathogens, the identified SNPs and candidate genes were compared with QTLs and candidate genes previously associated with other bovine diseases, longevity, and reproductive and health traits (http://www.animalgenome.org). In addition, the identified candidate genes were also compared with human candidate genes previously identified for CD, IBD, and colorectal cancer (http:// www.ebi.ac.uk/gwas).

Genomic estimated breeding values for MAP survival indexes within MDMs and validation of the genomic predictions
gEBVs for MAP survival indexes within MDMs were calculated for each animal in the reference population using the gBLUP model of GCTA 1.93.2 (30). Subsequently, gEBVs for MAP survival indexes within MDMs were predicted in a validation population. PBMCs were isolated from animals with high (N = 8) and low (N = 8) gEBVs, differentiated to MDMs, infected ex vivo with MAP, and at 2 h and 7 d p. i. the intracellular load was quantified in the BACTEC MGIT system as described in point 2.3. An unpaired student t-test with the Welch-Satterthwaite correction was used to compare the experimental MAP survival indexes calculated for the animals with high and low gEBVs (GraphPad Prism 8, San Diego, CA, US). Differences were considered significant when P-values were ≤ 0.05. Correlations between the gEBVs and experimental MAP survival indexes within infected MDMs were calculated with the Pearson correlation implemented in R 4.1.2.

Functional phenotyping of bovine MDMs in response to MAP infection
As indicators of bovine MDMs performance, intracellular MAP load within MDMs at 2 h and 7 d p. i. and NO − , EREG, Gal9, and C3 secretion were measured in MAP-infected MDMs isolated from 75 cows ( Figure 2). The NOsecretion was measured in the supernatants of MAP-infected MDMs at 2 h p.i. while EREG, Gal9, and C3 were measured at 7 d p. i. To estimate the technical reproducibility of each phenotype, the technical coefficients of variation (CVs) were calculated for each phenotype. The technical CV calculated using the duplicate measures within one experiment for C3, EREG, Gal9, NO -, and MAP load with macrophages at 2 h and 7 d p. i. were 6.47% (0% -17.54%), 4.76% (0% -10.92%), 4.08% (0% -9.48%), 14.07% (0.89% -76.69%), 2.00% (0% -7.12%), and 7.11% (0% -33.89%), respectively. To estimate the variation within the population, the CV were calculated for the measures of each phenotype within the population. The CV of Gal9, C3, EREG, NO -, and MAP load within MDMs at 2 h and 7 d p. i. were 7, 9, 14, 41, 12, and 31%, respectively. The biological reproducibility of the MAP load within macrophages at 2 h and 7 d p. i. was previously assessed by calculating the CVs of biological triplicates within three different experiments; 9.40 and 10.72%, respectively (24). The biological CV considers the precision of the quantification of MAP load within macrophages to correctly characterize the cow or the fidelity of the reading of the phenotype. The correlations between the six phenotypes were analyzed using Spearman's rank correlation ( Table 1). The strongest positive correlation (r = 0.80) was found between the intracellular MAP load at 2 h and 7 d p. i., while the highest negative correlation (r = -0.64) was observed between the secretion of NOmeasured at 2 h p. i. and the intracellular MAP load at 7 d p. i. The NOmeasured at 2 h p. i. also correlated negatively with the intracellular MAP load at 2 h p. i. (r = -0.29). The secretion of C3 was positively correlated with the MAP intracellular load at 2 h and 7 d p.i.; r = 0.34 and r= 0.48, respectively. On the other hand, negative correlations were found between NOand Gal9, and NOand EREG expression; r = -0.24 and r = -0.35, respectively. Positive correlations were found between Gal9 land EREG, and Gal9 and C3 expression; r = 0.28 and r = 0.44, respectively.

Variance components and heritability (h 2 ) estimates
The intracellular MAP load at 2 h and 7 d p. i. and the extracellular NO -, EREG, Gal9, and C3 levels were the quantitative phenotypes analyzed. The variance components, additive genetic variance (s 2 G ) and residual variance (s 2 e ), and h 2 explained by all the SNPs were calculated using the genome-wide complex trait analysis (GCTA) software 1.93.2. Table 2 summarizes the variance components and h 2 estimates for each phenotype. The highest h 2 estimates were obtained for the MAP load within MDMs at 2 h (h 2 = 0.87) and 7 days p. i. (h 2 = 0.83). The estimated h 2 for C3 (h 2 = 0.78), EREG (h 2 = 0.74), NO -(h 2 = 0.44) and Gal9 (h 2 = 0.13) levels were also calculated. While the h 2 estimates for intracellular MAP load, C3 (h 2 = 0.78), and EREG (h 2 = 0.74) were high; NO -(h 2 = 0.44) and Gal9 (h 2 = 0.13) estimates were moderate to low.

GWAS
To identify associations between genetic variants and phenotypes, a GWAS was performed for each phenotype. Only the bacterial load within MDMs at 2 h and 7 d p. i. showed SNPs with a statistically significant association. Six SNPs surpassed the P FDR correction for the MAP load at 2 h p. i. and two for the MAP load at 7 d p. i. Overlap was seen in two SNPs associated with MAP load within MDMs at 2 h and 7 d p. i. These two SNPs were located on BTA17 and BTA18. Manhattan plots showing -log 10 (P-values) of association between each SNP and MAP load within MDMs at 2 h and 7 d p. i. are presented in Figures 3A, B,    The plots showed a distribution close to the expected line; l median = 1.00 and l median = 1.01 for MAP load at 2 h and 7 d p. i., respectively. This indicates that significant values were not overestimated due to population stratification or cryptic relatedness.

SNPs associated with MAP load within MDMs and candidate genes identification
The SNPs surpassing the significance threshold (P FDR ≤ 0.05), Pvalues, and candidate genes are presented in Table 3. The six identified SNPs were located on BTA2, BTA17, BTA18, and BTA21. Three of the SNPs were in intergenic regions, two in introns, and one is a downstream variant. The most genome-wide significantly associated SNP (P = 3.28 × 10 -7 ) was located on BTA2. The regression coefficients (b-values) of the six identified SNPs were all negative suggesting that all the SNPs were associated with a low bacterial load within MDMs. None of the six SNPs identified by GWAS was in strong linkage disequilibrium. A total of five candidate genes and one miRNA were identified within 50,000 bp upstream and downstream of the identified SNPs including the Cysteine and Serine Rich Nuclear Protein 3 (CSRNP3), U6 spliceosomal RNA (U6), Oxysterol Binding Protein Like 6 (OSBPL6), Coiled-Coil Domain Containing 92 (CCDC92), ENSBTAG00000042858, and bta-mir-2365. Two of the six identified SNPs were in intergenic regions and showed no evidence of a candidate gene within 100,000 bp distance.

Comparison of the identified SNPs and candidate genes to previously reported QTLs and candidate genes
The SNPs and candidate genes that were identified in the current study had not been associated with bovine PTB,

Validation of the genomic predictions
The gBLUP model (30) was used to estimate gEBVs of each individual in the study population based on the effect of the six SNPs with evidence of association with the survival indexes of MAP within MDMs and their genomic relationships. Subsequently, the six SNPs with evidence of association with the phenotype were used to predict the gEBVs of the validation population. For the validation of the genomic predictions, 16 animals with high (N = 8) and low (N = 8) gEBVs were selected. The distribution of the gEBVs is presented in Figure 4A. MDMs of the selected animals were isolated and infected ex vivo with MAP as previously described in Materials and Methods. Intracellular bacterial load was quantified at 2 h and 7 d p. i. and the experimental MAP survival indexes were calculated for the animals with low (mean = 75.7) and high gEBVs (mean = 84.1). There were significant differences in the killing capacity of macrophages from cows with low and high gEBVs (P-value = 0.0006) and the data demonstrated a tighter grouping of the macrophages from cows with low gEBVs having superior killing capacity ( Figure 4B). Pearson correlation coefficient between the gEBVs and experimental MAP survival indexes was positive (0.78) (P-value = 0.0003). These results demonstrated a significant effect of the six identified SNPs on MAP survival within macrophages and confirmed the existence of a strong genetic effect on macrophages' performance in response to MAP infection.

Discussion
Macrophages are important antigen presenter cells that participate in the destruction of pathogens by producing antimicrobial components after phagocytosis (39)(40)(41)(42)(43)(44). Some animals might be able to induce a strong early innate immune response and limit MAP load within macrophages better than others. A considerable proportion of this inter-individual variability might be genetic (16). Defining the relevant phenotype is the main challenge in identifying resistant animals and the genetic  profile of resistance against MAP infection. Since ex vivo assays measure the capacity of a cell population to respond to a challenge, functional testing represents the best way to evaluate resistance and vaccine protection (45). MGIA assays using bovine MDMs were previously used to assess macrophages´s capacity to clear MAP strains isolated from different hosts (23, 24), to test novel therapeutic approaches for the treatment of MAP-infected animals (46), and to test vaccines' efficacy and protection (47).
Our study is the first in using a quantitative MGIA as an unbiased measure of the bovine MDMs´ability to control MAP survival in vitro and as a correlate of resistance to MAP infection. Since the MGIA models MAP infection, it reveals immune parameters induced by in vivo infection and contributes to understanding the control of intracellular MAP load. In addition, evidence exists that the ability of the macrophages to produce an antimicrobial profile early after infection can cause differences in disease susceptibility or resistance in connection to their cellular responses (48)(49)(50). Previously, next-generation RNA (RNA-Seq) sequencing was used to study the transcriptomic response to MAP infection of the macrophages from cows that were naturally infected and identified as positive for JD (JD (+); n = 22) or negative for JD (healthy/resistant, JD (−); n = 28) (51). In addition to identifying genetic variants from RNA-seq data, SNP variants were also identified using the Bovine SNP50 DNA chip. Significant variants from JD (+/-) macrophages were identified by a genome-wide association study in a case-control approach and revealed two novel eQTLs on BTA4 and 11 (P< 5 × 10-7 ). This study was conducted using MDMs from cattle that were naturally infected with MAP. In our study, however, we used MDMs obtained from cattle that were negative by PCR for the detection of MAP DNA in fecal samples before blood collection and in subsequent annual samplings. Our study is the first in using a quantitative MGIA as an unbiased measure of the bovine MDMs´ability to control MAP load and as a correlate of resistance to MAP infection. In our previous study (22), the integration of RNA-Seq data and genotypes from a cohort of cows naturally exposed to MAP allowed the identification of four cis-eQTLs regulating the expression of genes with important roles in the induction of the innate immune response early after infection; the C3, EREG, Gal9, and NOS1 genes. As shown in Supplementary Figure 1, none of the animals with the highest levels of expression of these genes showed diffuse lesions, which suggested that the measure of the levels of expression of these genes might be used to evaluate macrophages performance and for the identification of cattle able to limit MAP infection. The goal of the current study was to identify genetic loci associated with resistance to MAP infection by combining genotypes and phenotypes including the intracellular MAP load within infected MDMs and the extracellular expression of NO -, Gal9, C3, and EREG. NOS1 belongs to the family of NOsynthases which synthesize NOfrom L-arginine. NOis one of the major macrophage-killing mechanisms (52). It was previously demonstrated that MAP interferes with NOsynthase expression early after infection (53). EREG belongs to the Epidermal Growth factor (EGF) family together with the Heparin Binding EGF-like Growth factor (HBEGF); acting both as mitogenic stimulators via binding to EGF receptors (EGFRs). EREG overexpression leads to downstream signaling events including Mitogen-Activated Protein Kinase 8 (MAPK8) phosphorylation and activation which results in apoptosis induction. Using RNA-Seq transcriptomic profiling, we previously detected down-regulation of EREG in peripheral blood and ileocecal valve from cows with diffuse lesions versus controls which may favor MAP survival (54). Gal9 is an S-type lectin with antimicrobial activity in infected macrophages by causing macrophage activation, Interleukin 1ß (IL1ß) secretion, and pathogen clearance (55,56). Higher C3 levels may allow increased opsonophagocytosis and effective bacterial clearance of several microorganisms including mycobacteria (47,57). In our study, negative correlations were observed between the extracellular NOand the MAP load within infected MDMs at 2 h and 7 d p. i. The correlations between EREG and Gal9 and intracellular MAP load within MDMs were not statistically significant. On the other hand, positive correlations were found between MAP load within macrophages and C3 expression at 2 h and 7 d p. i. One of the critical factors that influence the power of detection of SNP (s)-associated with a trait is how much of the phenotypic variance is explained by the associated SNP (s). The highest the h 2 estimate, the chance of detecting SNPs increases. The h 2 estimates obtained in our study for the MAP load (log CFUs) within infected MDMs at 2 h and 7 d p. i. were 0.87 and 0.83, respectively. These findings suggest a strong influence of host genetics on MAP load within infected macrophages. To our knowledge, the h 2 estimate for the macrophages´performance assessed by measuring MAP load within macrophages had not been estimated before and is higher than h 2 estimates calculated for other immunocompetence traits. For instance, the analysis of MDMs phenotypic variation revealed an h 2 = 0.77 for the NOproduction by MDMs in response to E.coli (18). After performing the GWAS, we did not identify SNPs significantly associated with the NO -, Gal9, EREG, and C3 expression. However, six SNPs located on BTA2, BTA17, BTA18, and BTA21 were significantly associated with a lower MAP load within MDMs after 2 h p. i. This discordance can be explained because the MGIA is a functional assay that assesses the combined effects of a range of immune and metabolic mechanisms that influence mycobacteria growth in vitro, which is one of its advantages over measuring single predefined immune markers. It is important to highlight that for the assessment of viable intracellular MAP, we used the BACTEC MGIT 960 system, a fully automated solution for mycobacterial liquid culture that provides accurate measures of MAP load.
We identified six SNPs associated with low MAP load within MDMs in dairy cattle. Interestingly, the two identified SNPs located on BTA18 overlapped with QTLs associated with the somatic cell score (QTL:3556), length of productive life (QTL:3555), and reproductive traits such as stillbirth (QTL:11362), birth Index (QTL:30537), and dystocia (QTL:2707) highlighting the important role of this region (23. 5-31.4 Mbp) on chromosome 18 in health, reproduction, and longevity (31)(32)(33)(34)58). This suggests that an animal that can limit MAP infection has also improved udder health, fewer reproductive problems, and, consequently, a longer productive life. In addition, the SNP located on BTA21 was located within QTLs associated with IgG levels (QTL:66212), somatic cell score (QTL:5446 and QTL:10190), and calving ease (QTL:11126) (35-38). This finding suggests that there might be genetic loci that impact both the primary innate immune response and the secondary production of IgG which is activated later if innate mechanisms are not able to eliminate the pathogen. Finally, the identified SNP on BTA17 was located within a QTL associated with calving ease (36). The identified SNPs on BTA17, BTA18, and BTA21 overlapped with QTLs for longevity, udder health, and reproductive traits. The identified markers could be useful for marker-assisted selection in dairy cattle breeding programs for PTB resistance in combination with other strategies like management and biosecurity.
Five candidate genes and one miRNA were identified 50,000 bp upstream and downstream of some of the identified SNPs; Cysteine and Serine Rich Nuclear Protein 3 (CSRNP3), the Coiled-Coil Domain Containing 92 (CCDC92), Oxysterol Binding Protein like 6 ( O S B P L 6 ) , U 6 , a n o n -c h a r a c t e r i z e d p r o t e i n (ENSBTAG0000004285), and the bta-mir-2365. Specific functions of the four identified candidate genes are summarized in Supplementary Table 2. CSRNP3, also named TGF-Beta Induced Apoptosis Protein 2, is a transcriptional activator involved in positive apoptosis regulation. Apoptotic clearance of pathogens including MAP plays a critical function in the resolution of the infection by the removal of the infected macrophages along with the ingested organisms. CCDC92 is an interferon-g stimulated protein that plays an important role in the regulation of innate immunity and attenuation of microbial infections (59). The OSBPL6 is an intracellular lipid receptor that contributes to the maintenance of cholesterol homeostasis by regulating cellular cholesterol trafficking and efflux. Previously, genetic variants in candidate genes involved in cholesterol metabolism (bta04979) such as Angiopoietin-like 4 and 8 (ANGPTL4 and ANGPTL8), Low Lipropotein Receptor (LDRL), and Neutral cholesterol Ester Hydrolase 1 (NCEH1) appeared enriched in animals with PTBassociated diffuse lesions (22). Consequently, reduced levels of plasma cholesterol were detected in cattle with diffuse lesions. Moreover, we observed that the cholesterol route (bta04977) was dysregulated in the ileocecal valve of cows with diffuse lesions with four upregulated Apoliprotein genes (APOA1, APOC3, APOA4, APOB) matching this route (54). Recent data suggest that MAPinfected macrophages accumulate intracellular cholesterol droplets and depict a foam phenotype during infection providing an enriched environment for MAP survival (60,61). In addition to the observations that increased lipid biogenesis and transport to the place of the infection may favor MAP growth by the modulation of host immune response leading to an anti-inflammatory milieu, MAP may potentially mobilize lipid body content as a source of nutrients to enhance their survival and replication in host cells (62). In an opposite manner, our study suggests that animals with genetic variants in candidate genes favoring cholesterol homeostasis such as the OSBPL6 might be able to limit MAP infection by reducing the lipid body content within infected MDMs. The U6 spliceosomal RNA and the bta-miR-2365 were also identified as candidate genes in our study. The U6 is a small non-coding RNA involved in the post-transcriptional regulation of gene expression by affecting both the stability and translation of mRNAs. Using miRmap, we predicted target genes for the bta-miR-2365. When we used ClusterProfiler with the predicted target genes, an enrichment in genes associated with the endocytosis (bta04144) and the Hippo signaling pathway (bta04390) was found. The Hippo signaling pathway controls organ size by coordinately regulating cell growth, proliferation, and apoptosis (63). Emerging evidence has suggested that the innate immune response was extensively regulated by multiple core components of the Hippo cascade while the activation of innate immune signaling conversely led to substantial modulations of the Hippo pathway (64). In most cases, the Hippo signaling potentiates the innate immune response. In addition, the Hippo pathway-mediated processes are interconnected with those of other key signaling cascades, such as those mediated by TGF-beta and Wnt growth factors. Therefore, the bta-miR-2365 along with the other candidate genes identified in our study could modulate MAP's entry and survival within macrophages. The candidate genes identified in our study were also compared to candidate genes previously identified in CD, IBD, and colorectal cancer to determine if there was any overlap. We found that CSRNP3 was previously found to be associated with colorectal cancer (65). In clear cell renal cell carcinoma (ccRCC), the CSRNP3 may serve as a prognostic biomarker to predict the overall survival of patients (66). CSRNP3 might impact the immune environment of ccRCC through immunocyte infiltration of natural killer cells and plasmacytoid dendritic cells.
The sample size used in the current study could be considered small in comparison to traditional GWAS analysis using a larger population. However, the assessment of the macrophagesṕ erformance by measuring MAP load within infected MDMs is a functional and controlled trait that can be measured in the absence of environmental and physiological effects. Other phenotypes, such as the NOresponse of bovine MDMs in response to E.coli or S. aureus, used approximately 60 samples to identify associated cis-SNPs (18). Further studies are needed to assess the capacity of macrophages to kill MAP as a phenotype for resistance to MAP infection in other cattle populations and ruminant breeds.

Conclusions
We report a new phenotype, the control of MAP early survival within infected macrophages, for evaluating resistance to MAP infection in dairy cows. This resistance is genetically determined and is mediated at the level of macrophages´function. Considering how difficult it is to assess resistance to MAP infection, our approach represents a major advance in the field. The SNPs and candidate genes identified here revealed the association between MAP load within macrophages and genetics, which helps fill the knowledge gap between genetics, innate immune responses, and resistance to MAP infection in cattle. Our results define a heritable and distinct immunogenetic profile in MAP-infected macrophages associated with cholesterol homeostasis regulation and with the induction of apoptosis and IFNg production. A strong effect of host genetics in limiting bacteria load within MDMs very early after infection was observed (h 2 estimate = 0.8) which opens the possibility of ranking Holstein cows based on predicted MAP survival variation within infected MDMs. The identified SNPs might be used to develop genetic evaluations for immunocompetence in the Spanish breeding program which would allow producers to select cattle more resistant to MAP infection and likely to other intracellular pathogens as well; ultimately reducing the prevalence of diseases and the dependence on antimicrobials, preventing economic losses, increasing the length of cattle productive life, and improving food safety.

Data availability statement
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

Ethics statement
The animal study was reviewed and approved by the Ethics Committee of Animal Experimentation of NEIKER.

Author contributions
MC measured intracellular MAP load within MDMs, extracellular EREG, C3, Gal9, and NOlevels, and performed the GWAS. GB-B performed the gBLUP model and contributed to the statistical analysis of the data. MC and GB-B performed the validation of the genomic predictions. MA-H was the principal investigator of the project and participated in project management, experimental design, and data analysis. All authors contributed to the article and approved the submitted version.