Response to Oxidative Burst-Induced Hypoxia Is Associated With Macrophage Inflammatory Profiles as Revealed by Cellular Genome-Wide Association

Background In mammalian species, hypoxia is a prominent feature of inflammation. The role of hypoxia in regulating macrophage responses via alteration in metabolic pathways is well established. Recently, oxidative burst-induced hypoxia has been shown in murine macrophages after phagocytosis. Despite the available detailed information on the regulation of macrophage function at transcriptomic and epigenomic levels, the association of genetic polymorphism and macrophage function has been less explored. Previously, we have shown that host genetics controls approximately 80% of the variation in an oxidative burst as measured by nitric oxide (NO-). Further studies revealed two clusters of transcription factors (hypoxia-related and inflammatory-related) are under the genetic control that shapes macrophages’ pro-inflammatory characteristics. Material and Methods In the current study, the association between 43,066 autosomal Single Nucleic Polymorphism (SNPs) and the ability of MDMs in production of NO- in response to E. coli was evaluated in 58 Holstein cows. The positional candidate genes near significant SNPs were selected to perform functional analysis. In addition, the interaction between the positional candidate genes and differentially expressed genes from our previous study was investigated. Results Sixty SNPs on 22 chromosomes of the bovine genome were found to be significantly associated with NO- production of macrophages. The functional genomic analysis showed a significant interaction between positional candidate genes and mitochondria-related differentially expressed genes from the previous study. Further examination showed 7 SNPs located in the vicinity of genes with roles in response to hypoxia, shaping approximately 73% of the observed individual variation in NO- production by MDM. Regarding the normoxic condition of macrophage culture in this study, it was hypothesized that oxidative burst is responsible for causing hypoxia at the cellular level. Conclusion The results suggest that the genetic polymorphism via regulation of response to hypoxia is a candidate step that perhaps shapes macrophage functional characteristics in the pathway of phagocytosis leading to oxidative burst, hypoxia, cellular response to hypoxia and finally the pro-inflammatory responses. Since all cells in one individual carry the same alleles, the effect of genetic predisposition of sensitivity to hypoxia will likely be notable on the clinical outcome to a broad range of host-pathogen interactions.

Background: In mammalian species, hypoxia is a prominent feature of inflammation. The role of hypoxia in regulating macrophage responses via alteration in metabolic pathways is well established. Recently, oxidative burst-induced hypoxia has been shown in murine macrophages after phagocytosis. Despite the available detailed information on the regulation of macrophage function at transcriptomic and epigenomic levels, the association of genetic polymorphism and macrophage function has been less explored. Previously, we have shown that host genetics controls approximately 80% of the variation in an oxidative burst as measured by nitric oxide (NO -). Further studies revealed two clusters of transcription factors (hypoxia-related and inflammatory-related) are under the genetic control that shapes macrophages' pro-inflammatory characteristics.
Material and Methods: In the current study, the association between 43,066 autosomal Single Nucleic Polymorphism (SNPs) and the ability of MDMs in production of NOin response to E. coli was evaluated in 58 Holstein cows. The positional candidate genes near significant SNPs were selected to perform functional analysis. In addition, the interaction between the positional candidate genes and differentially expressed genes from our previous study was investigated.
Results: Sixty SNPs on 22 chromosomes of the bovine genome were found to be significantly associated with NOproduction of macrophages. The functional genomic analysis showed a significant interaction between positional candidate genes and mitochondria-related differentially expressed genes from the previous study. Further examination showed 7 SNPs located in the vicinity of genes with roles in response to hypoxia, shaping approximately 73% of the observed individual variation in NOproduction by MDM. Regarding the normoxic condition of macrophage culture in this study, it was hypothesized that oxidative burst is responsible for causing hypoxia at the cellular level.

INTRODUCTION
Macrophages are crucial players in both the innate and adaptive immune responses. These cells are equipped with mechanisms to destroy pathogens, regulate inflammation and present antigens to the cells of the adaptive immune system (1,2). In addition, macrophages are known as effector cells in antibody-mediated immune responses (3). During steady-state, macrophages in some organs (i.e. liver) are generated from embryonic precursor while in some other organs, they are constantly replenished by monocytes (e.g. intestine) (4)(5)(6). The difference in the origin of macrophages during infections does not vary between most of the tissues as macrophages are mainly monocyte-derived (5). The genetic and epigenetic mechanisms of the monocyte transformation to macrophages, and biological pathways that are employed by macrophages in various actions, have been reported in detail by various research groups (7)(8)(9)(10)(11)(12)(13). However, the genetic polymorphisms that are associated with the magnitude of the macrophage functional characteristics have remained less investigated (14,15). For instance, the paradigm of macrophage polarization between the two extremes of M1 and M2 in respect to stimuli combination, epigenetic mechanism and the transcriptomic profile have been reported in detail (16), but the effect of host genetic polymorphism on macrophage polarization, epigenetic mechanism and transcriptomic profile have been less investigated. Our group has recently used the magnitude of macrophage nitric oxide (NO -) production in response to Escherichia coli (E. coli) as an index to reveal the effect of host genetic polymorphism on functional characteristics of monocyte-derived macrophages (MDM) (17). This recent study of bovine MDM, as well as another report from a study using a mouse model, have shown that approximately 80% of the individual variation is described by the host genetics (17,18). The follow up study on the transcriptome of MDM that were classified based on NOlevel, showed a distinct inflammatory profile between high and low NOresponder groups. The bioinformatic analysis then revealed that activation of the Hypoxia-induced factor-1a (HIF1A) pathway was positively associated with NOproduction ability of MDM (19).
Induction of the HIF1A pathway in macrophages after exposure to pathogens or stimulation with a pathogen-derived component such as lipopolysaccharide has also been reported by other groups (20,21). From the cellular metabolic point of view, mitochondria are known as the regulators of cellular metabolism and the main intracellular organelle to regulate oxygen and energy metabolites in cellular function (22,23). The role of mitochondria in response to hypoxia is very well documented (24). In addition, mitochondria have recently gained attention as a master regulator of intracellular signalling with an important role in inflammation (25,26). Specifically, the role of mitochondria in the regulation of macrophage function and its association with intracellular NAD(P)H content have been reported (27)(28)(29). However, the genetic mechanism that might regulate macrophages function via mitochondria has remained less clear.
In the current study, it was hypothesized that the magnitude of a macrophage inflammatory response is associated with polymorphism in genes that interact with mitochondria. In order to test this hypothesis, NOresponse of bovine MDM against E. coli was measured as the phenotypic index of the inflammatory response. A genome-wide association analysis was employed to discover the genetic polymorphisms that are associated with the NOindex. As the last step, the interactions between the genes nearby significant SNPs and mitochondria-related genes were investigated. The results of this study revealed a significant network that was formed between the positional candidate genes and the differentially expressed genes in response to E. coli in MDMs that are annotated with "mitochondrion" term.

Animals
Fifty-eight lactating Holstein cows were selected for this study from the University of Guelph Dairy Innovation Centre herd. These animals were not diagnosed nor treated for any diseases in the lactation period when samples were collected. Blood samples were collected as part of a previous study in groups of 12 individuals per sampling visit (17). All the procedure and handling of the animals were approved by the animal care committee of the University of Guelph (AC4449).
Aldrich, St. Louis, MO) was loaded into the Sepmate tubes (STEMCELL Technologies, Vancouver, BC) and whole blood samples were overlaid on the top of Histopaque-1077. After centrifugation for 10 minutes at 1200 ×g, the layer of cells above the Histopaque-1077 was collected and washed 3 times with Phosphate Buffered Saline (PBS) to obtain the purified BMCs. Purified BMCs were cultured at the concentration of 1×10 6 cells per square centimetre of the culture flask (Nunc, Thermo Fisher Scientific Inc., Mississauga, ON) for 2 hours in Monocyte Attachment Medium (PromoCell, Heidelberg, Germany) at 37°C. Non-adherent cells were removed by washing using PBS and the medium was replaced with AIM V ® Medium (Thermo Fisher Scientific Inc., Mississauga, ON) containing 5 ng/ml recombinant bovine Granulocyte-Macrophage Colony-Stimulating Factor (GM-CSF, Kingfisher Biotech, St. Paul, MN) in the presence of 5% CO 2 . After 6 days of incubation, adherent cells were detached from the flask using TrypLE ™ Select Enzyme (Thermo Fisher Scientific Inc., Mississauga, ON). Phenotypic characteristics of macrophages (CD14 + , CD205and strong auto-fluorescence) were analyzed using flow cytometry to determine the proportion of macrophages among the harvested cells (30)(31)(32). Harvested cells were labeled with RPE conjugated mouse anti-bovine CD205 (Clone: IL-A114, Bio-Rad, Mississauga, Ontario) and ALEXA FLUOR ® 647 conjugated mouse anti-human CD14 (Clone: TÜK4, Bio-Rad, Mississauga, Ontario), separately. The labelled samples were analyzed using the BD Accuri C6 flow cytometer (Becton Dickinson, Franklin Lakes, NJ) and the data from the flow cytometer were analyzed by FlowJo V.10 (FlowJo LLC, Ashland, OR). The fluorescence emission of unlabeled harvested cells was compared with the emission of unlabeled fresh BMCs in 533/30 filter excited by the blue laser to test the auto-fluorescence as described previously (17).
The MDMs from each sample were seeded in 2 wells of a 48well plate at the density of 2×10 5 cells per well in AIM V ® medium containing 5 ng/ml GM-CSF and incubated overnight. One well was assigned as the control and the other well was assigned as the treatment group and was exposed to inactivated E. coli (MOI: 5) for 48 hours. Supernatant from each well was collected and the concentration of nitric oxide (NO -) was measured with the Measure-iT ™ High-Sensitivity Nitrite Assay Kit (Thermo Fisher Scientific Inc., Mississauga, ON). The concentration of NOin the challenge group was corrected by the control group. These data collected in the previous study were used in the analyses described below (17).

Genotyping and Quality Control
DNA was extracted from hair follicles collected from each individual and genotyping was performed with the Illumina Bovine SNP50 BeadChip by Zoetis Canada (Kirkland, Quebec, Canada). The initial dataset contained 45,187 SNP markers that are used in routine official genomic evaluation in Canada by the Lactanet (Guelph, ON, Canada). Details of quality control were explained in Wiggans et al. (33). The sporadic missing genotypes were imputed using 50,000 reference Holsteins from the Lactanet database by FImpute software (34). In the present study SNPs located on the X chromosome were not included and due to relatively small sample size, 1,716 SNPs with minor allele frequency < 1% were excluded resulting in 43,066 SNPs for further association analysis.

Statistical Analysis
The UNIVARIATE procedure of SAS (V 9.4) was used to test the distribution of the dependent variable (NOresponse to E. coli) for normal, Log-normal, Weibull, and Gamma distributions. The goodness-of-fit tests based on the empirical distribution function for Gamma distribution were >0.50. The p values for goodness-of-fit tests based on the empirical distribution function, Anderson-Darling, Cramer-von Mises, and Kolmogorov-Smirnov, for Gamma distribution, were >0.50, > 0.25 and >0.50, respectively. Therefore, the dependent variable was transformed based on the method described by Krishnamoorthy et al. (2008) to achieve normal distribution (35).
The GLM procedure of SAS (V. 9.4) was used to test the effect of age in month, days in milk (classes: 1 = 0-20 days, 2 = 21-105 days, 3 = 106-235 days, 4 = > 235 days) (36,37), days in pregnancy (classes: 1 = non pregnant, 2 = 1 -120 days, 3 = 121 -180 days, 4 = 181 -220 days) (37) and sampling group (classes: 1-5). The statistical model was as follow: Where; y ijkl is the vector of the phenotype [the cubic roots of nitrite concentration of culture supernatant 48 hours after treatment with E. coli (n = 58)]; m is the overall average of the response. b 1 is the linear coefficient of the fixed regression on age (a i ) (in months), g j is the fixed effect of j th class of sampling group, m k is the fixed effect of k th class of days in milk, p l is the fixed effect of l th class of days in pregnancy and e ijkl is the random residual effect.
Genome-wide association analysis was carried out by SNP & Variation Suite (ver. 8.8.1, Golden Helix Inc.) using single SNP regression analysis on numerically coded alleles and the cubic roots of NOconcentration as the continuous dependent variable. The statistical model was as follow: Where y is the cubic root of nitric oxide concentration in the culture supernatant 48 hours after treatment with E. coli, m is the overall average of the response; b is the gene substitution effect for the SNP, l is the vector of eigenvalues, e is the random residual effects; c is a vector of genotypes for the SNP coded as described, below.
The additive, dominance and recessive models were utilized, separately. For the additive genetic model, the SNP alleles were numerically coded as 0, 1, and 2 for "dd", "Dd", and "DD" genotypes, respectively. For the dominance genetic model, alleles were coded as 0 for "dd"; and 1 for "Dd" and "DD". For the recessive genetic model, alleles were coded as 0 for "dd" and "Dd"; and 1 for "DD", where d is the major allele and D is the minor allele (38,39). The principal component analysis was performed on the genotypic data and the top three eigenvalues were included in the analysis to account for the population stratification. The details of the estimators and assumptions are provided in the SNP & Variation Suite (ver. 8.8.1, Golden Helix Inc.) manual. To reduce the probability of false discovery in multiple comparisons, a permutation test was employed by randomizing the observations over the genotypes in 10,000 permutation samples. The permutation value was calculated based on the direct count of the number of samples that have resulted in equal or smaller p values compared to the unrandomized pvalue divided by the total number of samples. The permutation value of 0.0005 was set as the cut-off for the significant association. The linkage between the significant SNPs were evaluated using the Expectation-Maximization algorithm in SNP & Variation Suite (ver. 8.8.1, Golden Helix Inc.) (40). The SNPs with linkage disequilibrium R-Squared index of greater than 0.5 were considered linked and the second SNPs in each pair were removed for the quantitative analysis.
The significant SNPs were used in the K-fold genomic prediction algorithm in SNP & Variation Suite using Bayes C (p) method in 50,000 iterations in 6 folds of prediction (41).

Identifying Positional Candidate Genes
Candidate genes were selected based on two strategies; 1) Position-dependent; and 2) Comparative genomics (42). First, the coordinates of significant SNPs from UMD3.1 were converted to ARS-UCD1.2 by using BovineMine tools from Bovine Genome database. Genes (Ensemble Gene 97) that are located in the vicinity of the significant SNPs (up to 50K bp upstream and downstream of the significant SNPs) were considered as the positional candidate genes. Due to the fact that the annotation of the bovine genome is not complete, the comparative strategy was employed on the SNPs when the first strategy was not able to retrieve any genes. In the comparative strategy, 1Kb up and downstream of the significant SNPs were aligned against the human genome (GRCh38.p12, Ensemble Genes 97), and genes that were located in the identical region were selected as the positional candidate genes. The Gene Ontology (GO) terms for the positional candidate genes were extracted from Ensemble using 'biomaRt' (v. 2.41.7) package in R (v. 3.6) (43).

Network Analysis
Activation of HIF1A that was found in the previous study, as well as the known role of mitochondria in response to hypoxia led to investigating the interaction between the positional candidate genes and mitochondria, the genes that have been found to be differentially expressed (DE) between the high and low nitric oxide responder were selected from our previous study (GSE136722) (19). Then, GO terms for these genes were extracted from Ensemble using 'biomaRt' (v. 2.41.7) package in R (v. 3.6). Only, the DE gene that has been annotated with "mitochondrion" (GO:0005739) were selected for interaction analysis. The protein identifiers of these genes were combined with the list of all positional candidate genes from the current GWAS. The protein-protein interaction network in the final gene set (mitochondria-related DE genes and positional candidate genes from the GWAS) was analyzed using the Search Tool for Retrieval of Interacting Genes/Proteins (STRING, v. 11) (44).

RESULTS
The flow cytometry assay confirmed that more than 98% of the harvested cells after 6 days of incubation had characteristics of macrophages (CD14+, CD205-, and strong autofluorescence in (p-value = 0.25), days in milking (p-value = 0.47), or pregnancy (p-value = 0.73) on the production of NOin response to E. coli (Model p-value = 0.52, R-Squared = 0.21). Therefore, the cubic root of NOconcentrations without any adjustment for environmental effects were used in the association analysis. The association of 60 SNPs with the concentration of NOwas found to be statistically significant (permutation value < 0.0005). These SNPs were located on 22 autosomal chromosomes. The pattern of genetic effects for 22, 20 and 18 SNPs followed additive, dominant, and recessive models, respectively (Figures 2-4 and Supplementary File 1, Sheet Table 1). The linkage analysis showed strong LD among 9 pairs of SNPs ( Figure 5). The 51 informative SNPs were able to predict the NOconcentration with 85% correlation between the observed and the predicted NOconcentration with an R-Squared of 0.708 for the linear regression between the observed and the predicted NOconcentration ( Figure 6). This R-Squared was considered as the genomic heritability of the NOindex (45).
By employing the position-dependent strategy, 72 genes, including 14 novel genes, 8 non-coding genes and 2 micro-RNA were identified as positional candidate genes. By employing the comparative genomic approach, 8 coding genes were identified based on the similarity to the human genome. No genes were found in the vicinity of the 8 SNPs at the end of screening for the positional candidate genes (Supplementary File 1 -Sheet Table 1). The top 3 R-Squared for a single SNP were 0.413 (adjacent to THBS2), 0.332 (adjacent to PARPBP), and 0.306 (adjacent to ENSBTAG00000045919) (Figure 7). After removing all non-coding and novel genes and combining genes from all 3 steps of the screening, the final list of positional candidate genes consisted of 58 genes (Supplementary File 1 -Sheet Table 1). In this set of 58 genes, the "vesicle-mediated transport" GO term (GO:0016192) was in common between 8 genes (ARL1, STXBP6, PIK3C2A, SORL1, YIPF5, MCFD2, STX2, and VAMP4), calcium-related GO terms of "calcium ion binding" (GO:0005509), and "calcium:sodium antiporter activity" (GO:0005432), were in common between 6 genes (SNED1, THBS2, LC8A1, FGF14, CDH20, PDE1A, and LRP1B) and "response to hypoxia" (GO:0001666) were in common between 4 genes (APAF1, CD38, SLC18A1, and SMAD3).
From the previously data set of 839 DE genes between high and low NOresponders (Fold change > 1.5 and FDR p-value < 0.1), 44 genes were found to be annotated with the GO term of "mitochondrion" (GO:0005739) (Supplementary File 1 -Sheet Table 3). The interaction analysis using STRING (v. 11) on the combined list of genes (positional candidate genes and mitochondria-related DE genes = 102 genes) resulted in one network of 100 connected nodes ( Figure 8). The enrichment pvalue of this network was 6.87e-6 with the average local clustering coefficient of 0.345. Only two genes, SLC18B1 and CNBD1 (both were selected at the first step of screening) were not connected to the network. The functional enrichment analysis on the network showed that "response to oxygencontaining compound" (GO:1901700) and "response to FIGURE 2 | Distribution of -Log10 regression p values from additive genetic model analysis of SNPs in association with MDM nitric oxide production. In this analysis, the alphabetic codes of aa, ab, and bb for each SNP were transformed to 0, 1, and 2, respectively. The regression analysis was performed between normalized nitric oxide response of Monocyte-Derived Macrophages culture treated by Escherichia coli (n = 58) and numerically coded SNPs. To account for multiple comparison errors, permutation analysis (10,000 iterations) analysis for each SNP and the permutation p-value of 0.0005 were considered as the threshold for each SNPs. Significant SNPs are shown in red.  Table 4).

DISCUSSION
Identifying the genetic elements that shape the host response to pathogens has been gaining more attention as the treatment and prognosis of diseases are moving towards precision (veterinary) medicine. Individual differences in response to pathogens or treatments have been observed in outbred populations. However, identifying the genetic mechanisms that control this variation is extremely difficult. More than 8,000 genes in the human genome are known to contribute to immune responses (Ensemble, Genes 97) (46). In addition, the effective population size of the human population is very large (47). The combination of these two factors requires a large sample size to discover genetic control of immunity at the whole-organism level. Reductionist approaches are shown to be a successful alternative to overcome some of the barriers of studies on whole-organisms. Investigating the effect of genetic polymorphism on regulation of subsystems, instead of the whole organism, can overcome the limitation of the complexity of the trait. However, in the reductionist approach, FIGURE 4 | Distribution of -Log10 regression p values from recessive genetic model analysis of SNPs in association with MDM nitric oxide production. In this analysis, the alphabetic codes of aa, ab, and bb for each SNP were transformed to 0, 0, and 1, respectively. The regression analysis was performed between normalized nitric oxide response of Monocyte-Derived Macrophages culture treated by Escherichia coli (n = 58) and numerically coded SNPs. To account for multiple comparison errors, permutation analysis (10,000 iterations) analysis for each SNP and the permutation p-value of 0.0005 were considered as the threshold for each SNPs. Significant SNPs are shown in red. the effect of environment, interactions with other cells and the physiological condition of host are removed which may limit the application of the findings. The genetic control of the subsystems is much simpler than the whole organism (48). Hence, these studies require smaller sample size and the generated data is free from environmental effects (49). In 2013, the concept of cellular genome-wide association study was proposed by Dennis Ko et al. (50). Since then, this approach has been used to reveal the effect of genetic polymorphism on regulation of hostpathogen interaction (51)(52)(53). It has been shown that the genetic associations that are identified in this approach are relevant to resistance to disease based on clinical data (54,55). Genetic studies using a reductionist approach in livestock can further address the limitation of effective population size in comparison to studies on human. The heritability of NOproduction in bovine MDM in-vitro was estimated to be 0.78. This heritability is similar to the heritability of the NOproduction by mouse peritoneal macrophages (0.8), indicating the majority of variation is described by genetics in this in-vitro model (17,18). Further investigation of the cellular mechanisms shaping the NO-based phenotype showed that one network of transcription factors is controlled by the host genetics. This network was composed of two clusters: hypoxia-related genes including HIF1A and inflammatory-related clusters. The induction of HIF1A under the normoxic condition in macrophages and neutrophils has been reported previously (20,21). This activation is thought to be hypoxia-independent (20). However, the role of the oxidative burst in the activation of HIF1A was shown by abrogation of the HIF1A pathway when ROS generation was inhibited in macrophages (20,56). Furthermore, Lue et al. showed that oxidative burst can cause hypoxia at the cellular level and activate HIF1A (57).
In the current study, 60 SNPs were found to be significantly associated with NO-production. Although the sample size of this study seems to be small for a GWAS, sample size is a relative concept. In association studies, it depends on several factors, for instance: the level of genetic diversity in the population, the complexity of the phenotype and its heritability, environmental and physiological effects, and the population structure (58).  Escherichia coli based on the cow's genotype for significant SNPs nearby genes with role in response to hypoxia. The y-axis represents the nitric oxide concentration in µM, the animal genotypes are represented along the x-axis. "The lower and upper hinges correspond to the first and third quartiles. The upper whisker extends from the hinge to the largest value no further than 1.5 * IQR from the hinge (where IQR is the inter-quartile range, or distance between the first and third quartiles). The lower whisker extends from the hinge to the smallest value at most 1.5 * IQR of the hinge." (Summary of statistics for geom_boxplot function in 'ggplot2' R package).
The effective population size (N e ) of Holstein, based on linkage disequilibrium, has been estimated to be approximately 58 (59). The small N e indicates a high level of inbreeding, limited genetic diversity and long genetic linkage (59,60). The phenotype that was analyzed in the current GWAS is a relatively simple trait and highly heritable (h 2 : 0.77-0.99), measured in a fully controlled environment, and shown to be free from physiological effects (17). In addition, samples were collected from one sex and one breed to reduce the possibility of population stratification. Considering these factors, the sample population was sufficient to identify loci that are significantly associated with the trait. A similar design in a study on the human population has successfully identified loci that are associated with cytokine response by using 500 samples (51). It should be noted the effective population size of the aforementioned human population was approximately 150,000 which is 2,500 fold larger than effective population size of Canadian Holsteins (61). When the effective population size is smaller, much lower samples will be adequate to find QTLs (61). This concept in human studies is known as familial GWAS or GWAS on isolated population (62). In cattle, 50 samples and less were successfully used in a genome-wide eQTL mapping analysis and GWAS (63,64).
The similarity between the genomic heritability and the pedigree-based heritability indicates that the identified loci are the main genetic regulators of macrophage pro-inflammatory responses. The significant interaction network between the genes that are nearby significant SNPs and mitochondria-related genes that were associated with the NOindex, indicates the role of mitochondria in the regulation of the inflammatory profile of MDM. The functional annotation analysis of the network showed that "response to oxygen-containing compound" (GO: 1901700) was the most enriched biological process in the network. "Response to hypoxia" was one of the other enriched terms among 4 positional candidate genes. The SNPs nearby genes that are annotated with the response to hypoxia (APAF1, CD38, SMAD3, and SLC8A1) described a notable proportion of the observed variation (33.1%, 26.4%, 18.5%, and 16.9%, respectively). In addition, the SNPs that are nearby two novel genes, ENSBTAG00000025812 (an orthologue of CD38) and ENSBTAG00000045919 (has a role in oxidative stress-induced sensing, Reactome pathway: R-BTA-2559580), explained 32.4% and 30.6% of the observed variation, respectively. The SNP that is nearby PARPBP described 24.9% of the observed variation. Although this gene does not have a direct role in hypoxia, the product of this gene binds to PARP1 and enhance its function (65). PARP1 is an ADP-ribose polymerase which can be found in mitochondria with a role in the cellular response to oxidative stress (GO:0034599). Our previous study has shown that PARP1 is significantly downregulated in NO-high responder macrophages. The PARP1 gene has a positive role in mitochondrial depolarization (GO:0051901). Depolarization of the mitochondrial membrane is a mechanism that results in sensitivity to hypoxic condition and inability to cope with hypoxic conditions (66). By only including SNPs nearby genes that are associated with response to hypoxia (APAF1, CD38, SMAD3, SLC8A1, PARPBP, ENSBTAG00000025812, and ENSBTAG00000045919) in the general linear model (to include additive and non-additive genetic effects), 72.9% of the variation in NOresponse was described by these 7 SNPs. Therefore, the difference between the broad-sense heritability (0.80) of macrophage NOproduction and the R-Squared of the aforementioned model indicates that approximately 7% of the genetic regulation of macrophage response is controlled by the other 53 SNPs, and probably other loci that were not detected in this study. This large proportion of variation that is described by these 7 SNPs shows that response to hypoxia, probably via mitochondria, is the major genetic mechanism that controls the MDMs functional characteristics.
The results of our previous study showed that HIF1A was one of only 4 TFs found to be activated at both 3-and 18-hours postexposure to E. coli. In addition, the cluster containing the HIF1A was expanded at 18 hours in comparison to 3 hours. The synergic interaction between HIF1A and NFkB in macrophages can lead to increased production of NO - (67). We have previously shown a strong positive correlation between NOproduction and phagocytosis (17). It probably indicates that the pathway from phagocytosis to oxidative burst, hypoxia, HIF1A activation, pro-inflammatory response in macrophages is a dosedependent pathway.
The role of other candidate genes (i.e. PASK), biological process (i.e. "vesicle-mediated transport"), molecular function (i.e. "calcium ion binding", and "calcium:sodium antiporter activity") that were found in the current study should also be emphasized. PASK codes a protein that is a nutrient-sensing molecule with a role in regulating cellular energy balance via regulating lipid and glucose metabolism, mitochondrial respiration and phosphorylation (68,69). The role of Ca +depended exocytosis in macrophages is an example of the interaction between the "vesicle-mediated transport" and Ca ions that regulates phagocytosis which is the primary triggers of the oxidative burst-induced hypoxia (57,70).
In conclusion, the findings of the current study align with the previously published mechanisms in mice and humans regarding the role of mitochondria and response to hypoxia controlled by HIF1A of macrophages, thereby broadening the importance of this mechanism across species. Furthermore, these results show the effect of genetic polymorphism on regulation of response to oxidative burst-induced hypoxia is a major mechanism that controls macrophage pro-inflammatory responses such as NOproduction. It should be mentioned that, although this genetic predisposition is reported in MDM in the current study, the genes are shared among all leukocytes of the host. Hence, in an individual who carries unfavourable alleles, other cells, such as neutrophils and monocytes, that are capable of oxidative burst will likely also share an impaired response (71,72). Also, if the hypoxic condition is caused by other factors such as inflammation or solid tumours, then all the cells of the immune system in a hypoxic tissue (such as lymphocytes) may not be able to maintain their function and cope with the hypoxic condition (73). Therefore, the final consequence of genetic sensitivity to hypoxia will likely have a broader effect on the outcome of infection than just the effects that are directly or indirectly linked to macrophages. It should also worth noting that the results of the current study may not fully represent the function of macrophages in the host. Macrophages in the host are transformed from monocytes in a dynamic environment and work in collaboration with other cells of the immune system. These interactions are missing in the model that we used. Moreover, the phase of quantitative trait loci is different between breeds of cattle, and the results of the current study should be validated on other populations.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.

ETHICS STATEMENT
The animal study was reviewed and approved by The animal care committee of the University of Guelph.

AUTHOR CONTRIBUTIONS
ME and BM conceptualized the study and designed the experiments. ME and ST ran the experiments. ME and MS ran the statistical procedures. ME drafted the manuscript and ST, MS, and BM critically reviewed the manuscript. BM spearheaded the project as the principal investigator. All authors contributed to the article and approved the submitted version.

FUNDING
This research was funded by grants to BM from the Natural Sciences and Engineering Research Council (NSERC). This paper is also a contribution to the Food from Thought research program supported by the Canada First Research Excellence Fund. The funders did not have any role in the design, collection, analysis, and interpretation of data and in writing the manuscript. ME was supported by the Natural Sciences and Engineering Research Council via scholarship.