Expression Profiles and Interaction of MicroRNA and Transcripts in Response to Bovine Leukemia Virus Exposure

Bovine leukemia virus (BLV) infection in cattle is omnipresent, which causes significantly economical losses worldwide. The objective of this study was to determine microRNA (miRNA) and transcript profiles and to establish their relationship in response to exposure to the virus. Small noncoding and messenger RNA were extracted and sequenced from serum and white blood cells (WBCs) derived from seven BLV seropositive and seven seronegative cows. Transcriptomic profiles were generated by sequencing RNA libraries from WBC. Bta-miR-206 and bta-miR-133a-3p were differentially expressed in serum (P < 0.05). In WBC, bta-miR-335-3p, bta-miR-375, and bta-novel-miR76-3p were differentially expressed (P < 0.03). There were 64 differentially expressed transcripts (DETs). Gene ontology (GO) analysis of the DETs overexpressed in the seropositive group with GOs of response to stimulus and immune system process predicted that the DETs could potentially negatively regulate viral life cycle and viral entry or release from host cells. In addition, the DETs depleted in the seropositive group could play a role in the downregulation of antigen processing and presentation of endogenous peptide antigen via MHC class I. The differentially expressed miRNAs targeted 17 DETs, among which the expressions of bta-miR-133a-3p and bta-miR-335-3p were significantly negatively correlated with the expressions of ENSBTAT00000079143 and ENSBTAT00000066733, respectively. Under high prediction criteria, 90 targets of the differentially expressed miRNAs were all non-DETs. The most enriched biological process GO term of the targets was the RNA-dependent DNA biosynthetic process, which could be associated with virus replication. These results suggested that the differentially expressed miRNAs fine-tune most of the target genes in responding to BLV exposure. In addition, Bta-miR-206 interacted with BLV regulatory genes rex and tax by targeting their coding regions. A further study of the miRNAs and the genes may reveal the molecular mechanisms of BLV infection and uncover possible ways to prevent the infection.

Bovine leukemia virus (BLV) infection in cattle is omnipresent, which causes significantly economical losses worldwide. The objective of this study was to determine microRNA (miRNA) and transcript profiles and to establish their relationship in response to exposure to the virus. Small noncoding and messenger RNA were extracted and sequenced from serum and white blood cells (WBCs) derived from seven BLV seropositive and seven seronegative cows. Transcriptomic profiles were generated by sequencing RNA libraries from WBC. Bta-miR-206 and bta-miR-133a-3p were differentially expressed in serum (P < 0.05). In WBC, bta-miR-335-3p, bta-miR-375, and bta-novel-miR76-3p were differentially expressed (P < 0.03). There were 64 differentially expressed transcripts (DETs). Gene ontology (GO) analysis of the DETs overexpressed in the seropositive group with GOs of response to stimulus and immune system process predicted that the DETs could potentially negatively regulate viral life cycle and viral entry or release from host cells. In addition, the DETs depleted in the seropositive group could play a role in the downregulation of antigen processing and presentation of endogenous peptide antigen via MHC class I. The differentially expressed miRNAs targeted 17 DETs, among which the expressions of bta-miR-133a-3p and bta-miR-335-3p were significantly negatively correlated with the expressions of ENSBTAT00000079143 and ENSBTAT00000066733, respectively. Under high prediction criteria, 90 targets of the differentially expressed miRNAs were all non-DETs. The most enriched biological process GO term of the targets was the RNA-dependent DNA biosynthetic process, which could be associated with virus replication. These results suggested that the differentially expressed miRNAs fine-tune most of the target genes in responding to BLV exposure. In addition, Bta-miR-206 interacted with BLV regulatory genes rex and tax by targeting their coding regions. A further study of the miRNAs and the genes may reveal the molecular mechanisms of BLV infection and uncover possible ways to prevent the infection.
Bovine leukemia virus infection leads to significant economic losses by causing reduced milk production, premature culling, lymphosarcoma-related deaths, trade restrictions, and decreased resistance to infectious diseases in cattle (8,(10)(11)(12). Preventive and therapeutic strategies were developed against BLV infection including hematological, genomic, or serological methods to identify BLV-positive animals (13). Cattle with BLV are culled as well as corrective management and veterinary practices are indicated (13,14). There is currently no effective vaccine against BLV (12,13).
Attempting to find a solution to decrease the prevalence of BLV, an early study showed a considerable genetic influence on susceptibility to BLV infection (15). Recent studies on BLV-resistant cattle identified that DRB3 alleles correlate with a low proviral load profile (16)(17)(18). Based on the analysis of a US Holstein cattle population, the heritability estimate of leukosis incidence was approximately 8%, indicating that some genes may modulate leukosis (19). Furthermore, nine single nucleotide polymorphisms (SNPs) and the genes where these SNPs reside or neighbor have been associated with the incidence of BLV infection status (20). To further gain insight into the underlying mechanisms of the maintenance of latency, lymphoproliferation in the persistent lymphocytosis stage, and initiation of tumorogenesis, cattle cytokine imbalance has been analyzed during different stages of BLV infection. Results indicated that short-termed immune response of IL-12p40 and IFNγ followed by overexpression of IL-10 might modulate disease progression to persistent lymphocytosis (9,21). In addition, epigenetic modulation of BLV expression revealed that methylation in the viral promoter may be associated with latency (22).
MicroRNAs (miRNAs) are an abundant class of small, conserved, noncoding RNA molecules that posttranscriptionally regulate gene expression via base-pairing with complementary sequences within messenger RNAs (23). Their importance in development, physiology, and disease in many organisms was identified (23-26). Our initial analysis reported that certain miRNAs were potentially associated with BLV infection (27).
The objective of this study was to determine the association of miRNAs and expressed genes to establish putative candidate miRNAs that could potentially be used as biomarkers for the condition, based on the profiles of miRNA and transcript presented in BLV seropositive and seronegative samples.

Sample Test and Collection
Animal care and sample collection were performed according to the management protocol approved by the Institutional Animal Care and Use Committee of the National Animal Disease Center, in Ames, IA, United States.
The animals used in this study have been described previously (28). In brief, fourteen Holstein females were selected at the National Animal Disease Center, in Ames, IA, United States. Animals were randomly sampled from the dairy herd at the Center until seven positive and seven negative animals were identified. All animals included in this study were considered healthy according to the attending veterinarian. Samples came from 3 heifers and 11 cows with at least one calving and at midlactation. A jugular venipuncture with PAXgene tubes (PreAnalytiX GmbH, Hombrechtikon, Zurich, Switzerland) was used to obtain white blood cells from blood samples. Prior to RNA isolation, the tubes were incubated at room temperature for 2 h and placed in a refrigerator (4 • C) before being centrifuged at 3,000 g for 10 min to collect serum and white blood cells (WBCs). Sera were obtained from a second blood sample collected through jugular venipuncture in serum separator vacutainer tubes (SST, TM BD, Franklin Lakes, NJ, United States). The tubes were incubated at 37 • C for 30 min and then centrifuged at 1,250 g for 30 min. Isolated sera were stored at −80 • C until processed. Sera were used to establish IgG reactivity to BLV with a direct enzyme-linked immunosorbent assay (ELISA), using the IDEXX Leukosis Serum X2 Ab Kit (Idexx Laboratories, Westbrook, ME, United States). A sample was considered positive if the sampleto-positive ratio (S/P) was >115%. The S/P was determined by the ratio of the difference between the optical density 450 nm (OD450) of the sample minus the OD450 of the negative control, divided by the difference between the OD450 of the mean positive control minus the OD450 of the negative control. Based on the ELISA test, seven animals were identified as seronegative, and seven were seropositive.

RNA Isolation and Sequencing for Small RNA Analysis
The total RNA was extracted from serum, and WBC samples using the MagMAX TM mirVana TM Total RNA was eluted in 100 µl of RNase-free water. The concentration and quality of small RNAs in each sample were determined using a 10-40 nucleotide gate on an Agilent 2100 Bioanalyzer Small RNA chip (Agilent Technologies, Santa Clara, CA, United States).
The purified RNAs extracted from each sample were used to prepare individual libraries using the NEBNext Multiplex Small RNA Library Prep Kit (New England BioLabs, Ipswich, MA, United States). Library concentration and purification were performed using the QIAquick PCR purification kit (QIAGEN, Germantown, MD, United States). Each library was run on an Agilent 2100 Bioanalyzer High Sensitivity DNA chip (Agilent Technologies, Santa Clara, CA, United States) to determine the quality and quantity of the prepared library between 135 and 170 base pairs. Then, 30 ng of each library was pooled (14 libraries in the pool), and the size was selected using AMPure XP beads (Beckman Coulter, Indianapolis, IN, United States). Following the size selection, library pools were concentrated using the QIAquickPCR purification kit (QIAGEN, Germantown, MD, United States) and eluted in RNase-free water. The Agilent 2100 Bioanalyzer High Sensitivity DNA chip was used to determine the concentration of each library pool. The library pool was sequenced as single-end 50 base pair reads using the Illumina HiSeq 3000 System (Illumina, San Diego, CA, United States).

RNA Isolation and Sequencing for Transcriptome Analysis
The WBC RNA was extracted using the mirVana TM kit (Ambion, Carlsbad, CA). The quality of the mRNA was first checked using a Nanodrop (Thermo Fisher Scientific, Willington, DE). If samples had a measured λ260/λ280 > 1.9, they were run on Agilent 2000 bioanalyzer using the RNA 6000 Nano chip (Agilent, Santa Clara, CA). Samples with RIN # > 7.0 with 1-10 µg of total RNA were used to prepare libraries.
The libraries were constructed using the NEB Next Ultra RNA library kit for Illumina with the NEBNext PolyA Magnetic Isolation Module (NEB, Ipswich, MA). After final clean-up, 1 µl of the libraries were run on an Agilent 2000 bioanalyzer using the high sensitivity DNA chip. Based on the average size and concentration, the individual libraries were pooled to an equal molar concentration. The pools were sequenced using an Illumina HiSeq 3000 sequencer and were run at 2× 100 bp at the Iowa State University DNA Sequencing Facility (Ames, IA, USA).

miRNA Analysis
The bovine reference genome was downloaded from Ensembl Genes 97 (https://uswest.ensembl.org/info/data/ftp/index. html), and the miRNA precursor and mature sequences were downloaded from miRBase (http://www.mirbase.org, Release 22.1). FastQC was used to evaluate serum and WBC small noncoding RNA raw sequences (29), and then cutadapt was used to select high-quality sequences (30). The serum and WBC small noncoding RNA sequences were counted using miRDeep2 version 2.0.0.8 (31). The counted BLV miRNA raw reads were then used to identify differentially expressed miRNAs by running the DESeq2 package (32). The small RNA sequences are available on the NCBI SRA under BioProject accession number PRJNA378560.

Differentially Expressed Transcript
The WBC RNA sequencing reads were selected and evaluated in a similar way as small noncoding RNA, as previously described. Reads were mapped by STAR (33) against the bovine reference genome. Raw transcript sequences were sorted with samtools (34) and then counted with RSEM (35). Raw counts were analyzed using the DESeq2 package (32) to identify differentially expressed transcripts (DETs). Raw RNA-seq data were deposited in NCBI SRA under BioProject accession number PRJNA839936.

Target Gene Prediction
The bovine 5 ′ and 3 ′ UTR and coding sequences were obtained by running the R package "biomartr" (Version 0.9.0) (36). The BLV noncoding 5 ′ and 3 ′ UTRs and coding sequences were extracted from the BLV genome deposited in NCBI (https://www. ncbi.nlm.nih.gov/assembly/GCF_000853665.1). Target genes of the differentially expressed miRNAs were separately predicted by miRanda and PITA (37,38). The criteria to select target transcripts are maximum energy ≤-15 kcal/mol for miRanda and energetic score ≤-10 kcal/mol for PITA. The reported target genes were those that were predicted by both programs.

Statistical and Gene Functional Analyses
Gene set enrichment was tested using the Fisher's exact test. Gene ontology (GO) of interested transcript sets was analyzed using OmicsBox (www.biobam.com). Spearman correlation parameters between the differentially expressed miRNA and transcript across the 14 animals were estimated using R (3.6.1).

MicroRNA and Transcript Sequencing and Mapping
Small noncoding RNA sequencing from serum and WBC generated 568,377,432 clean reads ( Table 1). Among these, 143,755,444 small RNA sequences from serum were mapped to the bovine genome (ARS-UCD1.2), and 2,758,717 were mapped to mature bovine miRNA. From WBC, 362,041,990 sequences were mapped to the bovine genome, and 140,595,431 were mapped to mature bovine miRNA. Therefore, a higher number of small noncoding RNA sequences were derived from WBC than from serum. WBC transcriptome sequencing produced a total of 759,903,436 sequences, of which 598,380,902 were uniquely mapped to the bovine genome and 106,416,066 sequences were mapped to multiple bovine loci.

Differentially Expressed miRNAs
Small noncoding RNA sequences were mapped to 836 and 1,086 mature miRNAs in miRBase with at least one read in serum and WBC, respectively (Supplementary Table S1). Novel miRNAs were predicted by miRDeep2 with at least 200 reads presented in at least 6 samples (39). A total of 22 and 252 novel miRNAs were identified in serum and WBC, respectively (Supplementary Table S2). Among those miRNAs, the mature miRNAs bta-miR-206, bta-miR-133a-3p, bta-miR-335-3p, and bta-miR-37 and one novel miRNA bta-novel-miR76-3p were differentially expressed ( Table 2).

Differentially Expressed Messenger RNA
White blood cell RNA sequences were aligned to 22,509 bovine transcripts with a raw sequence count >10 for the 14 samples, among which 64 transcripts were differentially expressed (Figure 1). There were 43 DETs that were upregulated in the seropositive group and 21 DETs upregulated in the  In silico Functional Analysis of Differentially Expressed Messenger RNA Gene ontology analysis of the 43 DETs presented abundantly in the seropositive group revealed 17 associated biological processes at the second level of GOs (Figure 2). Transcripts ENSBTAT00000080234, ENSBTAT00000028059, ENSBTAT00000013838, and ENSBTAT00000083540 were associated with response to stimulus, and ENSBTAT00000022991 and ENSBTAT00000028059 were related to the immune system process. The specific biological processes of high-level GOs of those transcripts suggested that the animals might be actively responding to BLV infection through the regulation of gene expression. For example, the transcript ENSBTAT00000028059, encoded by a gene named TRIM8, is associated with both response to stimulus and immune system process. The specific biological processes of the transcript described by its highlevel GOs include not only negative regulation of viral release and entry into the host cell, defense response, and innate immune response but also negative regulation of viral gene expression, viral process, viral transcription, and viral life cycle (Supplementary Table S4). The 21 upregulated transcripts in the seronegative group share the same GO terms at the second level as the 43 transcripts from the seropositive group (Figure 2); however, the specificities of the GO terms are different from that of the 43 transcripts. For example, ENSBTAT00000048998, encoded by gene TAP2 with the 2nd level GO of the immune system process, has high-level GOs without negative regulating viral-related biological processes, but involving antigen processing and presentation, transport, antigen processing and presentation of endogenous peptide antigen, antigen processing and presentation of peptide antigen via MHC class I, as well as organic substance, nitrogen compound, amide, and peptide transport (Supplementary Table S5).

miRNA Target Analysis
With the mapped 22,509 cow transcripts as query, miRanda and PITA predicted 5,741 targets by the differentially expressed miRNAs (37,38), seventeen of which were the DETs identified in this study (Supplementary Table S6). Among the 17 DETs, ENSBTAT00000075709 and ENSBTAT00000076256 were targeted at their 3 ′ untranslated regions (UTR); the remaining DETs were targeted in the coding region by the miRNAs. The coding regions of GCC2-203, CENPF-201, KTN1-202, and OSBPL8-204 were separately targeted by two of the three differentially expressed miRNAs bta-miR-133a-3p, bta-miR-335-3p, and bta-miR-206 (Supplementary Table S6). Nine and eight targeted DETs were expressed abundantly in the seropositive and seronegative groups, respectively.
There are 5,724 of the target transcripts, which are non-DETs. When the criteria were increased to ≤-20 kcal/mol for miRanda and PITA, the number of non-DET targets decreased to 90 (Supplementary Table S7). The enrichment analysis of the 90 transcripts showed that all the significant GOs were overrepresented. RNA-dependent DNA biosynthetic process and RNA-directed DNA polymerase activity are the most significantly enriched GOs in the biological process and the molecular function categories, respectively ( Table 3).
The BLV genome consists of 8,714 nucleotides and encodes several structural and regulatory protein genes (40,41). Both miRanda and PITA predicted that the BLV regulatory genes rex and tax were targeted by bta-miR-206 (Table 4).

Correlation Between Differentially Expressed miRNA and Their Targets
Spearman correlations were estimated between differentially expressed miRNAs and the target DETs. Among the 21 pairs of miRNA and transcript predicted by miRanda and PITA (Supplementary Table S6), the expressions of bta-miR-133a-3p and bta-miR-335-3p were significantly negatively correlated with the expressions of ENSBTAT00000079143 and ENSBTAT00000066733, respectively (P < 0.05) ( Table 5), and the others presented trends without being significant. Both ENSBTAT00000079143 and ENSBTAT00000066733 were abundantly presented in the seropositive group.

DISCUSSION
The small RNA sequencing of the libraries identified two differentially expressed miRNAs in serum and three miRNAs in WBC. Those miRNAs are not overlapped, indicating that the cellfree circulating miRNAs in serum were not related to miRNAs in WBCs in the collected blood samples, but instead may be selectively released from different types of cells under varied pathological statuses (42). Among the five differentially expressed miRNAs, except for bta-novel-miR76-3p which have not been characterized, the other 4 miRNAs are all related to cancer.
The miRNA-206 expression has been reported to be significantly higher in human breast cancer tissue than that in paracancerous tissue (43). The human miRNA-206 inhibits both invasion and stemness of breast cancer and epithelial ovarian cancer (44,45). In this study, the expression of bta-miR-206 was observed to be relatively higher in seronegative serum than in seropositive serum, which may have similar functions as in human. In silico prediction of this miRNA targets to BLV genome identified that it interacted with BLV regulatory gene rex and tax by targeting their coding regions, which provided a possible explanation of how this miRNA regulates BLV invasion. Similar to miRNA-206, miRNA-133a was also a tumor suppressor in breast cancer (46). Consistent with the human study, bta-miR-133a-3p shared a similar expression pattern as bta-miR-206 in this study. One of the target DETs of the miRNA, ENSBTAT00000028059, encoded by TRIM8 sowed negative regulation of viral life cycle, as well as viral entry and release from host cells by GO analysis. This result may deliver further insight into the molecular mechanisms of cancer-suppressing functions of bta-miR-133a. The differential expressions of miR-206 and miR-133a were also reported in swine which was related to the infection of Aujeszky's disease virus (47).
MicroRNA-335 was reported to be functionally the same as miRNA-206, which inhibits the migration of breast cancer through different mechanisms (43,48). The expression of bta-miR-335-3p in WBC was significantly high in the seronegative animals and similar to bta-miR-206 and bta-miR-133a-3p in serum. miR-375 has been implicated in several types of cancers, including cervical cancer (49), skin cancer (50), and liver cancer (51). The miR-375 expression was high in Merkel Cell Carcinoma (MCC) cell lines and tissues compared with non-MCCs (50). A functional assay indicated that miR-375 can inhibit the proliferation and invasion of hepatocellular carcinoma cells (51). In this study, the expressions of bta-miR-375, as well as bta-novel-miR76-3p, were higher in seropositive WBC than that in seronegative WBC. In seropositive WBC, the expression of miRNA-375 was significantly high, while the expressions of miRNA-206 and miRNA-335 were significantly low, but these miRNAs were all reported to have cancer inhibition functions, which might indicate the dual functional roles of miRNAs (52)(53)(54).
From the normalized transcript counts, we observed a clear pattern between seropositive and seronegative groups (Figure 1). However, many transcripts showed that some individuals have similar reads in both groups.  This may be due to specific characteristics of BLV infection; seronegative animals may also carry the virus but are at an early stage of infection. The seropositive animals might be at later phases of infection. Therefore, it is important to develop a method to accurately identify seropositive or seronegative animals at different infection stages. Regulation of gene expressions in response to BLV infection was observed in our analysis. In silico functional analysis of the abundant DETs in the seropositive group showed that the expressions of some genes might encounter the viral infection through negative regulation of the viral life cycle and interaction with host cells (Supplementary Table S4). In contrast, the abundant DETs in the seronegative group suggested that some genes of the animals fight viral infection through antigen processing and presentation of endogenous peptide antigens (Supplementary Table S5); therefore, the animals probably initiated their different defense strategies under the different stages of viral infection by altering their immune-related gene expression (55).
MicroRNAs as key regulators in growth and development have been demonstrated in animals and human (56)(57)(58)(59). In response to pathogen infection, miRNAs fine-tune various physiological and pathological processes (60)(61)(62). The complex roles of miRNAs are not only to suppress gene expression but also to activate gene expression in particular contexts (63)(64)(65). Our Spearman correlation analysis results in two significantly negative correlations among the expression values of 21 pairs of miRNA and target transcript, which implicated a big impact of some miRNAs on BLV infection as reported in other organisms (57,58,66). In contrast, the majority of the correlations between miRNAs and their targets were not significant, which might indicate miRNA's small effects on a large number of target genes (67). It is also possible that different BLV infection stages of the studied animals may alter the effects of the miRNAs. Furthermore, our results suggested that the miRNAs were not always significantly associated with their targets, and in silico predictive programs may not guarantee to identify significantly associated target genes.
The differentially expressed miRNAs targeted 3 ′ UTR, 5 ′ UTR, and coding regions of the cow transcripts in this study. Only 17 of the target genes were differentially expressed, but 5,724 targets were not differentially expressed. This result demonstrated the fine-tune characteristics of the miRNAs to their targets. Gene set enrichment analysis of the 90 target transcripts under high selection criteria for miRanda and PITA indicated that the significantly enriched GOs were all overrepresented ( Table 3). The GOs of RNA-dependent DNA biosynthetic process and RNA-directed DNA polymerase activity are the most significantly enriched in the biological process and the molecular function categories, respectively ( Table 3). These results may indicate the augmentation of the RNA virus in infected cattle. It has been reported that Epstein Barr Virus encoded EBNA1 could hinder the presentation of MHC I-associated peptides (68,69). Similarly, our results suggested that BLV has employed some strategies to negatively regulate antigen processing and presentation of peptide antigen via MHC class I ( Table 3); therefore, it might overcome the recognition of the infected cell and enable viral propagation in the host. The genes in the related processes were not significantly expressed but targeted by some of the differentially expressed miRNAs.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found online repositories. The names of the repository/repositories and accession number(s) can be found below: https://www.ncbi.nlm. nih.gov/, PRJNA378560 and PRJNA839936.

ETHICS STATEMENT
The animal study was reviewed and approved by Institutional Animal Care and Use Committee, National Animal Disease Center.

AUTHOR CONTRIBUTIONS
EC and JL conceived and designed the experiment. HM analyzed the data and wrote the manuscript. HM and EC interpreted the results. HM, JL, and EC reviewed and edited the manuscript. All authors contributed to the article and approved the submitted version.

FUNDING
This research was part of an intramural project of the USDA/ARS National Animal Disease Center. The USDA had no role in study design, data collection, analysis, interpretation of results, or preparation of the manuscript.