Relationship between Allelic Heterozygosity in BoLA-DRB3 and Proviral Loads in Bovine Leukemia Virus-Infected Cattle

Simple Summary Bovine leukemia virus (BLV) caused a severe cattle neoplastic disease called enzootic bovine leukosis (EBL). EBL causes significant economic losses in farming by reducing milk production, reproductive performance, and fertility, and through cattle culling or death. The BLV proviral load (PVL) represents the quantity of BLV genome that has integrated into the host’s genome in BLV-infected cells. It has been reported that PVLs differ according to the genetic background of the host, and some studies of BLV-associated host factors have reported on polymorphisms within the bovine major histocompatibility complex (MHC), namely bovine MHC is bovine leukocyte antigen (BoLA-DRB3). However, there is a great diversity in the PVLs associated with carrying various combinations of these alleles, especially for heterozygous alleles. Therefore, this research investigated whether heterogeneity in BoLA-DRB3 allele combinations would affect PVLs during BLV infections in different ages and breeds of cattle in Japan. This is the first report where the association between heterozygous allelic combinations and BLV PVLs phenotypes (HPLs, LPLs) was analyzed. Our findings augment current understanding about the immunological role played by BoLA heterozygosity in BLV-associated PVLs and biocontrol in BLV infections. Abstract Enzootic bovine leukosis is a lethal neoplastic disease caused by bovine leukemia virus (BLV), belongs to family Retroviridae. The BLV proviral load (PVL) represents the quantity of BLV genome that has integrated into the host’s genome in BLV-infected cells. Bovine leukocyte antigen (BoLA) class II allelic polymorphisms are associated with PVLs in BLV-infected cattle. We sought to identify relationships between BoLA-DRB3 allelic heterozygosity and BLV PVLs among different cattle breeds. Blood samples from 598 BLV-infected cattle were quantified to determine their PVLs by real-time polymerase chain reaction. The results were confirmed by a BLV-enzyme-linked immunosorbent assay. Restriction fragment length polymorphism-polymerase chain reaction identified 22 BoLA-DRB3 alleles. Multivariate negative binomial regression modeling was used to test for associations between BLV PVLs and BoLA-DRB3 alleles. BoLA-DRB3.2*3, *7, *8, *11, *22, *24, and *28 alleles were significantly associated with low PVLs. BoLA-DRB3.2*10 was significantly associated with high PVLs. Some heterozygous allele combinations were associated with low PVLs (*3/*28, *7/*8, *8/*11, *10/*11, and *11/*16); others were associated with high PVLs (*1/*41, *10/*16, *10/*41, *16/*27, and *22/*27). Interestingly, the BoLA-DRB3.2*11 heterozygous allele was always strongly and independently associated with low PVLs. This is the first reported evidence of an association between heterozygous allelic combinations and BLV PVLs.


Introduction
Enzootic bovine leukosis (EBL), a lethal neoplastic disease caused by bovine leukemia virus (BLV), belongs to the Retroviridae family (genus Deltaretrovirus) [1]. BLV infection has three stages: the silent aleukemic stage; the persistent lymphocytosis (PL) stage, which is characterized by the polyclonal expansion of non-neoplastic CD5+ B-lymphocytes; and the EBL stage, as manifested by malignant CD5+ B-cell lymphoma [2]. While approximately 30% of BLV-infected cattle progress to the PL stage, less than 5% of them develop fatal B-cell lymphosarcoma after a long latent period [3]. EBL causes significant economic losses in farming by reducing milk production, reproductive performance, and fertility, and through cattle culling or death [4]. The proviral load (PVL) of BLV represents the quantity of its genome that has integrated into the genomic DNA of BLV-infected cells. BLV infected animals with high levels of PVL is most likely to develop malignant CD5+ B-cell lymphoma. Thus, PVL is an important index for estimating the stage of a BLV infection in a susceptible animal because it is associated with disease progression, making it a useful way of determining BLV transmission [5].
However, it has been reported that PVLs differ according to the genetic background of the host, and some studies of BLV-associated host factors have reported on polymorphisms within the bovine major histocompatibility complex (MHC) [6][7][8][9]. The bovine MHC, bovine leukocyte antigen (BoLA), which is located on chromosome 23, includes three DRB class II genes. BoLA-DRB3 is highly polymorphic with 130 described alleles and is the only active one of the three DRB genes [10]. Polymorphisms in BoLA-DRB3 are reportedly associated with susceptibility to several infectious diseases (e.g., BLV-induced lymphocytosis, mastitis, and dermatophilosis) [11]. However, there is a great diversity in the PVLs associated with carrying various combinations of these alleles, especially for heterozygous alleles. Therefore, this research investigated whether heterogeneity in BoLA-DRB3 allele combinations would affect BLV viral loads during BLV infections in different ages and breeds of cattle in Japan.

Blood Samples
Peripheral blood samples (n = 598) were collected in EDTA tubes from different cattle breeds from four commercial beef farms in Kyushu, Japan. These farms raised mainly Japanese Black cattle (523), Holstein cattle (68), Crossbreed (4) and Jersey (3) ( Table 1). All animals recruited in this study were apparently healthy without showing any clinical signs. The samples were collected in January to December 2019. DNA was extracted from 300 µL of each blood sample within a week of its collection using the Wizard ® Genomic DNA Purification Kit (Promega, Madison, WI, USA) according to the manufacturer's instructions. Purified DNA samples were each suspended in 50 µL of distilled water. Plasma was separated from 1000 µL of each blood sample by centrifugation (1500× g, 5 min, 4 • C) and then stored at −20 • C until analysis.

Real-Time PCR Quantification of BLV-Associated PVLs
The concentration and purity of each DNA sample was determined with a NanoDrop 8000 spectrophotometer (Thermo Fisher Scientific, Wilmington, DE, USA), and each sample was diluted to 20 ng/µL. Proviral loads were quantified on the Applied Biosystems 7300 Real-Time PCR System (Applied Biosystems, CA, USA). The reaction mixture contained 5 µL of 2× Cycleave PCR ® Reaction Mix SP (Takara-Bio Inc., Shiga, Japan), 0.2 µL of probe/primer mix for BLV, 0.1 µL of ROX™ Reference Dye II (Takara-Bio Inc., Shiga, Japan), 1 µL of DNA template (20 ng/µL), and PCR-grade water to a final 10 µL volume. To determine PVLs, calibration curves were generated from the measured concentrations of a dilution series of positive control DNA containing the BLV TAX gene (Takara-Bio Inc., Shiga, Japan) and the PVL was expressed as the number of provirus copies/50 ng of genomic DNA [12].

Diagnosing BLV Infections Using BLV-Enzyme-Linked Immunosorbent Assay (ELISA)
All plasma samples were analyzed using a BLV gp51 antibody detection ELISA kit (JNC, Tokyo, Japan) according to the manufacturer's instructions.

BoLA-DRB3 exon2 Cloning and Sequence Analysis
Ambiguous RFLP-PCR results were subjected to sequence analysis verification. Each amplicon extracted from an agarose gel using the QIAquick Gel Extraction Kit (QIAGEN, Hilden, Germany) was cloned into T-vector pMD20 (Takara-Bio Inc., Shiga, Japan), and sequenced. Sequencing was performed in both directions using M13 forward and reverse primers with the Big Dye ® Terminator v3.1/1.1 cycle sequencing kit (Applied Biosystems). The resulting data were analyzed by the Applied Biosystems 3730 DNA Analyzer.

Statistical Analyses
Hierarchical clustering in R (version 3.2.2, https://cran.ism.ac.jp/bin/windows/ base/old/3.2.2/ (accessed on 14 August 2015)) by Euclidean distance calculation was used to statistically group the BLV PVLs into the following three clusters: the high PVL group (HPL), the medium PVL group (MPL), and the low PVL group (LPL). BoLA-DRB3 allele frequencies were calculated. Spearman's rank correlation rho test was conducted on BLV PVLs versus cattle age. Pairwise Wilcoxon rank sum testing was conducted to compare the BLV PVLs among the cattle breeds. The mean PVL for each allele combination was calculated on a heatmap. Violin plots were used to show the BLV PVL distribution patterns for each BoLA-DRB3 allele combination. Associations between the different BoLA-DRB3 alleles and BLV PVLs were estimated using the Poisson regression model. The Poisson distribution assumes that the expected and variance values are equal, but this is not always true, causing dispersion of the data when the variance is higher than average. When overdispersion occurs, one of the ways to estimate its parameter is to use the negative binomial distribution [15], which is why we used this regression model in the (MASS) package (http://www.stats.ox.ac.uk/pub/MASS4/ (accessed on 5 May 2003)). The BLV PVL was the outcome, and BoLA-DRB3 alleles, cattle age and breed were the explanatory variables. Regression model is conducted into two steps: Univariate regression analysis which used for all variables (BoLA-DRB3 alleles, cattle age, and breed). Any variables significantly associated with BLV PVLs at the p < 0.20 level were subsequently selected for multivariate regression analysis. The final model was obtained with p levels for the remaining variables of <0.05 [16]. Odds ratios (ORs) were calculated using https: //www.medcalc.org/calc/odds_ratio.php (accessed on 1 April 2017) to define the BLV susceptible and resistant alleles in the cattle, and p-values of less than 0.05 were considered statistically significant. It was reported that MPL and LPL groups were a low-risk spreaders group and HPL was a high-risk spreaders group. They also assumed that if OR > 1 it will be as resistant allele, OR < 1 it will be as susceptible allele [17].

Quantifying BLV-Associated PVLs and Grouping Cattle According to These Loads
Plasma samples (598) from cattle were confirmed to be BLV-positive without showing of any clinical signs through anti-gp51 antibody detection using BLV-ELISA. The BLV PVLs in peripheral blood mononuclear cells from these cattle were quantified by real-time quantitative PCR, as described earlier. Real-time PCR produced negative results for the 119 samples deemed positive by BLV-ELISA. The PVLs for BLV in terms of copies/50 ng were as follows: minimum 5.0 (log 10 0.699), median 767.8 (log 10 2.885), and maximum 9107.0 (log 10 3.959) (Figure 1). The samples fell into three distinct groups based on the Euclidean distance calculation between each BLV PVL in R by hierarchal clustering. Table 2 shows their statistical grouping (278 HPL, 62 MPL, and 258 LPL).

BLV PVLs and Their Relationships with Cattle Age and Breed
The study cattle population comprised two main breeds (Holstein and Japanese Black), and various ages. Therefore, we investigated the relationships between PVL, and the different cattle breeds and ages. There was no significant correlation between BLV PVLs and cattle age (Figure 2A). When we conducted comparison between the BLV PVLs in the different cattle breeds, no significant difference between different cattle breeds and BLV PVLs was found ( Figure 2B).

Relationships among BLV and BoLA-DRB3 Alleles, Age, and Breed
We conducted RFLP-PCR typing of BoLA-DRB3 alleles to define their relationship with the different BLV PVL groups. Twenty-two different BoLA-DRB3 alleles were identified in the 598 cattle samples. Table 3 shows the defined allele frequencies and median values for log 10 BLV PVL versus each BoLA-DRB3 allele. Based on these results we showed each allele frequencies among the main two cattle breeds in this study (Japanese Black (JB) and Holstein (H)). BoLA-DRB3.2*3, *8, *11, *24, *28 were more frequent among JB) than H. While BoLA-DRB3.2*7 was only appeared as a heterozygous allele in JB cattle. However, BoLA-DRB3.2*22 was more frequent among H than JB (Table S1). All these previously mentioned alleles were categorized as belonging to LPL, whereas BoLA-DRB3.2*1, *10, *16, and *41 belonged to HPL (Figure 3). BoLA-DRB3.2*1 and BoLA-DRB3.2*41 were only appeared in JB cattle. However, BoLA-DRB3.2*10 and BoLA-DRB3.2*16 were more frequent among JB than H (Table S1). Overall, the inbreeding of the Japanese Black and Holstein herds may have a greater interfere with the prevalence of specific alleles.
Associations between PVL in BLV infections and different BoLA-DRB3 alleles, age, and breeds were estimated using a univariate negative binomial regression model. The calculated confidence intervals (95%) of the log 10 PVL median values of all the BoLA-DRB3 alleles showed the wide range and diversity of the BLV-associated PVLs for each allele from the different combinations, especially the heterozygous ones (Table 3). Multivariate negative binomial regression modeling and OR calculations were conducted for the BLV PVLs and the different BoLA-DRB3 alleles to define which cattle carried BLV susceptible (OR < 1) or resistant alleles (OR > 1). Table 4 shows that BoLA-DRB3.2*11, *3, *7, *8, *22, *28, *20, and *24 alleles were associated with LPLs and MPLs in BLV infections, whereas BoLA-DRB3.2*10 was only associated with HPLs in the infections. Statistically significant allele frequencies were calculated for the different cattle breeds.   The alleles were determined to be resistant when OR > 1 and susceptible when OR < 1, LPL: low proviral load, MPL: medium proviral load, HPL: high proviral load. Bold means significant p value.

Discussion
Our RFLP-PCR genotyping identified 22 previously reported different BoLA-DRB3 alleles. While our multivariate negative binomial regression modeling and OR calculations showed that BoLA-DRB3.2*11, *8, *28, *7, *3, *24, and *22 alleles were significantly associated with LPLs and MPLs in the BLV infections, the BoLA-DRB3.2*10 allele was significantly associated with HPLs. Nevertheless, the heat map of the log 10 mean values of the BLVassociated PVLs for each BoLA-DRB3 allele combination and the violin plots showed the possibility of various allele combinations affecting PVLs in the BLV infections.
BoLA-DRB3 is a highly polymorphic gene and several of its alleles are favorably or unfavorably associated with BLV-associated PVLs and disease progression. Here, we found that the BoLA-DRB3.2*7 allele was statistically associated with LPLs (p-value < 0.01, OR = 1.98), a result supported by a previous study where the BoLA-DRB3*002:01 (DRB3.2*7) allele was reportedly associated with BLV resistance and LPLs (OR > 1) [16]. Surprisingly, we found that BoLA-DRB3.2*7 in combination with BoLA-DRB3.2*10 was associated with HPLs. Our finding that BoLA-DRB3.2*10 is associated with HPLs (OR=0.69) is consistent with the previously reported findings that BoLA-DRB3*016:01 (DRB3.2*10) was associated with HPLs (OR < 1) [17,18]. It is possible that the BoLA-DRB3*016:01 (DRB3.2*10) allele has higher affinity than the BoLA-DRB3*002:01 (DRB3.2*7) allele for antigen presentation to T cell receptors, and as a susceptible allele this causes HPLs. It is also possible that the polymorphisms in the peptide-binding cleft in MHC class II molecules that are responsible for the affinity of the bound peptides influence the magnitude and extent of the elicited adaptive immune response and, therefore, the PVL [19]. High-dose pathogen exposure is also postulated to overcome the immune system, even in the most resistant animals [20]. Monroy et al. reported that the nature of the interaction between two heterozygous alleles combination of the BoLA-DRB3 gene and BLV resistance or susceptibility is not yet clear [21]. Indeed, although 70% of cattle carrying RR (resistant-resistant) alleles were not infected with BLV, only 33% of SS (susceptible-susceptible) cattle were infected [21]. Moreover, the BoLA-DRB3.2*7-BoLA-DRB3.2*16 combination is associated with a wide and variable PVL range between MPL and HPL. BoLA-DRB3*015:01 (DRB3.2*16) was also determined to be an allele for BLV susceptibility related to HPL (OR < 1) [17]. Furthermore, RR (resistant-resistant) genotypes and SS (susceptible-susceptible) genotypes were found to be associated with LPLs and HPLs, respectively, whereas RS (resistant-susceptible) genotypes covered a wide PVL range [8]. While our results show that the BoLA-DRB3.2*8 allele is statistically associated with LPLs (p-value < 0.001, OR = 1.40), elsewhere cattle carrying the BoLA-DRB*012:01 (DRB3.2*8) allele were reported to be at low risk of BLV infection [22]. Thus, we found that BoLA-DRB3.2*7 allele combinations had synergistic effects and were mostly associated with LPLs. Consistent with this finding, a previous report showed that the RR (resistant-resistant) allele combination carried by all the BLV-infected study animals was associated with LPLs [8]. This may also reflect the high degree of plasticity related to the function of different genes and variants of the MHC system together with other immune defense genes.
Interestingly, highly variable PVL ranges have been reported when different allele combinations occur in the presence of heterozygous susceptible BoLA-DRB3 alleles. BoLA-DRB3.2*10 and *16 are reported to be the susceptible alleles associated with HPLs [17,18,22,23]. Although widely variable PVLs were associated with other described allele combinations (e.g., BoLA-DRB3.2*1, *23, *24, *27, among others) in our violin plots, the homozygous susceptible alleles were mostly associated with HPLs. Heterozygous individuals potentially have a wider variety of effective T lymphocytes than homozygous individuals, which generates a more diverse T-cell repertoire [24]. The susceptible BoLA-DRB3*006:01 (DRB3.2*10) allele, which encodes Glu-Lys (EK) at positions 70 and 71 in the BoLA-DRβ chain, may be associated with HPLs [18]. It has been shown that TNF-α expression levels increased significantly in sheep that were experimentally infected with BLV, and that this was responsible for virus elimination in the early infection stages [25]. When associated with a polymorphism in the MHC class II β chain, polymorphism in the TNF-α promoter gene might influence the host's response to BLV infection, something that might also occur as a consequence of linkage disequilibrium [8]. In BLV infections, it has been reported that IL-2 production is associated with asymptomatic animals, and that IL-10 production increases in animals with lymphocytosis [26]. These cytokines are mainly produced by activated antigen specific CD4+T lymphocytes, and their activation depends on the presentation of antigenic epitopes by MHC class II molecules. Therefore, it is reasonable to find that certain MHC class II alleles can confer strong protective immunity in BLV-infected animals, whereas others may fail [9].
Another conceivable mechanism controlling PVLs may be related to previously identified haplotypes via the strong linkage between BoLA class II DR and DQ genes [18]. The class II BoLA molecule repertoire is not limited to the intra-haplotype combination of α/β chains; instead, it can be generated by the inter-haplotype combination of α/β chains, as seen among highly polymorphic BoLA DQA (more than two loci) and DQB (more than two loci) genes [27], because these inter-haplotype α/β chain combinations exist with heterozygous BoLA class II alleles in cattle. We hypothesize that the results from the above-mentioned studies indicate that the great diversity seen with BLV-associated PVLs may result from inter-haplotype pairing between the DQA and DQB genes of heterozygous allele combinations. A more detailed study of the linkage between DR and DQ genes and its effect on BLV PVLs is underway. Further studies in this area are important because it has been reported that many of the genetic and epigenetic factors that are involved in BLV infection exist in addition to BoLA-DRB3 alleles [8].
We identified various BoLA-DRB3 allelic types that are associated with MPL or LPL (BoLA-DRB3.2*11, *7, *28, *8, *3, *24, and *22, but not limited to BoLA-DRB3.2*11). The selection of individuals to be parents of the next generation of cattle on the basis of protective alleles is the next logical step in the direction of creating more disease-resistant stock; however, complete information with regards to the different components of BoLA such as the haplotype combination of DQ and DR genes is necessary. It was reported that PVLs vary according to the genetic background of the host, despite having the same clinical status [7]. Thus, in addition to improving production traits, the future application of our findings could also be simultaneously applied to genetic selection as a sort of biocontrol in an integrated strategy. The limitation of this study is that the sample size was not large enough to show significant differences in the BLV PVL values for each heterozygous allele combination. This is the first report where the association between heterozygous allelic combinations and BLV PVLs (HPLs, LPLs) was analyzed. Our findings show that the heterozygous BoLA-DRB3.2*11 allele is always strongly and independently associated with very low PVLs, whatever the allele combinations. Moreover, allelic divergence and polymorphism in heterozygous allele combinations were found to greatly affect PVLs in BLV infections.
Our findings enhance current understanding of the immunological role played by BoLA heterozygosity in BLV-associated PVLs and the biocontrol of BLV infections.
Supplementary Materials: The following are available online at https://www.mdpi.com/2076-261 5/11/3/647/s1, Table S1: BoLA-DRB3 allele frequencies among the main two cattle breeds in this study Japanese Black (JB) and Holstein (H). Figure  allele combinations. The green dot represents the median, the thick light blue bar in the center represents the interquartile range, and the thin light blue line represents the rest of the distribution, except for the red points determined to be "outliers" using a method that is a function of the interquartile range. The colored area on each side of the light blue line is a kernel density estimation to show the distribution shape of the data.

Institutional Review Board Statement:
The authors confirm that the ethical policies of the journal, as noted on the journal's author guidelines, were adhered to and that the study received approval via appropriate ethical review by the Animal Ethics Committee of the Faculty of Agriculture at the University of Miyazaki (University of Miyazaki, 1-1 Gakuen-Kibanadai-Nishi, Miyazaki 889-2192, Japan) (number 2018-03-29-Z20).

Data Availability Statement:
The data presented in this study are available in supplementary martials.