Transcriptomic analysis of bovine monocytes in response to non-cytopathic bovine viral diarrhea virus infection

Background: Monocytes are signicant players in the detection of invading pathogens, particularly in pathogen defense. Bovine Viral Diarrhea Virus (BVDV) can cause a persistent infection and immune suppression if animals are infected with an non-cytopathic (ncp) biotype. However, its exact role in ncp BVDV-infected bovine monocytes remains poorly understood. To explore the immune suppression mechanisms of ncp BVDV, we used a transcriptomics approach to nd genes with differential expression patterns in monocytes during infection with ncp BVDV over time. Results: Bovine monocytes were sampled at 2 and 24 h post-infection (hpi) to represent the early and late stages of an ncp biotype strain of bovine viral diarrhea virus infection. Compared with the non-infected cells, 9959 and 7977 differentially expressed gene (DEGs) were identied at 2 and 24 h hpi, respectively. These DEGs were associated with signal transduction, immune response, apoptotic process, cellular process , binding and cellular component. The differential expression proles of select the type I interferon signaling pathway , interferon (IFN)-stimulated genes (ISGs), and genes involved in the innate immune response, including IRF7, DDX3X, TLR13, DDX58(RIG-I), MVAS, TLR9, TRAF6, IRF1, IFIT1, STAT1, ISG20, TRIM25, MX1,NLRX1, CYLD, SIKE1 and ZAP70 were conrmed by real-time quantitative PCR and consistent with the RNA-seq data. These results indicated that infection with ncp BVDV could activate type I interferon signaling pathway in bovine monocytes and induces weak ISGs responses, which extends our present understanding how the virus modulates the immune response and leads to better understanding behind the immunopathogenesis of ncp BVDV. Conclusion: Our transciptome anslysis provides useful initial data towards better understanding of the infection mechanisms used by ncp BVDV, while highlighting the potential molecular relationships occurring between the virus and the host’s immune response. TRIF: cGAS: Cyclic GMP-AMP synthase; STING: Stimulator of interferon genes ; MAVS: Mitochondrial antiviral-signaling; PBMC: peripheral blood mononuclear cells; TCID50: The 50% tissue culture infectious dose; mock: Uninfected control; FPKM: fragments per kilobase per million mapped reads; DEG: differentially expressed gene; GO: Gene ontology; NF-κB: nuclear factor kappa-B; GAPDH: glyceraldehyde–3-phosphate dehydrogenase; qRT-PCR: real-time quantitative reverse transcription-polymerase chain reaction;

variety of pattern recognition receptors (PRRs). Currently identi ed PRRs can be divided into the Toll-like receptors (TLRs), the retinoic acid-inducible gene I (RIG-I) receptor-like receptors (RLRs), and the NOD-like receptors (NLR). All three PRRSs recognize the viral PAMPs [7].Each of the three receptors cooperates with their adaptor proteins to deliver information in the presence of viral PAMPs and interferons. The Tolllike receptors affect the Total Recordable Injury Frequency (TRIF) adaptor protein, the Cyclic GMP-AMP synthase (cGAS) receptor effects the Stimulator of interferon genes (STING) adaptor protein, and the RIG-I receptor interacts with the Mitochondrial antiviral-signalingMAVS) adaptor protein. These connecting molecules recruit and activate downstream signaling pathways to produce corresponding in ammatory cytokines, chemokines, and Type I (IFN-α/β) and III (IFN-) interferons [8,9]. Much research has demonstrated that evading innate immunity is vital to viral persistence because it evades host immunity [10].
Previous studies have revealed many crucial discrepancies in the results of ncp and cp biotypes in variety of cell types. ncp BVDV does not induce IFN synthesis in cultured cells in vitro, suggesting that this could be a defense mechanism to evade innate host immunity that might be critical for establishing persistent infections [11,12,13]. However, the molecular mechanisms by which BVDV evades host cell immunity remains unclear. RNA-seq has recently been developed for transcriptomics analysis of BVDV infection in MDBK cells. Research was performed with microarrays that showed gene expression changes in host cells infected with ncp BVDV compared with cp BVDV [14,15]. However, it was reported that a large number of genes related to BVDV invasion and replication mechanisms still need to be characterized to better understand the interaction between BVDV pathogenetic mechanisms and host immune defenses.
Although ncp BVDV is thought to play an important role in BVDV-infected bovine monocytes [16][17][18][19], little is known about how BVDV disturbs host cell gene expression pro les. During this research, we showed that BVDV infected bovine monocyte gene expression pro les, which provided meaningful information from genetic analyses that could not be discerned using biochemical methods.

Animals
Nine conventionally reared, healthy BVDV-free cows from the Shihezi University Animal Experimental Station were used. All animal experiments were approved by the Shihezi State University Institutional Animal Care and Use Committee. All animals were con rmed to be free of BVDV, infectious bovine rhinotracheitis virus, bovine parain uenza virus, and Mycoplasma bovis infection using enzyme linked immunosorbent assay (ELISA) kits (Idexx Labs Inc., USA) and by reverse transcription PCR (RT-PCR) or PCR for viral nucleic acid detection [18]. As we expected, all animals were BVDV mRNA negative.

BVDV stocks and infection
The ncp BVDV strain was named BVDV-LC (GenBank accession NO. MK102095, identi ed and kept in our lab), and was classi ed as genotype 1q and non-cytopathic (ncp). For titration of virus stock, four replicates of 10-fold serially diluted virus (starting from 1/10) were inoculated on MDBK cell monolayer in 96-well culture plates. After 48 h incubation, the culture plates were xed at 4 °C for 30 min with ice cold absolute ethyl alcohol and subjected to immuno uorescence staining with uorescein isothiocyanate (FITC)-conjugated polyclonal anti-BVDV (VMRD, USA) antibody (diluted 1:1000). The uorescence signals were observed undera uorescence microscopy (ZEISS) and viral titers was expressed as the 50% tissue culture infective dose (TCID50)/mL by Reed-Muench method.
For infection of bovine monocytes, virus dilutions were made in DMEM 4Mm L-glutamine, 4.5 g/l glucose, 1.5 g/l sodium bicarbonate and 10% horse serum. 5×10 6 monocytes were added to each well of a 6 well tissue culture plate and adsorbed with ncp BVDV biotypes at the same MOI of 0.1 for 2 h. After infection 2 h and 24 h, at least 10 7 cells were pooled in one tube. Uninfected control (mock) cells were treated similarly but were dosed with media and no viral inoculum.All data were determined using triplicate monocyte cultures.

RNA extraction
Cells were harvested and prepared for RNA extraction 2 hours and 24 hours after infection. The monocytes were rinsed three times with ice-cold PBS, and whole RNA was extracted by DNase digestion using the RNeasy mini kit following the manufacturer's instructions (Qiagen). RNA was quanti ed using a NanoDrop spectrophotometer (Thermo Fisher) and stored at -80℃. For RNA-seq analysis, the RNA quality was assessed with an Agilent Bioanalyzer (Agilent Technologies).

RNA sequencing and analysis
The RNA elimination was performed using Ribo-Zero rRNA removal kit (Epicenter, USA) following the manufacturer's instructions. Illumina sequencing libraries were prepared using the TruSeq Stranded mRNA Sample Preparation Kit (Illumina, USA) according to the manufacturer's protocol. RNA sequencing was performed on the Illumina HiSeq 4000 platform at the Novogene Bioinformatics Institute (Beijing, China).
The sequence data of the reference genome was checked against the NCBI database. The ltered mass reads were aligned with the reference genomic sequence using Bowtie 2. Relative transcript abundance was metered as fragments per kilobase per million mapped reads (FPKM). Genes were clustered into functionally relevant groups using the eggNOG (Genetic Evolution Genealogy: Unsupervised Orthologous Group) database and the metabolic pathways were analyzed using the KEGG (Kyoto Encyclopedia of Genes and Genomes) database. Visualization of mapping outcomes and differentially expressed gene (DEG) analysis was performed using the CLRNASeqTM program (ChunLab, Korea).
Infection and control time points were performed in triplicate using puri ed monocytes, and three biological replicates were generated. ncp BVDV expression and infected monocyte gene expression were contrasted with gene expression from the inoculum and uninfected monocytes. The fold change in gene expression was calculated from three different sample types.
Real-Time Quantitative Reverse Transcription PCR (qRT-PCR) cDNA was synthesized from 2 μg of total RNA using a reverse transcription kit (Fisher Scienti c, USA) and by following the manufacturer's instructions. Gene expression was quanti ed using a LightCycler 480 PCR platform (Roche Applied Science, Indianapolis, IN, USA). The qRT-PCR procedure began with an initial step of 10 minutes at 95 °C followed by 45 cycles at 95 °C for 30 seconds, 58 °C for 30 seconds, and 72 °C for 30 seconds. Gene transcription levels were quanti ed relative to GAPDH gene expression using the relative cycle threshold (ΔCT) method. Primers were made using Primer 5.0 software (Table 1).

Statistical analysis
All data are shown as means ± SE. Differences were evaluated using ANOVA followed by the Student's ttest. Statistical signi cance was de ned as P < 0.05.

RNA sequence analysis of ncp BVDV-infected monocytes
The results show that more differential expression, either up-or down-regulated was seen at 2 h postinfection ( hpi), the expression of 9959 genes was signi cantly changed compared with those of the controls, of which 4968 genes were upregulated, and 4991 genes were down-regulated (Fig.1A, Additional le 1: Table S1). At 24 hpi, the expression of 7977 genes was signi cantly different from that of the controls, with 4184 genes up-regulated and 3793 genes down-regulated (Fig.1A, Additional le 2: Table  S2). These signi cantly altered genes were involved in signal transduction, immune response, apoptotic process, cellular process, binding and cellular component, with most of the differential genes being involved in signal transduction.
The DEGs were ltered to determine which genes were present at both of the time points (2 hpi and 24 hpi). At 2 hpi time point, 9959 genes responded to ncp BVDV infection, and at 24 hpi time point, 7977 genes responded; however, only 5709 genes were common at both time points (Figure. 1A). Two libraries were packaged and analyzed against the non-redundant NCBI database using the BLAST program. Next, we mapped and annotated the genetic libraries using the Gene ontology (GO) database. To identify pathways involved in immune activation after ncp BVDV infection of bovine monocytes, transcript records were further analyzed. The GO terms classify the function of BVDV-infected cell transcripts, offering 4014 transcripts involved in molecular functions, 4476 involved in biological processes, both and 3341 involved in discretely cellular processes ( Figure. 1B). Analysis revealed that differentially expressed genes are involved in a variety of biological processes, including immune responses and immune system procedures, cell death, macromolecular localization, apoptosis, cell death regulators, antigen processing and presentation, protein localization, and others ( Figure. 1C, Additional le 3: Table S3, Additional le 4: Table S4).
Through careful examination of the information set ( Figure. 2B and Table 2), several ISGs [20] mainly produced by type I interferons, were differentially expressed in ncp BVDV-infected cells. Figure.2A and B show strati ed heat maps of all of the up or down regulated genes in ncp BVDV-infected cells, and they were all associated with immune responses and type 1 IFN signaling. Overall, the RNA-seq analysis indicated that ncp BVDV infection resulted in different expression of genes related to inherited immunity, type I IFNs, stimulatory cytokines, and cell signaling pathways, leading to the creation of antiviral responses in monocytes.
Partial validation of the DEGs involved in the type 1 IFN pathway by qRT-PCR Sequences found using RNA-seq analysis were validated with qRT-PCR to assess if mRNA expression was upregulated or downregulated in relation to the genes and to the immune responses. The genes were identi ed if they were recognized as being involved in type I IFN responses, were reported to be related to innate antiviral immune responses [21,22], and if they were common to both infected groups. Using GAPDH as an internal control, the transcription levels of 31 DEGs were measured (between the 2 and 24 hour time points) with qRT-PCR with the primers listed in Table 1. In conclusion, the RNA-seq and qRT-PCR results from bovine monocytes at 2 and 24 hpi were fairly consistent ( Figure.

Discussion
The complex and unique nature of BVDV continues to challenge infectious disease researchers,veterinarians,and the cattle industry. The different genotypes and biotypes of BVDV ability to induce persistently infection, as well as its ability to interfere both innate and adaptive immunity of the host, make it di cult for prevention and control [5]. Although Mais[18] evaluated the effect of cp and ncp BVDV infection of bovine monocytes to determine their role in viral immune suppression and uncontrolled in ammation by using proteomics, this study was the rst to analyze transcriptionally diversi ed gene expression after in vitro ncp BVDV infection of bovine monocytes. In this study, using RNA-seq, we demonstrated that most of the DEGs between ncp BVDV infected monocytes for 2 and 24 hours, involved immune responses, suggesting that ncp BVDV infection plays an important role in bovine monocytes. We observed increased levels in a variety of antiviral genes such as IFI27, IFI30, MX1, TRIM25, DDX41, ISG20, STAT1, TRAF6, CD14, CD40, IL1B, IL6, MAPK12, MAPK13, TNFSF13B, and TBKBP1 (Table 2), It seems to play a fundamental role in ncp BVDV infection in monocytes. The expression levels of the IFI27, IFI30, ISG20, STAT1, TRAF6 and TRIM25 genes obtained from the RNA-seq data were con rmed by the realtime PCR mRNA expression data (Table 1 and Figure. 3). These proteins might be upregulated in bovine monocytes because of the ncp BVDV infection. However, we observed that type I IFN genes were reduced such as IRF7, DDX3X, DDX58, TLR13, TLR9, STING, and IFIT1, which was possibly caused by the ncp BVDV viral particles [12,13]. On the other hand, upregulated DEGs, such as NLRX1, CYLD, SIKE1, ATG12, and RNF125, were shown to interfere with host immune responses in previous documents [23][24][25][26].
The DEAD-box proteins are characterized by the conserved sequence Asp-Glu-Ala-Asp (DEAD) is a presumptive RNA unwindase involved in many cellular processes, such as RNA combination and RNA secondary framework changes. This gene codes proteins involved in the RNA helicase-DEAD box protein motif and the caspase employment domain (CARD). It is related to the viral double-stranded (ds) RNA identi cation and coordination of immune responses [27,28,29].
Our research showed that DHX16, DHX34, DHX37, DHX58 (Table 2) were not signi cantly DEGs after ncp BVDV infection, but DHX58 expression was increased ( Figure. 3). In other studies, viral dsRNA stimulation did not induce high expression of DEAD-box proteins, and viral dsRNA was shown to be degraded by some proteases, such as BVDV Erns protein [30,31]. We found that DHX16, DHX34, DHX37, DHX58 DDX3X, DDX10, DDX51, DDX54, and DDX55 were differentially expressed after the ncp BVDV infection of monocytes. Post-downregulation was not signi cant ( Table 2), but DHX41 expression was signi cantly increased, suggesting that DDX41 plays an important role in RNA binding and secondary structure changes [32,33].
Toll-like receptors (TLRs) are one of the most essential innate immune receptors that identify many types of PAMPs from different pathogens and activate downstream immune responses. TLRs are grouped into two groups: TLR1, TLR2, TLR4, TLR5, TLR6, and TLR11, which are concerned with membrane proteins that detect microbial membrane components, while TLR3, TLR7, TLR8, and TLR9 are concerned with nucleic acid recognition of proteins, cells, intracytoplasmic molecules, and microorganisms [34,35]. Nevertheless, studies showed that infections and administrative mechanisms restrict the TLR13 signaling pathway in mollusks. In this study, TLR8, TLR9, and TLR13 expression in BVDV infected bovine monocytes were downregulated (Table 2 and Figure. 3), which may be the reason for the decreased type 1 IFN expression that was seen [36][37][38][39][40].
This study found common groups of genes related to important cellular processes at 2 hpi. For example, inhibition of NF-κB activation was induced by upregulated IKKE and SIKE1, negative regulators of the NF-κB signaling pathway, and by downregulated p65/RelA-NF-κB [26]. These inhibitory actions could lead to blocked or delayed antiviral responses (Table 2). These results are interesting because they show that infected living cells exhibit strong immunosuppressive phenotypes. In fact, the results also show downregulation of vital immune response genes, for example, the PIDD1, IKBKG, IKBKB, INFGR2, PRKCQ and EIF2AK2 genes (Table 2). Moreover, the expression of IKKA and SIKE1 transcripts indicates the blocking of NF-κB in infected cells, which were support by the qRT-PCR results ( Figure. 3 and Table 2).
Our research demonstrates that there is antiviral regulation into bovine monocytes upon activation, as shown by upregulation of the NLRX1, CYLD, SIKE1, ATG12, and RNF125 genes (Table 2 and Figure. 3). These ndings are vital for the understanding of the molecular mechanisms by which ncp BVDV evades innate immunity. NLRX1 is involved in antiviral signaling. By inhibiting the virus-induced RLH (RIG-like helicase)-MAVS interaction, NLRX is a negative regulator of the MAVS-mediated antiviral response [23]. RNF125 mediates DDX58/RIG-I ubiquitination inducing DDX58/RIG-I degradation of Lys-181 [24]. SIKE1, IKK-epsilon, and TBK1 genes are physiological inhibitors that act on viral and TLR3-triggered IRF3 to inhibit TLR3-mediated interferon-stimulated response elements (ISRE) and IFN-β promoters. Activation of the ISRE and IFN-β promoters can be achieved by disrupting the interaction between IKBKE or TBK1 and TICAM1/TRIF, IRF3, and DDX58/RIG-I, which does not inhibit the NF-κB activation pathway [26].
Type I IFNs are key mediators of host immune responses to viral infections. These IFNs induce the expression of many ISGs of which some could be antiviral. Besides, they regulate innate and adaptive immunity by activating immature DCs, strengthening NK cell responses, and boosting the survival and effector functions of T and B lymphocytes [41]. In this research, a signi cant decrease in ZAP70, TNFSF11, V-TCR, and CARD11 was regulated in ncp BVDV-infected bovine monocytes. ZAP70 regulates T cell activation by modulating TCR expression on T cell surfaces [42,43]. CARD11 is involved in the costimulatory signals necessary for T cell receptor (TCR)-mediated T cell activation. CARD11 connects with the DPP4 and guides T cell multiplication and NF-κB activation in a T cell receptor/CD3-dependent manner [44]. TNFSF11 is a necessary regulator of T cell and dendritic cell interactions and might regulate T cell-dependent immune responses [45]. Most serine/threonine kinase-expressing genes are increased upon ncp BVDV infection of monocytes, which is quite interesting because only MAP2K6 is signi cantly reduced. MAP2K6 plays a vital role in regulating cytokine to cytokine responses, particularly as it relates to the induction of IL-6 production. IL-6 is increased 2 hpi and is decreased by 24 hpi (Table 2 and Figure. 2) [46, 47].