Genome-wide scan for commons SNPs affecting bovine leukemia virus infection level in dairy cattle

Background Bovine leukemia virus (BLV) infection is omnipresent in dairy herds causing direct economic losses due to trade restrictions and lymphosarcoma-related deaths. Milk production drops and increase in the culling rate are also relevant and usually neglected. The BLV provirus persists throughout a lifetime and an inter-individual variation is observed in the level of infection (LI) in vivo. High LI is strongly correlated with disease progression and BLV transmission among herd mates. In a context of high prevalence, classical control strategies are economically prohibitive. Alternatively, host genomics studies aiming to dissect loci associated with LI are potentially useful tools for genetic selection programs tending to abrogate the viral spreading. The LI was measured through the proviral load (PVL) set–point and white blood cells (WBC) counts. The goals of this work were to gain insight into the contribution of SNPs (bovine 50KSNP panel) on LI variability and to identify genomics regions underlying this trait. Results We quantified anti–p24 response and total leukocytes count in peripheral blood from 1800 cows and used these to select 800 individuals with extreme phenotypes in WBCs and PVL. Two case-control genomic association studies using linear mixed models (LMMs) considering population stratification were performed. The proportion of the variance captured by all QC-passed SNPs represented 0.63 (SE ± 0.14) of the phenotypic variance for PVL and 0.56 (SE ± 0.15) for WBCs. Overall, significant associations (Bonferroni’s corrected -log10p > 5.94) were shared for both phenotypes by 24 SNPs within the Bovine MHC. Founder haplotypes were used to measure the linkage disequilibrium (LD) extent (r2 = 0.22 ± 0.27 at inter-SNP distance of 25−50 kb). The SNPs and LD blocks indicated genes potentially associated with LI in infected cows: i.e. relevant immune response related genes (DQA1, DRB3, BOLA-A, LTA, LTB, TNF, IER3, GRP111, CRISP1), several genes involved in cell cytoskeletal reorganization (CD2AP, PKHD1, FLOT1, TUBB5) and modelling of the extracellular matrix (TRAM2, TNXB). Host transcription factors (TFs) were also highlighted (TFAP2D; ABT1, GCM1, PRRC2A). Conclusions Data obtained represent a step forward to understand the biology of BLV–bovine interaction, and provide genetic information potentially applicable to selective breeding programs. Electronic supplementary material The online version of this article (10.1186/s12864-018-4523-2) contains supplementary material, which is available to authorized users.


Background
The Bovine leukemia virus (BLV) is an oncogenic δretrovirus that persistently infects cattle [1]. BLV infection elapses without noticeable symptoms in most of the animals while 30-40% of them evolve into a sub-clinical persistent lymphocytosis (PL) characterized by an expansion of proviral infected B-lymphocytes in peripheral blood increasing white blood cells (WBC) counts (leukocytosis). Less than 10% of infected animals will develop the severe form of the disease leading to lethal lymphomas or lymphosarcoma in various organs [2]. The BLV is transmitted among cows through direct exposure to blood, body fluids or consumption of contaminated colostrum/milk with proviral DNA [3,4]. In the absence of effective measures controlling BLV transmission, the virus is massively propagated reaching high levels of endemicity. BLV is world-wide distributed. In North America BLV's prevalence reached 70-85% at farm level [5][6][7], and in Argentina more than 80% of the lactating cows in the main dairy production area are infected [8]. International trade restrictions of livestock products from affected herds [9] impact on the regional economy, but the total economic loss is significantly higher if we consider milk production dropping and faster culling of asymptomatic BLV carriers compared to BLV free herds [10,11]. These conditions are emphasised in lymphocytotics animals and also in those with high anti-BLV humoral response [12][13][14][15][16]. Traditional control strategy implies segregating/eliminating infected cows after BLV testing, but this is not economically affordable when the regional prevalence is high [17]. A suitable long term approach to contain endemic livestock diseases, especially those where the therapeutics or preventives treatments failed, is DNA-based breeding for host resistance to infection [18]. Mastitis incidence is currently being incorporated in genomic evaluations for dairy cattle ( [19,20]).
A key component to set up host genetic studies is properly defining and recording a consistent phenotype [21,22]. The high PVL has been associated with progression of the BLV pathogenesis [23,24], and with enhanced viral dissemination [25]. Previous evidence pointed to DNA polymorphisms in candidates genes (in particular, alleles from the bovine MHC class II DRB gene) as contributors to the inter-individual variation in the PVL and PL [26][27][28][29][30][31][32][33][34]. However, genetic variation for resistance/susceptibility to diseases in cattle is usually polygenic, i.e. susceptibility to paratuberculosis [35,36], mastitis [37,38] bovine respiratory disease complex (BRDC) [39] and recently observed in a genome-wide mapping for genetics determinants of BLV incidence [40,41].
With the final objective in mind of developing (implementing) genomics selection programs to control BLV dissemination it is essential to understand the biological mechanism implicated. In an attempt to address this questions we utilized antibodies anti-p24 as cost-effective predictors of the level of PVL in blood [8] on over 1800 animals, characterized with WBC counts, in order to recruit a population of Holstein and Holstein x Jersey crossbred dairy cows representing extreme phenotypes of level of infection (LI) in vivo, and therefore, with different individual risk for viral transmission. Afterwards, we integrate these phenotypes in a case-control genome-wide association study (GWAS) using the Illumina BovineSNP50v2 BeadChip to investigate genomic regions associated with BLV infection.

Study population
Blood samples were obtained from 1800 Holstein and Holstein x Jersey cows representing 25 half-sisters families from 16 dairy herds located in the central region of Argentina. Farm owners' consent was obtained before animal sampling. All animals have full pedigree and were over 3 years old at the time of sampling and between 1 st and 9 th lactation (78.7% concentrated between 2 nd and 4 th lactation). The procedures followed for extraction and handling of samples were approved by the Institutional Committee for Care and Use of Experimental Animals of the National Institute of Agricultural Technology (CICUAE-INTA) under protocol number 35/2010, and followed the guidelines described in the institutional Manual. Blood samples were collected in heparinized tubes. Plasma and fresh blood were frozen at − 20°C until used.
Quality and concentration of genomic DNA extracted from whole blood samples using a commercial kit (Blood Genomic DNA AxyPrep™, Axygen Biosciences, Union City, USA) was assessed using a micro−volume spectrophotometer (NanoDrop™ Technologies, Inc. Wilmington, USA). DNA samples were normalized to a minimum concentration of 20 ng/μl.

Phenotyping
The anti-p24 ELISA assay was used to discriminate BLV infected animals [42], followed by total WBC counts. In a previous study we showed that the individual anti-BLV humoral response is a cost effective indirect measure of infection in vivo [8]. Therefore, animals were first screened by anti-p24 ELISA and then divided into two groups of contrasting phenotypes on the basis of their antibodies level. The LI in selected animals was then corroborated by determining blood PVL set-point.
Serology: The recombinant viral core-p24 protein was used to detect plasma anti-BLV antibodies for each sample as described in [42] Briefly, a ratio sample to positive (S/P), called % of reactivity (% R) is obtained: S/P = [(OD Sample -OD WS-)/(OD WS + 1/7 -OD WS-)] × 100, where OD = optical density, WS-= negative serum, WS + 1/7 = serum international standard weak positive control (diluted 1:7 in negative serum) used as reference value. The % R of all samples tested was indicated as "low" (LR), for records between a cut-off value (25%) [42] and the weak positive serum value (100%); and as "high" (HR) in those samples with measurements higher than 200%.
PVL relative quantification: PVL determination was performed over LR and HR animals by real-time PCR using TaqMan technology as described in [42,43]. A fragment of BLV pol gene (target) was amplified simultaneously with the reference gene 18S. Fetal Lamb Kidney cells (FLK) containing 4 copies/cell of BLV proviral DNA, in a final concentration of 1% in BLV-free peripheral blood mononuclear cells (PBMCs) were used as internal control. This proportion of BLV-infected cells corresponds to a low level of natural infection [42], proved to be 5% in aleukemics animals [44]. The relative PVL was defined according to the equation described by Pfaffl et al. [45]: where E = reaction efficiency; Ct = threshold cycle. The control PVL was indicated as 1 and all samples tested were referred to this. The detection limit was 1 infected cell /2000 uninfected cells. The PVL was qualified as "No Detectable" (ND) if a Ct value could not be obtained (below the detection limit), as "low" if R < 1 (LPVL) and as "high" if R > 1 (HPVL). WBC counts: Total leukocyte count was performed manually using a haemocytometer. Fresh blood was diluted 1:20 in a 1% Gentian Violet and 2% acetic acid solution. Normal WBC counts in BLV-seronegative adult cattle are in the range of 2100 to 11,500 cells/μl (mean: 6800 ± 2 SD) [46].  [47]. SNPs were mapped to the bovine assembly UMD 3.1 [48]. The SNPs identification was reassigned to its corresponding "rs" catalogued in dbSNP [49] using the SNPchiMp v.3 tool [50].
Quality control (QC) of the genotypic data was carried out using a set of tools and routines provided by PLINK v1.07 [51]. Individuals with a genotype call rate (CR IND ) < 90% were excluded from further analysis. Mendelian inheritance discordances were checked and set as missing. SNPs with a call rate (CR SNPs ) < 90%, a deviation from HWE p < 1.10 − 08 and Minor Allele Frequency (MAF) < 0.01 were removed from the study. After the QC, the data set was composed of 44,174 SNPs and 925 individuals (898 cows and 27 bulls).

Case-control definitions
The case-control phenotypes to conduct GWAS for LI in vivo considered 396 cases (HPVL) and 373 controls (ND/LPVL). The alternate phenotype was defined by the number of circulating white blood cells (WBCs), which is indicative of the presence of PL, a sub-clinical symptom of BLV infection [52], clustering 283 cases (WBC ≥ 11,500 cells/μl) and 611 control cows (WBC < 11,500 cells/μl).

Linkage disequilibrium
We performed chromosomal haplotypic phase assignment using the software DualPHASEv1.1 included in the PHA-SEBOOK 2.3 package [53] using the default settings. To calculate the extent of linkage disequilibrium (LD) only maternal founders chromosomes were used. Phasing algorithms allowed the imputation of not assigned genotypes [54] so that the missing genotypes or those set as missing due to Mendelian inconsistencies were imputed. Measure of LD extent was addressed by calculating the correlation of alleles in two loci (r 2 and D') for each chromosome with Haploview v4.2 [55]. The calculations were performed on all pairs of syntenic SNPs with an inter-marker distance < 1 Mb. The r 2 measure was then averaged at intervals of 10 kb and plotted versus the inter-SNPs distance.

Population stratification
We inspected the genealogy of bulls in the sampled families (http://www.sra.org.ar/rrgg/) evidencing a strong co-ancestry among them (Additional file 1: Figure S1). In order to prevent spurious associations due to population substructure we performed a multi-dimensional scaling (MDS) analysis on identity by state (IBS) pairwise based distances matrix as implemented in PLINKv1.07 [51]: 1 − (IBS2 + 0.5IBS1)/N, where IBS2 is the number of SNPs in which two individuals share two IBS alleles; IBS1 is the number of SNPs in which only one allele is shared IBS and N is the number of markers used. Finally, the distance matrix was partitioned into a number of uncorrelated variables (or Principal Components, PC). The first two PCs were used to visualize differences in ancestry among individuals (Additional file 2: Figure S2).

Association analysis
The GWAS were performed using LMMs as implemented in GCTA v1.24 [56]. The general model in the observed scale can be written as: where y is a vector n × 1 that contains binaries values for PVL or WBC; b is a vector that contains the fixed effects of the SNP i to be tested for association (considering an additive genetic effect) and the co-variables Age (A j ) (j = 3 to 12); Herd (H k ) (k = 1 to 16); Lactation Number (L l ) (l = 1 to 9); Percentage of Holstein breed (PH m ) (m = 0, 0.125, 0.25, 0.50, 0.75 and 1); Bull (B n ) (n = 1 to 25) and PC derived from the co-ancestry analysis (PC1 -PC8). The random effects included the additive polygenic effect u of the animal; and the vector e of residuals effects, X and Z stand for the incidence matrix that relate the values in y to the fixed effects in b and the random effects in u, respectively. A liability threshold model l = Xb + g + e was utilized in order to relate observations on the observed scale to liabilities on the unobserved continuous scale. Liability phenotypes are distributed as N(0,1), On the liability scale, the random effects included the additive polygenic effect g of the animal, distributed N(0, σ g 2G), where G represents the genetic relationship matrix (GRM) calculated from the proportion of shared alleles among individuals for all QC-passed SNPs; and the vector e of residuals effects, distributed~N(0, σ e 2I )), where I is the identity matrix. Other terms are the same as in the linear model but on the liability scale [57,58]. The heritability on the scale observed Þ was transformed to that on the liability scale [58]: l is the heritability on the underlying scale; h 2 o is the heritability on the observed scale; K is the disease prevalence and zis the height of the density function of the standard normal probability at the threshold value (K).
This model was utilized to calculate the genetic correlation between PVL and WBC. According to Lee et al. [59] the genetic correlation was defined as covðg 1 ; g 2 Þ1 z 1 where, gis the genetic value on the observed scale, g * is the genetic value on the underlying scale, z is the height of the density function of the standard normal probability at the threshold value (K), P is the proportion of cases in the sample and K is the prevalence for features 1 (PVL) and 2 (WBC), respectively.
The association analysis implemented in GCTA v1.24 uses the restricted maximum likelihood method (REML) to solve the equations of LMMs [56,57].
The genetic variance for PVL and WBCs traits was calculated on the 29 individual autosomes and on the X bovine sexual chromosome. The genetic variance was decomposed into separate chromosomal segments following: where, g i is a random genetic effects vector corresponding to the complete genome or a particular chromosome. The phenotypic variance was partitioned so that varðyÞ i is the variance of the i-th genetic factor with its related GRM.
The family-wise error rate was controlled by Bonferroni's correction α′ ≈ α/n, where α is the significance level of 5% and n is the number of SNPs tested. The α′value represents the threshold that the p-value must reach for each test to be considered significant.

Gene content of the genomic regions identified in the association studies
Gene annotation for BTA 23 (UMD3.1 bovine assembly, [48]) was retrieved using BioMart tools [60]. Firstly, SNPs were assigned to a particular gene if they were located on exons, introns or untranslated regions (UTRs). Secondly, SNPs lying within a 20 kb window upstream and downstream from gene boundaries, aiming to capture promoter/ expression regulatory regions [61]. Finally, remaining SNPs were assigned to the closest gene from 20 to 30 kb, within the range of useful average LD values found (r 2 > 0.2) [62][63][64][65]. Gene to SNPs mapping was performed using BedTools package [66]. Gene functional classification based in Gene Ontology terms (GO) was performed using the PANTHER classification system website (http://www.pantherdb.org/) [67].

Results
Genome-wide distributed SNPs were tested for association with both PVL and WBCs traits using logistic regressions assuming an additive model that incorporated as covariates 8 PCs and other confounding factors Age (A); Herd (H); Lactation number (L); Percentage of Holstein (PH) and Bull (B). The observed statistical values of association were ordered from lowest to highest and plotted against their expected values and genomic inflation factor (λ) was calculated using PLINK. Both for PVL (Additional file 3: Figure S3, λ = 1.21) and WBCs (Additional file 4: Figure S4, λ = 1.30), the data was best fitted with models considering PCs1-8 + A_H_L_PH_B. Given that kinship and/or cryptic relationship might be contributing sources of distortion in the case-control associations we decided to use LMMs as implemented in GCTA v1.24 [56] that incorporates each of QC-passed SNP as a fixed effects and the polygenic animal effect as random. QQ plots showed a proper fit for stratification and a pronounced deviation to the right, which is suggestive of true associations (Fig. 1).
The Manhattan plots revealed significantly associated SNPs with HPVL in peripheral blood (Fig. 2) and with a high number of circulating WBCs (Fig. 3).
All SNPs exceeding the significance threshold after Bonferroni's correction (−log 10 p > 5.94) are located on BTA 23. Notably, SNPs mapping on BTA 23 largely captured most of the genetic variation for both traits, influencing BLV levels of infection. Overall, a group of 24 SNPs were significantly associated with high PVL, whilst 11 SNPs were associated to high WBCs counts. Of these, only 1 SNP (rs110277740) is not shared with those associated with PVL (Table 1)   is explained by SNPs located outside this MHC region, so we decided to repeat the GWAS for PVL but excluding the SNPs located between~20-36 Mb of BTA 23, which could be masking the SNPs effects in non-MHC regions. As result, we did not find any significantly associated SNP (Additional file 5: Table S1). Unless a significant SNP is itself the causal variant, it will capture the effects of other variants in LD responsible for the phenotype. The clustering of genes with related functions in the linkage disequilibrium pattern in the bovine MHC region makes it difficult to detect effects of individual loci. We characterized the LD extent (Additional file 6: Figure S5) in our Holstein and Holstein x Jersey population. At distances ≤30 kb, the average correlation value r 2 was 0.26 ± 0.30 (34% of pairwise SNPs showed r 2 > 0.25). This value rapidly decays to 0.14 ± 0.21 when the inter-markers distance was between 50−100 kb. About 60% of the inter-SNPs distances in this work is in the range of 25-50 kb with an average r 2 = 0.22 ± 0.27, indicating that causal genetics variants effects could be retrieved from SNPs hits.
We found five significantly associated SNPs pairs having a strong allelic correlation (r 2 > 0.  SNPs were assigned to putative biological function if they fell within a coding, intronic, or regulatory region. For those SNPs mapping outside genic regions we scrutinised the closest gene located at a distance where a moderate LD was kept. A total of 44 genes were identified using the criteria described above. Table 2 Table 2). The two remaining SNPs, rs41587216 and rs110473048, were at 39.3 kb of CRISP1 (ENSBTA G00000008085) and 21.7 kb GCM-1 (ENSBTAG0000000 8148), respectively.
The PVL phenotype association study included all significant SNPs for the GWAS considering WBCs -except the SNP assigned to Activator of Transcription Basal 1 gene (ABT1, rs110277740) -and other 14 significant SNPs. The functional annotation resulted in a synonymous substitution in TNXB (ENSBTAG00000001444) and a SNP located in the 3'-UTR region of SCL44A4 (ENSBT AG00000005675). Also, there are two SNPs located in intronic regions of NEU1 and TRAM2 (ENSBTAG000 00009112 and ENSBTAG00000005674, respectively). The   Odorant receptors (Olfactory receptor). Finally, rs11 0579760, rs110794231 and rs41587536 polymorphisms are at 30.7, 20.3 and 22.3 kb of BoLA-DRB3 (ENSBTAG00000013919) BoLA-A (ENSBTAG00000 022590) and SFTA2 (ENSBTAG00000038810), respectively. Genes represented by the SNPs were assigned to Gene Ontology (GO) functional terms and the PANTHER (Protein Class, PC) (http://pantherdb.org/) classification system. The analysis was performed for all genes represented by significantly associated SNPs with any of the two traits. A total of 11 PC "parents" terms were identified being the three most common categories "protein defense/immunity", (PC00090) (21,7%), "signalling molecule" (PC00207) (21,7%) and cytoskeletal protein (PC00085) (8,7%). In Additional file 7: Table S2 a description of each PC term is displayed along with their "child" terms providing specificity to the functional description. Notably, a single gene might be contained in more than one PC class. Afterwards, we classified genes under the category GO "Biological Processes" (BP) wherein each gene identification is assigned according to their function within a network of proteins that collectively carry out a particular process within the cell. Figure 5 presents the distribution of 12 "parent" terms found according to GO BP classification, these represent general cellular processes being the top five terms "immune system process" (16%), "response to stimulus" (16%), "metabolic process" (14%), "cellular process" (12%) and "localization" (11%). The search for GO terms and/or KEGG pathways overrepresented on the gene set identified vs those present in whole bovine gene set using FatiGO method [68] did not produce significant results (data not shown).

Discussion
The humoral immune response plays a crucial role in the control of viral infections in cattle and humans, and the level of response varies from individual to individual, being this variability highly heritable [69,70]. In cattle, we showed that the humoral response to BLV infections indeed varies among animals and, interestingly, it can be used to infer the circulating PVL in an infected animal [8]. In this study, we use the % of Reactivity (anti-p24) to screen a population of BLV naturally infected animals aiming to select cows most likely to have extreme values of PVL set-point, which was subsequently confirmed by qPCR. The predictive ability of serological markers to detect high PVL conditions was good being the ROC area under the curve > 0.80 (data not shown). We then performed a case-control study setting levels of infection (LI) as phenotype aiming to estimate (a) the proportion of phenotypic variation that can be attributed to genetic components and (b) the cattle genome mapping positions of genetic polymorphisms underlying this variation. A positive genetic correlation between PVL and WBCs of magnitude~0.81 (SE ± 0.10) supports the idea of shared or genomic regions influencing both traits. We used kinship captured from QC pruning variants to assess the additive genetic variance of LI produced by BLV in the population. The percentage of variation reached was 63% (h 2 GRM_PVL = 0.63) and 56% (h 2 GRM_WBC = 0.56) on the risk scale or liability. On this scale, the phenotypic variances do not depend on the disease prevalence, so h 2 lia values allow comparisons between populations and environments where incidence varies. When partitioned, Fig. 5 Distribution of the GO "Biological Process" (BP) categories for genes represented by the significantly associated SNPs in GWAS for PVL and WBC counts the proportion of the phenotypic variance explained by markers for both phenotypes was distributed on different chromosomes. In a pedigree-based study by Abdalla et al. [71], the contribution of genetic components to the incidence of BLV in Holsteins and Jersey herds (approximately 30% of herd level prevalence) was estimated to be 8%, showing that environmental factors could play a major role over genetics in susceptibility to BLV occurrence. Paradoxically, the implementation of good veterinary practices and proper management in order to reduce iatrogenic issues, considered one of the main factors of infection among cows within herds [44], failed to reduce the prevalence ([7, 43]). However, we focused on phenotypes representing the BLV level of infection in cows. The LI is directly related to the cows' viral transmission ability to BLV free herd mates. Methods incorporating GRM as used in this study, have been successfully tested in selective breeding of animals and plants for complex traits [72,73] and proposed to predict human disease phenotypes [74][75][76]. The high heritability values obtained for LI have remarkable practical implications in genetic improvement programs aiming to reduce the dissemination of BLV in heavily infected populations. Sources of non-additive genetic variance, such as dominance or epistasis among markers, were not considered in this study; however, the evidence indicates that the additive genetic variance is the major contributor to the phenotypic variance in complex characters [77][78][79]. To our knowledge, this is the first report that estimates the level of phenotypic variation in PVL and WBC explained by SNPs markers in BLV infected animals.
The molecular basis influencing the inter-individual difference in LI in vivo during aleukemic phase remains unexplained. The mechanisms of BLV infection involve an infectious cycle after a primary infection, where the virus replicates and spreads within the host [80] by cellto-cell virion transference, but this process is promptly abrogated by the anti-BLV host immune response [81]. Subsequently, the PVL is maintained by polyclonal expansion of surviving infected B-cells elapsed by mitotic division where viral factors hijack host cellular physiological pathways.
To evaluate which genes might be influencing the set-point of PVL and WBC counts, we identified genomic regions harbouring significantly associated SNPs (p < 1.13 × 10 − 6 ). These SNPs were located in the MHC region on BTA 23, the most relevant gene cluster of the immune sub-genome with high gene density and exceptional diversity.
Polymorphisms in the MHC loci −featured by an intricate structure of linkage disequilibrium [82][83][84][85]− are associated with viral infections outcomes, such as HIV-1, HBV, HCV and HPV [86][87][88]. Fellay et al. [89,90] confirmed through GWAS that polymorphisms in HLA are leading host genetic determinants on the viral load set-point in HIV−1 infected individuals. To investigate whether effects in others loci were masked by the BoLA region, we have also performed the same genome-wide analysis with the PVL level but excluding SNPs within the bovine MHC. No tested SNP reached significance level (Additional file 5: Table S1). However, when BLV incidence (based on ELISA test) was considered as phenotype for GWAS in US Holstein and Poland Holstein-Frisian populations, many genes outside of the MHC appears to be modulating this trait [40,41]. The phenotype evaluated in these studies (BLV-seronegatives vs BLV-seropositives cows) indicated novel potentially associated genomic regions, however they do not corroborate previous finding suggesting the MHC as the main genomic region harbouring genes influencing BLV infection [28,31,33].
Genomic positions highlighted by GWAS emphasised the potential influence of polymorphic variants in MHC class I and II genes related to their central role in cellmediated immunity. Four SNPs shared the functional term antigen processing and presentation (GO:0019882). The SNP with the lowest p-value for both conditions studied was rs110525467 (p PVL = 5.32 × 10 − 16 and p WBC = 5.90 × 10 − 10 ) producing a synonymous substitution in the MHC class II gene BoLA-DQA1. The rs110579760 (p PVL = 3.66 × 10 − 09 ) located 30.7 kb downstream BoLA-DRB3 was significantly associated with PVL, but not with WBC (Table 1). After infection, a strong cytotoxic T-cell response is developed in peripheral blood against specific viral epitopes [91]. In this regard, two significantly associated SNPs in strong LD (r 2 = 0.96) were functionally related with BoLA−A and SFTA2 (Surfactant Associated 2) genes. The magnitude of the immune response to viruses exerted by CD8+ T cells and its association with polymorphisms in HLA class I molecules is widely described in the literature [92][93][94]. In the case of BLV infections in cattle, previous studies have identified associations between polymorphisms in MHC class I genes and resistance to PL in different breeds [26,27,95,96]. The closest SNP (rs110794231) to BoLA-A is 35.2 Kb away corresponding to the SFTA2 gene. The SFTA2 gene encodes a secreted hydrophilic protein with surfactant properties with a role in host defence pathogen recognition patterns (bacteria, viruses and fungi) [97,98].
A particular haplotype block located in the MHC Class III region, encompasses a number of genes with a plethora of biological effects that may explain the associations detected (Fig. 4). These genes play relevant roles in the immune response, specifically in the regulation of apoptosis and proteasome degradation (LTA, TNF, BAG6) [107], HIV−1 progression (BAT1-5, LY6) [86], T cells response to interferon and viral shedding prevention during innate response to infections (LY6) [108,109].
Tumor necrosis factor (TNF) and LTA are proinflammatory cytokines classified to the GO BP "apoptotic process" (GO:0006915). Polymorphisms in TNF have been reported to contribute partially to the development of lymphosarcoma and PVL in BLV infected cattle [30,110] and early elimination of BLV in experimentally infected sheep [111].
Several of the targeted genomic regions in the GWAS might be involved in viral intra-host spreading after primary infection according to mechanisms proposed in [81]. For example, in the case of cells-cells fusions the SNP rs110473048 (p PVL = 1.78 × 10 − 12 and p WBC = 9.97 × 10 − 08 ) is at 7.21 kb from the Transcription Factor (TF) GCM1 (Glial Cells Missing 1) gene that regulates the expression of the ENV gene (located in the endogenous retrovirus HERV-W) [112] encoding for a fusogenic protein called Syncytin-1. Recently, it has been postulated that infection by an exogenous virus can trans activate ENV ectopically by increasing the expression of GCM1 [113] and promoting syncytia formation in breast cancer and endometrial carcinoma [114,115].
Two significantly associated SNPs mapped within and close to TRAM2 (Translocation-Associated Membrane Protein 2) gene, involved in the biosynthesis and excretion of type I collagen [116] that constitutes the extracellular matrix components [117] where virions are stored. The SNP rs109343703 causes an intronic substitution and the rs110499907 SNP is located 24.7 kb upstream of it. Both SNPs were associated with PVL and not with WBCs, and the nearest SNP in BoLA-DRB3 was in strong LD (r 2 > 0.85) with them. Another associated polymorphism (rs110836188) produces a synonymous substitution in the TNXB gene (Tenascin-X), which encodes an extracellular matrix protein involved in matrixcell adhesion (GO:0007160) and organization of collagen fibres (GO:0030199). Ultimately, analogous to "immunological synapses" some viruses subvert cellular processes to disseminate in the host through "virological synapses" [118,119]. In this process, cytoskeletal proteins are required for polarization of cells involved in the synapses [120,121]. Crucially, the SNP rs110155623 associated with both phenotypes is located in the putative regulatory region of CD2AP (CD2-associated protein). CD2AP organizes the cytoskeleton in polarization sites (GO:2,000,249) interacting directly with the Actin protein [122]. Two other GWAS hit variants could capture the effects of genes related to regulation and polymerization of the cytoskeleton structures. A substitution in an intron in the polycystic kidney and hepatic disease 1 (PKHD1) gene (rs41641297) and a SNP (rs110742604) positioned in regulatory regions of TUBB5 and FLOT1 genes. The TUBB5 gene encodes β5-tubulin which belongs to the β-tubulin family.
Finally, four SNPs were assigned to TFs ( Table 2) regulating transcriptional processes in promoters depending on RNA polymerase II (GO:0006357). TF genes are candidates for host determinants controlling BLV intra−host spreading, due to a tight mechanism of transcriptional regulation exerted by the trans-activating protein Tax, modulating several host signalling pathways such as AP-1, NF-κB and CREB to induce oncogenic transformation [123].
Altogether, in addition to classical immune response related genes, in this study we pointed novel genes/mechanisms as possibly involved in PVL control and hence BLV dissemination.
Despite our results do not conclusively identify causal genomic variants; we have successfully contracted the landscape of genes potentially related to BLV infections to a short list of very specific candidate genes suitable for experimental validations by gene-expression assays. These results also encourage fine mapping and validation studies. Furthermore, functional studies integrating multi-omics data will allow discerning the underlying molecular mechanisms of the BLV−host interaction.

Conclusions
Results obtained in this study revealed relevant genetic variants in the bovine MHC strongly influencing the control of BLV infection in cattle. Functional inspection of genes harbouring statistically significant associated SNPs brought preliminary evidence about possible biological pathways underlying the LI in infected cows. The proviral BLV load is associated with the stage of the disease progression in the infected cows and directly determines its ability to transmit BLV to healthy animals. Therefore, genome wide identification of genetics variants associated with low proviral load BLV infections would be useful not only to design therapeutics and preventive treatments but also as alternative control programs based on selective breeding of animals through genomic selection in dairy cattle.

Additional files
Additional file 1: Figure S1. Pedigree of the Holstein and Jersey crosses cattle population under study. Red lines connect the individual with his father. Blue lines connect the individual with his mother. (DOCX 252 kb).