Transcriptome analysis reveals differential immune related genes expression in bovine viral diarrhea virus-2 infected goat peripheral blood mononuclear cells (PBMCs)

Bovine viral diarrhea virus (BVDV) is an economically important viral pathogen of domestic and wild ruminants. Apart from cattle, small ruminants (goats and sheep) are also the susceptible hosts for BVDV. BVDV infection could interfere both of the innate and adaptive immunity of the host, while the genes and mechanisms responsible for these effects have not yet been fully understood. Peripheral blood mononuclear cells (PBMCs) play a pivotal role in the immune responses to viral infection, and these cells were the target of BVDV infection. In the present study, the transcriptome of goat peripheral blood mononuclear cells (PBMCs) infected with BVDV-2 was explored by using RNA-Seq technology. Goat PBMCs were successfully infected by BVDV-2, as determined by RT-PCR and quantitative real-time RT-PCR (qRT-PCR). RNA-Seq analysis results at 12 h post-infection (hpi) revealed 499 differentially expressed genes (DEGs, fold-change ≥ ± 2, p < 0.05) between infected and mock-infected PBMCs. Of these genes, 97 were up-regulated and the remaining 352 genes were down-regulated. The identified DEGs were found to be significantly enriched for locomotion/ localization, immune response, inflammatory response, defense response, regulation of cytokine production, etc., under GO enrichment analysis. Cytokine-cytokine receptor interaction, TNF signaling pathway, chemokine signaling pathway, etc., were found to be significantly enriched in KEGG pathway database. Protein-protein interaction (PPI) network analysis indicated most of the DEGs related to innate or adaptive immune responses, inflammatory response, and cytokine/chemokine-mediated signaling pathway. TNF, IL-6, IL-10, IL-12B, GM-CSF, ICAM1, EDN1, CCL5, CCL20, CXCL10, CCL2, MAPK11, MAPK13, CSF1R and LRRK1 were located in the core of the network and highly connected with other DGEs. BVDV-2 infection of goat PBMCs causes the transcription changes of a series of DEGs related to host immune responses, including inflammation, defense response, cell locomotion, cytokine/chemokine-mediated signaling, etc. The results will be useful for exploring and further understanding the host responses to BVDV-2 infection in goats.


Background
Bovine viral diarrhea virus (BVDV) is the prototypic member of the genus Pestivirus in the family Flaviviridae, and two main different BVDV species have been recognized as BVDV-1 and BVDV-2 [1]. BVDV infection decreases productive performance and causes considerable economic losses in cattle industry worldwide [2]. Infections with both species of BVDV can induce similar diseases, from subclinical infections to severe clinical diseases including acute diarrhea, respiratory diseases, reproductive failures, congenital defects, and increased mortality due to immunosuppression [2][3][4]. Persistent infection is the common type of infection in cattle and the persistently infected (PI) animals are considered the main source of BVDV transmission [2,5]. PI animals had also been detected in heterologous species, which amplify and facilitate the reservoirs for BVDV [6,7]. Evidence of BVDV infection exists in 7 families (over 50 species) of Artiodactyla including Antilocapridae, Bovidae, Camelidae, Cervidae, Giraffidae, Suidae, and Tragulidae [6]. The circulation of BVDV-1 and BVDV-2 in cattle, pigs had been identified in China [8][9][10]. Our previous study has identified the prevalence of BVDV-1 in Chinese goat herds [11]. The infection of BVDV-2 in goat or sheep has been confirmed in India, Korea [12][13][14]. Therefore, the risk and prevalence of BVDV-1 and BVDV-2 in goat/sheep herds needs urgent attention.
Recently, high-throughput RNA-Sequencing (RNA-Seq) technologies give the opportunity to produce large numbers of sequence data in non-model organisms, and this method is better than the traditional microarray analysis [15,16]. It provides a thorough understanding of the host defense mechanisms and immune evasion strategies of viral infection [17]. The transcriptional landscape in the host upon virus infection facilitates the understanding of host immune responses and defense mechanisms upon the pathogenic microorganism infection at whole mRNA level, and provides new approaches to the potential control of virus infections.
Two biotypes of BVDV are recognized: cytopathic (cp) and non-cytopathic (ncp) strains. In experimental infected calves, BVDV-specific antibody is first detected shortly after viral clearance for both biotypes. T cell proliferative responses are detectable by 3-4 weeks postinfection with cpBVDV; while delayed to about 6-8 weeks after ncpBVDV infection [18]. The primary site of BVDV replication is immune tissue, viral replication results in altered cell function or cell death in different lymphoid populations. The resulting immune suppression occurs in all acute BVDV infections [19]. Peripheral blood mononuclear cells (PBMCs), including lymphocytes, monocytes and macrophages, play a pivotal role in the host innate or adaptive immune responses to viral infection. PBMCs were the main target of BVDV infection, infection of lymphocytes and monocytes by BVDV resulted in lymphoid depletion of B cells, T helper cells, cytotoxic T cells and γ-δ T cells [20,21]. PBMCs have been proved to be a suitable model for characterizing the host immune responses to virus infection and have been utilized for the evaluation of immune responses to animal viruses [17,22,23]. Global transcriptome analysis has been employed to explore the molecular events of host interaction with BVDV in bovine originated cells [24][25][26][27]. However, no report on BVDV-goat interactome is available to date. In this study, Illumina sequencing method was used to identify the transcriptome changes in BVDV-2 infected goat PBMCs. For the first time, we obtained the differentially expressed transcriptome profile in the goat PBMCs during BVDV-2 infection. The results will be helpful for better understanding the host responses to BVDV-2 infection and its relationship to viral pathogenesis in goats.

Determination of BVDV-2 replication in goat PBMCs
To confirm the replication of BVDV-2 in goat PBMCs, RT-PCR and qRT-PCR were performed. As shown in Fig.  1a, 5'-UTR fragment with~290 bp was amplified in infected goat PBMCs from 6 to 24 h post-infection (hpi). qRT-PCR detection showed similar results, the BVDV genome copy numbers increased from 6hpi and reached high levels at 12 and 24 hpi (Fig. 1b). These results confirmed the BVDV-2 infection in goat PBMCs. To explore the effect of early BVDV replication on gene expression, 12 hpi was selected for sampling and RNA-Seq. In addition, BVDV nucleotide was not detected in the mock infected PBMCs at any of the experimental time points.

DEGs analysis and functional annotation
After the gene mapping and the Cuffdiff analyses in terms of FRKM, a total of 449 genes were identified as significantly differentially expressed for infected group (B2), when comparing with the mock-infected control group (N, fold change (FC) ≥ ± 2, p < 0.05). Among the 449 genes, 97 were up-regulated and 352 genes were down-regulated (Fig. 2, Additional file 4: Figure S1, Additional file 1: Table S1).

Functional and PPI analysis of DEGs
After GO enrichment analysis, DEGs were enriched into different GO terms. For immune related DEGs, significant enrichment was observed in regulation of locomotion/ localization, immune response, inflammatory response, immune system process, defense response, regulation of cytokine production of BP group and in cytokine activity, chemokine activity, receptor binding of MF group. (Fig. 4 and Additional file 2: Table S2).

Partial validation of RNA-Seq data
To further validate the RNA-Seq data, DEGs with annotations associated with immune responses from RNAseq were selected for qRT-PCR analysis. As shown in Table 3, eighteen selected genes exhibited a concordant direction both in RNA-Seq and qRT-PCR analysis. The correlation coefficient between RNA-Seq and qRT-PCR results was high (R 2 = 0.91). Some of the DEGs were further determined by Western blot or ELISA, as shown in Fig. 7, the expression of Annexin A2 decreased obviously ( Fig. 7a) and the expression of TNF-α, GM-CSF and IL-6 were increased significantly (Fig. 7b). In addition, Viperin (used as a control) expression showed no change between infected and mock-infected group, which was consistent with the RNA-Seq result. These results confirmed that the differential expression genes identified by RNA-Seq is reliable.

Discussions
BVDV is one of the most important viral diseases of various species of domestic animals. Regardless of clinical presentation, all BVDV infections result in significant loss of immune tissue, the resulting immune suppression and increased severity of subsequent infections [28]. The different genotypes and biotypes of BVDV, wide spectrum of susceptible host, ability to induce persistently infection, as well as its ability to interfere both innate and adaptive immunity of the host, make it difficult for prevention and control. Microarray and RNA-Seq analysis are the two main techniques for transcriptome analysis. There were many reports about transcriptome changes upon animal virus infection by microarray analysis [29][30][31]. As a revolutionary tool, RNA-Seq has been widely used in recent years, such as BTV, PRRSV, PPRV and NDV [22,[32][33][34]. Several studies had explored the effect of BVDV infection on mRNA expression changes in different bovine cells [25,27,35], but understanding of BVDV-host interaction is still far from complete. Comparing to the  studies in cattle little is known regarding the cellular impact of BVDV-2 infection in goats. Herein, in the present study, goat PBMCs were infected with a ncp BVDV-2 strain, and the transcriptome was evaluated at 12 hpi to identify the global gene expression changes, to understand and delineate the mechanism of early responses induced by virus infection. Transcriptional changes of the DEGs involved in immunological processes were analyzed specially. The results from GO, KEGG and PPI analysis indicated that various numbers of DEGs were involved in different biological processes of host immune responses. The innate immunity is the first defense line against infectious disease. The PBMCs were the gatekeeper of virus spread during systemic viral infections and include several subpopulations that may cooperate in the activation processes [36]. These cells play important role in innate and adaptive immune responses. The pathogenassociated molecular pattern receptors (PAMPs) such as Toll-like receptors (TLRs) or RIG-I-like receptors (RLRs) recognize pathogens, then released a series of related antiviral cytokines. In this study, under BVDV-2 infection, the expression of TLRs, RLRs, interferon (IFN) and the ISGs were not significant induced. Contrarily, the mRNA levels of lysozyme, beta-defense, competent of complement system (C1q, C1r, C1s, CFD) and IFITM3 (one of the ISGs) was found to be down-regulated. A previous study showed that TLR3, type I IFN gene was up-regulated in ncp BVDV-infected monocytes at 1hpi but not at 24hpi [37]. It has also been reported that ncpBVDV dose not induce type I IFN in vitro and block the induction of type I IFN by dsRNA or other viruses and interferon tau-stimulated ISGs expression [18,38]. The interaction of ncp BVDV with its host cells impairs both of the innate and adaptive immunity [39]. During BVDV infection, viral RNA firstly triggers IFN synthesis, while the viral RNase E rns protein inhibits IFN expression; in addition, N pro promotes the degradation of IRF-3, which effectively blocks IFN expression in BVDV-infected cells [39][40][41]. So, the inhibition of immune genes mentioned above might result from the interference of BVDV-2 infection on the related immune pathways. The inhibition of IFN synthesis and other competent of innate immunity might play an important role in escaping innate immunity and the establishment of effective infection for BVDV in host cells.
The complement system, which consists of both soluble factors and cell surface receptors, is one of the major innate defense systems. The main role of complement system is to protect against infections, it also links the innate and adaptive immune responses [42]. In the present study, C1q, C1r, C1s, CFD were down-regulated, while C3 was up-regulated; imply that BVDV-2 infection in goat PBMCs may inhibit the classical or alternative pathways of the complement system activation.
Inflammation was one of the important anti-microorganism responses upon infection. Once exposed to infectious agents, host cells produced certain cytokines, such as TNF-α, IL-1 and IL-6, resulting in the development of inflammation. TNF-α is an essential mediator of inflammation and also facilitates the transition from innate to adaptive immunity. IL-6 affects both inflammation and adaptive immunity. It promotes some aspects of inflammation, especially in response to tissue damage and severe infections since it is a major mediator of the acute phase reaction and of septic shock. Many viruses, including BVDV, could induce inflammation after infection. Studies in cattle and sheep have shown that the pathology of BVDV associated with the replication of virus and the production of pro-inflammatory cytokines and induced inflammatory responses [43][44][45]. One study examined the mRNA expression in tracheo-bronchial lymph nodes of BVDV infected beef calves and found that high virulence BVDV-2 strain induced pro-inflammatory (TNF-α, IL-12, IL-1β, IL-2, IFN-γ) and anti-inflammatory (IL-4 and IL-10) cytokines, while low virulence BVDV-1a strain only up-regulated IL-12 and IL-15 gene expression [46]. In an other study, both cp and ncp BVDV biotypes suppressed pro-inflammatory cytokines TNF-α, IL-1β, and IL-6, but did not change IL-12 and INF-γ gene expression in bovine PBMCs [37]. As a member of pestivirus, virulent CSFV infection in pigs resulted in the secretion of cytokines associated with inflammation or apoptosis such as TNF-α, IL-2, IL-4, IL-6, and IL-10 [47]. The BVDV-2 strain used in this study induced fever, viremia and lymphopenia in experimentally infected goats (unpublished data). The clinical signs, together with the pattern of increased TNFα, IL-6, IL-12, IL-10 and IL-17 mRNA level observed in BVDV2-infected goat PBMCs in this study, suggested that BVDV-2 induced an acute inflammatory response in the early stage of infection. The enriched pathways related to cytokine/chemokine signaling provide obvious evidence that specific immune responses are activated during that acute phase of BVDV infection. The serum amyloid A (SAA) family comprises a number of apolipoproteins, which include acute-phase SAAs (A-SAAs) and constitutive SAAs (C-SAAs). A-SAAs have been identified in all vertebrates and could be induced as much as 1000-fold during inflammation [48]. Transcription of SAAs members LOC102168428 (SAA like) and LOC100860781 (SAA3) were found to be upregulated in this study. Members of A-SAAs have been identified to be activated upon infection and proven useful as an inflammatory marker for many viral diseases of animals [49][50][51]. Similarly, Ganheim C. et al. reported a significant acute phase response with elevated values of SAA and other acute phase proteins in calves experimentally infected with BVDV [52]. So, SAA might be considered as a diagnostic marker for BVDV induced inflammation.
Viral infection usually induces the recruitment of inflammatory cells to the infection site, and this activity is regulated by various cytokines and chemokines. Chemokines are a great superfamily of at least 50 small (8-to 10-kDa) structurally related chemoattractant proteins, which play a key role in initiating the innate and subsequently adaptive immune responses. Chemokines have many roles in the regulation of leukocyte development, angiogenesis, tumor growth and metastasis [53]. They coordinate the migration of leukocytes and hence dictate the course of many inflammatory and immune responses by recruiting different immune cells toward sites of infection [54]. Monocytes/macrophages are part of the first line of defense and have been shown to release a variety of chemokines in response to infection [31]. Studies have determined chemokines expression upon the infection of several viruses, such as PRRSV, PEDV and PCV2 [30,32,55]. Study by Helal et al. showed that low expression of CXCR4 and high expression of IL-10 is associated with the production of PI calves in a herd level [56]. Trzeciak-Ryczek A. et al. reported that ncpBVDV infection increases CXCR4, CXCL12 mRNA expression in bovine PBMCs [57]. In another study, detection of chemokine profile in cp and ncp BVDV infected bovine monocytes/macrophages demonstrated up-regulation of several key chemokines of the CCL and CXCL families to cpBVDV, but not ncpBVDV [31]. In our study, BVDV-2 infected goat PBMCs showed increased transcription of CCL3, CCL4, CCL5, CCL20, CXCL10 and decrease of CCL2, suggesting that such chemokines may contribute to the recruitment and regulation of macrophages or other inflammatory cells to the infected sites after infection.
GM-CSF was another up-regulated DEG found in our study. It is not only an inducer of differentiation and proliferation of granulocytes/macrophages but also involved in a wide range of biological processes in both innate and adaptive immunity [58]. GM-CSF has been widely used as adjuvant for vaccines and has been shown powerful as an important therapeutic target in several autoimmune and inflammatory diseases [58,59]. The induction of GM-CSF expression has been identified in several viruses and shown to alter pathogenic "M1-like" macrophage inflammation after influenza A virus infection [60][61][62]. The high expression of GM-CSF in goat PBMCs indicated the regulation role of it in host immune responses after BVDV-2 infection.

Conclusions
In conclusion, the changes of a series of immunity related genes in BVDV-2 infected goat PBMCs were observed by RNA-Seq analysis. The results will be useful for exploring and further understanding host responses to BVDV-2 infection in goats.

Goat, PBMCs culture and virus infection
Goats used for blood collection were purchased from a goat farm (private farm, 210 goats) located in Jurong county, Jiangsu province. The goats were vaccinated against goat pox virus, foot and mouth disease virus, peste des petits ruminants virus and Mycoplasma mycoides subsp.capri. The animal was further screened for BVDV and viral specific neutralizing antibodies using RT-PCR and virus neutralization test, respectively. EDTA anticoagulant whole blood were collected from the goats (n = 3) and PBMCs were separated using Histopaque-1077 (Sigma) by density gradient centrifugation at 500×g for 20 min, then washed three times with RPMI-1640 medium at 500×g for 10 min. Cells from each animal were suspended to 1× 10 6 cells/mL with complete RPMI 1640 medium (RPMI 1640 containing 10% FBS) and seeded in six well plates (n = 4, 8 wells/animal). Two plates were infected with BVDV2 (2 × 10 5 TCID 50 /mL, MOI = 0.1) and the other two plates were served as mockinfected control. Samples from one infected plate and one mock-infected plate were collected at 6, 12 and 24 hpi for virus replication detection. At 12 hpi, medium was discarded and the cells of the left two plates (infected and mock-infected) were washed twice with PBS, mixed and pelleted for RNA-Seq.
RNA-Seq transcriptome libraries were prepared by TruSeqTM RNA sample preparation Kit from Illumina (San Diego, CA), using 1 μg of total RNA. Shortly, mRNA was isolated with polyA selection by magnetic oligo (dT) beads and fragmented into small pieces using fragmentation buffer. cDNA synthesis, end repair, Abase addition and ligation of the Illumina-indexed adaptors were performed according to Illumina's protocol. Libraries were then selected for cDNA target fragments of 200-300 bp on 2% Low Range Ultra Agarose followed by PCR amplified for 15 cycles using Phusion DNA polymerase (NEB). After quantified by TBS380, paired-end libraries prepared for both the infected and mockinfected samples (B2 and N) were sequenced with the Illumina HiSeq PE 2X151bp read length. RNA-Seq was performed by Shanghai Biozeron Biotech Co. Ltd.

Reads quality control and mapping
The raw paired end reads were trimmed and quality controlled by Trimmomatic with default parameters (http://www.usadellab.org/cms/index.php?page=trimmomatic). Clean reads were then separately aligned to the reference caprine genome (https://www.ncbi.nlm.nih. gov/assembly/GCF_001704415.1) with orientation mode using TopHat software (http://tophat.cbcb.umd.edu/), which can align RNA-Seq reads to a genome in order to identify gene expression and exon-exon splice junctions. It is built on the ultrafast short read mapping program Bowtie2 to map with default parameters.

Differential expression genes (DEGs) analysis and annotation
To identify DEGs between the two samples, the expression level for each transcript was calculated using the fragments per kilobase of exon per million mapped reads (FPKM) method. Cuffdiff (http://cufflinks.cbcb.umd.edu/) was used for differential expression analysis and the DEGs were selected using the following criteria: the logarithmic of fold change was greater than 2 and the p-fdr should be less than 0.05. All DEGs were subjected to Gene Ontology (GO) annotation based on GO database.

Functional and protein-protein interaction (PPI) analysis of DEGs
To understand the functions of the DEGs, GO functional enrichment and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis were carried out by Goatools (https://pypi.org/project/goatools/) and KOBAS (http://kobas.cbi.pku.edu.cn/). DEGs were significantly enriched in GO terms and KEGG pathways when their p-value was less than 0.05.
PPI network among the DEGs was analyzed using the STRING (http://string-db.org/) database, which included direct and indirect associations of proteins. After analyzing the result from STRING analysis and expression change information for each DEG, the network figure was drawn for the selected DEGs (connected with one or more DEGs) by using Cytoscape software.

RT-PCR and quantitative real-time RT-PCR (qRT-PCR)
Total RNAs were extracted from infected and mockinfected PBMCs using TransZol UP reagent (Transgen, Bio, Inc.,China) according to the manufacturer's instruction.
The RT-PCR was performed to identify the replication of BVDV and carried out with EasyScript onestep RT-PCR supermix (Transgen, Bio, Inc., China) in a 20 μl reaction mixture containing 10 μl of 2 × R-Mix, 20 pM of each primer (F: 5'-ATGCCCWTAGTAGGAC-TAGCA-3', R: 5'-TCAACTCCATGTGCCATGTAC-3'), 0.4 μl of E-Mix, 2 μl extracted RNA and 6.6 μl ddH 2 O. The reaction was run in a thermocycler (Mjmini, BIO-RAD) with the following program: reverse transcription at 45°C for 30 min; denaturation at 94°C for 5 min, 35 cycles composed of denaturation at 94°C for 30 s, annealing for 30 s at 54°C and extension at 72°C for 30 s; and was terminated with a final extension of 10 min at 72°C. Amplification products were detected by electrophoresis in 1.2% agarose gels.
The qRT-PCR amplification was carried out with TransScript one-step qRT-PCR supermix (Transgen, Bio, Inc., China) in a 20 μl reaction mixture containing 10 μl of 2 × Supermix, 20 pM of each primer (Table 4), 0.5 μl of E-Mix, 0.4 μl of passive reference Dye and 2 μl extracted RNA. The reaction was run in ABI Step One instrument as the following procedure: samples were incubated at 45°C for 5 min firstly; then heated at 94°C for 30 s and a two-step cycle (5 s at 94°C, 30 s at 60°C) was repeated for 40 cycles. GAPDH was used as the internal control and relative quantification of target gene expression was the target transcript in infected group to that of mock-infected group and expressed as -ΔΔCt.

Virus titration
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 fixed at 4°C for 30 min with ice cold absolute ethyl alcohol and subjected to immunofluorescence staining with BVDV-specific mouse monoclonal antibody Mix (RAE2020, AHVLA, UK; 1:200 diluted in PBS) and FITC-conjugated goat anti-mouse IgG (BOSTER, Wuhan, China; 1:200 diluted in PBS). The fluorescence signals were observed undera fluorescence microscopy (ZEISS) and viral titers was expressed as the 50% tissue culture infective dose (TCID 50 )/mL by Reed-Muench method.