Genetic and non-genetic factors distinctly shape the variation of the immune response in cattle

Domain: biology
Published: 2026-01-15
Editorial outcome: major_revision

Abstract
Increasing agricultural production while reducing the reliance on synthetic inputs such as antibiotics is an important challenge. For cattle breeding, this implies better understanding the genetics underlying meat production and the immune response. Here, we use systems immunology to investigate the genetic and environmental drivers of immune variation in Belgian White Blue male cattle, a breed historically bred for meat production. While seasonality and other non-genetic factors account for much of the immune variation observed, genome-wide association studies identify loci with major effects on specific immunophenotypes. Genetics also emerges as the primary driver of cytokine production. Finally, we develop a predictive model linking genetic data to cytokine responses. Our findings support the selection of cattle with improved immunity and advance our understanding of mammalian immune variation.

Introduction
Given the profound impact of immune responses on health, it is necessary to understand their variability in a given population and how this variability relates to susceptibility to infections and immune-mediated diseases1–5. In humans and livestock species, this variation is determined by intrinsic (such as sex, age, and genetics), as well as extrinsic factors (such as diet, smoking, physical activity, microbial exposure)5. Despite recent advances, our understanding of the precise mechanisms and the extent to which these diverse factors contribute to the heterogeneity of immune responses among individuals and populations remains in its nascent stages6.

There is a growing consensus that between 20 and 40% of variation in immune traits is inherited in humans2,4,7. However, the majority of the underlying genetic variants remain largely unidentified7,8. Until recently, the study of genetic factors affecting immune responses has mainly been performed through the study of “in-born errors of immunity”, i.e. monogenic immune deficiencies9. While these approaches have successfully revealed critical immune pathways, they are less effective at elucidating how genetic variation acts in more subtle ways to regulate transcription or translation of key immune mediators4. In the past decade, strategies and technologies, collectively termed “systems immunology”, have emerged to explore the immune system with unprecedented precision7,8,10–16.

In general, systems immunology aims to identify the key components of immune networks, such as the immune cells that interact within specific systems and the molecular mediators they employ, and to determine how these components are altered when the system is perturbed12,17. To this end, these approaches were usually coupled with high-throughput “omics” technologies that measure mRNA, proteins, metabolites and epigenetic modifications that orchestrate gene transcription18,19. Systems immunology thus holds the potential to deepen our understanding of the mechanisms underlying disease and to provide insights into host-pathogen co-evolution18,19. Some recent projects, such as the Human Immunology Project Consortium (HIPC, https://immunespace.org), the Human Functional Genomics Project (HFGP, http://www.humanfunctionalgenomics.org), and the Milieu Interieur Project (https://www.milieuinterieur.fr), collectively provide a comprehensive grasp of interindividual variation in human immune responses. Interestingly, when integrated with the exploration of the genomic landscape governing these responses, these studies, identified some of the most important genetic variants associated with immune cell phenotypes and cytokines production7,10,15,20–27.

Strikingly, while systems immunology approaches have demonstrated their value in various contexts, they remain underutilized in livestock species. Although several pioneering studies have begun to explore systems-level immune responses in production animals, including pigs, cattle, and chicken28–36, these remain relatively rare and focus primarily on transcriptomic data or examine a small number of specific immunophenotypes. This is particularly surprising given that resistance or susceptibility to infectious diseases is a major determinant of production efficiency37. Moreover, the increased susceptibility of animals to infectious diseases is a threat to humans as it contributes to the emergence of zoonotic diseases and the propagation of antimicrobial resistance38. The scarcity of systems immunology studies in animals becomes even more perplexing when considering the notable advancements witnessed in the past decade in animal genetics. Genomic selection, for instance, has been introduced in several major livestock species, leading to an impressive doubling of genetic advancements in some cases39. These developments highlight the potential for a better understanding of immune responses in these populations through a formalized and integrative systems immunology framework.

The Belgian White Blue (BWB) cattle breed is the most muscular breed in the world and is attracting considerable interest to increase meat production in purebreds and industrial crossbreeding. Until now, selection programs in the BWB breed have focused on production traits and have not included health-related traits, particularly immune resistance to pathogens. Moreover, the respective effects of host and environmental factors, as well as the precise nature of the genetic variants that affect heterogeneity in the immune system in bovines, are largely unknown. The objectives of this study were to investigate the immune response in BWB cattle in a comprehensive manner by integrating both environmental and genetic factors that contribute to immune variation. Specifically, we aimed to identify genetic variants associated with immune traits that could be used to improve animal health through marker-assisted selection. Finally, we sought to demonstrate the potential of this genetic information to predict immune responses, thereby opening avenues for sustainable livestock production.

Results
Measurement of immune parameters in BWB cattle
Over the course of an 18-month period, we examined 246 BWB calves enrolled in the RESISTOMICS cohort, a bovine population cohort gathered to understand immune variation in cattle and raised in standardized conditions at the Center for Bovine Selection (CBS) in Ciney, Belgium. These calves, aged between 1 and 7 months upon entry in CBS, were evaluated for immune parameters at two distinct time points (28 days interval) by a broad immunophenotyping approach (Fig. 1a, b). Briefly, five unique flow cytometry staining panels were applied to immune cells, enabling the collection of 21 immune cell population percentages for each blood sampling (Fig. 1a and Supplementary Fig. 1–4, Supplementary Tables 1–3). The concentrations of 10 cytokines (IL-1β, IL-6, IL-4, IFN-γ, IL-17, TNF-α, IP-10, IL-8, MCP1 and IL-10) and of three antibody isotypes (IgM, IgG1, and IgG2) in serum were also measured (Supplementary Fig. 2a). In parallel, we measured the ex vivo production of 6 cytokines (IL-1β, IL-4, IL-6, IL-10, IL-17 and IFN-γ) by peripheral blood mononuclear cells (PBMCs) and whole white blood cells (WWBCs) stimulated by 15 different stimuli (Supplementary Fig. 2b). 66 cytokine-stimulus measurement pairs were excluded from further analysis because they showed no difference in comparison to the non-stimulated states (Wilcoxon signed rank test, p > 0.01) or had mostly (over 50%) values below the detection limit (Supplementary Fig. 2b). In total, the quantitative values of 148 immune parameters were compared for each animal at each time point.

Correlations between immune compartments
We first analyzed correlations between immunophenotypes to reveal some underlying patterns of variation at the population level. Those analyses were performed while correcting for age and season effects as several studies highlighted the importance of these factors on immune parameters40,41. At steady state in blood, concentrations of immunoglobulin isotypes were positively correlated with each other, while concentrations of IP-10 (CXCL10) was positively correlated with those of IL-10 and MCP-1 (CCL2) (Fig. 2a) (FDR < 0.05). Secondly, we explored the correlation between concentrations of circulating immunoglobulins and cytokines and the proportions of the different immune cell populations in blood (Fig. 2b). Notably, these analyses revealed a strong positive correlation between the percentage of a specific monocyte population circulating in the blood (Live CD172+CD16+CD14+), known as intermediate monocytes in human42, and the blood concentration of IL-10 (r = 0.419; FDR = 5.8e−08) (Fig. 2b). Other less strong and less significant correlations were also found, such as a positive correlation between CD21-IgM+ B lymphocytes and circulating IgM levels (r = 0.315, FDR = 0.000546).

We then performed unsupervised clustering of the concentrations of cytokines produced by immune cells in response to ex vivo stimulations (Fig. 2c). Four main clusters were observed. Cluster 1 encompassed IL-10 production in PBMCs and WWBCs. Cluster 2 mainly grouped monocyte-derived cytokines (IL-6 and IL-1β) in WWBCs. Cluster 3 comprised IFN-γ and IL-1β responses in PBMCs, while cluster 4 exhibited IL-17 and IL-6 responses in PBMCs. Although a distinct response to PMA/ionomycin stimulation was identified, it was surprising to note that the cytokine response primarily followed a cell-type (WWBCs vs PBMCs) and a cytokine-dependent pattern (Fig. 2c), not a stimulation defined one (Supplementary Fig. 5). These cell-type (WWBCs vs PBMCs) and cytokine-dependent patterns were confirmed when these parameters were integrated using Non-metric Multi-Dimensional Scaling (NMDS) as significant differences in cytokine productions were noted between PBMCs and WWBCs in both the first (FDR = 4.87e−09) and second (FDR = 2.27e−16) NMDS dimensions (Fig. 2d), suggesting the major influence of specific whole-blood components, such as neutrophils, on the variation in cytokine production. When comparisons were based on the cytokines, significant differences were also highlighted in both NMDS dimensions, especially when IL-10 production was compared to IFN-γ, IL-17, and IL-1β productions (Fig. 2d).

Finally, to comprehensively understand the interaction between immune cell populations in blood and ex vivo cytokine production, their relationship was evaluated across the 15 stimulations (Fig. 2e). It first identified a major negative correlation between the percentages of neutrophils and the production of most of the measured cytokines by WWBCs. Strikingly, percentages of circulating CD14+ monocytes exhibited positive correlations with the ex vivo production of monocyte-derived cytokines (IL-1β, IL-6 and particularly IL-10) in PBMCs and especially in WWBCs. In parallel, NK, T, and B cell percentages in blood displayed strong positive correlations with IFN-γ and IL-17 production after A. fumigatus stimulation of WWBCs (FDR < 0.05) (Fig. 2e). Conversely, negative correlations emerged between some monocyte-derived cytokine productions and lymphocyte populations. Specifically, CD3+γδ-, CD3+γδ+ and B cell subpopulations were inversely linked with HKLM- and HKMT-stimulated IL-6 production by PBMCs (FDR < 0.005). Similarly, those populations of lymphocytes negatively correlated with the production of IL-10 from PBMCs after various stimulations (Fig. 2e).

Contribution of inter- and intra-individual factors to immune variation
Out of the 246 animals, 211 underwent two separate samplings (which enabled an investigation into the evolution of immunophenotypes over time). Initial sampling (T1) occurred when the animals entered the CBS, followed by standardized vaccines or treatments (Supplementary Table 1). Re-sampling took place 28 days after entry (T2) (Fig. 1a). When we compared correlations between data obtained at T1 or between data obtained at T2, we observed similar patterns (Supplementary Fig. 6a, b). Interestingly, when we correlated immunophenotypes measured at T1 with those measured at T2, the cytokine-dependent patterns described above were still observed suggesting that animals display functional immune properties that are conserved over time (Supplementary Fig. 6c). To assess how the immune landscape of calves evolved over time, a longitudinal ANOVA model was developed, incorporating two independent variables: individual calf ID and sample time point (Fig. 3a and Supplementary Fig. 7)43. Each immune parameter was tested separately as response variable (Supplementary Fig. 7). The majority of the variance in each parameter was explained by a model depicting stable intraindividual immune profiles across longitudinal sampling (Fig. 3a). The median R2 across the parameters was 0.66, ranging from 0.50 to 0.89 (Supplementary Fig. 7). Notably, 65.97% of models maintained statistical significance (Bonferroni-corrected p < 0.05) even after correction for multiple testing (Supplementary Fig. 7a, b). Variation between repeated samples of the same individuals contributed minimally to the observed total variation, with a median proportion of R2 of 0.0091 (range 1.21e−05–0.067). In contrast, interindividual variation explained a much larger effect, with a median proportion of R2 of 0.663 (range 0.523–0.891) (Fig. 3a and Supplementary Fig. 7c, d), with some exceptions where sampling time (intraindividual variation) was significant (Supplementary Fig. 8). Consistent with this result, moderate but significant (FDR < 0.05) correlations between time points were found for most immunophenotypes, with a median correlation coefficient of 0.323 (range −0.0295 to 0.796). Interestingly, the temporal correlation between immunophenotypes measured in serum was significantly higher than for the production of cytokines after stimulation (FDR < 0.05; Fig. 3b). Furthermore, cytokine production following stimulation in WWBCs exhibited slightly higher interindividual variability, but lower intraindividual variability compared to PBMCs (p < 0.05) (Supplementary Fig. 7e). To further explore whether cumulative minor parameter changes could anyway lead to intraindividual variation over time, NMDS on paired samples was conducted. As expected, samples did not primarily cluster by sampling point (Supplementary Fig. 9a). However, NMDS2 revealed a difference between T1 and T2 (p = 0.0034) (Supplementary Fig. 9a). This difference originated from immune cells proportions (p = 0.005) and from the production of cytokines after stimulations (p = 0.00041) (Supplementary Fig. 9a). To analyze this sampling point variation further, an integrative analysis method called multiomics factor analysis (MOFA) was used to identify shared sources of variation and correlation between T1 and T2 (Supplementary Fig. 9b)44–46. This method enabled the detection of hidden parameters (namely latent factors (LFs)) that explain variance across the whole dataset and their potential relationship with differences between sampling points. For the immunoglobulins and cytokines in serum, percentages of immune cells and production of cytokines after stimulation, the MOFA model explained 2.11%, 66.05%, and 54.27% of the variance, respectively. Seven LFs were identified to explain the variance in the combined dataset. LF1 explained most of the variance observed across the percentages of immune cells (Supplementary Fig. 9b) but did not explain differences between T1 and T2 (Supplementary Fig. 9c). LF1 was mainly impacted by the variation in the percentages of neutrophils and T cells among samples (Supplementary Fig. 9d). LF6, explaining 8.99% of the total variance in the model, was the most significant LF that could distinguish between T1 and T2 (Supplementary Fig. 9c) (Bonferroni-corrected p = 1.37e−08). We therefore assessed the underlying factors contributing to LF6 to better understand the biology behind these sampling time-point differences. It revealed that IL-17 production and percentages of γδ+ and γδ- T cells were the main constituents of LF6 reflecting the fact that T2 samples displayed lower percentages of γδ+ T cells at steady state and an increased capacity of producing IL-17 upon stimulation (Supplementary Fig. 9e). Collectively, these findings suggest that, as seen in human, individual specificities are the main source of variation of immunophenotypes and the repeatability of the measures between sampling points is high. Anyway, some immunophenotypes display variation explained at least partially by sampling point and could be used in the future to monitor longitudinal variation.

Effect of age and season on the interindividual variation of immune parameters
Age significantly influences immune profiles10,16,43,44. We therefore investigated age effect on the data collected at T1 using a Spearman correlation test after correcting for the effect of season. Our analysis revealed both negative and positive correlations between immune parameters and age (Fig. 4a). As the variation of ages is limited in the cohort (Fig. 4b), we conducted permutation tests to underscore the robustness of the age-related impact on immune traits. 95.8% of immunophenotypes yielded results consistent with the actual estimates (Supplementary Data 1). Specifically, age was linked to an increase in CD21+IgM- B cells, accounting for 8.9% of the variance, and elevated serum IgM levels (FDR < 0.05) (Fig. 4a, b). Age was also related to increased serum TNF-α levels (FDR < 0.05). In terms of cytokine production induced by stimulations, the levels of IL-10 induced by several stimulations in PBMCs decreased with age (FDR < 0.05) (Fig. 4a, b and Supplementary Data 1). Despite the presence of some correlations, age exerted only a moderate influence (average adjusted R2 = 0.013, range from 2.35e−06 to 0.089) on immune profile variation within our cohort of healthy, young calves, which showed limited age variability.

Seasonality is another important factor that has previously been shown to influence immune cell populations and gene expression10,16,47. Similarly, we observed a seasonal effect on immune cell subpopulations and cytokine productions (Fig. 4c, d). Importantly, this was not related to some batch effect due to sampling point as these differences were observed from data collected on a weekly basis over two years and shown by month. During spring, B and T cell subsets were more abundant. Conversely, intermediate (CD14+CD16+ in CD172+) and non-classical monocytes (CD14− CD16+ in CD172+) were found at higher levels during winter and classical monocytes (CD14+CD16− in CD172+) during summer. Neutrophils levels displayed a peak during autumn (Fig. 4c, d and Supplementary Data 1). In line with immune cell patterns, cytokine production was similarly influenced by season. Interestingly, similar to classical monocytes, the levels of IL-6 (monocyte-derived cytokines) induced by stimulation of PBMCs were higher in summer (Fig. 4c, d and Supplementary Data 1). Through a resampling approach, 96% of immunophenotypes showed consistent results when compared to the original full dataset in more than 70% of iterations (Supplementary Fig. 10). Overall, these findings underscore the significant role of seasonal environmental factors in modulating variations in immune cell populations and cytokine production.

Identification of genetic variations that impact immunophenotypes
Based on these results, we conducted genome-wide association studies (GWAS) for each immunophenotype using a mixed linear model correcting for the effect of age and season to identify common genetic variants (minor allele frequency > 0.05) influencing immunophenotypes. This analysis unveiled ten independent genome-wide significant quantitative trait loci (QTLs) (p < 1.3e−07) from the immunophenotypes collected at T1 (Fig. 5a and Supplementary Data 2–4). When we repeated those analyses for immunophenotypes collected at T2, 6 regions could again be identified as significant at T2 through the identification of significant loci (p < 1.3e–07) in high linkage disequilibrium (LD) with loci identified at T1 (Supplementary Fig. 11a and Supplementary Data 2). Finally, we re-sampled an independent cohort of 239 calves to test the replication of the signals obtained from the cytokine measurements performed directly on serum (IL-8 and IL-10). In both cases, we again identified significant loci (p < 1.3e−07), each of which was in strong LD with the loci identified in the first cohort at both T1 and T2 (Supplementary Fig. 11b and Supplementary Data 2). We then conducted conditional GWAS by adjusting those immunophenotypes on the 10 leading associated variants and found no independent loci that reached genome-wide significance. Finally, we looked at the effect of these independent QTLs on all the immunophenotypes in order to reveal some patterns of effects (Fig. 5b). Most of these QTLs displayed trend, nominal or suggestive signals with other immunophenotypes suggesting their general importance in bovine immune responses. We then specifically investigated these 10 different loci.

Genetic variants associated with cytokines production. We first identified a QTL on chromosome 6 that governs serum IL-8 levels and accounts for 31.6% of the phenotypic variation. This QTL was associated with genetic variation in the vicinity of the IL-8 gene (CXCL8, Fig. 5c). These variants were located either in an intronic region of the gene or in 5’ or 3’ untranslated regions (Supplementary Data 3). The lead SNP identified from T1 dataset (6:88812260:G:T) was associated with significant differences in IL-8 concentrations in serum with the reference TT allele producing less detectable IL-8 (Fig. 5d). Then, based on transcriptome-wide gene expression profiling of resting whole blood cells from the same animals (Dr C. Charlier, ENA accession number: PRJEB94322), we analyzed expression quantitative trait loci (eQTL) governing the expression of IL-8 gene (Fig. 5c) and colocalized them with genetic variants associated with variation of IL-8 detection in serum. A strong colocalization evidence was revealed between the cis-pQTL and cis-eQTL of IL-8 (PH4 = 0.944) suggesting that cis-acting variants regulate IL-8 expression at steady state. Finally, in order to identify the possible consequences of this genetic variation on animal health, we tested the colocalization of those SNPs with loci recently associated with clinical traits in the Holstein breed48. These pQTL and eQTL colocalized with a region associated with cow livability, a health-related herd management trait that refers to the ability of a cow to remain alive and stay in the herd during a specific period (PH4 = 0.949 and 0.95, respectively), suggesting that it could have a significant impact on the health of those animals (Fig. 5c, e and Supplementary Data 4).

Similarly, the region associated with variation in IL-10 production was located on chromosome 16 and encompassed the IL-10 gene itself (Supplementary Fig. 12). Alleles associated with lower IL-10 production bear 4 coding SNPs (LD with lead SNP > 0.8, Supplementary Data 3). However, these coding mutations, potentially impeding protein recognition, do not appear to be the cause of the difference observed, as detection of the protein remained possible (albeit at lower levels) in animals homozygous for the mutant allele (Supplementary Fig. 12). Unexpectedly, this pQTL region is quite extensive, encompassing numerous SNPs that could account for the variability in IL-10 production. When we tried to colocalize it with cis-eQTL for IL-10 expression identified as described above (transcriptomic data from resting whole blood cells from the same animals), we did not detect any colocalization, suggesting that this difference in protein detection was not associated with gene expression differences. Interestingly, one of the mutations identified is located at a site potentially involved in the splicing of the IL-10 mRNA. In the future, additional analyses will be necessary to pinpoint the molecular mechanisms responsible for the observed defect in IL-10 production within specific genotypes and to decipher the potential functional consequences on the health of those animals.

Two additional loci functioned mainly as trans-pQTLs (> 1 Mb from the gene encoding the protein; Supplementary Fig. 12). The first, chr5:53539628:C:T, was associated with IL-1β production in PBMCs following imiquimod stimulation, explaining 11% of phenotypic variance. This intergenic SNP is located near SLC16A7 (MCT2), and individuals with the TT genotypes showed higher cytokine levels than CC homozygotes. The second locus, chr7:99347326:C:CA, was linked to IL-4 production after PMA/Ionomycin stimulation. The lead SNP lies upstream of FAM174A and ST8SIA4, and the minor allele was associated with increased IL-4 concentrations (p = 4.2e−08). These two signals were, however, not replicated with T2 dataset (Supplementary Data 2).

Genetic variants associated with proportions of immune cells. Six independent genetic loci were identified as being significantly associated with variations in the percentages of certain immune cell populations. Firstly, genetic polymorphism in the region of WC1 genes, containing 13 genes in cattle49, was associated with the percentage of WC1.2+ in γδ+ T cells (Supplementary Fig. 12). The results point towards a dominant allele negatively controlling the expression of WC1.2 antigens on γδ+ WC1+ T cells. It will be important to analyze the functional consequences of this variation in the future, given that the differential expression of these genes determines the type of γδ+ T cell subpopulations and that each of these subpopulations has specific functional properties against pathogens frequently encountered by cattle49.

A second signal linked to variation in the percentage of WC1.2⁺ cells among γδ⁺ WC1⁺ T cells was identified near the KCNJ2/KCNJ16 genes, which encode potassium channel components (Supplementary Fig. 12). Colocalization analysis showed moderate evidence (PH4 = 0.58) between the KCNJ2 cis-eQTL (identified as above from resting whole blood cell transcriptomes) and this trans-pQTL, suggesting that KCNJ2 expression may regulate this T cell subset and play a role in cattle immunity. This signal was, however, not replicated with T2 dataset (Supplementary Data 2).

Thirdly, genetic variation on chromosome 11 upstream of the GNLY gene (distance = 143,003 bp) was found to be linked with the percentage of CD8+ in NK cells (p = 2.2e−08) (Supplementary Fig. 12). GNLY is an essential effector molecule of killer lymphocytes involved in immune response defense against bacterial and viral infection50.

Fourthly, genetic variations around the MHC locus on chromosome 23 were associated with variations of the percentage of CD4+ in CD3+ cells (distance = 105,595 bp, p = 1.16e−07) and of the percentage of CD8+ in CD3+ cells (distance = 421,160 bp, p = 2.49e−09) (Supplementary Fig. 12). The LD between the lead SNPs linked with the percentage of CD4+ in CD3+ cells and the percentage of CD8+ in CD3+ cells is minimal (r2 = 0.556076), suggesting that the percentages of CD4+ and CD8+ T cells were determined by distinct variations, even if the quality of genetic information at that location is poor.

Finally, genetic polymorphism on chromosome 16 (chr16:18662313:C:T) was associated with variation in the different monocyte populations among CD172+ cells (Supplementary Fig. 12). However, no SNP hit was highlighted at T2.

In addition, we identified 93 genome-wide suggestive signals (p < 2.6 × 10−6; Supplementary Data 5) warranting further investigation.

Determinants of interindividual variation of cytokine profiles at baseline or after stimulation ex vivo
We then systematically quantified the effects of genetic and non-genetic factors on the immunophenotypes. To assess the relative contributions of age, season, and genetics (based on significant and suggestive SNPs) to interindividual variation in immune cell proportions, we used PERMANOVA to estimate the variance explained by each factor (Fig. 6a and Supplementary Data 6)51,52. Notably, without controlling for covariates, season (21.5%) emerged as the primary contributor to this interindividual variation, overshadowing the effect of the GWAS loci (12.2%). Using a similar approach, we calculated the proportion of variance explained by age, season, immune cells, immunoglobulins, and the same genetic loci for the entire profile of cytokines measured in serum or induced by stimulations (Fig. 6a). In contrast to immune cell profiles, genetic factors accounted for the largest proportion of variance in cytokine traits, consistent with human studies53. By comparison, immune cell proportions and season had moderate effects (5.3% and 7.5%, respectively) across most cytokine–stimulus combinations, while circulating immunoglobulin levels and age had minimal impact. These factors cumulatively explained 21.6% and 36.2% of the variance in cytokine production at steady (serum) and activated states, respectively (Fig. 6a). Subsequently, we calculated the proportion of variance explained by genetic (using here all genetic markers from SNP array (38,828 SNPs across the genome)) and non-genetic factors for each independent immunophenotype (Fig. 6b, c and Supplementary Data 7). Interestingly, and in line with previous study10, we found that genetic factors explained a significantly greater proportion of variance in T-cell immune traits compared to B-cell traits (Student’s t test, p = 0.01) (Fig. 6b). In contrast, unlike previously reported in human7, we did not highlight differences in variance explained by genetic factors between innate and adaptive immune cells (Student’s t test, p = 0.582). This result could be attributed to the small size of the cohort and to the fact that we used the percentage of immune cell populations as parameters instead of Mean Fluorescence Intensity (MFI) measurements of individual receptor expression.

In addition, we determined which was the dominant factor (i.e. the factor accounting for more variance than the others) for each immunophenotype. For proportion of immune cells, season was the dominant factor in 57% of the cases and genetics was in 43% (Fig. 6b and Supplementary Data 6). In contrast, interindividual differences in stimulated cytokine production were found to be primarily explained by genetics in 77% of cases. In 10.5% of the cases, season was the dominant factor and for 11.5% of them, it was immune cells (Fig. 6c and Supplementary Data 7). To further gauge the significance of the estimated contribution of genetics (h2 SNP) to the observed total variance, we performed 1000 permutations of the animal IDs in our dataset and then applied the analysis pipeline to the permuted data to obtain the empirical distribution of the estimates of h2 SNP (null distribution). The estimate of h2 SNP from the original data was then compared to the estimate of h2 SNP from the permuted data. Among the 21 immune cell phenotypes and the 104 stimulated cytokines immunophenotypes (Figs. 6c), 7 and 35 exhibited significant h2 SNP estimates, respectively (p < 0.05, Fig. 6b, c). In particular, 12 of these belonged to IL-10 production induced by stimuli in WWBCs and PBMCs.

Finally, we compared data on stimulus-dependent cytokine production from the two different types of cells (PBMCs and WWBCs) from the same individuals used in our study. This highlighted that genetic factors, proportions of immune cell populations, and the concentration of immunoglobulin levels had a greater effect on stimulus-dependent cytokine production in WWBCs than in PBMCs (Mann–Whitney U test, p < 0.05) (Supplementary Fig. 13a). When we compared the effect of these different factors in function of the stimulations used, we did not identify any major effect of any of these factors in function of a particular type of stimulus (Supplementary Fig. 13b). In contrast, while genetics was still the major factor explaining variation for the production of each cytokine, it appeared that season explains more variance for IL-1β production (and to a lesser extent for IL-6) than for the other cytokines (IFNγ, IL-10, IL-17 and IL-4) (Supplementary Fig. 13c). For example, season explained 18.7% and 16.3% of the variance of IL-1β production by PBMCs following P. multocida and M. haemolytica-stimulation, respectively. These results also show that, although genetics is the most important factor contributing to interindividual variation in cytokine production (mean heritability = 0.241, ranging from 1e−06 to 0.756) (Fig. 6c), the majority of the variation in our population could still be explained by non-genetic factors, which is consistent with human studies8,53.

Cytokine production in response to stimulations can be partly predicted on the basis of the immune cell composition and genetic factors
Previous studies have proposed that cytokine responses to ex vivo stimulation could be predicted from steady state measurements (i.e. genetics and blood cell composition) through machine learning approaches with various degrees of success53–56. We tested therefore the possibility of predicting cytokine production using a cross-validation strategy that involved randomly splitting of the data into training (70% of the animals) and validation (30% of the animals) sets multiple times to predict cytokine productions. We initially explored whether immune cells proportions could be used to predict the cytokine response to ex vivo stimulations as previously shown in human (Fig. 7a, b and Supplementary Data 8)55. The average correlation coefficient between predicted and measured stimulus-induced cytokine productions was 0.214 (range −0.083 to 0.542). Interestingly, in comparison to PBMCs, cytokines induced by stimulation in WWBCs are more likely to be predicted by immune cell composition (Wilcoxon signed-rank test, p = 1.574e−09). Notably, five cytokine-stimulus pairs exhibited prediction accuracy exceeding 0.45 in WWBCs, including A. fumigatus-induced IL10 production (mean 0.542), HKMT-induced IL-10 production (mean 0.537), A. fumigatus-induced IL-17 production (mean 0.495), A. fumigatus-induced IFNγ production (mean 0.478), and HKMT-induced IL-6 production in WWBCs (mean 0.454) (Fig. 7b). Secondly, genetic factors (corresponding here to SNPs with a p-value < 5e−05 in the GWAS analysis) were tested to predict cytokine production in response to stimulation (see Methods)53. When we used genetic information from GWAS analyses performed on the entire cohort, correlations of approximately 0.75 were observed between predicted and measured values (Supplementary Fig. 14). However, in this setup, the test animals contribute to the identification of GWAS signals and thus to model construction. Therefore, in a second approach, to ensure the models were completely independent from the test cohorts, we repeated the GWAS analyses for each model using only the data from the corresponding training cohort (Fig. 7a, b). In this setting, the average correlation between predicted and measured stimulus-induced cytokine production was 0.096 among the cytokine-stimulus pairs (range −0.102 to 0.421). However, some, such as the prediction of IL-10 production, consistently showed a correlation of the order of 0.3. Finally, the inclusion of both immune cell proportions and genetic factors in the model yielded to notable enhancement in predictive accuracy compared to genetics only (FDR < 0.05, Student’s t test) with an average correlation of 0.151 between predictions and measures (range −0.066 to 0.422) (Fig. 7b). To further validate the accuracy of prediction, we trained the model on 70% of the cohort with cell proportions measurements from T1 sampling and predicted cytokine production for the other 30% of the cohort at T2 (Fig. 7c, d and Supplementary Data 9). The predictions made using immune cell proportions, genetics, or genetics plus immune cells resulted in comparable accuracies at T2 than the ones obtained for T1, with mean values of 0.204, 0.070, and 0.205, respectively. In general, models based on cell population percentages gave better results than those based on genetics alone (Fig. 7). However, a clear improvement was observed for IL10 production estimates, which could be expected given the prior identification of a highly significant pQTL region for production of this cytokine (Fig. 5 and Supplementary Fig. 12).

Methods
Ethics
All procedures were conducted in accordance with the European Convention for the Protection of Vertebrate Animals used for Experimental and Other Scientific Purposes (CETS 123) and were approved by the University of Liège Animal Ethics Committee (protocol number 2509).

Population cohorts
246 BWB male calves housed at the center for bovine selection (CBS) in Ciney have been enrolled in the first cohort of this study. All calves were sampled on the first week after entry (time point T1) and 211 individuals were sampled 4 weeks later (time point T2). The 35 animals not included in the second sampling (T2) were excluded primarily due to logistical constraints. Importantly, no mortality occurred during the study period. On the week after T1, all calves were vaccinated and treated with anti-parasite drugs (Supplementary Table 1). To validate certain GWAS hits in an independent cohort of animals, an additional 239 male calves were sampled one month after entry at the CSB.

Cell Processing
Blood was sampled at the jugular vein in K2 EDTA tubes with a Vacutainer system. For WWBCs purification, blood was incubated in red blood cell (RBC) lysis buffer according to the manufacturer’s instructions (Thermo Fisher), with some minor modifications. Briefly, blood (5 mL) was first incubated for 7 minutes with lysis buffer (10 mL) before the addition of 30 mL PBS containing 2% fetal calf serum (FCS). Cells were spun down (335 g for 5 minutes at 4°C) and once again treated with 500 µL of RBC lysis buffer for 3 min and spun down at 335 g for 5 min at 4°C after addition of 30 mL PBS containing 2% FCS. Cells were resuspended in 10 mL PBS 2% FCS and counted.

For PBMCs isolation, 5 mL of blood was mixed with the same volume of PBS 2% FCS and added on top of 5 mL of Ficoll (Ficoll-Paque PLUS, Fischer Scientific) before centrifugation at 697 g for 10 min at RT, no brake. Interphase containing purified PBMCs was transferred to a new tube and washed with 10 mL PBS 2% FCS and centrifuged at 335 g for 5 min at 4°C. Cell pellets were treated with 500 µL of RBC lysis buffer for 3 minutes, washed with 45 mL PBS 2% FCS and spun down at 335 g for 5 min at 4°C. Cells were resuspended in 10 mL PBS 2% FCS and counted.

Flow cytometry analysis
Before staining, most unlabeled antibodies were labeled with lightning-link kits (Abcam) following the manufacturer’s instructions. 106 WWBCs were stained with standard panels of immunophenotyping antibodies (Supplementary Table 2 for a list of antibodies, clones, fluorochromes, and manufacturers) and incubated for 15 minutes at 4°C with a volume of 50 µL per well. Cells were washed with 200 µL fluorescence-activated cell sorting (FACS) buffer (PBS containing 2 mM EDTA, 0.5% bovine serum albumin (BSA), and 0.1% NaN3) and centrifuged at 335 g for 5 min at 4°C. Next, cells were incubated with second mix in FACS buffer for 15 min at 4°C and washed with PBS with a volume of 200 µL per well. Finally, cells were stained with Zombie aqua (BioLegend) according to the manufacturer’s instructions and suspended in 200 µL FACS buffer. Data were acquired with a BD LSR Fortessa X-20 flow cytometer equipped with 50-mW violet 405-nm, 50-mW blue 488-nm, 50-mW yellow-green 561-nm and 40-mW red 633-nm lasers and a ND 1.0 filter in front of the forward scatter (FSC) photodiode, using BD FACSDiva software (BD Bioscience). Data was analyzed using FlowJo software (TreeStar).

Each sample was analyzed with 5 multicolor antibody panels: (i) NK cells, (ii) γδ T cells, (iii) CD8+ CD4+ T cells, (iv) Monocytes, and (v) B-neutrophils illustrated in Supplementary Fig. 1. For each panel, the live leukocytes subpopulation (P1) was discriminated using side-scattered area (SSC-A) and forward scatter area (FSC-A) to select cells, and forward scatter height (FSC-H) and FSC-A to select single-cell population, and finally Zombie aqua staining to select live cells. In panel 1, within the lymphocytes selected by FSC-A/SSC-A, γδ T cells populations (γδ+, CD8+ in γδ+, WC1+ in γδ+, WC1.2+ in WC1+ γδ+) were identified. In panel 2, T cells populations (CD3+ γδ-, CD8+ in CD3+ γδ-, CD4+ in CD3+ γδ-, CD25+ in CD4+ CD3+ γδ-) were identified. In panel 3, NK cells (CD335+, CD2+ in CD335+, and CD8+ in CD335+) were characterized. In panel 4, granulocytes and lymphocytes were discriminated by FSC-A and SSC-A and further characterized by staining with the CH138A antibody (identifying granulocyte and therefore mostly neutrophils in blood) and for IgM and CD21 expression to identify B cell subsets. In panel 5, subpopulations within monocytes were separated as CD172+, CD14+ in CD172+, and CD16+ in CD172+.

Stimuli
Bacteria. HKMT (Heat-killed Mycobacterium tuberculosis), HKLM (Heat-killed Listeria monocytogenes) and Flagellin (Recombinant flagellin from S. typhimurium) were purchased from InvivoGen and used at a final concentration of 100 μg/mL (HKMT), 108 cell/mL (HKLM) or 1 μg/mL (Flagellin). Mannheimia haemolytica and Pasteurella multocida were produced and heat inactivated by the Bacteriology Laboratory from the veterinary faculty (ULiège) and used at 4 × 105 cfu/mL (Mannheimia haemolytica), and 2 × 103 cfu/mL (Pasteurella multocida).

Fungi. Aspergillus fumigatus extract was purchased at the UK National Institute for Biological Standards and Control (NIBSC) and used at 1% dilution.

Virus. Bovine herpesvirus 1 (BoHV-1) and BoHV-4 were amplified on Madin-Darby bovine kidney cells (ATCC CCL-22) and purified on sucrose gradient. The purified productions were heat inactivated at 56°C for 60 min and used at 105 pfu/mL (BoHV-1) and 5.104 pfu/mL (BoHV-4). Parainfluenza 3 (PI3) was amplified on bovine turbinate (BT) cells and concentrated by ultracentrifugation (100,000 g, 2 h at 4°C) after removal of the cell debris by low-speed centrifugation (1000 g, 10 min at 4°C).

TLR ligands and non-microbial stimuli. Pam3CSK4 (Pam3Cys-SerLys4), PolyI:C (High molecular weight), LPS (E. coli serotype 055:B5), Imiquimod (R837), and ssRNA were purchased from InvivoGen and used at a final concentration of 1 μg/mL (Pam3CSK4, LPS, Imiquimod, and ssRNA), 10 μg/mL (PolyI:C). Phorbol 12-myristate 13-acetate (PMA) and Ionomycin were purchased from Sigma and used at 50 ng/mL and 1 μg/mL, respectively.

Ex vivo stimulations
2.5 × 105 WWBCs or PBMCs were added to a 96-well plate and cultured with stimuli in 250 μL of IMDM + glutamax (Gibco) 10% FCS, 2% penicillin-streptomycin, 1% amino acid, 50 μM β-mercaptoethanol for 48 h at 37°C and 5% CO2 except for PMA and Ionomycin which were incubated at 37°C for 5 h. Supernatants were collected and stored at −20°C until used for ELISA.

ELISA analysis of cytokines and immunoglobulin levels
Concentrations of bovine IL-1β, IL-6, IL-4, IFN-γ, IL-17, TNF-α, IP-10 (CXCL-10), IL-8, and IL-10 as well as serum levels of IgG1, IgG2, and IgM were determined by ELISA following manufacturer’s recommendations (see Supplementary Table 3 for a list of antibodies, clones, and manufacturers). Briefly, the 96-well microtiter plates were coated with capture antibody in PBS overnight at 4°C and then were blocked for 2 h at 20°C with 200 μL of 0.1% PBS-Tween 20 for anti-Ig ELISA or 5% bovine serum albumin (BSA) for cytokines ELISA. 50 μL of serum dilution or two times diluted supernatant were added into the wells and incubated for 2 h at 20°C, then the detection antibody was added into wells after washing four times with wash buffer (PBS containing 0.1% tween 20) and incubated for 1 hour at room temperature (RT). Finally, streptavidin-HRP was incubated for 30 minutes at RT. After each incubation, the plates were washed four times with wash buffer. The plates were developed using 3,3’,5,5’-tetramethylbenzidine. The reaction was stopped after suitable color development by the addition of 50 μL of 0.2 M sulfuric acid and the absorbance at 450 nm was read using an ELISA plate reader. Detection ranges were 31.25 pg/mL-2000 pg/mL (IL-6), 15.625 pg/mL-1000 pg/mL (IL-1β and IL-4), 25 pg/mL-6250 pg/mL (IFN-γ), 62.5 pg/mL-4000 pg/mL (IL-10), 0.6 pg/mL-632 pg/mL (IL-17), 125 pg/mL-8000 pg/mL (TNF-α), 0.123 ng/mL-30 ng/mL (IP-10), 0.39 ng/mL-25 ng/mL (IL-8), 0.048 ng/mL-200 ng/mL (IgG1 and IgG2) or 0.024 ng/mL-100 ng/mL (IgM).

Genomic DNA extraction, SNP array genotyping, and 57,979 SNPs
Genomic DNA was extracted from whole blood (EDTA) on a KingFisher™ instrument following standard procedures. SNP genotyping for all the animals was conducted at the GIGA Genomics platform (ULiège) using two custom bovine medium density (MD) BeadChip arrays (Illumina) following the manufacturer’s instructions. The two arrays, namely EuroGenomics MDv1 and MDv2, interrogate, respectively, 50,809 and 57,979 genomic positions. For each animal, final genome-wide genotypes were obtained with the GenomeStudio Genotyping (GT) Module from Illumina. After intersecting the two array positions and excluding markers on sex chromosomes, a total of 43,173 genetic markers from SNP array genotyping was kept for subsequent analysis. In parallel, a total of 156 animals were selected for whole-genome sequence by random selection. The sequencing was performed by Illumina HiSeq 2000 instrument in the GIGA Genomics platform. The bioinformatic analysis to call variants from the sequence data is fully described in Gualdron Duarte et al.85 and is summarized here. Sequence reads were mapped to the bovine reference genome ARS-UCD1.2 using BWA-MEM software (version 0.7.5)86. After sorting and marking duplicates, variant calling was performed on individual recalibrated BAM files using HaplotypeCaller (GATK4 v4.1.7.0)87 to obtain individual genomic variant call format (GVCF) files. Joint genotyping of all these GVCF files was then done using GenotypeGVCFs (GATK4 v4.1.7.0)87. After applying the variant quality recalibration procedure with (GATK4 v4.1.7.0)87, variants with a quality score below the 97.5 threshold, minor allele frequency (MAF) < 0.01 and multi-allelic sites were filtered out. Genotype imputation to the sequence level was then performed for the individuals genotyped on an MD array using ShapeIT4.2 in combination with Minimac485, using sequenced Belgian Blue individuals as a reference. Finally, a total of 11,486,822 SNPs were included in the genome-wide association studies, with 246 calves from T1 sampling. As a validation cohort, 211 animals from T2 sampling were used.

Quantification and statistical analysis
Data pretreatment. The Shapiro-Wilk normality test was used to test the normality of both raw (Supplementary Fig. 3) and log-transformed data (Supplementary Fig. 4) data. A threshold of p > 0.05 was employed to define a normal distribution. Percentages of the various immune cell populations, were normalized using an inverse rank transformation (IRT) algorithm10. A log2 transformation was used to transform immunoglobulin levels as well as all cytokine levels (Supplementary Fig. 4). In cases where the response fell outside of the ELISA measurement range, the corresponding lower or upper value of the ELISA plate detection range was used as the measured value. Raw cytokine production followed a non-Gaussian distribution (Supplementary Fig. 3) while some log2 transformed cytokine concentrations followed a normal distribution (Supplementary Fig. 4). The Wilcoxon signed rank test was used to compare the significance of cytokine responses after different types of stimulation to the unstimulated condition (IMDM medium). Considering the fact that some cytokines were not normally distributed, the data were further normalized using the IRT algorithm88, as described above.

Correlations among immune parameters. To determine the correlation between immune parameters including cytokines concentrations, immunoglobulin isotypes concentrations, and immune cells percentages, we first corrected the normalized immunophenotypes for age and seasonal effects using a linear model. Pairwise Spearman correlation analysis was used to calculate associations using the normalized and corrected immunophenotypes, and an unsupervised hierarchical clustering approach (hclust) was used using these correlations as distance. The same methodology was used to calculate the relationship between immune cells percentages and immunoglobulins concentrations, as well as the relationship between immune cell percentages and cytokine production. After correction for multiple testing (foreach category of dataset), the significance of the correlations was considered (FDR < 0.05)89. As part of the preprocessing for multidimensional scaling analysis (MDS) using monoMDS from the vegan package or hierarchical clustering with heatmap in stats package, Euclidean distances were calculated from correlation matrices (using Spearman).

Genome-wide association analyses. Prior to the GWAS, values of immunophenotype data that were more than 1.5 times the interquartile range above the upper quartile and less than 1.5 times the lower quartile were removed (removal of outliers). The association analysis was used to identify common genetic variants that influence interindividual variation in immunophenotypes using leave-one-chromosome-out (LOCO) mixed linear model, as implemented in GCTA (v1.26.0) software90, with age and season as covariables. The approach used by Druet et al.91 was applied to estimate genome-wide significance at p < 0.05 using Bonferroni corrections for multiple-testing. By comparing uncorrected and corrected p-values, the number of independent tests was set as 382,000. As a result, the p-value threshold for genome-wide significance association was 1.3e−07 (0.05/382,000). Suggestive association threshold was set to p < 2.6 × 10−6 (1/382,000). Inflation or deflation in p-values was assessed by genomic inflation factor (λ) and also visually inspected by quantile–quantile (Q–Q) plot. λ is calculated as the median of the χ2 test statistics (1 degree of freedom) divided by its theoretical median under the null distribution. To further evaluate whether the genetic loci of immunophenotype contain several independent signals, we conducted conditional and joint analysis utilizing the stepwise model selection approach in GCTA-COJO. The proportion of variance in phenotype explained by a given SNP (PVE) was calculated92. Subsequently, Variant Effect Predictor (VEP) (http://www.ensembl.org/Tools/VEP) was used for the identification of the candidate genes located in significant loci and predicted functional consequences of identified mutations93.

Colocalization analysis. Immune-related QTLs identified in this study were colocalized with eQTLs derived from transcriptomic data from resting whole blood cells from the same animals (ENA accession number: PRJEB94322) or with genomic regions associated with health-related phenotypes in cattle using the R package ‘coloc’ (cran.r-project.org/web/packages/coloc). For genomic regions associated with health-related phenotypes, summary statistics from large-scale single-marker GWAS for seven cow health-related variables in Holstein cattle, including cow livability, mastitis, ketosis, hypocalcemia, displaced abomasum, metritis, and retained placenta were extracted from a recent study48. The chromosomal locations were initially transformed from UMD3.1 to ARS-UCD1.2 using liftOver. Bayesian test was then used to assess the likelihood of QTL and clinical traits sharing the same underlying causative variants94. The Bayesian co-localization produces posterior probabilities for five hypotheses related to a shared variant between the two traits: (H0) no association with either trait; (H1) association with expression only; (H2) association with disease only; (H3) association with both traits but independent causal variants; and (H4) association with both traits driven by a shared causal variant. Strong evidence of colocalization was defined as a posterior probability of ≥ 0.8 for H4.

Estimation of explained variance. Feature selection based on permutational multivariate analysis of variance using distance matrices (PERMANOVA) was firstly used to estimate the contributions of different factors to interindividual variations of each dataset (e.g., cytokines production and immune cells)52 (Fig. 6a). In brief, interindividual variations for each dataset was estimated using vegan’s (version 2.5.5) adonis function with 1000 times permutation95. Only response variables and candidate explanatory variables with a permutational FDR < 0.05 were kept. Only suggestive and significant SNPs from SNP-genotyping array were retained as genetic variants in this analysis. To address collinearity of the selected factors, the set of candidate explanatory variables was correlated to itself. When candidate explanatory factors displayed Spearman correlation greater than 0.4, the factor displaying the lowest correlation with the immunophenotype was removed from the model53. All independent explanatory factors were further included in PERMANOVA to estimate the combined contribution to interindividual variation. In order to estimate the variance explained by the different factors on individual immunophenotypes (Fig. 6b, c), a correlation-based factor-selection approach was used53. In brief, if a factor (age, season, antibody concentration and frequency of immune cell population) had a significant association (Spearman FDR < 0.05) with the immunophenotype, it was added to the set of candidate explanatory variables. To address collinearity of the selected factors, the set of candidate explanatory variables was selected using the same procedure stated above. The adjusted R2 was used as the measure of explained variance, accounting for the inflation caused by adding predictors. The heritability was estimated using genome-based restricted maximum likelihood (GREML) implemented in GCTA software, with age and season as covariates. Meanwhile, we permuted the immunophenotypes and covariate datasets 1000 times randomly before running the GWAS with these datasets to generate 1000 random profiles for each immunophenotype. Using the pipeline described above, the heritability for the permuted immunophenotypes dataset for each run was calculated. This analysis produced a distribution of 1000 heritability estimates for each immunophenotype. A measured estimate was considered significant if it was in the top 5% of the permuted distribution of immunophenotype estimates.

Cytokine level prediction. For machine learning, Elastic Net was used to predict cytokine levels96. 206 individuals sampled at both T1 and T2 were selected (5 animals out of 211 were excluded based on missing values). These animals were firstly randomly divided into training (70% of the dataset) and test sets (30% of the dataset) to validate the prediction. Predictions (always based on measured made at T1) were either made for cytokine production measured at T1 (Fig. 7a, b) or at T2 (Fig. 7c, d). The training dataset consisted of either percentages of immune cells alone, or genetic information alone, or immune cell percentages + genetic information. Importantly, in order to define the training genetic information, we used for prediction either SNPs with a p-value < 5 × 10−5 obtained from the whole cohort (Supplementary Fig. 14) or SNPs with a p-value < 5 × 10−5 obtained from GWAS analyses performed on the training animals only (Fig. 7) in order to avoid that the same animals would contribute to the GWAS results while being part of the test set. Briefly, for the training dataset, a 5×5-fold cross-validation approach was applied and the trained model with the lowest MSE (Mean Squared Error) for each run was chosen for the prediction of test set. All the procedure was repeated 5 times and the prediction accuracy was assessed by calculating Spearman correlations between the measured cytokine levels and the predictions. A t test was then used to investigate whether there was a significant improvement in power for cytokine production prediction when immune cell percentages and genetics were used together in the model.

Other statistical analyses. To investigate the effect of sampling time on immune variation, linear regression modeling was performed using the base R function (lm), with R2 proportions calculated using the relaimpo package43. The majority of the modeling shown uses a categorical variable (such as subject sampling point or calf number) to model a continuous variable (immune parameter). lm provides an ANOVA model in this case (for example, immune parameter ~ sampling point + calf’s ID). The Bonferroni correction was used for multiple testing, as implemented in base R. Wilcoxon signed rank test or dependent t-test were used to compare the different immunophenotypes from the same individual at different time points (Supplementary Fig. 7). To test for an association of immunophenotypes with age, we normalized the immune traits as described above, then used a linear model to correct for the effects of season. The significance and direction of the correlations were determined using Spearman correlation analysis. A cosinor model was fitted with the month of the year, in a one-year cycle period, using age as covariates to determine the effect of season on the immune traits measured. The fitted model was compared to an alternative regression model that did not include these terms, and significance (multiple correction testing; FDR < 0.05) was determined using a variance test between the full model and the alternative regression model16,47. After controlling for multiple testing (FDR < 0.05), significance was declared. MOFA was performed in the R statistical environment using the MOFA+ package (version 1.0)97, and latent factors were determined using MOFA’s TopWeights function.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "‘A systems immunology approach reveals distinct roles of genetic and non-genetic factors in shaping variation of immune responses in cattle. Shifang Li et al\nIt is now well accepted that human immune systems vary significantly, with some individuals being more susceptible to infection and immune mediated pathology than others. While it is known that genetic variation contributes significantly to these differences in human immune system function, genetics are responsible for only a proportion of the variation in immune activity. Age, environment, stress and nutrition are also critical influencers. It should be no surprise that other species display similar immune variation and susceptibility to infection which should have implications for breeding programmes most of which have not previously taken the opportunity to develop breeds and strains that might be more resistant to infection, therefore reducing the need for antibiotic usage. To quote the paper under review ‘the selection process for these animals has prioritized production traits, without considering health and disease resistance.’ This paper is one of the first to systematically search for genetic contribution to immune variation in cattle. It describes significant analysis of large amounts of data collected on two occasions from 246 Belgian White Blue calves.\nIt is important that such comprehensive systems -immunology studies be carried out on economically relevant species in order to understand the variation in immune activity of livestock, given the likelihood of their involvement in the emergence of pathogens that could cause future epidemics and pandemics. This present study is an important, interesting and timely investigation that will be of interest to diverse groups of scientists and researchers not exclusively immunologists. Overall, the article is well written, with high quality figures and statistical analysis and incisive and relevant discussion.\nHowever a number of issues need to be addressed:\n1. The choice of breed: Belgian White Blues are heavily selected. Are they more inbred than other breeds of cattle? Would they therefore have less variation in their genetics and immune responses? Are there other limitations to the use of this breed? Even ‘though only one gene is responsible for the double muscling that is characteristic of the breed’, this physiological feature is likely to have significant metabolic implications which may themselves impact on immune activity. This issue should be addressed.\n2. Some basic analysis of the genomic structure of Belgian Blues should be presented to compare with other breeds, eg CNV status has been demonstrated to be important in human & mouse studies; comment on CNV in immune genes of BWB?\n3. Are the immune cell populations in BWB similar to what has been published for other breeds?\n4. The abstract should be rewritten – over half is background and no data are presented\n5. The technological approaches used in the study are not introduced appropriately nor their limitations discussed\n6. Age has been demonstrated to significantly impact human immune function – this cannot be addressed here. This limitation should be highlighted.\n7. In Lines 101-103, we read that ‘the respective effects of host and environmental factors, as well as the precise nature of the genetic variants that affect heterogeneity in the immune system in bovine, are largely unknown’ However this issue is not really addressed in the paper; some comment needs to be made in the discussion.\n8. The Results sections are very long, and include conclusions/discussion points/summarizations that could be moved to ‘Discussion’. Concise reporting of results will make this section easier to read.\n9. The ‘identification of genetic variations that impact immunophenotypes’ section could certainly be separated out into smaller sections.\n10. The authors go back and forth describing immune cell/genetic associations and cytokine stimulation/genetic associations in the ‘identification of genetic variations that impact immunophenotypes’ section; this section structure could be improved – maybe a table would help. For example, associations are found with gamma delta populations but differential responses to cytokine stimulations between them are described in a different section.\n11. The lack of stimulating specific responses found should be discussed? Is this species or breed specific?\n12. Given the amount of immune data collected and analysed, 10 independent QTLs correlated in the GWAS seems low. Is the figure lower than anticipated? Does it reflect the interbreeding of BWBs?\n13. Were genetic associations carried out with PBMC stimulations or whole blood stimulations? This is not very clear. A comment should be included in the discussion that ‘whole blood’ systems provide clearer data for reproducible analysis https://pubmed.ncbi.nlm.nih.gov/28943400/\n14. Where do the health trait data come from? From the cohort of this study? or is this a separate analysis? This section is unclear.\n15. Other references that should be included and would help interpretation and discussion – especially around non genetic effects that might help explain immune variation\na. Variation in cytokine responses in cattle: https://pubmed.ncbi.nlm.nih.gov/34446770/\nb. Stimulant specific responses in cattle: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC9169910/\nc. Seasonal and vitamin D effects on cytokine responses and immune cell populations\nhttps://www.nature.com/articles/s41598-023-37427-z https://www.nature.com/articles/s41598-021-98343-8.\nd. https://pubmed.ncbi.nlm.nih.gov/37407588/\n16. Why was the IL-1B association with SLC16A7 placed in supplementary? This seems an important and interesting finding – particularly in light of the comments above about possible immunometabolic influences of the BWB specific loss-of-function mutation in the myostatin (MSTN) gene.\n17. ‘In addition to significant associations (p<1.3e-07), we identified 93 genome-wide suggestive signals at p<2. 417 6e-06 (Supplementary Fig. 10 and Supplementary Table 8) 418 that warrant further analysis’. This section is confusing. What do the authors mean by ‘suggestive signals’? If they did not meet p value cut-offs, maybe less could be written about them as this section is already quite long.\nMinor Issues\nThere are many grammatical and linguistic errors throughout the manuscript - some examples are presented below:\nL63 ‘limited to decipher’…\nL67’ also directly in humans ‘ … The word ‘directly is redundant here\nL77’ rewrite: ‘prospect to refine our understanding of …’\nL85 ‘they remain notably absent in animals’ should be ‘remain notably absent from animals’ L90 ‘This becomes…’ what becomes?\nL298 ‘genetic variations that impacts’…\nL 605 What is ‘livibility’?\nL 614-616 needs a reference.", "reviewer_2": "Li et al present a systems immunology study of cattle. These approaches have been previously applied to human populations and mouse models so it is welcome to see this being extended to another species. Extensive immune phenotyping was performed for circulating immune cells, baseline cytokines and antibodies, and induced cytokines after microbial stimulation. They also performed genome wide SNP analysis and integrated the different data sets through different statistical methods. The study is generally well performed and the results will be of interest to the field. However there are some concerns that need to be addressed to improve the confidence of the reader in the results presented and in their interpretation by the authors.\nWhat was the justification/rational for having two distinct cellular populations isolated from whole blood; PBMCs and WWBC ? Why not just use WWBC for consistency and avoid differences due to cell composition that will be a result of the cellular preparation? This could be the reason driving the two main clusters in Fig 1e & f. This is also supported by the correlations in Fig 1g where it is mostly negative correlations in PBMCs and positive correlations in WWBC with the main difference due to neutrophils, which normally are removed during PBMC preparation. It could also explain why the samples do not cluster more by stimuli as is mostly observed in other studies as acknowledged by the authors. What if you analyze the two cell preparation types separately, do you then see stronger stimulus specific clusters ?\nAlso you cannot call either of these cell types as “tissues” given that PBMCs don’t exist as a natural population anywhere in vivo, never mind constituting a tissue.\nDid the authors measure the same cytokines in the FCS batches used for cell preparation? This is less of a concern in human or mouse studies but may well cause interference in the measures of cattle cytokines some of which are surely present in FCS. This could potentially explain why there were no real stimulus specific clusters which is an unexpected result.\nThe longitudinal analysis is interesting as this is often missing from human studies. It makes sense that IgG levels are among the most stable phenotypes as they would not be expected to change so much over such a time period, and not surprising that induced cytokines are the most variable but the authors fail to mention the technical variability that will be strongly impacting the induced cytokine responses. Could they address this by comparing the intra-individual variability between the WWBC and PBMC responses, this would be a justified use of these two different sample types.\nWhich vaccines and anti-parasite drugs did the cattle receive after T1 and would this be expected to effect the immune phenotypes measured at T2 ?\nThe seasonal results are also interesting, can the others say whether there were known housing, feeding or antibiotic treatment changes of the cattle associated with these different seasonal periods? Similarly can they exclude any change of reagents, lots, experimental operators to exclude underlying experimental factors that be inadvertently aligned with season.\nHow do the authors explain the failure to replicate half of the GWAS hits at the second timepoint given the previously described low intra-inter individual variance, this is a concern.\nThe genetic associations with IL8 and its link with good health, does it mean low or high IL8 is good for cattle health ?\nCould the lack of an eQTL for IL10 just be due to differences in kinetics of gene expression that are not the same for the protein. This is impossible to interpret as no details of the transcriptomic data are given beyond “based on transcriptome-wide gene expression profiling of resting PBMCs from the same animals (Dr C. Charlier, unpublished data)”. I think we have moved beyond the referencing of unpublished data sets, and I believe if the data is to be cited it has to included or at the very least referenced with a publicly available reference (point for editor to decide).\nThe decomposition of variance explained is quite limited by the very narrow age range (6 months) and focus on a single sex (males) thereby emitting real life variability due to these two important factors, this limitation needs to be acknowledged. In fact it could also be presented as a strength of the study that by excluding variability from these two factors the authors are better powered to detect and quantify variability from genetics and seasonality.\nThe part on “prediction” of immune phenotypes is not very convincing, a mean correlation value between observed and predicted phenotypes of 0.2 or 0.15 is a long way from predicting immune responses. Again no evidence of any correcting for false positives so just by chance there will be some correlation values higher than the average.\nA lof of immune phenotypes are analyzed considering the sample size so the use of FDR is very important. More details should be provided on exactly how and when FDR was applied to statistical testing which was not so clear.", "reviewer_3": "• What are the noteworthy results?\nThis study uses a system immunology approach coupled with genome wide association studies to examine non-genetic and genetic factors affecting immune variation in Belgian White Blue cattle. The study carries out an impressive set of different types of statistical analyses applied to an extensive set of immune parameters collected at 2 time points (28 days apart) plus ex-vivo cytokine responses to 15 different stimuli, of 246 BWB calves, coupled with genomic data. Noteworthy results included the clustering of immune parameters and the relationship between different immune compartments identified by NMDS, the finding that inter-individual variation in immune variation exceeded within-individual variation, genetic variation in immune phenotypes and identification of genetic variants associated with immunophenotypes (although I question the validity – see below), and that non-genetic factors have a higher impact on immune phenotypes than host genetic factors.\n• Will the work be of significance to the field and related fields? How does it compare to the established literature? If the work is not original, please provide relevant references.\nThe data collected and systems approach (consisting of a large variety of statistical tests) used to analyse the data, and the results certainly offer new insights into immune system variation in livestock. But I believe the authors grossly over-state the novelty and significance. They claim that systems immunology is a novel field in livestock, however I don’t believe this is the case. There are several studies in livestock that use similar approaches to those applied in this study and warrant citation in my view (see some examples below). The authors don’t cite a single livestock study, which suggests to me that they haven’t done their homework and investigate the state of the art in livestock systems immunology.\nhttps://www.nature.com/articles/s41541-018-0078-0\nhttps://link.springer.com/article/10.1186/1753-6561-5-S4-S32\nhttps://policycommons.net/artifacts/15817350/systems-immunology-approaches-to-decipher-individual-variability-of-immune-capacity/16708209/\nhttps://www.cabidigitallibrary.org/doi/full/10.5555/20123233872\nhttps://www.sciencedirect.com/science/article/abs/pii/S0301622601003207\nhttps://www.cambridge.org/core/journals/genetics-research/article/on-the-need-for-combining-complementary-analyses-to-assess-the-effect-of-a-candidate-gene-and-the-evolution-of-its-polymorphism-the-example-of-the-major-histocompatibility-complex-in-chicken/149BA4746F376F7A0241F745EB8884F0\n• Does the work support the conclusions and claims, or is additional evidence needed?\nIn my view the study lacks the statistical power and rigour in the analyses to back up the strong statements of the authors. In particular 246 individuals is an extremely small number of individuals for a GWAS. Evidence and theoretical expectations indicate that such small GWAS are prone to result in false positives, which may explain the strong GWAS hits observed in this study. Furthermore, I was missing information about the relatedness structure of these calves (e.g. how many full- or half-sibs) This information needs to be provided and included in the GWAS, as it can confound SNP effect estimates. It seems to me that the estimated SNP effects could be grossly inflated as relatedness, the influence of polygenic effects and other relevant factors (e.g. cell type / stimulus used, time of measurements, pen, indoor / outdoor rearing ) that may contribute to variation in the immune parameters are not accounted for in the models. I have similar concerns regarding the use of machine learning approaches for such a small, and standardized dataset, especially for predicting continuous (rather than discrete) response. In fact, the predictive ability of the models wasn’t very high (r=0.2 is not good predictive power in my mind), even if the authors claim that cytokine response can be predicted. The low correlation, and also the fact that the time points are only 1 month apart suggest to me that this is a gross over-statement. It would be interesting to see if this claim holds after known immune challenge.\n• Are there any flaws in the data analysis, interpretation and conclusions? Do these prohibit publication or require revision? Is the methodology sound? Does the work meet the expected standards in your field?\nSee my concerns above. I believe these would need to be addressed first before the study can be published. I believe the dataset is too small to provide conclusive results, in particular regarding the GWAS results. However, the study may serve as a proof of concept study.\n• Is there enough detail provided in the methods for the work to be reproduced?\nI’m not qualified to comment on the description of the molecular methods. The statistical methods are adequately described in my view. The authors provide extensive supplementary information.", "reviewer_4": "General comments\nThe authors report very interesting results on the individual variability of a wide range of immune traits in Belgium White Blue cattle selected for muscle development and being acknowledged as the most muscular breed worldwide. They claim to be first reporting large scale studies on the variability of immunity traits in livestock (here in cattle), using methods referred to as systems immunology. It is actually disappointing that the authors do not refer to other publications in livestock (cattle, pigs) on this topic with similar approaches but only to reports in humans and in mice. In addition, references reporting the influence of the host genetics on immunity traits would be expected.\nThe immune traits comprise a set of parameters measured in the serum (cytokines, non-specific immunoglobulins), a set of parameters for the counts of blood cell sub-populations by flow cytometry, and a set of parameters that measure cytokine levels after ex-vivo stimulations of either whole white blood cells or PBMCs by various inducing factors. This collection of immune parameters collected at two time points for 211 animals is impressive (148 parameters per animal, twice). The results provide data on genetic and non-genetic factors (age, season) found involved in the observed individual variability of the measured immune traits, and the most striking results are the candidate genomic regions found significantly associated with some immune traits that include response traits after ex-vivo stimulation. These results suggest that genome information might predict immune response variabilities for further use in breeding programs. The effect of age is reported as limited. All animals are young and the two sampling points quite close suggesting that immune phenotypes are quite stable at least in the young age. The authors claim that immunity trait levels at T2 are very well predicted by data collected at T1. However, the level and accuracy of the prediction are not very high for all traits and predicting immune traits one month later has limited impact.\nThe authors explain that they have adopted a systems immunology approach that usually aims to integrate multi-omics data with immune phenotypes in order to address the complexity of biological systems in a holistic vision with the aim to identify properties of the systems and further find clues on how the systems may react to perturbations. Indeed, in this report, the authors estimated correlations between immune phenotypes and between two time points and searched for associations between immune trait variations and genome polymorphisms, with successful results.\nMajor concerns\nA total of 246 animals born in various farms were housed in the collecting farm at difference ages (1 to 6 months) and at various seasonal periods over 18 month spanning January to December with pics in March, May and June, September and October. All animals were sampled the first week after entry (T1) and a subset of 211 animals were collected 4 weeks later (T2). The authors have to explain why 35 animals had to be removed from the experiment (mortality, infectious disease?) and how the health status of the animals was assessed for the experiment. Have the authors included the farm origin of the animals in the model?\nIt is not always very clear whether the data used for the analyses come from T1 or T2 sampling.\nIt would be expected that immune parameters are categorized by their main contribution either to innate or to adaptive immunity.\nThe animals were all vaccinated between T1 and T2 but we do not know at what time post-vaccination the T2 sampling has been done. Is it the same time lapse between T1 and T2 for all animals? Have the authors looked at possible influence of vaccination on parameters measured at T2? What vaccines were injected? Have the authors measured specific antibodies? Interesting to monitor individual variability to the response to vaccines.\nThe order of the result presentation needs to be improved in order to help the understanding. The correction for age and season is important for the correlation estimates but it seems that these corrections were done by relying on publications and not on the analysis of the effect of these two factors on individual variability. This needs to be explained. Actually, it would be interesting to know whether age and season influence all immune parameters or groups of parameters. Thus, it would be appropriate to have the section on the influence of age and season (non-genetic parameters) before the section on correlations, and to end with the genetic analysis and the prediction studies. It would help to know about non-genetic effects that are corrected all over the analyses.\nThe authors have produced many figures that are insufficiently self-explanatory even with the associated legend. They must select the main figures that best highlight the results presented in the main text.\nIn some result sections, the authors have included comparisons with other published data that should feed the discussion, highlighting similarities and differences with other published results.\nLine 197: the title refers to findings that have been reported for long. It would be expected to have data on the value ranges and the coefficient of variation per parameter.\nThe genetic data that show candidate regions are very interesting. The number of animals being quite small for conducting GWAS, have the authors checked all possible biases? To what extend do you think that your GWAS results might be biased by your low number of animals? You mention the MAF only in the materials and methods part but it may be worth mentioning it earlier and maybe even providing the number of animals that carry the minor frequency allele for your highest associated SNP.\nSince phenotypic data at T1 and T2 are highly correlated at intra-individual level, how do the authors explain the differences between the number of QTLs found significantly associated? Have the authors the capacity to estimate heritability of each parameter?\nThe discussion with possible outcomes is limited. Among topics to be discussed: what about immune parameter values to select and what associations with production traits and the sensitivity to disease over animal life course; based on correlations ranging from 0.5 to 0.9 between values from T1 and T2, and also data on latent variable stability between T1 and T2, what recommendations for time sampling; based on immune phenotype correlation estimates, what strategy could be used to simplify the phenotyping and propose methods for large scale immune phenotyping. The discussion on animal age showing differences with other studies should be discussed regarding the young age of the animals and the close sampling time points. Discussion on how prediction of response using basal information directly obtained from blood data together with knowledge on genetic determinism of the variability of the basal parameters could be jointly exploited for animal breeding is missing.\nOther comments\nThe authors use alternatively re-stimulation or stimulation in the text: suggestion to only refer to stimulation.\nIntroduction:\n- L.100: the Achilles heel of this breed is more the calving ease than the sensitivity to diseases\n- L103-108: It looks more like an abstract than a description of what is coming. Please rephrase. And do not hesitate to add explicitly the objectives of the study.\nResults:\n- L114: Why is there such an age difference at the entrance? Especially considering the low number of animals below 5 months, is it worth keeping them?\n- L198: What happened to the calves not in T2? Are they those with health issues? Did they die? If so, did they display some specificities at T1?\n- L299: provide the information that those data are based on whole genome sequence association studies\n- L306: to what calculation does the suggestive thresholds correspond to?\n- L440 and L524: Why do you use on the one hand only significant and suggestive SNP and on the other hand a new threshold (where does it come from by the way)? Is it impossible to use the entire genome for the calculation available at L440?\nMaterials and methods:\n- L920: double spacing before “genotyping”?\n- L1040: Why Spearman and not Pearson here?\nFigures:\n- Not enough self-explanatory. Abbreviations should be explained\n- Fig1a and 1b should be separated from the rest of the figure since 1a represents the sampling design and the overall protocol, and 1b provides information on the number of animals per age and sampling month. These two figures do not fit with the general title that refers to results on correlations.\n- Fig. 2c: Meaning of an immunological distance?", "reviewer_5": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts."}, "author_response": "Response: We sincerely thank the reviewer for the positive and encouraging comments. We fully concur that immune variation in livestock species has been historically underexplored, despite its clear relevance to disease resistance, antibiotic usage, and overall animal health, productivity and welfare.\nWe particularly appreciate the reviewer’s emphasis on the importance of applying comprehensive systems immunology approaches to economically significant livestock species, especially given their potential role in global health and zoonotic disease emergence.\nHowever a number of issues need to be addressed:\n1. The choice of breed: Belgian White Blues are heavily selected. Are they more inbred than other breeds of cattle? Would they therefore have less variation in their genetics and immune responses? Are there other limitations to the use of this breed? Even ‘though only one gene is responsible for the double muscling that is characteristic of the breed’, this physiological feature is likely to have significant metabolic implications which may themselves impact on immune activity. This issue should be addressed.\nResponse: We thank the reviewer for raising this insightful point concerning the choice of breed and its potential implications for genetic diversity and immune variation.\nIt is true that the Belgian Blue cattle breed has been intensively selected, including the fixation of a few major genes, but selection has been performed at levels similar to those of other typical cattle breeds.\nConsequently, the Effective population size (Ne) is low to moderate at around 100, similar to that of many other breeds (see, for example, Druet et al., 2013, Animal Genetics 44(6)758-762, for these aspects). Inbreeding levels are high at around 0.10 over the last 50 generations, which is similar to levels observed in other breeds (see, for example, Solé et al., 2017, Genet Sel Evol 49, 92).\nOverall, diversity is high enough to perform selection and association studies. Thus, following rigorous quality filtering, we identified over 15 million variants (SNPs and indels) in a sample of 250 individuals (the genome size being 2.7 Gb). Many successful association studies have been performed in different cattle breeds, demonstrating that there is sufficient genetic variation, and the Belgian Blue cattle breed has similar levels of variation (Charlier et al., 2016, Genome Res 26(10):1333-1341). In agreement with this, we recently identified several large-effect variants with a GWAS in that breed (Gualdron-Duarte et al., 2023, Genet Sel Evol 55, 83).\nIn the present study, the limiting factor is likely to be the sample size as mentioned in the discussion (Lines 597-599). The identification of several significant signals in the present study demonstrates that there is sufficient variation and that large-effect variants are segregating.\nIn conclusion, to respond to the reviewer's comment, Belgian Blues have been selected in a similar way to other cattle breeds, they have similar levels of inbreeding and genetic diversity to other cattle breeds and this variation is sufficient to identify large effect genetic variants by GWAS with our sample size.\nWith respect to the double-muscling phenotype, which is primarily driven by mutations in the MSTN (myostatin) gene, we agree with the reviewer that downstream metabolic and physiological consequences-such as alterations in energy balance, muscle metabolism, and tissue oxygenation-may modulate immune responses. We have incorporated a statement into the Discussion (Lines 623-626) to acknowledge this possibility, emphasizing that metabolic characteristics may interact with immune phenotypes and should be taken into consideration when interpreting our results and extrapolating findings to other breeds.\n2. Some basic analysis of the genomic structure of Belgian Blues should be presented to compare with other breeds, eg CNV status has been demonstrated to be important in human & mouse studies; comment on CNV in immune genes of BWB?\nResponse: We thank the reviewer for this comment regarding the genomic structure of Belgian Blue cattle and the potential contribution of copy number variation (CNV), particularly within immune-related genes.\nWe fully agree that CNVs are increasingly recognized as important modulators of immune function in both humans and mice, and we acknowledge their likely relevance in livestock species. Although the present study focused on SNP-based genome-wide association analyses using imputed sequence data, we recognize that structural variants, including CNVs, represent a significant and complementary layer of genomic variation that may underlie immune phenotypic diversity.\nWhile CNV analysis was not performed in the current study, we have now included a statement in the Discussion (Lines 638-640) to highlight this limitation and to propose that future investigations incorporating CNV profiling, particularly in immune-related loci, may yield additional insights into the genetic architecture of immune variation in cattle.\nIn addition, outside the scope of this study, we have very recently received funding to pursue our research in this area, and we plan to sequence 50 BWB genomes using “long reads” technology, specifically to address the CNV issue, as suggested by the reviewer.\n3. Are the immune cell populations in BWB similar to what has been published for other breeds?\nResponse: To the best of our knowledge, panels similar to the one used in our study to determine the proportion of leukocytes have not been reported in other cattle breeds. However, regarding peripheral blood mononuclear cells (PBMCs), the frequencies and distributions of major immune cell subsets in BWB calves, including T cells, B cells, monocytes, and natural killer (NK) cells, were broadly consistent with those previously reported in Holstein cows, such as 30-day-old Holstein calves and Holstein cows at 260 days of gestation and at calving (Casaro et al., Vet Immunol Immunopathol 2022). Although some differences were observed in specific subpopulations, these were relatively minor and may reflect differences in genetic background, selection pressures, environmental influences, or management practices. The reviewer comment underscores the importance of conducting similarly comprehensive immunophenotyping studies across diverse cattle populations. This has now been mentioned in the Discussion (Lines 626-628).\n4. The abstract should be rewritten – over half is background and no data are presented\nResponse: In response to this comment, we have revised the Abstract to provide a more balanced and informative summary of the study. The updated version significantly reduces the length of the background section and places greater emphasis on the core findings, including the extent of inter-individual immune variation, the relative contributions of non-genetic and genetic factors, and the identification of genomic regions associated with specific immune traits. We have also tried to comply as far as possible with the length constraints imposed by the journal (150 words).\n5. The technological approaches used in the study are not introduced appropriately nor their limitations discussed\nResponse: As requested by the reviewer, we tried to more clearly introduce and contextualize the technological approaches employed in our study and we discussed more explicitly their limitations.\nFirstly, although the results section appears long, we have kept the first paragraph “Measurement of immune parameters in BWB cattle” in its entirety to give the reader a detailed overview of all the experiments carried out.\nSecondly, we have expanded the Discussion to include some dedicated sentences that outline additional limitations associated with these methodologies (Lines 597-641). Specifically, we discussed: (1) the statistical and biological constraints of GWAS conducted in moderate-sized cohorts, including limited power to detect variants with small effect sizes and the potential inflation of effect estimates (Lines 597-599); (2) the interpretive challenges associated with ex vivo cytokine stimulation assays, including susceptibility to variability introduced by cell viability, batch effects, and differences in stimulation intensity. Moreover, we discussed how whole blood stimulation could improve our readouts (Lines 599-607); (3) the benefits that could be gained from conducting a longitudinal study (Lines 607-617); (4) the limitation of using only blood samples (Lines 617-620); (5) the absence of multi-omics integration (e.g., transcriptomics, proteomics, metabolomic), which restricts our ability to infer underlying regulatory mechanisms and limits the depth of systems-level insight (Lines 620-623); (6) the breed-specific context (Lines 623-632); (7) the exclusion of outliers, which could prevent us from identifying very interesting profiles (Lines 632-638); (8) the fact that we are not currently addressing structural variants (Lines 638-640).\nOverall, these revisions provide a more comprehensive perspective on the strengths and limitations of our approach and help position the study within the broader context of immunogenomic research.\n6. Age has been demonstrated to significantly impact human immune function – this cannot be addressed here. This limitation should be highlighted.\nResponse: We agree that age is a well-established determinant of immune function in humans and likely exerts a similar influence in cattle. However, the relatively narrow age range of our study population, primarily between 4 and 6 months, limits our capacity to comprehensively assess age-related variation in immune responses.\nIn response, we have explicitly acknowledged this limitation in the Discussion section (Lines 609-615). We note that the restricted age distribution of animals in this cohort precludes a thorough investigation of immune development across life stages or the onset of immunosenescence. We therefore recommend that future longitudinal studies incorporating a broader age range will be necessary to fully elucidate the role of age in shaping immune function in cattle.\n7. In Lines 101-103, we read that ‘the respective effects of host and environmental factors, as well as the precise nature of the genetic variants that affect heterogeneity in the immune system in bovine, are largely unknown’ However this issue is not really addressed in the paper; some comment needs to be made in the discussion.\nResponse: We partly agree with the reviewer. Indeed, although our study is carried out on a limited cohort and needs to be confirmed by further studies, the questions mentioned by the reviewer seem to us to be addressed. Indeed, (1) we highlight that immune variation for the parameters measured is primarily of inter-individual origin (Figure 3). (2) Then, and in particular in Figure 6, we assess the proportion of variance in the various immune parameters explained respectively by age (although our cohort is not well designed for this), season and genetics. Season appears to be a major determinant, which is also shown in figure 4. (3) Finally, with regard to genetics and in particular SNPs, we show that this is the major factor explaining cytokine production at steady state in serum, but especially after ex vivo restimulation. We have also identified 10 genetic loci that are strongly associated with variations in certain immunophenotypes (Figure 5).\nHowever, we agree with the reviewer that the Discussion could have been more elaborated. In response, we have revised the Discussion to explicitly acknowledge these limitations and to clarify the extent to which our findings inform these questions. Specifically, we discussed (Lines 597-641): (1) the current limitations in our ability to functionally characterize the identified genetic variants, which remain candidate loci requiring further validation through fine mapping and functional studies; (2) the inherent complexity of disentangling host genetic influences from environmental factors, even within a well-controlled experimental design, and the need for integrative analytical frameworks and larger cohorts to enable more precise decomposition of variance components. In particular, we stressed that future studies should take more non-genetic factors into account, notably, BMI, sex, viral infection history, vaccine history, parasitic infections, metabolism, and epigenetic regulation.\nFinally, we emphasize that this study represents an important initial step toward elucidating the respective contributions of genetic and non-genetic factors to immune heterogeneity in livestock. We highlight the importance of future research to fully resolve the mechanistic underpinnings of immune variation in cattle for selection of more resistant animals but also for a better understanding of mammalian immune responses.\n8. The Results sections are very long, and include conclusions/discussion points/summarizations that could be moved to ‘Discussion’. Concise reporting of results will make this section easier to read.\nResponse: We agree with the reviewer that this section can appear long and that many results are reported. Moreover, some sentences were clearly discussion elements. We have tried to improve it and make it more concise.\nIn particular, we have tried to reposition the elements of discussion within the discussion itself. This has been done for elements that could lead to a general discussion. For some more specific aspects, comparison with specific results in the literature has been kept within the results section.\nSpecifically, the following elements were repositioned in the discussion:\n(1) The discussion on the correlations between cytokine production is now in the Discussion section (Lines 519-539).\n(2) In line with the reviewer’s suggestion, we have thoroughly revised the part of “Identification of genetic variations that impacts immunophenotypes” in the Results section to enhance its clarity and focus (Lines 278-389). Specifically, we have substantially shortened the result section regarding suggestive significant loci (Lines 388-389).\n(3) The discussion on the prediction of immune responses (Lines 512-518).\nIn the end, the length of the results section has been reduced, and we believe we've also improved its organization, as requested by the reviewer.\n9. The ‘identification of genetic variations that impact immunophenotypes’ section could certainly be separated out into smaller sections.\nResponse: We agree that the section detailing the identification of genetic variants associated with immunophenotypes was overly dense and would benefit from greater structural clarity through subdivision into focused subsections.\nIn response, we have reorganized this portion of the Results into two delineated subsections, each addressing a distinct aspect of the analysis: (1) genetic variants associated with cytokines production; (2) genetic variants associated with immune cells.\nMoreover, we have removed results regarding suggestive significant loci.\nWe hope that this restructuring enhances the readability and logical flow of the section.\n10. The authors go back and forth describing immune cell/genetic associations and cytokine stimulation/genetic associations in the ‘identification of genetic variations that impact immunophenotypes’ section; this section structure could be improved – maybe a table would help. For example, associations are found with gamma delta populations but differential responses to cytokine stimulations between them are described in a different section.\nResponse: We agree that this is a difficult part to write in order to keep clarity and make all the connections. In response, as stated above, we have restructured the relevant section of the Results into clearly labeled subsections that address:(1) genetic variants associated with cytokines production; (2) genetic variants associated with proportions of immune cells.\nWe hope that this revised structure enhances the clarity of the presentation and facilitates a more intuitive understanding of the scope and implications of the genetic findings.\n11. The lack of stimulating specific responses found should be discussed? Is this species or breed specific?\nResponse: We agree that the limited detection of strong stimulation-specific genetic associations is a noteworthy observation that warrants further discussion.\nThis point is addressed in the modified Discussion section (Lines 519-539). We acknowledge that this pattern may reflect species-specific or breed-specific immunogenetic regulation. In fact, similar findings have been reported in human functional genomics studies, including in healthy African populations (Boahen et al., Am J Hum Genet 2022) and in European cohorts with Lyme borreliosis (Botey-Bataller et al., Nat Commun 2024), where limited stimulus-specific genetic effects were also observed.\nConversely, other studies - such as those conducted in Dutch healthy populations - have reported functional clustering according to the type of stimulus applied (Li et al., Cell 2016). These divergent findings suggest that multiple factors may shape the degree of genetic control over stimulation-specific immune responses.\nWe acknowledge that it remains unresolved whether the limited number of stimulation-specific associations reflects biological characteristics of the species or breed, or methodological limitations such as sample size, timing of measurements, or assay sensitivity (discussed in detail in the Discussion section, Lines 597-641). We now highlight this issue as a priority for future studies involving diverse cattle breeds, expanded cohorts, or additional time points following immune stimulation.\n12. Given the amount of immune data collected and analysed, 10 independent QTLs correlated in the GWAS seems low. Is the figure lower than anticipated? Does it reflect the interbreeding of BWBs?\nResponse: We agree that this outcome is indeed somewhat lower than anticipated, and we believe several factors likely contributed to this: (1) sample size and statistical power: with a cohort of 246 animals, the study was likely underpowered to detect associations with small to moderate effect sizes, particularly given the complex, polygenic nature of immune traits. It is therefore plausible that additional true associations exist but fall below the threshold of statistical significance; (2) stringent multiple testing correction: we applied conservative genome-wide significance thresholds to minimize the risk of false positives. While this approach increases the robustness of detected signals, it also reduces the number of loci that reach statistical significance, (3) we think that a methodological variation (e.g. working with purified cells and not whole blood as mentioned by the reviewer) may have decreased our ability to identify loci of interest.\nIn contrast, as indicated in the response to the reviewer's first comment, we believe that, although in general the level of inbreeding can affect the number of QTLs by GWAS analysis, this does not seem to be the case here compared to other cattle breeds that display the same level of inbreeding and similar levels of genetic variation.\nWe have expanded the Discussion section to explicitly address some of these points (Lines 597-641) and to emphasize that, although the number of significant QTLs is limited, the findings provide an important foundational reference. Moreover, we highlight the potential for future studies in larger or more genetically diverse cohorts to build upon these results and uncover additional loci relevant to immune variation in cattle.\n13. Were genetic associations carried out with PBMC stimulations or whole blood stimulations? This is not very clear. A comment should be included in the discussion that ‘whole blood’ systems provide clearer data for reproducible analysis https://pubmed.ncbi.nlm.nih.gov/28943400/\nResponse: We confirm that cytokine stimulation assays were conducted using both peripheral blood mononuclear cells (PBMCs) and whole white blood cells (WWBCs). Genetic association analyses were performed on cytokine levels measured in both sample types following stimulation with 15 different agents. We have ensured that this methodological detail is conveyed consistently throughout the manuscript.\nMoreover, in response to the reviewer’s suggestion, we have also expanded the Discussion to highlight the advantages of using whole blood stimulation systems (Lines 599-607). As noted in the referenced study, whole blood assays offer greater physiological relevance and improved reproducibility across studies, owing to minimal sample processing and the preservation of native cellular interactions. We concur that these features enhance the robustness of the cytokine data and further support its application in genetic association analyses.\n14. Where do the health trait data come from? From the cohort of this study? or is this a separate analysis? This section is unclear.\nResponse: Unfortunately, the health trait data referenced in the manuscript do not originate from the same cohort used for immunophenotyping. Rather, they were obtained from external public datasets comprising a larger commercial population of Holstein cattle (Lines 999-1003). These data were utilized to explore potential overlap or relevance between the immune-associated genomic regions identified in our study and loci associated with health traits in an independent population.\nTo address potential ambiguity, we checked that the relevant sections of the Results (Lines 316-318) and Methods (Lines 999-1003) clearly distinguish the origin of the health trait data and prevent confusion regarding cohort composition. Additionally, we have included a brief discussion of the limitations inherent in cross-cohort comparisons and emphasized the importance of future studies that integrate immune phenotyping and health trait data within the same individuals to enable more robust and biologically informative analyses (Lines 588-596).\n15. Other references that should be included and would help interpretation and discussion – especially around non genetic effects that might help explain immune variation\na. Variation in cytokine responses in cattle: https://pubmed.ncbi.nlm.nih.gov/34446770/\nb. Stimulant specific responses in cattle: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC9169910/\nc. Seasonal and vitamin D effects on cytokine responses and immune cell populations\nhttps://www.nature.com/articles/s41598-023-37427-z https://www.nature.com/articles/s41598-021-98343-8.\nd. https://pubmed.ncbi.nlm.nih.gov/37407588/\nResponse: We thank the reviewer for recommending these references. We agree that they provide important context for interpreting our findings, particularly in relation to non-genetic influences on immune variation in cattle. In response, we have incorporated the following references into the Discussion to enhance and substantiate our interpretations:\nThe study on variation in cytokine responses in cattle (Reid. C et al., Sci. Rep. 2021) is now cited in the context of ex vivo stimulation using whole blood rather than purified cells (Line 607).\nThe reference on stimulant-specific immune responses (Lesueur J. et al., Front Immunol, 2022) use a very nice methodological approach and has been cited for that (Line 607). However, it was more difficult for us to directly compare the conclusions obtained because the analyses were carried out for different purposes and the results are therefore not directly comparable.\nStudies addressing seasonal and vitamin D–related effects on immune cell populations and cytokine responses (Reid. C et al., Sci. Rep. 2021; Flores-Villalva S. et al., Sci. Rep. 2021), have been cited in our discussion of seasonal influences on immune traits, supporting the biological plausibility of the non-genetic effects detected in our dataset (Line 622).\n16. Why was the IL-1B association with SLC16A7 placed in supplementary? This seems an important and interesting finding – particularly in light of the comments above about possible immunometabolic influences of the BWB specific loss-of-function mutation in the myostatin (MSTN) gene.\nResponse: We thank the reviewer for drawing attention to the association between IL-1B response and SLC16A7, and for emphasizing its potential relevance within the broader context of immunometabolic regulation and the known MSTN loss-of-function mutation in Belgian White Blue (BWB) cattle.\nWe fully agree that this represents an intriguing and biologically meaningful finding. Our initial decision to report it in the Supplementary Information was primarily due to space limitations and the absence of strong colocalization evidence (PH4<0.8) between the protein QTL (pQTL) for SLC16A7 and its corresponding cis-eQTL. This may suggest a context-specific regulatory effect on SLC16A7 expression. Moreover, the identification of this locus was not replicated when we used T2 data (Supplementary table 5). Therefore, due to the lack of strong colocalization evidence, the absence of replication at T2 and the space limitations of the Figure, we present this result in the supplementary figure.\n17. ‘In addition to significant associations (p<1.3e-07), we identified 93 genome-wide suggestive signals at p<2. 417 6e-06 (Supplementary Fig. 10 and Supplementary Table 8) 418 that warrant further analysis’. This section is confusing. What do the authors mean by ‘suggestive signals’? If they did not meet p value cut-offs, maybe less could be written about them as this section is already quite long.\nResponse: We agree with the reviewer. We have therefore replaced this section with a shortened version that mentions these results without going into specific details (Lines 388-389).\nMinor Issues\nThere are many grammatical and linguistic errors throughout the manuscript - some examples are presented below:\nL63 ‘limited to decipher’…\nL67’ also directly in humans ‘ … The word ‘directly is redundant here\nL77’ rewrite: ‘prospect to refine our understanding of …’\nL85 ‘they remain notably absent in animals’ should be ‘remain notably absent from animals’\nL90 ‘This becomes…’ what becomes?\nL298 ‘genetic variations that impacts’…\nL 605 What is ‘livibility’?\nL 614-616 needs a reference.\nResponse: We thank the reviewer for careful reading of the manuscript and for highlighting several grammatical and linguistic issues.\nThe specific issues noted by the reviewer have been addressed as follows:\nLine 63: ‘limited to decipher’ has been replaced by “less effective at elucidating” (Lines 57-58).\nLine 67: The sentence has been removed as likely non necessary (Line 61).\nLine 77: The whole sentence has been rewritten and is now “Systems immunology thus holds the potential to deepen our understanding of the mechanisms underlying disease and to provide novel insights into host–pathogen co-evolution.” (Lines 68-70).\nLine 85: Revised to “remain underutilized in livestock species.” (Line 79).\nLine 90: The sentence has been rephrased to increase clarity and is now : “The scarcity of systems immunology studies in animals becomes even more perplexing when considering the notable advancements witnessed in the past decade in animal genetics.” (Lines 87-89).\nLine 298: Corrected to “genetic variations that impact…”(Line 278).\nLine 605: The typographical error “livibility” has been corrected to “livability”. This term is the one used in the referenced publication mentioning it. We have given a definition of it to clarify it (Lines 319-321).\nLines 614–616: Supporting references have been added and the sentence has been rephrased (Lines 573-576).\nIn addition, we have performed a final language and style revision to ensure the manuscript reads clearly and concisely throughout."}, {"round": 2, "reviewer_comments": {"reviewer_1": "The response to reviewers' comments is very impressive and this manuscript is enormously improved as a result. It is now ready for publication", "reviewer_2": "I thank the authors for their thorough responses, and the effort they have made to address the concerns raised which have improved the manuscript overall.\nI still disagree with some of their opinions, in particular related to the challenging comparisons of PBMC v WWBC. For example they highlight a strength of including stimulated PBMCs in that it may “capture more nuanced immune variability” by showing greater variability than in WWBC, but there is no evidence shown that this higher variabilty is biological and not techncial. In fact the new analyses presented in Fig S7 indeed suggests that these specific cytokines in PBMCs do show higher technical variability than in WWBC.\nThe measurement of cytokines in FCS is interesting as it demonstrates variable levels of IL10 at phsyiological levels, it would be helpful to note this in the manuscript so other teams are aware of this potential confounding factor.\nHowever overall I commend and congratulate the authors on their work.", "reviewer_3": "In my view the manuscript has improved considerably by focusing primarily on the statistical significant results and providing additional context and a more honest picture of the limitations of the study. I'd like to thank the authors for addressing my concerns and providing convincing evidence for the rigour of the analyses and validity of the results. I fully agree, even if the sample size was somewhat limited, the novelty of the research approach in livestock and the findings are of great interest in the animal science community\nI only have a few minor remaining comments:\nL292-294: Could you please provide info how you ‘adjusted those immunophenotypes on the 10 leading associated variants”. And how the lack of identifying additional significant loci indicates that “single major QTLs affect these immunophenotypes”. Do you mean those 10 QTL?\nL427: suggest to change “significance of the estimates of genetics explained variance (h2_SNP)” to “significance of the estimated contribution of genetics (h2_SNP) to the observed total variance”.\nIn the methods section (L989-990) the authors state “The proportion of variance in phenotype explained by a given SNP (PVE) was calculated as described previously.90” Please check if this is the right reference as it does not seem to provide this information.\nPlease check consistency in spelling of interindividual / inter-individual throughout\nL454-455: The section header seems a bit of an overstatement given the mean correlation between predicted and measured cytokine production was <0.1, and the highest one was 0.4. Maybe replace with ‘partly predicted’?\nL475: ‘using a described methodology’ with a reference is a bit unsatisfying for the reader. Could you provide some more detailed information about the methodology, even if only referring to the methods section? I did not see the info in the methods section.\nL513: Is this not a bit of an over-statement given the at best moderate correlations between observed and predicted measures? Suggest to tone the statement down.\nL484. Add full-stop", "reviewer_4": "The authors have thoroughly revised the manuscript and have nicely answered to most questions. Congratulations.\nMain comment\n1) It is still mandatory to better acknowledge previous published work on the influence of host genetics on immune response variability in farmed animals. It is true that the number of scientific reports relying on systems immunology approach is limited in livestock but it is not accurate to assimilate this to a lack of results or even interest on links between genetics and immunity traits in the same animal species. There have been significant publications in this field for long time (list of “old” references on pigs below). And of course, it is highly expected that implementing systems immunology approaches more widely in farmed animals will contribute to better take into account the role of the genetics to improve animal health. This must be clearly stated in the abstract and in the text.\n- Abstract, lines 32-33: replace “However, recent human studies show that genetics strongly influence immune responses” by “However, studies on humans and livestock also show that genetics strongly influences immune responses, …”\n- Introduction line 46: not only “In humans”, add references for livestock species also\n- Discussion line 508: “echoing patterns seen in humans”. This also echoes patterns seen in pigs for instance. Other references must be added such as Flori et al., 2011, (ref #28)\nFor the record, examples of pioneer work in the pig:\n- Jensen, P. T., and K. Christensen. 1980. Genetic studies on the in vitro PHA transformation of porcine lymphocytes. Vet. Immunol. Immunopathol. 2:133\n- Rothschild, M. F., H. L. Chen, L. L. Christian, W. R. Lie, L. Venier, M. Cooper, C. Briggs, and C. M. Warner. 1984. Breed and swine lymphocyte antigen haplotype differences in agglutination titers following vaccination with B. bronchiseptica. J. Anim. Sci. 59:643.\n- Rothschild, M. F., H. T. Hill, L. L. Christian, and C. M. Warner. 1984. Genetic differences in serum-neutralization titers of pigs after vaccination with pseudorabies modified live-virus vaccine. Am. J. Vet. Res. 45:1216.\n- Edfors-Lilja, I., B. Gahne, and H. Petersson. 1985. Genetic influence on antibody response to two Escherichia coli antigens in pigs. II. Difference in response between paternal half-sibs. Z. Tierz. Zu¨chtungsbiol. 102:308.\n- Edfors-Lilja et al., 1998 https://pubmed.ncbi.nlm.nih.gov/9670960/\n- Mallard, B. A., B. N. Wilkie, B. W. Kennedy, and M. Quinton. 1992. Use of estimated breeding values in a selection index to breed Yorkshire pigs for high and low immune and innate resistance factors. Anim. Biotechnol. 3:257.\nIn cattle, the authors could refer to Martin et al. 2022. Identification of genome regions and promising candidate genes linked to innate immune capacity on young Holstein calves. Proceedings of 12th World Congress on Genetics Applied to Livestock Production (WCGALP), R.F. Veerkamp and Y. de Haas (eds) DOI: 10.3920/978-90-8686-940-4 _94. This reference could also be used in the discussion regarding identified genome regions.\nOther comments:\n2) mentioning antibiotics as synthetic intrants is questionable because antibiotics are delivered upon prescription and not as permanent intrants (at least they should not). Even if it is important to use them with parsimony, their status is actually different from the status of synthetic fertilizers for plants.\n3) Line 512: replace “can” by “could” (important to be cautious because the prediction capacity in other populations or breeds in other environments may not be guaranteed)\n4) Be careful with abbreviations. LD for Linkage Disequilibrium (line 286) is not defined for example and LD stands for Lethal Dose in the journal abbreviations list.", "reviewer_5": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts."}, "author_response": "Response: We sincerely thank the reviewer for the thoughtful and constructive feedback, and we greatly appreciate the recognition of our efforts to improve the manuscript.\nRegarding the reviewer’s comments about the comparison between PBMCs and WWBCs, we agree that distinguishing biological from technical variability is challenging. WWBCs are certainly the cellular preparation that most closely reflects the natural conditions found in blood. It is also clear that the manipulations required to isolate PBMCs are more likely to introduce technical variation. Nevertheless, we consider it appropriate to present both types of cellular preparations, since PBMCs, being depleted of neutrophils, for instance, allow for a clearer detection of genuine immune variations within other cell populations such as monocytes and lymphocytes. This approach has also been employed in previous studies (Li et al., Cell 2016; Boahen et al., Am J Hum Genet 2022; ter Horst et al., Cell 2016, etc). In order to address the reviewer’s comment, we have further emphasized the relevance of WWBCs in our discussion (Lines 587-595).\nAs suggested, we have also added a note in the manuscript acknowledging that variable cytokine levels in FCS may represent a potential confounding factor and should be considered in future analyses (Lines 595-597).\nResponse: We thank the reviewer for the positive and encouraging comments. We are pleased that the revisions have clarified the scope and limitations of the study and strengthened the presentation of our key findings.\nI only have a few minor remaining comments:\nL292-294: Could you please provide info how you ‘adjusted those immunophenotypes on the 10 leading associated variants”. And how the lack of identifying additional significant loci indicates that “single major QTLs affect these immunophenotypes”. Do you mean those 10 QTL?\nResponse: Among the measured immunophenotype, we performed GWAS for each immunophenotype and 10 QTLs were found for the corresponding immunophenotype. In order to test whether additional genetic variants independently influence the phenotype after accounting for the effect of the major (lead) QTL, we included the lead SNP as a covariate in the model. This procedure (condition analysis) allowed to reveal whether other nearby or distant variants contribute to the trait beyond the primary association signal. Here, after conducting condition analysis for those 10 immunophenotypes, no other significant loci were identified, indicating that these 10 major QTLs exert strong genetic effects on the specific immunophenotype.\nWe have slightly modified the sentence in the manuscript to improve clarity (Lines 284-285).\nL427: suggest to change “significance of the estimates of genetics explained variance (h2_SNP)” to “significance of the estimated contribution of genetics (h2_SNP) to the observed total variance”.\nResponse: We have revised the sentence accordingly (Lines 415-416).\nIn the methods section (L989-990) the authors state “The proportion of variance in phenotype explained by a given SNP (PVE) was calculated as described previously.90” Please check if this is the right reference as it does not seem to provide this information.\nResponse: The formula and method for calculating PVE are provided in the supplementary information in that reference.\nPlease check consistency in spelling of interindividual / inter-individual throughout\nResponse: We have harmonized the spelling to use inter-individual consistently throughout the manuscript.\nL454-455: The section header seems a bit of an overstatement given the mean correlation between predicted and measured cytokine production was <0.1, and the highest one was 0.4. Maybe replace with ‘partly predicted’?\nResponse: We agree with the reviewer’s observation and have changed the section title to “Partly predicted cytokine responses” to better reflect the predictive accuracy.\nL475: ‘using a described methodology’ with a reference is a bit unsatisfying for the reader. Could you provide some more detailed information about the methodology, even if only referring to the methods section? I did not see the info in the methods section.\nResponse: The detailed description of the methodology is provided in the Methods with a specific section (Lines 905-924).\nL513: Is this not a bit of an over-statement given the at best moderate correlations between observed and predicted measures? Suggest to tone the statement down.\nResponse: We have softened the statement to more accurately reflect the moderate correlations observed between the predicted and measured values (Lines 500-501).\nL484. Add full-stop\nResponse: We think that the full-stop is there. Sorry if we have not seen where the problem is.\nResponse: We sincerely thank the reviewer for the positive feedback and kind words. We are pleased that the revisions have satisfactorily addressed the previous comments and improved the overall quality of the manuscript.\nMain comment\n1) It is still mandatory to better acknowledge previous published work on the influence of host genetics on immune response variability in farmed animals. It is true that the number of scientific reports relying on systems immunology approach is limited in livestock but it is not accurate to assimilate this to a lack of results or even interest on links between genetics and immunity traits in the same animal species. There have been significant publications in this field for long time (list of “old” references on pigs below). And of course, it is highly expected that implementing systems immunology approaches more widely in farmed animals will contribute to better take into account the role of the genetics to improve animal health. This must be clearly stated in the abstract and in the text.\n- Abstract, lines 32-33: replace “However, recent human studies show that genetics strongly influence immune responses” by “However, studies on humans and livestock also show that genetics strongly influences immune responses, …”\n- Introduction line 46: not only “In humans”, add references for livestock species also\n- Discussion line 508: “echoing patterns seen in humans”. This also echoes patterns seen in pigs for instance. Other references must be added such as Flori et al., 2011, (ref #28)\nResponse: We thank the reviewer for this important comment and fully agree that some important work had been previously done on the genetic basis of immune response variability in livestock. We intended to indicate that, to date, no large systems immunology studies have been conducted in livestock species. However, we acknowledge that this statement may have inadvertently downplayed the valuable work that has already been carried out in this field. We have therefore done our best to revise the manuscript accordingly in response to the reviewer’s comment.\nAbstract (lines 32–33): We agree with the comment. However, based on the suggestion from editor, we take the edited abstract from the editor to comply fully with the journal’s requirements.\nIntroduction (line 46): We have included some additional relevant studies in livestock species as suggested. They are provided as references 28-35.\nDiscussion (line 508): We have added the reference of Flori et al., 2011. (now line 500)\nFor the record, examples of pioneer work in the pig:\n- Jensen, P. T., and K. Christensen. 1980. Genetic studies on the in vitro PHA transformation of porcine lymphocytes. Vet. Immunol. Immunopathol. 2:133\n- Rothschild, M. F., H. L. Chen, L. L. Christian, W. R. Lie, L. Venier, M. Cooper, C. Briggs, and C. M. Warner. 1984. Breed and swine lymphocyte antigen haplotype differences in agglutination titers following vaccination with B. bronchiseptica. J. Anim. Sci. 59:643.\n- Rothschild, M. F., H. T. Hill, L. L. Christian, and C. M. Warner. 1984. Genetic differences in serum-neutralization titers of pigs after vaccination with pseudorabies modified live-virus vaccine. Am. J. Vet. Res. 45:1216.\n- Edfors-Lilja, I., B. Gahne, and H. Petersson. 1985. Genetic influence on antibody response to two Escherichia coli antigens in pigs. II. Difference in response between paternal half-sibs. Z. Tierz. Zu¨chtungsbiol. 102:308.\n- Edfors-Lilja et al., 1998 https://pubmed.ncbi.nlm.nih.gov/9670960/\n- Mallard, B. A., B. N. Wilkie, B. W. Kennedy, and M. Quinton. 1992. Use of estimated breeding values in a selection index to breed Yorkshire pigs for high and low immune and innate resistance factors. Anim. Biotechnol. 3:257.\nIn cattle, the authors could refer to Martin et al. 2022. Identification of genome regions and promising candidate genes linked to innate immune capacity on young Holstein calves. Proceedings of 12th World Congress on Genetics Applied to Livestock Production (WCGALP), R.F. Veerkamp and Y. de Haas (eds) DOI: 10.3920/978-90-8686-940-4 _94. This reference could also be used in the discussion regarding identified genome regions.\nResponse: We thank the reviewer for this suggestion and for drawing our attention to the proceeding by Martin et al. (2022), which identified genomic regions and candidate genes linked to innate immune capacity in Holstein calves.\nAccordingly, we have now cited Martin et al. (2022) in both the Introduction, to highlight previous research in livestock species, and in the Discussion, when discussing the validation in other breeds.\nOther comments:\n2) mentioning antibiotics as synthetic intrants is questionable because antibiotics are delivered upon prescription and not as permanent intrants (at least they should not). Even if it is important to use them with parsimony, their status is actually different from the status of synthetic fertilizers for plants.\nResponse: We thank the reviewer for this comment. We agree that the usage of antibiotics is increasingly regulated in high-income countries. However, there are still in use in some countries, particularly where regulatory frameworks are less strict, or enforcement is limited. For example, this is discussed in the following reference and in some publications of FAO.\nMatheou, A. et al. Antibiotic Use in Livestock Farming: A Driver of Multidrug Resistance? Microorganisms 13, 1–19 (2025).\n3) Line 512: replace “can” by “could” (important to be cautious because the prediction capacity in other populations or breeds in other environments may not be guaranteed)\nResponse: We agree with the reviewer and have replaced “can” by “could” (line 500).\n4) Be careful with abbreviations. LD for Linkage Disequilibrium (line 286) is not defined for example and LD stands for Lethal Dose in the journal abbreviations list.\nResponse: We thank the reviewer for pointing out the issue with abbreviations. We have now defined LD (Linkage Disequilibrium) at its first occurrence (line 279) and verified all abbreviations throughout the manuscript to ensure consistency and clarity.\nResponse: We thank the reviewer for sharing this information and for the contribution to evaluating our manuscript."}]

Decision letter