Microarray analysis of differential gene expression profiles in blood cells of naturally BLV-infected and uninfected Holstein–Friesian cows

The aim of the present study was to examine gene expression changes in response to bovine leukemia virus (BLV) infection, in an effort to determine genes that take a part in molecular events leading to persistent lymphocytosis (PL), and to better define genes involved in host response to BLV infection. Using bovine 70-mer oligonucleotide spotted microarrays (BLOPlus) and qRT-PCR validation, we studied global gene expression profiles in blood cells in vivo of 12 naturally BLV-infected Polish Holstein cows, and 12 BLV non-infected controls of the same breed and reared in herds with high BLV seroprevalence. With an arbitrary cut-off value of 1.5-fold change in gene expression, we identified the down-regulation of 212 genes (M value ≤−0.585) and the up-regulation of 158 genes (M value of ≥0.585) at 1% false discovery rate in BLV-positive animals in comparison to the BLV-negative group. The gene set enrichment analysis demonstrated that the differentially expressed (DE) genes could be classified to diverse biological processes, including immune response of host blood cells. Interestingly, our data indicated the potential involvement of the innate immunity, including complement system activation, NK-cell cytotoxicity and TREM-1 signaling, during the BLV-induced pathogenesis. We showed the occurrence of numerous regulatory processes that are targeted by BLV-infection. We also suggest that a complex network of interrelated pathways is disturbed, causing the interruption of the control of B-cell proliferation and programmed cell death. Electronic supplementary material The online version of this article (doi:10.1007/s11033-016-4088-6) contains supplementary material, which is available to authorized users.


Introduction
Enzootic bovine leukosis (EBL) is the most prevalent, lethal neoplastic disease of cattle, causing considerable economic losses worldwide. The etiological agent of EBL is a B-lymphocytotropic deltaretrovirus, named BLV (Bovine Leukemia Virus), which infects cattle and has distinct clinical outcomes. Most cattle infected by BLV remain clinically silent in an aluekemic state, approximately 30% of infected animals develop a persistent lymphocytosis (PL) which is a benign polyclonal proliferation of B cells, and finally, lethal lymphoma or lymphosarcoma appears in less than 5% of infected animals, after a long latency period [40]. After infection, BLV distributes in the host as integrated proviral DNA, spreading via the cell divisions of the infected leukocytes [20]. The course of leukemogenesis caused by BLV, likewise the mechanisms underlying the aspect of host resistance/susceptibility to BLV infection and disease progression, are intricate and obscured. Although it is not entirely clear, it seems likely that the subtle balance between viral gene expression and efficient host immune response plays a decisive role in the outcome of BLV infection. Furthermore, tumor development in EBL is preceded by an accumulation of chromosomal aberrations and mutations in oncogenes e.g. the p53 gene [8,18].
An apparent puzzle of BLV-induced leukemogenesis is the fact that in vivo BLV is quiescent in the majority Abstract The aim of the present study was to examine gene expression changes in response to bovine leukemia virus (BLV) infection, in an effort to determine genes that take a part in molecular events leading to persistent lymphocytosis (PL), and to better define genes involved in host response to BLV infection. Using bovine 70-mer oligonucleotide spotted microarrays (BLOPlus) and qRT-PCR validation, we studied global gene expression profiles in blood cells in vivo of 12 naturally BLV-infected Polish Holstein cows, and 12 BLV non-infected controls of the same breed and reared in herds with high BLV seroprevalence. With an arbitrary cut-off value of 1.5-fold change in gene expression, we identified the down-regulation of 212 genes (M value ≤−0.585) and the up-regulation of 158 genes (M value of ≥0.585) at 1% false discovery rate in BLV-positive animals in comparison to the BLV-negative group. The gene set enrichment analysis demonstrated that the differentially expressed (DE) genes could be classified to diverse biological processes, including immune response of host blood cells. Interestingly, our data indicated the potential involvement of the innate immunity, including complement system activation, NK-cell cytotoxicity and TREM-1 signaling, during the BLV-induced pathogenesis. We showed the occurrence of numerous regulatory processes that are targeted by BLV-infection. We also suggest that a complex Electronic supplementary material The online version of this article (doi:10.1007/s11033-016-4088-6) contains supplementary material, which is available to authorized users. 1 3 of infected cells, however, vigorous humoral immune response has been elicited in BLV-infected cattle, indicating that BLV is not completely silent. Furthermore, the latent state is immediately reversed when infected lymphocytes are subjected to transient ex vivo cell culture. The potential repression mechanism of viral expression involves a poorly defined plasma blocking factor (PBF), related to fibronectin and inhibited by a platelet lysate [11,13,50].
A number of reports suggest that the multifunctional viral accessory protein Tax is the major transforming protein of BLV (reviewed in [13]). The ability of this viral transactivator to both transcriptionally regulate cellular gene expression and directly interact with cellular proteins provides the basis for leukemogenesis [22]. Microarray transcription profiling of ovine immortalized Clone2 B-cell line and peripheral blood-derived ovine B cells cultured in vitro [22] and human HeLa cells [4] transfected with Tax BLV constructs revealed that in vitro Tax was an efficient modulator of host gene transcription. A vast list of Tax-responsive genes contained, among others, genes encoding for cytokines, immune modulators, transcription factors, genes associated with cell cycle, DNA repair and apoptosis, signaling factors and adhesion molecules which could be involved in BLV-induced pathogenesis. It is not yet known whether similar changes in host gene expression occur in vivo during EBL progression in naturally infected cattle.
Evidence to date suggests that host genetic factors may play a role at different stages of infection with BLV: from initial infection to the likelihood of developing PL, through to the fatal lymphoma [14]. It was proved that the resistance and susceptibility to PL development was linked with the MHC class II BoLA-DRB3 gene [26]. The resistance occurs to be dependent upon the presence of polar motif ER at positions 70-71 within a highly polymorphic region encoding the putative peptide binding domain [52]. Recently, one of the alleles containing the ER motif, designated as BoLA-DRB3.2*0902, was shown to be significantly associated with genetic resistance to PL and with a low proviral load profile (LPL) [20]. It was proposed that animals presenting LPL profile and harboring the BoLA-DRB3.2*0902 allele might be unable to transmit BLV in herd conditions, and this suggestion underlay a markerassisted selection scheme for the simultaneous eradication of BLV and improvement of milk production traits implemented in Argentina [11]. However, the extent of variation in LPL profile not associated with the favorable BoLA-DRB3.2*0902 allele [11] implies that genetic resistance to BLV and the outcome of the disease could be under the control of multiple as yet unknown or unidentified genes, each contributing slightly to the phenotype. In addition, polymorphism in the predicted enhancer region of the bovine TNF-α (tumor necrosis factor, allele −824G associated with low transcription activity) was shown to contribute in the progression of BLV-induced lymphoma [23].
As variation in gene expression could be a causative mechanism underlying BLV resistance/susceptibility, the application of gene expression microarray technology enabling the large scale transcriptome profiling of host immune response promises the identification of networks of genes involved in disease pathogenesis, providing a key for understanding and, perhaps in the future, intervening against successful BLV evasion of host biodefenses. The aim of the present study was to examine gene expression changes in response to BLV infection, in an effort to determine genes that take a part in molecular events leading to PL, and to better define genes involved in host response to BLV infection. We focused on fresh whole blood cells, with their transcriptome freezed and captured immediately after bleeding using RNA stabilization solution, without culturing them because most of these cells in vivo usually do not express viral proteins including Tax, in contrast to either cultured PBMCs from BLV-positive animals or BLVinfected cell lines.

Animals, sample collection and BLV diagnostics
All procedures were approved by the local ethics committee of University of Warmia and Mazury in Olsztyn (no. 13/2008/N/T). Blood samples were collected from two dairy herds, kept in a free stall barn system and located in north-eastern and north-western parts of Poland respectively. BLV seroprevalence rates in both herds exceeded 70%. All sampled animals belonged to the Polish Holstein-Friesian breed and were reared and maintained in the infected herd at least 3 years. The blood was taken from jugular vein of lactating cows using Vacuette® Blood Collection Set and evacuated blood collection tubes: Vacu-ette® Serum, Vacuette® EDTA (Greiner Bio-One) and Vacuette® Tempus™ Blood RNA Tube (Applied Biosystems, Greiner Bio-One). In other to minimize seasonal effects, all the samples were collected within a single month (November, 2012). BLV serological diagnostics was carried out by ELISA technique using Pourquier® ELISA Bovine Leukosis Screening Kit (Institute Pourquier) according to the manufacture's protocol. Moreover, the diagnosis was validated with a molecular approach based on BLV proviral DNA detection using nested-PCR method [29]. Automated blood cell counting was performed with SF3000 automated hematology analyzer (Sysmex) according to the manufacture's instructions. Based on serological (ELISA), molecular (nested-PCR proviral detection) and 1 3 hematological assays, the experimental group was composed of 12 BLV-infected individuals with apparent hematological alterations characteristic for the lymphoproliferative stage of enzootic bovine leukosis (average lymphocyte count × 10 3 /μl, 15.16 ± 10.28) and the control group was composed of 12 BLV free animals with normal hematological blood indices (average lymphocyte count × 10 3 /μl, 3.43 ± 1.22). The difference in average age of the animals in the experimental group (5 years 2 months, SD = 0.81) and of the control group (4 years 8 months, SD = 0.58) was not significant.

RNA extraction
Blood samples of 3 ml volume were collected into Vacu-ette® Tempus™ Blood RNA Tube (Applied Biosystems, Greiner Bio-One) containing 6 ml RNA stabilization solution, mixed vigorously, transported on ice within a few hours to the laboratory and then frozen in −20 °C until RNA isolation. To minimize handling time during the RNA extraction procedure, samples were processed in small batches of 4. After thawing for 30 min. at room temperature each sample was poured into 50 ml Falcon conical tube (BD) and 3 ml of 95% molecular biology grade ethanol (Applichem) was added. The mixture was intensively vortexed on high speed for 5 min. and then centrifuged at 5200×g for 60 min at 0 °C. The foaming and the supernatant were poured off and the rim of the inverted Falcon tube was blotted on clean absorbent paper for 1 min. 500 µl of lysis solution supplemented with 5 µl of TCEP solution from PerfectPure™ RNA Cell and Tissue (5 Prime) was added and RNA pellet was dissolved by vortexing for 1 min. The lysate was pipetted onto a purification column (5 Prime). All other procedures including successive centrifugations, washing steps, DNase treatment and RNA elution were according to PerfectPure™ RNA Cell & Tissue manual except prolonged to 45 min. On column DNase digestion at room temperature. Purified RNA concentration, purity, and integrity were determined by A 260 , A 280 and A 230 measurements using a NanoDrop ND1000 Spectrophotometer (NanoDrop Technologies) and Agilent 2100 Bioanalyzer with Agilent RNA 6000 Nano Kit (Agilent Technologies). The remaining samples were aliquoted and stored at −80 °C until use.

Oligonucleotide microarray, aRNA synthesis and labelling, microarray experimental design, hybridization, and image acquisition
The bovine microarray platform used was Bovine Long Oligo Extension (BLOPlus), a 10 K spotted 70-mer oligonucleotide microarray (GEO platform accession: GPL9176), manufactured and purchased from Center for Animal Functional Genomics, Michigan State University. Cy3/Cy5 dye (GE Healthcare) labelled, linearly amplified, antisense aRNA was generated from 1 µg of each sample total RNA using Amino Allyl MessageAmp TM II aRNA Amplification kit (Ambion, Life Technologies) according to the manufactures' instructions. The quantity and the quality of the resulting aRNA was analysed by spectrophotometry and and Agilent 2100 Bioanalyzer with Agilent RNA 6000 Nano Kit (Agilent Technologies). The amount of fluorescent dye incorporated into aRNA was measured as its absorption at 550 nm (Cy3) and 650 nm (Cy5) using a NanoDrop ND1000 Spectrophotometer MicroArray builtin module (NanoDrop Technologies). Preliminary microarray dyE−swap experiments were performed in order to optimize the conditions of microarray hybridization, washing and scanning protocol. Afterwards a common reference design was employed for microarray hybridizations, thus all the aRNAs were labelled with Cy5 and co-hybridized with Cy3 labelled common reference pool aRNA. The common reference aRNA pool was set up by mixing equal amounts of total RNA from all samples in the study and pooled RNA was subsequently used for aRNA linear amplification. Twenty-four arrays were hybridized, representing 12 BLV-infected and 12 non-infected cows. For each microarray slide 60 pmol of Cy5 labelled aRNA was combined with 60 pmol of Cy3 labelled common reference aRNA (5-10 μl of total volume) diluted with SlideHyb#1 glass array hybridization buffer (Ambion, Life Technologies) up to 130 μl and denatured at 70 °C for 5 min. Microarray hybridizations were carried out using a Tecan HS400 Pro hybridization station (Tecan) according to following proto-

Microarray data processing, normalization and analysis
The linear models for microarray data (LIMMA, [45]) software was used to identify differential gene expression. Background correction was performed and within microarray print-tip LOWESS (locally weighted scatterplot smoothing) normalization was employed. No between array normalization was needed (data not shown). For each gene spot, t-statistic, B-statistic, M value (log 2 [Cy5/Cy3]), A value (1/2 [log 2 Cy5 + log 2 Cy3]) were calculated. Probability p values were corrected for multiple testing using the false discovery rate (FDR) [5]. An arbitrary fold change of FC ≥1.5 (represented by −0.585 ≤ M value ≥ 0.585) and FDR-adjusted (adj. p ≤ 0.01) significantly differentially expressed (DE) genes were subjected to meta-analysis using a web-based GeneGo MetaCore™ analytical suite (ver. 5.3, GeneGo, Thomson Reuters). The uploaded list of DE gene IDs was submitted to gene ontology (GO) over-representation analysis including functional categories: GO Processes, GO Molecular Functions, GO Localizations, GeneGo Process Networks and GeneGo Diseases. The in silico predicted functional protein networks were generated using the Analyze Networks (AN) buildin algorithm. Additionally, a list of transcription factors related to the expression regulation of DE genes was conceived by AN Transcription Regulation algorithm.

Validation of microarray results by real-time qRT-PCR
The validation of gene expression data was performed by qRT-PCR for 14 genes (up-regulated, down-regulated and with no regulation according to microarray data) using a LightCycler® LC 480 II system (Roche). For each target gene a total of 24 experimental cDNA samples in duplicate (from 12 BLV+ and 12 BLV− individuals) and a triplicate of no template control (NTC, molecular grade water) were amplified. All procedures including cDNA synthesis, primers design, PCR reaction set up, cycling conditions, melting curve analysis, amplicon's size and specificity control, standard curve preparation, PCR efficiency (E) and Error value calculation and quantification cycle (Cq) determination were performed as in [6].The primers and the target amplicons information are presented in Table 4. For qRT-PCR data normalization, a pair of RPLP0 and UCHL5 was employed, as the most stable reference genes for this experiment, based on a previous determination from a panel of 10 putative reference gene candidates [6]. The gene expression differences between BLV-infected and BLV uninfected groups were estimated with Relative Expression Software Tool REST 2009 V2.0.13 [37].

Microarray analysis of differentially expressed genes in blood cells of naturally BLV-infected and uninfected Holstein cattle
To better understand the role of the genes involved in the host response to BLV infection and BLV-induced lymphoproliferation and malignant transformation, we performed gene expression profiling of BLV-infected and uninfected Holstein cattle. The applied microarray platform comprised the two-color, 70-mer oligonucleotide spotted Bovine Long Oligo Plus (BLOPlus, 10 K) array and the obtained data was deposited in NCBI Gene Expression Omnibus with accession no. XXXXXX. An analysis of mRNA abundance in blood samples from 12 BLV infected cattle and 12 non-infected controls was performed. The BLV infected group comprised animals ca. 5 years old, seropositive with apparent hematological abnormalities (Table 1). Due to the strict BLV eradication policy in Poland, we were able to blood sample the BLV-positive individuals only once, so we could not officially categorize them as in the persistent lymphocytosis (PL) stage of disease. However, the high total lymphocyte counts (15.16 ± 10.28 × 10 3 /μl, Table 1), high total white blood cell counts (WBC), high titers of antibodies against p24 and gp51 and qualitative estimation of BLV proviral loads by PCR products densitometry analysis [29], (data not shown) indicated that this certainly was the case. The BLV-uninfected control group comprised cows of the same breed, similar age (Table 1), diagnosed as BLV-negative by both serological and molecular approaches. It is of interest to emphasize the fact that all the studied BLV-negative animals were housed for at least 3 years in the same 2 BLV-positive herds as their BLV-infected counterparts. Nevertheless they did not become infected. As resistance and susceptibility are opposite sides of the same coin, it is tempting to treat the BLV-negative group not only as the control group, but also to reverse the microarray scheme and to suspect the BLV-negative animals as constituting some degree of innate resistance to BLV infection. Because of the critical dependence of gene expression measurements on starting RNA quality, we included in our analysis only samples which were marked by high RNA quality indices, and without any significant differences between the compared groups (Table 1). A total of 24 microarrays were hybridized, scanned and analyzed using a common reference experimental design. With an arbitrary cut-off value of 1.5 fold change (FC) in gene expression, we identified the down-regulation of 212 genes (M value ≤−0.585) and the up-regulation of 158 genes (M value of ≥0.585) at a 1% false discovery rate in BLV-positive animals in comparison to the BLV-negative group. The maximum value of downregulation (FC = −3.21) was noted for two genes S100A4 (S100 Calcium-Binding Protein A4, formerly known as MTS1 or Metastasin) and CFD (Complement Factor D, formerly known as Adipsin). The most up-regulated gene (FC = 3.18) in blood cells from BLV-positive in comparison to BLV-negative cows was ADRA2A encoding Adrenoceptor Alpha 2A, a member of the G protein-coupled receptor superfamily. Lists of the top 25 down-and up-regulated genes in BLV-infected vs BLV-uninfected groups classified according to decreasing statistical significance of gene expression estimation (B-statistics) are presented in Tables 2 and 3 respectively. The complete lists of 212 and 158 DE genes are presented in the supplemental files (Tables 2S and 3S). Hierarchical clustering HCL heatmap for 370 differentially expressed genes and 24 samples (12 BLV-infected vs 12 BLV noninfected) analyzed in the study is presented in the supplementary Fig. 6S. Assuming gene expression fold changes FC ≥2 and FC ≥1.75 at p value (FDR-adjusted) <0.01 we were able to identify 33 and 70 up-regulated genes and almost twice as many down-regulated genes (64 and 130, respectively) in BLV-infected cattle in comparison to the control group. These values indicated that low (~1.5 fold) and medium (≥1.75 fold) gene expression changes were the most abundant classes, characteristic for BLVinduced progression to the lymphocytotic stage of the disease. The higher number of down-regulated genes than up-regulated ones may reflect the viral latency state, typical for the majority of BLV-infected cells in vivo and/or the evasion of the host's immune surveillance.

Validation of the microarray results by real-time qRT-PCR method
The accuracy of gene expression estimations from the microarray analysis was validated by the real-time qRT-PCR method. The analysis was performed on blood samples from the same 24 animals used in microarray analysis, although the technical replicate of RNA isolation was used for the reverse transcription reaction. A total of 14 genes were studied including 8 genes down-regulated according to microarray results (i.e.: CFD, ITCH, S100A, LGALS, TGFBI, CXCL8, DAP12, and C/EPB alpha, Table 2) and 4 up-regulated genes (i.e. MS4A1, HIF1A, MSH2 and ADAM9, Table 3) in the BLV-infected group in comparison to the BLVuninfected animals. Additionally, two genes with no differential expression between the studied groups were randomly selected, i.e. LTB encoding lymphotoxin beta (fold change, FC = −1.02, p value FDR adjusted = 0.929) and BNIPL encoding B-cell lymphoma (BCL2)/adenovirus E1B 19kD interacting protein (fold change, FC = 1.2, p value FDR adjusted = 0.36). For qRT-PCR data normalization, a pair of RPLP0 and UCHL5 gene was employed, based on stable expression in BLVinfected leucocytes which was determined previously [6]. The target amplicons information is presented in Table 4, and their specificity was confirmed by DNA sequencing (data not shown). The relative gene expression measurements between the BLV-infected and the BLV-uninfected groups were estimated with REST 2009 software (Fig. 1a) and indicated significant differences in gene expression levels in all analyzed genes except  Fig. 1b). The high value of the correlation coefficient holds out the hope that the rest of the gene expression measurements on the microarray are also valid.

Gene set enrichment analysis of DE gene lists generated from the microarray experiment
To improve the amount and the quality of the ontological information, a total of 370 DE genes identified with the BLOPlus microarray were updated in annotation terms and gene IDs using probe sequence information and BLAST analysis. The list of updated gene IDs and the expression data were input into online meta-analysis First, public gene ontology (GO) categories, i.e. GO Processes (biological process), GO Molecular Functions and GO Localizations (cellular component) were assessed, by that disclosing the biological motifs associated with the expression differences observed. The results indicated that the genes with differential expression (DE) could be classified to the numerous gene ontology terms, although with a substantial level of overlapping and redundancy (Fig. 2a,  c, d). The enrichment analysis reports, including the top-100 GO terms with calculated p value, FDR and the list of genes involved in each GO term, are presented in the supplementary MS Excel files (Tables 5S, 6S, 7S). Excluding some overlapping terms, we were able to select diverse biological process categories, significantly affected by BLVinfection and the host response, such as: immune system process, response to stress and wounding, regulation of immune response, cell activation, innate immune response, coagulation, regulation of programmed cell death, regulation of metabolic processes, regulation of cell communication, regulation of signal transduction, regulation of cell proliferation, and DNA repair. In addition, the more detailed GeneGo Process Networks Ontology, based on a manual curation database (MetaCore™ ver. 5.3, GeneGo, Thomson Reuters), further delineated the most affected processes and functional networks, among others: inflammation in particular NK-cell cytotoxicity, complement system, IL-2, TREM-1 and protein C signaling; cell adhesion and chemotaxies; immune response in a particular phagosome at antigen presentation, Th-17 derived cytokines, TCR signaling; apoptosis and signal transduction (TGFbeta, GDF, activin and leptin signaling) ( Fig. 2b and supplemental MS Excel Table 8S). Among GO molecular functions the most abundant category was represented by protein binding (Fig. 2c). On the other hand, GO Localizations enrichment analysis pointed to the importance of the extracellular region in BLV-induced pathogenesis (Fig. 2d).
Furthermore, the gene set enrichment analysis by GeneGo Disease Biomarkers Ontology linked the genes down-regulated in BLV-positive individuals with, among others: connective tissue diseases, numerous autoimmune diseases, viral diseases and lymphoma. In contrast, the genes with up-regulated expression in PL animals were connected with, as expected lymphoma, non-Hodgkin's lymphoma, numerous neoplasms, and, unexpectedly schizophrenia and mental disorders (Fig. 3a, b). Although we used in our comparisons animals similar to PL, which is thought to be a benign, preneoplastic stage of the disease, a total of 87 genes from our DE data set was linked with lymphoma and an even greater number of genes were linked with biomarkers of the neoplastic transformation of different organs. It may be suggested that deregulation of not a few but rather multiple genes established the transcriptional foundation of the initial transformation events induced by BLV towards a leukemogenesis. In addition, the observed enrichment of the genes associated with autoimmune diseases, sheds light on a new possible mechanism contributing to a failure of the host immune response against BLV.
We used the Analyzed Networks (AN) Transcription Regulation algorithm from MetaCore™ to identify the transcription factor (TF) pathways and networks that regulate the processes associated with BLV-induced pathogenesis. The AN Transcription Regulation algorithm revealed a list of 30 TFs potentially involved in the expression of differences observed in our study, including: CREB1, c-MYC, SP1, NF-κB, GCR-alpha, c-JUN, ETS1, C/EBP, p53, STAT1, STAT3 and HIF1A (supplementary MS Excel Table 9S). The three most connected factors: CREB1, c-MYC and SP1 regulate more than 50% of the genes in our DE data set. The CREB1 network (Fig. 4) illustrates CREB1 with its 159 targets (proteins regulated by CREB1) functioning in different areas of the cell. It should be emphasized that the BLV transactivator protein Tax binds and cooperates with the CREB transcription factor, and the former is one of the main targets of the Tax protein. Of note, transcription factors HIF1A (Fig. 5), C/EBP alpha, c-JUN (AP-1) which regulate the expression of many genes in our data, are simultaneously the targets for gene expression deregulation induced by BLV infection and the disease progression to the PL stage.

Discussion
In this study we used oligonucleotide bovine specific BLO-Plus microarrays and qRT-PCR analysis to examine gene expression changes in blood cells in response to retroviral BLV infection and disease progression to the lymphocytotic stage. Subsequently, the identified differentially expressed genes were used in the gene enrichment metaanalysis approach to reveal biological patterns associated with the expression differences observed. Before discussing the major items revealed by our study, we would like to emphasize the particular context of the obtained results. (1) The analysis of mRNA abundance between the naturally BLV-infected cattle and non-infected controls was performed on whole blood transcriptomes, using a dedicated tube system, which allowed for immediate cell lysis after   The various symbols used in the network have been described in detail in MetaCore Quick Reference Guide file publicly available at: https://portal.genego.com/help/MC_legend.pdf bleeding, and the stabilization of received mRNAs with RNA stabilization solution. In contrast, most of the gene expression data, analyzed in the context of BLV studies and reported previously [22,23,50], were generated with peripheral blood mononuclear cells (PBMC) separated by density gradient centrifugation. As we also collected blood samples into tubes with EDTA, we endeavored to isolate total RNA also from buffy coats and PBMCs prepared with the standard Ficoll-Hypaque method. Probably due to the prolonged time required for blood collection and transport from dairy farms to the laboratory, exceeding 24 h, these methods resulted in partially degraded RNA with highly variable indices of RNA quality (RIN values ranging from 3.2 to 7.8, data not shown), which was not suitable for reliable gene expression analysis. In contrast, Tempus™ Blood RNA Tube technology allowed for RNA isolation with yields sufficient for microarray analysis and the RNA quality parameters indicating negligible RNA degradation between the groups in comparison (Table 1). In addition to RNA degradation, the reliability of gene expression estimations of many immunologically important genes (including those mentioned in the context of BLV response, like TNFα, IL6, IL10, IL12, IFNγ, NFκB, IκB) is significantly affected by unintended gene regulation caused by phlebotomy [35], sample handling, [47], and uncontrolled activation of coagulation [39]. Furthermore, in vivo BLV is latent in the majority of infected cells, however, expression of BLV in samples of whole blood from BLVinfected animals is activated immediately upon incubation at 37 °C in the absence of any exogenous factors, except for anticoagulants or the removal the blood cells from plasma [13,46]. Consequently, an exaggerated immune response, based on differential mechanisms, could happen ex vivo in comparison to in vivo processes. Preservation of the blood transcriptomes in Tempus Blood RNA tubes or similar PAXgene tubes restricts ex vivo gene expression, allowing meaningful RNA assays and yielding transcript concentrations that are much closer to in vivo responses than can be obtained by other methods [39]. Nevertheless, the drawback of the whole blood transcriptome preservation approach is the fact that it is not possible at this time to distinguish which genes are deregulated because of their primal, direct participation in the host-pathogen interactions and which genes are affected as the result of secondary effects of BLV infection, such as the hematological profile alteration. Further investigations of selected genes reported in this study and involved in particular biological process of interest are needed to determine which molecular pathways result in BLV infection and disease progression and which are a result of BLV infection and disease progression. (2) The absence of clinically diagnosed lymphoma individuals that have been monitored in Poland in recent years and the comparatively high prevalence of animals with persistent lymphocytosis implied that a benign nonmalignant PL could be a convenient disease stage for the microarray analysis of host genetic factors that affect the outcome of BLV-infection. Lately, a strict BLV eradication policy has been implemented in Poland, and has resulted in the almost complete eradication of this previously highly endemic infectious agent. Newly diagnosed BLV seropositive individuals are eliminated from the herds immediately, and sometimes all animals within the herd (both seropositive and seronegative) are subjected to the compulsory slaughter. That is why we were able to collect blood samples from the BLV-positive individuals only once, so we could not examine the persistency of hematological aberrations and officially categorize them as the persistent lymphocytosis stage. However, the high total lymphocyte counts (mean 15.16 ± 10.28 × 10 3 /μl, Table 1), high total WBC counts, high titers of antibodies against p24 and gp51 and qualitative estimates of the proviral load (data not shown) [1,16,29] characterized all the animals within the BLV-positive group and indirectly indicated that all of them were PL.
The global gene expression data related to BLV-induced pathogenesis are relatively limited. To our knowledge, four research groups employed microarray technology to address potential gene expression alterations in host cells induced by BLV [4,12,22,50]. In general, the lists of differentially expressed genes reported in these studies coincide sparsely with each other. A similar trend is also apparent when we compare the DE genes identified in this study with the previous reports. There are only a few DE genes in common between the lists (e.g. MHC class II-DOB, HAT1 encoding histone acetyltransferase 1, APEX1 encoding multifunctional DNA repair enzyme, CCNG1 encoding cyclin G1), a few DE genes with invert expression (e.g. JUN oncogene, MAFB encoding v-maf musculoaponeurotic fibrosarcoma oncogene homolog B), and a number of differentially expressed genes encoding members of closely related gene families with overlapping molecular functions (e.g. ATF-activating transcription factor family of basic leucine zipper protein members, LMO-LIM domain only family members of transcription regulators, DUSP-dual specificity phosphatase family members). The relatively small coincidence between the lists of DE genes among the analyzed reports should not be surprising, taking into account that different experimental designs were applied in each case. In three studies, transformed and immortalized cell lines were used, including human epithelial HeLa cells [4] and ovine jejunal Peyer's patch-derived B-cell clone [22], which were subsequently transfected with plasmid vectors containing Tax BLV, in order to identify the spectrum of host genes regulated by Tax. In the latter, PBMC derived sheep B-cells cultured in vitro were also analyzed. In the third study [12] and its follow-up with siRNA knock-down of BLV Tax [34], the transcript profiling of a transformed but uninfected bovine lymphoblastoid cell line BL3 o was performed in comparison to its in vitro BLV-infected derivative cell line BL3*. However, no details about the identity of DE genes were presented in either report. On the other hand, the genes up-regulated in PBMC from 3 PL cows, cultured with purified plasma blocking factor, as compared with those cultured in medium were shown [50]. Furthermore, different microarray platforms were applied, including singlE−color human specific Affymetrix chips [4], twocolor human specific cross species microarray approaches [22,50], and the two-color bovine specific 7 K cDNA microarray [12]. In contrast, this study presents a comparison of the blood transcriptomes between naturally BLVinfected cattle and uninfected controls, captured immediately after bleeding without cell culturing and performed with the bovine specific microarray platform.
In our data we noted the interesting phenomenon that BLV infection and disease progression to PL is associated with higher numbers of down-regulated genes than upregulated ones (58% with FC ≥1.5 increasing to 65% when only FC ≥1.75 were considered). We hypothesize that, this may reflect the viral latency state, typical for the majority of BLV-infected cells in vivo and/or it may characterize BLV-induced immunosuppression, leading to the failure of an efficient immune response and underlying the host's evasion mechanism developed by BLV during evolution. The opposite result, with the dominance of up-regulated genes (exceeding 70-80% of total DE genes), was reported in in vitro cultured cells induced by the BLV Tax protein [4,22]. Such a redeployment towards gene induction is typical for many transcriptional profiles drawn from experiments utilizing in vitro cultures and purified factors in the study, dosed frequently with non-physiological ranges [9,19,50]. On the other hand, the results of the transcriptome profiling of blood cells from cattle naturally infected with pathogens like Mycobacterium bovis, Mycobacterium avium ssp. paratuberculosis, or Trypanosoma congolense, associated with chronic diseases, clearly demonstrated the importance of the gene repression mechanisms during the pathogenesis in vivo, which was reflected by the higher number of the down-regulated genes [21,31,33,51]. The range of gene expression changes identified in our data set, with a maximum down-and up-regulation fold change factor of around 3.2, indicated that BLV-induced progression to the PL stage of the disease was associated mainly with low (FC ~1.5) and moderate (FC ≥1.75) gene expression changes, with only rare examples of high levels of differential expression (FC ≥3). A similar picture of primarily small and moderate gene expression changes emerged from the analysis of ovine B-cells transfected with the Tax containing vector [22]. Likewise, among 122 genes deregulated by the wildtype Tax protein only 10 showed differential expression exceeding fourfold, with the rest of the gene expression changes within the range of 2-3 fold [4]. Apart from genes with high levels of differential expression, those exhibiting low differences could be of peculiar concern, taking into account that small changes in gene expression are supposed to be adequate to disturb cell homeostasis [22].
Our expression data support the importance of DNA repair processes during BLV-induced leukemogenesis. Philpott and Buehring reported that in the immortalized and transformed cell lines infected by the HTLV/BLV group of retroviruses [38], Tax inhibits basE−excision DNA repair of oxidative damage, thereby potentially inducing genomic instability and increasing the accumulation of mutations in cellular DNA. Of note, we found an increased expression of APEX1 (also known asAPE1/Ref1) mRNA in the BLV-infected group (+2.05-fold, p = 1.58E−07). APEX1 encodes a multifunctional protein engaged in the DNA base excision repair (BER) pathway of oxidative DNA damage (apurinic/apyrimidinic-endonuclease activity), as well as in transcriptional regulation. The latter is achieved as a redox coactivator of many transcription factors (among others: AP-1, Egr-1, NF-κB, p53, and HIF), and by its ability to indirectly bind to the negative calcium response elements (nCaRE) of some promoters acting as a transcriptional repressor [2,48]. The up-regulation of APEX1 mRNA was also reported in sheep B lymphocytes transfected with Tax-containing vector [22], and in tumor tissues and cancer cells of diverse origin [48], and its overexpression is associated with tumor cells' resistance to various anticancer drugs [44]. The necessity of APEX1 for cellular survival and its frequent overexpression in tumor cells strongly suggest a fundamental role of this protein in preventing cell death and in controlling cellular proliferation [48], probably through its regulatory role in the expression of cyclin-dependent kinase inhibitor p21 (CDKN1A) [44] and in the DNA binding activity of the p53 protein [7]. Furthermore, we found MSH2 (MutS homolog 2) encoding an important component of the post-replicative DNA mismatch repair system (MMR) as the most significantly upregulated gene in the BLV-infected group (+2.25-fold, p = 7.54E−09, Table 3). MutSα may also play a role in DNA homologous recombination repair, and may modulate cell cycle regulation and apoptosis. Defects of MSH2 are implicated in the development of lymphoblastic lymphomas in humans and mice, and are associated with the aberrant expression of the LMO2 (LIM domain only 2) gene [28], encoding a transcription factor which has a central and crucial role in hematopoietic development and neoplastic transformation. It is of interest that we also observed an increase in the expression of LMO2 gene (+2.16-fold, p = 4.43E−07) in BLV-infected animals in comparison to the BLV-negative control group. In vitro studies indicated that genetic instability associated with the downregulating of MutSα alpha is induced by the HIF-1A transcription 1 3 factor (Hypoxia inducible factor-1α) [24]. Hypoxia-oxygen deficiency in tissues is a key factor in the tumor microenvironment that has been tightly associated with tumor progression, metastasis and resistance to therapy [24] due to the increased rate of mutations in hypoxic conditions. We detected a significant over-expression of HIF-1A (+2.56fold, p = 1.58E−07) in the BLV-positive group and our meta-analysis identifies an HIF1A network (Fig. 5) with the HIF1A protein as a hub with multiple biochemical connections leading potentially to EBL associated pathogenesis. The activation and involvement of hypoxia-inducible factor 1 in human T-cell leukemia virus type 1 (HTLV-1) infected cell lines and primary adult T-cell leukemia (ATLL) cells was reported previously by [49]. It is of particular interest because HTLV-1 and BLV are characterized by similar genomic organization, similar strategies for gene expression, and similar pathologies. Figure 5 also shows direct and mutual protein interactions between HIF1A and AP1 (cJUN) and C/EBPα. It could be of interest that both CEBPA and c-JUN genes were downregulated in BLV infected cows in comparison to control group (−2.32-fold, p = 9.3E−06 and −1.85-fold, p = 5.53E−03, respectively), increasing the effect of hypoxic factors [25,43,53] in BLVassociated disease progression.
In addition, the results of our transcription profiling indicted some aspects of innate immunity which could be impaired in BLV-infected animals. We noticed the differences in the mRNA expression levels of genes involved in the activation of the complement system, NK-cell cytotixicity and the TREM-1 signaling pathway. Namely, in the BLV-infected group we identified a decrease in the level of mRNAs for C1qA and C1qC constituents of the complement subcomponent C1q (−2.37-fold, p = 3.41E−03 and −1.9, p = 1.49E−03, respectively). This is of particular interest, because it has been reported previously that, the interaction and binding of human complement subcomponent C1q with the gp-21 protein of HTLV-1 inhibits its infectivity [17]. Moreover, we also found in BLV-infected cattle an increased level of C1qBP (C1q-binding protein) mRNA (+1.99-fold, p = 3.03E−06), encoding a multifunctional protein known to inhibit complement component C1 activation by binding to globular heads of C1q molecules. C1qBP is also involved in cell proliferation, migration and apoptosis, and has been reported to be upregulated in many human cancers [30]. Finally, we observed a decreased expression of CFD (adipsin, complement factor D), CFB (complement factor B) and CFP (properdin, complement factor P) genes (−3.21-fold, p = 1.74E−10; −1.93-fold, p = 4.25E−03 and −1.66-fold, p = 3.8E−03, respectively), which were involved in the alternative pathway of the complement system activation. The role of the alternative pathway in the response against retroviral agents is unclear, but some clues indicate their contribution to leukemogenesis.
The CFD gene encodes a serine protease, which cleaves the complement factor B, releasing a small fragment Ba and a larger fragment Bb. The Ba fragment of complement factor B inhibits human B-lymphocyte proliferation [3], while the Bb fragment acts as B-cell growth factor [36] and is further complexed with C3b and properdin to form C5 convertase, participating in the complement cascade. Both adipsin and properdin genes were reported among 50 genes useful for leukemia class prediction, allowing the distinguishing of acute lymphoblastic leukemia (ALL) from acute myeloid leukemia (AML) patients [15].
The differences in leukocyte cell subpopulations found among the BLV-infected and non-infected cattle groups may be partially responsible for some of the identified gene expression changes. Nevertheless, the data showed in this report also reinforce the hypothesis that a host-or pathogen-or both-driven mechanism of innate immune gene repression might be of crucial importance for the outcome of BLV infection under natural conditions. It could be of interest that the cattle in the non-infected control group but exposed to BLV for at least 3 years in farm conditions were marked with the increased levels of mRNA transcripts associated with innate immunity. Disruption or deficiency of a suitable innate immune response may thus be crucial in BLV dissemination and progression to disease. Due to the paucity of information on innate antiviral response to BLV infection and disease progression, further studies are clearly needed.

Conclusion
Compared to earlier studies, we identified a significant number of new genes that have altered gene expression in BLV-infected blood cells. A higher number of down-regulated genes than up-regulated ones that have been noticed may reflect the viral latency state, typical for the majority of BLV-infected cells in vivo and/ or it may characterize BLV-induced immunosupression. The latter could lead to the failure of efficient immune response and underlie the host evasion mechanism developed by BLV during evolution. Taking into account their function, the differentially expressed genes like: MSH2, APEX1, HIF1A, LMO2, CEBPA, and ITCH cast a new light on the mechanism of leukemogenesis associated with BLV and should be of particular interest in further studies. Furthermore, we have identified several new response pathways important for BLV-induced pathogenesis involving innate immunity including for the first time the complement system and TREM-1 signaling pathway, that are deregulated in BLV-infected cells. It provides a key for better understanding of the role of innate immunity against retroviral agents and gives a promise for successful intervening against them. However it should be emphasized that there are numerous regulatory processes that are targeted by BLV-infection and the complex network of interrelated pathways is instead disturbed, causing the interruption of the control of B-cell proliferation and programmed cell death.