Skip to main content

RNA-Seq based transcriptome analysis during bovine viral diarrhoea virus (BVDV) infection

Abstract

Background

Bovine viral diarrhoea virus (BVDV) is the member of the genus Pestivirus within the Flaviviridae family and responsible for severe economic losses in the cattle industry. BVDV can employ ‘infect-and-persist’ strategy and ‘hit-and-run’ strategy to remain associated with hosts and thus contributes to BVDV circulation in cattle herds. BVDV have also evolved various strategies to evade the innate immunity of host. To further understand the mechanisms by which BVDV overcomes the host cell innate immune response and provide more clues for further understanding the BVDV-host interaction, in this descriptive study, we conducted a investigation of differentially expressed genes (DEGs) of the host during BVDV infection by RNA-Seq analysis.

Results

Our analysis identified 1297, 1732, 3072, and 1877 DEGs in the comparison groups mock vs. MDBK cells infected with BVDV post 2 h (MBV2h), mock vs. MBV6h, mock vs. MBV12h, and mock vs. MBV24h, respectively. The reproducibility and repeatability of the results were validated by RT-qPCR. Enrichment analyses of GO annotations and KEGG pathways revealed the host DEGs that are potentially induced by BVDV infection and may participate in BVDV-host interactions. Protein-protein interaction (PPI) network analyses identified the potential interactions among the DEGs. Our findings suggested that BVDV infection induced the upregulation of genes involved in lipid metabolism. The expression of genes that have antiviral roles, including ISG15, Mx1, OSA1Y, were found to be downregulated and are thus potentially associated with the inhibition of host innate immune system during BVDV infection. The expression levels of F3, C1R, KNG1, CLU, C3, FB, SERPINA5, SERPINE1, C1S, F2RL2, and C2, which belong to the complement and coagulation signalling cascades, were downregulated during BVDV infection, which suggested that the complement system might play a crucial role during BVDV infection.

Conclusion

In this descriptive study, our findings revealed the changes in the host transcriptome expression profile during BVDV infection and suggested that BVDV-infection induced altering the host’s metabolic network, the inhibition of the expression of antiviral proteins and genes within the complement system might be contributed to BVDV proliferation. The above findings provided unique insights for further studies on the mechanisms underlying BVDV-host interactions.

Background

Bovine viral diarrhoea viruses (BVDV), with positive, non-segment, single stranded RNA genome, comprise a heterogeneous group of viruses which belong to the genus Pestivirus, the Flaviviridae family [1]. BVDV is a globally-distributed virus and its infection has been detected not only in domesticated ruminants, but also in wild ruminants and wild boars [2,3,4,5,6,7]. It is responsible for numerous clinical disease syndromes in cattle and severe economic losses in the cattle industry.

BVDV is divided into two species, namely, Pestivirus A and Pestivirus B, which are formerly known as BVDV 1 and BVDV 2, respectively [8]. Given the genetic heterogeneity of BVDV, both species are divided into multiple genotypes. To date, Pestivirus A is divided into at least 17 genotypes, while Pestivirus B is divided into 4 genotypes [9]. Both BVDV species have two biotypes, including the non-cytopathogenic biotype (NCP) and the cytopathogenic biotype (CP) [1]. NCP BVDV can establish persistent infection by crossing the placenta and infecting the foetus, which is born as a persistently infected (PI) calf after 30 to 150 days of pregnancy [10]. PI cattle are the most important sources of BVDV infection and are the key factors affecting BVDV circulation in farms. PI cattle can carry and shed viruses throughout their lifetimes, which in turn produces more PI calves from PI cattle. Therefore, the detection and elimination of PI cattle are of great significance for the prevention and control of BVD/MD in cattle herds.

BVDV replication, occurring in the cellular cytoplasm, is a multistage process. The binding and entry of BVDV begin with attachment or interaction of the virion with specific host cell receptors, followed by internalization and pH dependent fusion of the viral envelope and cell membrane [11, 12]. After entry into the host cell is complete, viral RNA is released into the host cell cytoplasm and RNA translation begins. Initiation of the translation process is mediated by the IRES which binds specifically to the 40S ribosomal subunit in the absence of any additional translation initiation factors [13]. Then, the 80S complex was formed at the start codon AUG and initiate the translation of polyprotein which are cleaved into structural and non-structural proteins [14]. Viral non-structural proteins assemble into a functional replicase complex to catalyze transcription of positive-sense RNA into full-length complementary strand negative-sense RNA which provide template for the replicase complex to synthesize additional positive-sense RNA molecules [15]. The replication process begins with a positive-strand replicase complex comprised of viral and cellular components formed at the 3′ terminus of the genome. Progression from initiation to elongation occurs after the synthesis of nascent RNAs 8–10 nucleotides in length [16]. Elongation mediated by viral-polymerase displaces the positive strand from the RI template, allowing recycling of the template while elongation of the prior nascent strand continues. Then, capsid protein binds to viral genomic RNA and identifies the cytoplasmic region containing the envelope protein that budding in the endoplasmic reticulum cavity [17,18,19]. BVDV virions appear to mature in intracellular vesicles at the Golgi apparatus or endoplasmic reticulum. The intact virions are released by exocytosis with detection reported as early as 8 h post infection [20].

Viruses employ two different strategies to remain associated with hosts, including ‘infect-and-persist’ strategy and ‘hit-and-run’ strategy [21]. Viruses that use the ‘infect-and-persist’ strategy rather than the ‘hit-and-run’ strategy are capable of efficiently evading the immune defence systems of the susceptible hosts and have evolved elaborate mechanisms of evasion [22, 23]. BVDV can employ both strategies to remain associated with hosts and thus contributes to BVDV circulation in cattle herds [23].

The innate immunity serves as the first line of defence against foreign microbial pathogens. Considering that full induction of antiviral innate immunity severely limits viral replication, viruses have evolved diverse strategies to evade the innate immunity [24]. Current studies have shown that BVDV can evade the innate immune system via multiple strategies, such as inhibition of IFN synthesis, inhibition of IFN protein activity, and interference with expression or activity of the IFN-induced antiviral effector proteins [24]. Degradation of IRF-3 induced by BVDV Npro has been shown to decrease the production of interferon [25, 26]. BVDV inhibits interferon production and lymphocyte apoptosis through the degradation of viral RNA mediated by Erns [27, 28]. The mechanisms underlying BVDV-host interactions and the pathogenesis of BVDV infection are complex. Successful evasion from the host immune defence system is the basis for persistent BVDV infection. Although BVDV employs several strategies to evade the immune system, especially the innate immune system, these strategies do not reveal the complete mechanisms underlying BVDV-host interactions during infection. The present study aimed to investigate the changes in gene expression and identify the active or suppressed signal pathways induced by BVDV during the early stage of virus infection, which could be associated with the BVDV-host interactions. This work provided a theoretical basis for further studies on BVDV-host interactions during the early stage of viral infection.

Results

Growth curves of BVDV with different MOIs

To determine the period of intracellular replication, the infection dynamics of BVDV on MDBK cells was defined using one-step growth curve method. The growth curves of BVDV at 0.1, 1, and 10 MOI were determined by measuring the changes in the number of copies of the BVDV genome in the cell supernatant by RT-qPCR (Fig. 1). The growth curves indicated that the number of viruses remained more or less stable in the cell supernatant within 12 h-post infection, which was followed by a rapid increase in the viral load in cell supernatant. This observation suggested that it is the intracellular replication period of BVDV within 12 h post-infection. After 12 h post-infection, the virus began to be released from cells in large quantities. Due to all stages of BVDV life cycle were observed within 24 h post infection, the samples used for transcriptome sequencing were collected at the 2, 6, 12, and 24 h time points.

Fig. 1
figure1

Growth curves of BVDV in MDBK cells. BVDV showed similar proliferation trend after MDBK cell were inoculated with virus of different MOI. Each data point represented the average of three replicates. The error bar indicates standard deviation

Evaluation of transcriptome sequencing data

At least 6.0 Gb of clean data from each sample was obtained from transcriptome sequencing from Illumina NovaSeq 6000 platform and were available for further expression level analysis after quality control. Q30 percentages of clean data for all samples were higher than 94.36%, and the GC contents of the clean data for all samples ranged between 54.3 to 55.32% (Table 1). For further analysis, the high-quality clean reads were mapped to the reference Bos taurus genome (UMD3.1). Approximately 95.97 to 96.32% of clean reads were successfully mapped to the reference Bos taurus genome (UMD3.1). In addition, 91.58-93.44% of clean reads were uniquely mapped to the reference Bos taurus genome (UMD3.1). A total of 53,929 transcripts were reference-based assembled from the mapped data using StringTie software. Length distribution statistics of transcripts indicated that the length of more than half of the transcripts were greater than 1800 bp (Additional file 1: Table S1). After removing low-quality and redundant sequences, a total of 53,391 transcripts were obtained, corresponding to 26,740 reference transcripts. Simultaneously, 26,975 genes were assembled from the mapped reads, out of which 24,616 genes successfully mapped to the reference Bos taurus genome (UMD3.1). A total of 26,651 novel transcripts and 2359 of novel genes were predicted by gffcompare software and were subsequently annotated by comparison with the NR, Swiss-Prot, and Pfam databases (Additional file 1: Table S2). Further studies were performed based on the reference genes.

Table 1 Summary statistics for sequence quality control and mapped data of samples

Functional annotation and classification

To generally described the functions and pathways of genes obtained from RNA-Seq, functional annotation and classification were performed by comparing the sequences with the GO, COG, and KEGG databases. These analyses could provide general information for function and classification of genes or transcripts, which can only annotated up to level 2. GO is a gene function classification system that is often used to describe gene properties. Annotation results for the three databases were shown in Additional file 2: Figure S1. In the present study, GO terms were extracted using Blast2GO software. Genes annotated by GO database were classified into the three GO categories, namely, biological progress (BP), cellular component (CC), and molecular function (MF). The genes were assigned to a total of 63 GO terms, which were predominantly comprised of cellular process, single-organism process, biological regulation, regulation of biological process, metabolic process, cell, cell part, organelle, binding, catalytic activity, signal transducer activity. The assignments of the top 20 level 2 GO terms are shown in Additional file 2: Figure S2. Genes annotated with the COG database were best matched to four gene functional types (information storage and processing, cellular processes and signalling, metabolism, poorly characterized) and assigned to 24 functional categories. The top three functional terms were post-translational modification, protein turnover, chaperones, translation, ribosomal structure and biogenesis, and signal transduction mechanisms (Additional file 2: Figure S3).

The KEGG database contains abundant pathway information that can help understand the biological functions of genes at the systems level. In the present study, gene annotations using KEGG pathway database was performed using KOBAS (Additional file 2: Figure S4). Based on KEGG analysis, the genes were grouped into six first categories, corresponding to 43 s categories. The top five second-category KEGG pathways included signal transduction, sensory system, cancers overview, immune system, endocrine system, and transport and catabolism.

Differential gene expression analysis

To investigate changes in gene expression profiles induced by BVDV infection, FPKM expression values of the genes were calculated based on the read counts using featureCounts software. The FC values of each gene at different time post infection compared with control group were also calculated using DESeq2 software. The threshold values FDR ≤ 0.05 and |Log2FC| ≥ 1 were used to identify DEGs between groups at different time points post-infection relative to the control group. Differential expression analysis identified 1297, 1732, 3072, and 1877 DEGs in different comparison groups Mock vs. MDBK cells infected with BVDV post 2 h(MBV2h), Mock vs. MBV6h, Mock vs. MBV12h, Mock vs. MBV24h post-BVDV infection (Fig. 2, Additional file 2: Figure S5 and Additional files 3, 4, 5, 6, 7, 8, 9). Differential expression analysis for the comparison groups, Mock vs. MBV2h, Mock vs. MBV6h, Mock vs. MBV 12 h, and Mock vs. MBV24h, identified 743, 1068, 1863, and 1035 upregulated genes and 554, 665, 1342, and 853 downregulated genes, respectively (Fig. 2).

Fig. 2
figure2

Volcano plot of global DEGs in different comparison groups, Mock VS MBV2h (a), Mock VS MBV6h (b), Mock VS MBV12h (c), and Mock VS MBV24h (d). Red dots (Up) represent significantly up-regulated genes (P < 0.01, fold change ≥2); Yellow dots (Up) represent extremely significantly up-regulated genes (P < 0.05, fold change ≥2); Mazarine dots (Down) represent significantly down-regulated genes (P < 0.01, fold change ≤0.5); Wathet dots (Down) represent extremely significantly downregulated genes (P < 0.05, fold change ≤0.5); black dots (nosig) represent insignificantly differential expressed genes

During BVDV replication, double-stranded RNA (dsRNA) are generated to facilitate multiplication of progeny virus. dsRNA is an important pathogen-associated molecular pattern (PAMP) of BVDV that activates the host innate immunity. In our previous study, we detected the dynamic changes that occur in the BVDV negative strand in MDBK cells infected with 10 MOI BVDV using strand specific real-time PCR (Additional file 2: Figure S6). The negative strand of BVDV RNA in the host cells was detected at 2 h post-infection, and the amount of the negative strand of BVDV increased from 2 h to 6 h after BVDV infection. However, the amount of the negative strand of BVDV decreased at 12 h and gradually increased after 12 h. In addition, all biological processes in the life cycle of BVDV could be observed during the logarithmic growth phase. The DEGs identified in the comparison groups Mock vs. MBV2h, Mock vs. MBV6h, Mock vs. MBV12h, and Mock vs. MBV24h can help understand the mechanisms underlying the host response to BVDV infection. A total of 24 DEGs were identified, including ACLY, HMGCS1, INSIG1, LPIN1, HMGCR, and ACSS2 in group MBV2h vs. MBV6h (Additional file 7). The genes ACLY, HMGCS1, INSIG1, LPIN1, HMGCR, and ACSS2 were found to be significantly upregulated in comparison groups Mock vs. MBV6h, Mock vs. MBV 12 h, and Mock vs. MBV24h groups.

Enrichment analysis of GO terms and KEGG pathway

Compared with general description for properties of gene or transcripts with functional annotation, the lowest level of gene function and KEGG pathway can be annotated by enrichment analyses. In addition, through fisher test, significantly enriched gene function and KEGG pathway of DEGs will be determined. The most detail information on gene function and KEGG pathway provided by enrichment analyses which facilitate us to screen the unique insights on the response to BVDV-infection. Enrichment analyses of GO terms and KEGG pathways were performed for DEGs using GOatools and KOBAS, respectively. GO terms and KEGG pathways that satisfy the corrected p-value ≤0.05 were considered significantly enriched. The enriched GO terms were ordered based on the three categories, including biological processes (BP), cellular components (CC), and molecular function (MF) (Fig. 3). DEGs in comparison group Mock vs. MBV2h were primarily associated with GO terms related to stimulation, G protein-coupled receptor signalling pathway, transmembrane signalling transduction, such as sensory perception of smell, detection of stimulus, sensory perception, G-protein coupled receptor signalling pathway, G-protein coupled receptor activity, transmembrane signalling receptor activity (Fig. 3a). The majority of DEGs in comparison group Mock vs. MBV6h were associated with stimulation-related biological processes, cellular processes, and signal transduction-related molecular functions, such as detection of stimulus, G-protein coupled receptor signalling pathway, positive regulation of biological process/cellular process, negative regulation of biological process/cellular process, G-protein coupled receptor activity, signal transducer activity, signalling receptor activity, and transmembrane signalling receptor activity(Fig. 3b). The majority of DEGs in comparison group Mock vs. MBV12h were assigned to biological processes associated with stimulus, biosynthesis and metabolism, cellular components (cytoskeletal part, membrane part, cytoplasm), and molecular functions, such as transmembrane receptor activity, transmembrane signalling receptor activity, RNA binding, binding, ATP binding, protein binding (Fig. 3c and Additional file 10). The majority of DEGs in comparison group Mock vs. MBV24h were associated with the regulation of biological processes, cellular processes, and signal transduction, such as small molecule biosynthetic process, sensory perception of smell, regulation of signalling, positive regulation of biological process, positive regulation of cellular process, regulation of response to stimulus, G-protein coupled receptor activity, transmembrane receptor activity, signal transducer activity, protein binding, and small molecule binding (Fig. 3d and Additional file 11). DEGs identified from comparison group MBV2h vs. MBV6h were mainly assigned to biological processes associated with lipid synthesis and metabolism, such as lipid biosynthetic process, cholesterol biosynthetic process, secondary alcohol biosynthetic process, and lipid metabolic process (Additional file 2: Figure S7).

Fig. 3
figure3

GO enrichment analysis of DEGs in different comparison groups, Mock VS MBV2h (a), Mock VS MBV6h (b), Mock VS MBV12h (c), and Mock VS MBV24h (d). GO terms are on the x axis. Enrichment ratio of genes shown as GO terms for BP, CC, and MF. * means GO categories with significant enrichment

The genes in an organism act in a coordinated manner to perform their biological functions. Pathway analysis can help better understanding the biological functions of genes. Enrichment analysis of KEGG pathways associated with the DEGs can be used to determine the biochemical metabolic and signal transduction pathways. The analysis identified 24, 16, 27, and 44 KEGG pathways that were significantly enriched in DEGs in comparison groups Mock vs. MBV2h, Mock vs. MBV6h, Mock vs. MBV12h, and Mock vs. MBV24h, respectively (Fig. 4). The significantly enriched pathways in DEGs from comparison group Mock vs. MBV2h mainly included the complement and coagulation cascades, TGF-beta signalling pathway, FoxO signalling pathway, ErbB signalling pathway, PI3K-Akt signalling pathway, and signalling pathways associated with pathogen-infection. The pathways significantly enriched in DEGs from comparison group Mock vs. MBV6h mainly included complement and coagulation cascades, steroid biosynthesis, TGF-beta signalling pathway, and fatty acid biosynthesis. The significantly enriched pathways in DEGs from comparison group Mock VS MBV12h mainly included steroid biosynthesis and p53 signalling pathway (Additional file 12). The pathways that were significantly enriched in the DEGs from comparison group Mock vs. MBV24h mainly included amino acid metabolism, MAPK signalling pathway, HIF-1 signalling pathway, TNF signalling pathway, and mTOR signalling pathway (Additional file 13). KEGG pathways which were significantly enriched in all comparison groups included complement and coagulation cascades, TGF-beta, FoxO, ErbB, Hippo, and AGE-RAGE signalling pathways in diabetic complications. The PI3K-Akt, NOD-like receptor, HIF-1, MAPK, TNF, and mTOR signalling pathways play important roles in the recognition of pathogen-associated molecular patterns, signalling transduction, and the regulation of host immune system. These results provided important insights on the interactions between BVDV and host.

Fig. 4
figure4

KEGG Pathway enrichment analysis of DEGs in different comparison groups, Mock VS MBV2h (a), Mock VS MBV6h (b), Mock VS MBV12h (c), and Mock VS MBV24h (d). The name of KEGG pathway are on the x axis. Enrichment ratio of genes shown as the name of KEGG pathway for seven categories, environmental information processing (EIP), genetic information processing (GIP), cellular processing (CP), organismal systems (OS), drug development (DD), Human diseases (HM), Metabolism (M). * means KEGG pathway with significant enrichment

Protein-protein interaction network (PPI) analysis of DEGs

To further understand the biological relevance of DEGs, PPI network analyses of DEGs was performed using STRING database. The PPI protein-protein network was visualized using NetworkX package in Python.

PPI network analyses of DEGs in groups Mock vs. MBV2h and Mock vs. MBV24h were performed based on DEGs with FDR ≤ 0.05 and |Log2FC| ≥ 2. PPI network analyses identified a total of 233 interaction relationships among 89 genes within the DEGs from comparison group Mock vs. MBV2h (Additional file 2: Figure S9). The DEGs OAS1Y, FOS, EGR1, and PAI2 were found to play important roles in maintaining the tight connection of the whole network. The analysis additionally identified a total of 172 interaction relationships among 46 genes within the DEGs from comparison group Mock vs. MBV24h (Additional file 2: Figure S10). The DEGs OAS1Y, CDKN1A, FOSB, NOTCH3, FLT3, HDAC5, RELB, and CsDKN2B were found to play important roles in maintaining the tight connection of the whole network. Furthermore, 41 interaction relationships were identified among 14 genes in DEGs comparison group MBV2h vs. MBV 6 h (Fig. 5). The DEGs ACLY, HMGCS1, FDFT1, HMGCR, and ACSS2 were found to play important roles in maintaining the tight connection of the whole network. The network centrality coefficients of the genes are shown in Additional file 1: Table S3-S5. The DEGs OAS1Y, FOS/FOSB, EGR1, PAI2, CDKN1A, NOTCH3, FLT3, HDAC5, RELB, CDKN2B play key roles in multiple cellular activities (including cell cycle, apoptotic process, regulation of gene expression, and disease development) and signal transduction. The above results indicated that BVDV interferes with the expression of these genes to promote its own replication in host cells. The DEGs ACLY, HMGCS1, FDFT1, HMGCR, and ACSS2 were associated with lipid synthesis and metabolism, thereby suggesting that BVDV infection induced lipid synthesis and metabolism activity in host cells.

Fig. 5
figure5

Protein-protein interaction network analysis of DEGs in group MBV2h VS MBV6h. The STRING database (http://string-db.org/) was used to analyse the protein-protein interaction network based on the proteins corresponding to selected DEGs. Protein interaction relationship of selected DEGs existing in the database were extract for the construction of network. Then, the visualization of the network is carried out using network X under Python. The size of the node is proportional to node degree. ACLY, HMGCS1, FDFT1, were the three genes with the most node degree in the analysis, labelled in red blot, green blot and yellow blot, respectively

Quantitative real-time PCR (RT-qPCR) analysis

To validated the reproducibility and repeatability of DEGs identified from transcriptome sequencing, we randomly selected ten genes, namely, C8orf4, PSPH, ISG15, EGR1, SIGLEC10, FABP3, GRIP2, GALNT18, GATM, and IFITM1, for RT-qPCR analysis (Fig. 6). The results showed that these genes were significantly differentially expressed and were consistently upregulated or downregulated with the gene expression changes based on RNA-Seq, thereby suggesting that the DEGs obtained from transcriptome sequencing were reliable. The genes ACLY, ACSS2, INSIG1, HMGCR, HMGCS1, and LPIN1 were also validated by RT-qPCR (Fig. 7). The trends in differential expression of these genes confirmed by RT-qPCR in the Mock vs. MBV6h, Mock vs. MBV 12 h, and Mock vs. MBV24h groups were consistent with those of the transcriptome sequencing results. However, the genes ACSS2, INSIG1, HMGCR, HMGCS1, and LPIN1 were downregulated in the Mock vs. MBV group at 12 h, which was inconsistent with the transcriptome analysis. Correlation was measured by scatterploting log2 fold changes between RNA-Seq and RT-qPCR. As shown in Fig. 8, high correlation was observed between the results from RT-qPCR and RNA-Seq, with the correlation coefficient (R2) as high as 0.8004.

Fig. 6
figure6

Expression level of ten randomly selected gene validated by RT-qPCR. β-actin gene was used as an internal control and relative quantity of gene expression (fold change) of each gene was calculated with the comparative 2-ΔΔCT method. Values (RT-qPCR) shown were mean with SD

Fig. 7
figure7

Expression level of genes, ACLY, ACSS2, INSIG1, HMGCR, HMGCS1, LPIN1, validated by RT-qPCR. β-actin gene was used as an internal control and relative quantity of gene expression (fold change)of each gene was calculated with the comparative 2-ΔΔCT method. Values (RT-qPCR) shown were mean with SD

Fig. 8
figure8

Correlation of fold change analysed by data obtained using RT-qPCR (x axis) with RNA-Seq platform (y axis). Correlation analysis was performed using GraphPad software 6.0 (San Diego, CA)

Discussion

BVDV is responsible for the severe economic losses because of decreased production performance, slow foetal growth, diarrhoea, respiratory symptoms, reproductive dysfunctions, immunosuppression, concurrent infections, and persistently infected (PI) cattle [29, 30]. BVDV employed various strategies to evade the innate immune system. More in-depth studies on BVDV-host interactions will be of great importance in elucidating the mechanisms underlying the pathogenesis and the progression of persistent BVDV infection.

RNA-Seq is a recently developed high-throughput approach for transcriptome profiling that achieves considerably more precise measurement of the expression levels of transcripts and their isoforms compared to other methods [31]. The mechanisms responsible for the pathogenesis and the formation of persistent BVDV infection are complex, and RNA-Seq provides an effective means for the systematic study of host response to viral infection and for investigating BVDV and host interactions. In the present study, transcriptomes of MDBK cells inoculated with BVDV at the early stage of infection were sequenced using the NovaSeq 6000 platform. DEGs were identified in MDBK cells infected with BVDV at 2, 6, 12, and 24 h post-infection. The expression changes of randomly selected genes, which was analysed by RT-qPCR, showed a strong correlation with that identified with RNA-Seq. This supported the reliability of data from RNA-Seq analysis.

Metabolism provides the energy and materials required for all biological processes [32, 33]. Viruses are obligate parasites that derive all energy and materials necessary for infection and replication from the host cell [32, 34, 35]. Some viruses generate a suitable intracellular microenvironment for the viral life cycle by altering the host’s cellular metabolic network [36]. In the present study, analysis of the DEGs revealed that several key genes that were upregulated during the BVDV life cycle; these genes included ACLY, HMGCS1, HMGCR, ACSS2, INSIG1, and LPIN1, which are related to lipid metabolism. The LPIN1 gene, which encodes the lipin-1 protein, plays critical roles in adipocyte differentiation and lipid metabolism [37, 38]. Lipin-1 can bidirectionally regulate fat metabolism and exerts significant effect on fat deposition, which are potentially correlated with animal carcass traits and meat quality traits [39,40,41,42,43]. The INSIG1 gene primarily affects lipid metabolism by regulating sterol regulatory element binding protein (SREBP) and 3-hydroxy-3-methylglutaric coenzyme A reductase (HMGCR) and is correlated with economic traits of livestock [44,45,46,47,48]. The ACLY gene is related to intramuscular fat deposition in cattle, which is of great significance for beef quality improvement [49, 50]. Other genes, such as HMGCS1, HMGCR, and ACSS2, were also associated with growth performance and meat quality improvement in cattle [51,52,53,54,55,56]. The dysregulated expression of these genes can also be associated with the development of disease, such as cancers, viral infection, and diabetes mellitus [34, 35, 57,58,59,60]. Previous studies confirmed that positive-strand RNA viruses and enveloped viruses rewire lipid metabolism for replication [34, 61, 62]. In the present study, the genes that were upregulated during BVDV infection were related to lipid metabolism, thereby suggesting that BVDV might alter host lipid metabolism to facilitate itself survival. This phenomenon was similar to that of hepatitis C virus (HCV) infection in hepatocytes, in which tumour-like glutamine metabolism is induced to create an environment that is favourable for viral replication [61]. The above findings may provide a reasonable explanation for the low performance of cattle with persistent infections. However, further in vivo and in vitro validation are required.

Innate immunity is the first line of defence against pathogen invasion. Interferon-stimulated genes (ISGs) are important components of the host innate immune system that are instrumental in the defence mechanism against viral infection [63]. ISGs that were found to be differentially expressed during BVDV infection are interesting. OAS1Y, Mx1, and ISG15, which have direct antiviral activity, were found to be downregulated during BVDV infection. 2′-5’oligoadenylate synthase (OAS), a member of the double-stranded RNA family, is an important antiviral protein induced by IFN and is involved in cellular innate immune defence [64]. OAS protein can be activated by replication intermediate dsRNA, and activated OAS can catalyse the synthesis of 2′-5′ oligoadenylate (2′-5’A) [65]. Next, RNAse L is activated by 2′-5’A. Viral RNA is degraded by RNAse L. In a previous study, dynamic changes in the negative-strand RNA of BVDV during replication in MDBK cells were detected by strand-specific reverse transcription real-time quantitative PCR (ssRT-qPCR). The findings showed that the abundance of the negative-strand RNA of BVDV initially increased between the 2 h and 6 h time points, decreased between the 6 h and 12 h time points, and then gradually increased between the 12 h and 36 h points. Peak levels of the negative-strand RNA of BVDV were observed between 36 h and 48 h after infection of MDBK cells with BVDV (Additional file 2: Figure S6). The observed downregulation of the negative-strand RNA levels between the 6 h and 12 h time points may be associated with the degradation of dsRNA by the machinery of the host innate immune defence system, such as the OAS/RNAse L system. The ISG15 protein, which is encoded by IFN-stimulated gene 15, is a ubiquitin-like modified protein. ISG15 protein and its ubiquitin-like modification system participate in the innate immune response and are crucial to facilitate the antiviral effects of interferons [66]. Mx proteins, interferon (IFN)-induced GTPases, which exert potent antiviral activity against various single-stranded RNA viruses, are members of dynamin superfamily. Mx proteins inhibit viral replication through diverse strategies. Mx protein can inhibit virus replication by lysing ribonucleoprotein complexes, isolating viral nucleocapsid proteins in large membrane-associated perinuclear complexes, and inhibiting viral primary transcription [67,68,69,70]. The downregulated expression of virus resistance genes could be part of the mechanisms by which BVDV escapes from the host immune defence system, which is mediated by inhibiting the production of IFN in host cells.

The recognition of PAMPs and the activation of signalling transduction pathways are the basis for the activation of host immune defence system. The molecular components that participate in signal transduction pathways that activate the host immune defence system are blocked by pathogens. Enrichment analysis of KEGG pathways of DEGs will help understand the signal transduction pathways that are activated or inhibited by BVDV infection and provide clues into the interactions between BVDV and host. KEGG pathway enrichment analysis revealed that the DEGs from all treatment groups were enriched in pathways involved in complement and coagulation cascades, TGF-beta, FoxO, and AGE-RAGE signalling pathways in diabetic complications, and ErbB and Hippo signalling pathway [71,72,73]. The complement system, because of its dual effector and priming functions, is a mediator of innate immunity and acts as a nonspecific defence mechanism against pathogens. The complement system is composed of more than 35 soluble plasma proteins that play essential roles in innate and adaptive immunity [74]. A total of 11 DEGs were identified in this study, namely, F3, C1R, KNG1, CLU, C3, FB, SERPINA5, SERPINE1, C1S, F2RL2, C, which were enriched in pathways related to the complement and coagulation cascades and were downregulated during BVDV infection. C1, which is composed of C1q, C1r, and C1s, is the protein complex that initiates the classical pathway of complement activation [75]. Complement factor B (FB) participates in the formation of C3 convertase [76, 77]. The downregulated expression of the genes C1 and FB can inhibit the activation of complement system. The above findings suggested that BVDV inhibits the complement system and plays an important role during the early stages of BVDV infection. The inhibition of complement activation during viral infection has been confirmed in previous studies. Complement activation is known to be inhibited by HCV NS3/4A protease [78]. In the acute stage of rotavirus virus (RV) infection, C3 expression is significantly downregulated; C3 levels are restored to normal levels during the convalescent stage, suggesting that C3 participates in inflammatory immunity in the early stage of RV infection [79]. Signalling pathways, including TGF-beta, FoxO, and AGE-RAGE signalling pathways in diabetic complications, and ErbB and the Hippo signalling pathways play important roles in the regulation of cell proliferation, differentiation, apoptosis, immune function, and the development of disease [80,81,82,83,84,85,86,87]. KEGG pathways, such as PI3K-Akt and NOD-like receptor signalling pathways, were significantly enriched in DEGs in the Mock vs. MBV2h after BVDV infection, as well as the pathways HIF-1, MAPK, TNF, and mTOR signalling pathway, which were significantly enriched in DEGs in Mock vs. MBV24h after BVDV infection, also regulate cell progression, including the recognition of pathogen-associated molecular patterns (PAMP), signalling transduction, the regulation of host immune system regulation [88,89,90,91,92]. Our results indicated that these signalling pathways are closely related to BVDV infection. Further investigation of these enriched signalling pathways will help further understand BVDV infection and pathogenesis, identify targets for drug development, and contribute to the prevention and treatment of bovine viral diarrhoea.

Conclusion

In present study, we characterized the transcriptome profile of MDBK cells infected with BVDV. The dynamic changes of differentially expressed genes of BVDV infection contributed to understanding the molecular mechanisms of BVDV-host interaction. Our findings suggested that BVDV-infection induced altering the host’s metabolic network, the inhibition of the expression of antiviral proteins and genes within the complement system might be contributed to BVDV proliferation. The above findings provided unique insights for further studies on the mechanisms underlying BVDV-host interactions and could serve as the basis for the prevention and treatment of bovine viral diarrhoea.

Methods

Ethic statement

No animal experiment was involved in this study.

Virus and cells

Madin-Darby bovine kidney (MDBK) cells were cultured in Dulbecco’s modified eagle medium (DMEM) (Gibco) supplemented with 8% fetal bovine serum (FBS) (Gibco). MDBK cells were maintained at 37 °C in an atmosphere with 5% CO2. BVDV BJ-2016 strain is a NCP BVDV isolated from the commercial bovine fetal serum. The virus titer was detected by indirect immunofluorescence assay (IFA).

Growth curves of BVDV

The growth curves of BVDV at different multiplicity of infection (MOI) were done. Briefly, MDBK cells (2 × 106 cells/well) were plated into 6-well cell culture plate and cultured in DMEM with 8% FBS. After 12 h, MDBK cells were inoculated with 0.1, 1, 10 MOI BVDV, respectively. Three biological repeats were set for each titer. The cell supernatants were collected at 2 h, 6 h, 12 h, 24 h, 36 h, 48 h post-infection, separately. The viral genome was quantified by reverse transcription real-time quantitative PCR (RT-qPCR) [93]. Number of viral genomes per milliliter of cell supernatant was calculated according to the standard curve.

Sample collection and RNA extraction

MDBK cells were plated onto T25 cell culture flasks and grown in DMEM containing 8% FBS. After 12 h, MDBK cells were inoculated with 10 MOI BVDV and maintained at 37 °C with 5% CO2 atmosphere for 1 h. Afterwards, cells were washed thrice with PBS, and grown in DMEM without serum. Samples were collected at 2, 6, 12, and 24 h post-infection based on the BVDV growth curves. The experiment was performed with three biological replicates for error reduction.

Total RNA was extracted using TRIzol® Reagent (Invitrogen) according to the manufacturer’s instructions. Then, the quality and integrity of the total RNA were assessed using an Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA) and 1.2% agarose gel electrophoresis; the RNA concentration was measured using a NanoDrop 2000 instrument (NanoDrop Technologies, Technologies, Wilmington, DE). High-quality RNA samples with OD260/280 ratios ranging from 1.8-2.2 and OD260/230 ≥ 2.0, RNA integrity number (RIN) ≥8, and total RNA concentration ≥ 50 ng/μL were used for library preparation.

Library preparation and Illumina NovaSeq 6000 sequencing

RNA-Seq libraries were prepared using the Illumina TruSeq™ RNA sample preparation Kit (Illumina, San Diego, CA). Briefly, oligo (dT) magnetic beads were used to enrich for mRNAs, which contain poly A tails, from the total RNA. The enriched mRNAs were cleaved into 300-bp fragments using fragmentation buffer. The cleaved mRNA fragments served as the templates for double-stranded cDNA (dscDNA) synthesis using the SuperScript dscDNA synthesis kit (Invitrogen, CA) following the manufacturer’s instructions. The synthesized dscDNA was subjected to end-repair and A-tailed. The indexed adapters were ligated to the A-tailed dscDNA. The dscDNA with indexed adapters were purified and enriched by PCR for 15 cycles. The quantity and quality of enriched dscDNA were assessed using TBS380 and NanoDrop 2000 spectrophotometer and Agilent 2100 Bioanalyzer. The libraries were then used for paired-end (PE) sequencing with the Illumina NovaSeq 6000 platform.

Quality control and read mapping

Quality control of the reads generated from RNA-Seq was performed using SeqPrep (https://github.com/jstjohn/SeqPrep) and Sickle (https://github.com/najoshi/sickle). Adapter and primer sequences were removed, and sequences with lengths below 20 bp were discarded. After removing adapter and primer sequences, low-quality bases were trimmed from the 3′ end of the reads. After trimming low-quality bases, sequences with quality values less than 10 were discarded. Sequences with N ratios higher than 10% were also removed. The error rate (%), Q20 and Q30 values, GC-content (%), and sequence duplication levels of the resulting high-quality clean reads were then evaluated.

The high-quality reads were mapped to the reference Bos taurus genome (UMD3.1) (http://www.ensembl.org/Bos_taurus/Info/Index) using Hisat2 (http://ccb.jhu.edu/software/hisat2/index.shtml) and the reference-based assembly of transcripts was performed using Stringtie (http://ccb.jhu.edu/software/stringtie/). By comparing with the reference Bos taurus genome, on the one hand, genes with annotation information can be identified in assembled genome. On the other hand, transcripts without annotated information can be obtained, which are defined as new transcripts. Novel transcripts and genes were predicted by gffcompare software. And then, the gene prediction was used for calculation of gene expression values.

Functional annotation and classification

For functional annotation and classification, all transcripts and their corresponding genes were compared with the Clusters of Orthologous Groups of proteins (COG, http://www.ncbi.nlm.nih.gov/COG/), Gene Ontology (GO, http://www.geneontology.org) and Kyoto Encyclopedia of Genes and Genomes (KEGG, http://www.genome.jp/kegg/) databases. GO analysis was conducted using the BLAST2GO software with default parameters. COG functional classification was conducted using Blastx software in the STRING database. KEGG pathway annotation was performed using KOBAS (http://kobas.cbi.pku.edu.cn/).

Differential expression analysis

Read counts of each transcript or gene were calculated using featureCounts software. Gene expression values were expressed as reads per kilobase of exon per million fragments mapped (FPKM) using featureCounts software. To identified true differentially expressed genes (DEGs), the false discovery rate (FDR) was used for the rectification of the p-values. The fold change (FC) of each transcript or gene in different groups was calculated. Transcripts or genes with FDR ≤ 0.05 and |Log2FC| ≥ 1 were considered as significant DEGs. DESeq2 was used to perform differential expression analysis of genes or transcripts.

The DEGs were subjected to enrichment analyses of GO and KEGG pathways. Fisher’s exact test and multiple correction method, including Bonferroni, Holm, Sidak, and false discovery rate, were used for the enrichment analyses to correct p-values and the false positive rate. Functional enrichment analysis was performed using GOatools. GO terms with corrected p-value ≤0.05 were considered as significantly enriched. Enrichment analysis of KEGG pathways was conducted using KOBAS software. KEGG pathways with corrected p-value ≤0.05 were considered as significantly enriched.

Protein-to-protein interaction network analyses of DEGs was performed using STRING database (http://string-db.org/) and visualized the protein-protein interaction network relationship using NetworkX within Python.

Validation of transcriptome sequencing results

Reverse transcription quantitative real-time PCR (RT-qPCR) was performed to validated the DEGs identified from transcriptome sequencing. Ten randomly DEGs, namely, C8orf4, PSPH, ISG15, EGR1, SIGLEC10, FABP3, GRIP2, GALNT18, GATM, IFITM1, were selected for RT-qPCR validation. The gees ACLY, ACSS2, INSIG1, HMGCR, HMGCS1, LPIN1, were additionally quantified by RT-qPCR for analyses of genes expression. Total RNA extraction was performed using TRIzol reagent following the manufacturer’s instructions. cDNA synthesis was conducted using the TransScript One-Step gDNA Removal and cDNA Synthesis SuperMix (TransGen Biotech, Beijing, China). Sequence-specific primers of the randomly selected genes were designed using Primer Express software (Table 2). ChamQ™ Universal SYBR qPCR Master Mix (Vazyme Biotech Co., Ltd) was used to perform RT-qPCR following the manufacturer’s instructions. RT-qPCR was performed in 20-μL reaction volumes containing 10 μL of SYBR qPCR Master Mix, 0.8 μL of upstream and downstream primers (10 μM), 1 μL of cDNA template, and 7.4 μL of ddH2O. The following reaction profile was used: 95 °C for 30 s, followed by 40 cycles of 95 °C for 5 s and 60 °C for 30 s; melting curve analysis was performed to validate specific amplification. β-actin gene was used as an endogenous reference gene. RT-qPCR was performed in a 96-well plate on an ABI QuantStudio 7 Flex real-time system (Applied Biosystems, Foster City, CA). The detection was performed in triplicate for each biological replicate. The relative expression values of selected genes were calculated using the 2-ΔΔCt method and normalized against the expression levels of the β-actin gene.

Table 2 Primers for RT-qPCR in this study

Availability of data and materials

The datasets generated and/or analyzed during the current study are available at NCBI project PRJNA562724 (https://www.ncbi.nlm.nih.gov/sra/?term=PRJNA562724) with accession number of 15 objects (SRR10031502, SRR10031501, SRR10031495, SRR10031494, SRR10031493, SRR10031492, SRR10031491, SRR10031490, SRR10031489, SRR10031488, SRR10031500, SRR10031499, SRR10031498, SRR10031497, SRR10031496). Any reasonable requests are available from the corresponding author.

Abbreviations

BVDV:

Bovine viral diarrhea virus

COG:

Clusters of Orthologous Groups of proteins

CP:

Cytopathogenic biotype

DEGs:

Differentially expressed genes

DMEM:

Dulbecco’s modified eagle medium

FBS:

Fetal bovine serum

FC:

Fold change

FPKM:

The reads per kilobase of exon per million fragments mapped

GO:

Gene ontology

IFA:

Indirect immunofluorescence assay

IFN:

Interferon

KEGG:

Kyoto Encyclopedia of Genes and Genomes

MDBK:

Madin-Darby bovine kidney cells

MOI:

Multiplicity of infection

NCP:

Non-cytopathogenic biotype

PAMP:

Pathogen-associated molecular pattern

PBS:

Phosphate buffer saline

PPI:

Protein to protein interaction network

RT-qPCR:

Reverse transcription real-time quantitative PCR

References

  1. 1.

    Ridpath JF. Bovine Viral Diarrhea Virus. Encyclopedia Virol. 2008;24:374–80.

    Article  Google Scholar 

  2. 2.

    Harasawa R, Giangaspero M, Ibata G, Paton DJ. Giraffe strain of pestivirus: its taxonomic status based on the 5′-untranslated region. Microbiol Immunol. 2000;44(11):915–21.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  3. 3.

    Becher P, Orlich M, Shannon AD, Horner G, König M, Thiel HJ. Phylogenetic analysis of pestiviruses from domestic and wild ruminants. J Gen Virol. 1997;78(Pt6):1357–66.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  4. 4.

    Deng ML, Ji SK, Fei WT, Raza S, He CF, Chen YY, Chen HC, Guo AZ. Prevalence study and genetic typing of bovine viral diarrhea virus (BVDV) in four bovine species in China. PLoS One. 2015;10(4):16.

    Google Scholar 

  5. 5.

    Deng Y, Sun CQ, Cao SJ, Lin T, Yuan SS, Zhang HB, Zhai SL, Huang L, Shan TL, Zheng H, et al. High prevalence of bovine viral diarrhea virus 1 in Chinese swine herds. Vet Microbiol. 2012;159(3-4):490–3.

    PubMed  Article  PubMed Central  Google Scholar 

  6. 6.

    Gao YG, Wang SJ, Du R, Wang QK, Sun CJ, Wang N, Zhang PJ, Zhang LX. Isolation and identification of a bovine viral diarrhea virus from sika deer in China. Virol J. 2011;8:6.

    Article  CAS  Google Scholar 

  7. 7.

    Milicevic V, Maksimovic-Zoric J, Veljovic L, Kureljusic B, Savic B, Cvetojevic D, Jezdimirovic N, Radosavljevic V. Bovine viral diarrhea virus infection in wild boar. Res Vet Sci. 2018;119:76–8.

    PubMed  Article  PubMed Central  Google Scholar 

  8. 8.

    Smith DB, Meyers G, Bukh J, Gould EA, Monath T, Scott Muerhoff A, Pletnev A, Rico-Hesse R, Stapleton JT, Simmonds P, et al. Proposed revision to the taxonomy of the genus Pestivirus, family Flaviviridae. J Gen Virol. 2017;98(8):2106–12.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  9. 9.

    Yesilbag K, Alpay G, Becher P. Variability and Global Distribution of Subgenotypes of Bovine Viral Diarrhea Virus. Viruses. 2017;9(6).

  10. 10.

    Brownlie J, Clarke MC, Howard CJ. Experimental infection of cattle in early pregnancy with a cytopathic strain of bovine virus diarrhoea virus. Res Vet Sci. 1989;46(3):307–11.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  11. 11.

    Hulst MM, Moorman RJ. Inhibition of pestivirus infection in cell culture by envelope proteins E(rns) and E2 of classical swine fever virus: E(rns) and E2 interact with different receptors. J Gen Virol. 1997;78:2779–87.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  12. 12.

    Iqbal M, Flick-Smith H, McCauley JW. Interactions of bovine viral diarrhea virus glycoprotein E(rns) with cell surface glycosaminoglycans. J Gen Virol. 2000;81:451–9.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  13. 13.

    Pavlovic D, Neville DC, Argaud O, Blumberg B, Dwek RA, Fischer WB, Zitzmann N. The hepatitis C virus p7 protein forms an ion channel that is inhibited by long alkyl-chain iminosugar derivatives. Proc Natn Acad Sci. 2003;100:6104–8.

    CAS  Article  Google Scholar 

  14. 14.

    Meyers G, Thiel H-J. Molecular characterization of pestiviruses. Adv Virus Res. 1996;47:53–118.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  15. 15.

    Warrilow D, Lott WB, Greive S, Gowans EJ. Properties of the bovine diarrhoea virus replicase in extracts of infected MDBK cells. Arch Virol. 2000;145:2163–71.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  16. 16.

    Kao SJ-H, CC. Characterization of RNA products associated with or aborted by a viral RNA-dependent RNA polymerase. Virology. 1997;236:348–53.

    PubMed  Article  PubMed Central  Google Scholar 

  17. 17.

    Grummer B, Beer M, Liebler-Tenorio I. Greiser-Wilke1. Localization of viral proteins in cells infected with bovine viral diarrhoea virus. J Gen Virol. 2001;82:2597–605.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  18. 18.

    Gray EW, Nettleton PF. The ultrastructure of cell cultures infected with border disease and bovine virus diarrhoea viruses. J Gen Virol. 1987;68:2339–46.

    PubMed  Article  PubMed Central  Google Scholar 

  19. 19.

    Bielefeldt-Ohmann H, Bloch B. Electron microscopic studies of bovine viral diarrhea virus in tissues of diseased calves and in cell cultures. Arch Virol. 1982;71:57–74.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  20. 20.

    Nuttall PA. Growth characteristics of two strains of bovine virus diarrhea virus. Arch Virol. 1980;66:365–9.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  21. 21.

    Peterhans E, Schweizer M. BVDV: a pestivirus inducing tolerance of the innate immune response. Biol. 2013;41(1):39–51.

    CAS  Article  Google Scholar 

  22. 22.

    Tortorella D, Gewurz BE, Furman MH, Schust DJ, Ploegh HL. Viral subversion of the immune system. Annu Rev Immunol. 2000;18:861–926.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  23. 23.

    Peterhans E, Jungi TW, Schweizer M. BVDV and innate immunity. Biol. 2003;31(2):107–12.

    CAS  Article  Google Scholar 

  24. 24.

    Peterhans E, Schweizer M. Pestiviruses: how to outmaneuver your hosts. Vet Microbiol. 2010;142(1-2):18–25.

    PubMed  Article  PubMed Central  Google Scholar 

  25. 25.

    Gil LH, Ansari IH, Vassilev V, Liang D, Lai VC, Zhong W, Hong Z, Dubovi EJ, Donis RO. The amino-terminal domain of bovine viral diarrhea virus Npro protein is necessary for alpha/beta interferon antagonism. J Virol. 2006;80(2):900–11.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  26. 26.

    Hilton L, Moganeradj K, Zhang G, Chen YH, Randall RE, McCauley JW, Goodbourn S. The NPro product of bovine viral diarrhea virus inhibits DNA binding by interferon regulatory factor 3 and targets it for proteasomal degradation. J Virol. 2006;80(23):11723–32.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  27. 27.

    Bendfeldt S, Ridpath JF, Neill JD. Activation of cell signaling pathways is dependant on the biotype of bovine viral diarrhea viruses type 2. Virus Res. 2007;126(1-2):96–105.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  28. 28.

    Luo X, Pan R, Wan C, Liu X, Wu J, Pan Z. Glycosylation of classical swine fever virus E(rns) is essential for binding double-stranded RNA and preventing interferon-beta induction. Virus Res. 2009;146(1-2):135–9.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  29. 29.

    Houe H. Epidemiological features and economical importance of bovine virus diarrhoea virus (BVDV) infections. Vet Microbiol. 1999;64(2–3):89–107.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  30. 30.

    Ridpath JF. Practical significance of heterogeneity among BVDV strains: impact of biotype and genotype on US control programs. Prev Vet Med. 2005;72(1-2):17–30.

    PubMed  Article  PubMed Central  Google Scholar 

  31. 31.

    Wang Z, Gerstein M, Snyder M. RNA-Seq: a revolutionary tool for transcriptomics. Nat Rev Genet. 2009;10(1):57–63.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  32. 32.

    Goodwin CM, Xu S, Munger J. Stealing the keys to the kitchen. Viral manipulation of the host cell metabolic network. Trends Microbiol. 2015;23(12):789–98.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  33. 33.

    Sanchez EL, Lagunoff M. Viral activation of cellular metabolism. Virol. 2015;479-480:609–18.

    CAS  Article  Google Scholar 

  34. 34.

    Shin J, MacCarthy T. Potential for evolution of complex defense strategies in a multi-scale model of virus-host coevolution. BMC Evol Biol. 2016;16(1):233.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  35. 35.

    Cvirkaite-Krupovic V, Carballido-Lopez R, Tavares P. Virus evolution toward limited dependence on nonessential functions of the host: the case of bacteriophage SPP1. J Virol. 2015;89(5):2875–83.

    CAS  PubMed  Article  Google Scholar 

  36. 36.

    Strating JR, van Kuppeveld FJ. Viral rewiring of cellular lipid metabolism to create membranous replication compartments. Curr Opin Cell Biol. 2017;47:24–33.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  37. 37.

    Peterfy M, Phan J, Xu P, Reue K. Lipodystrophy in the fld mouse results from mutation of a new gene encoding a nuclear protein, lipin. Nat Genet. 2001;27(1):121–4.

    CAS  PubMed  Article  Google Scholar 

  38. 38.

    Nadra K, de Preux Charles AS, Medard JJ, Hendriks WT, Han GS, Gres S, Carman GM, Saulnier-Blache JS, Verheijen MH, Chrast R. Phosphatidic acid mediates demyelination in Lpin1 mutant mice. Genes Dev. 2008;22(12):1647–61.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  39. 39.

    Suviolahti E, Reue K, Cantor RM, Phan J, Gentile M, Naukkarinen J, Soro-Paavonen A, Oksanen L, Kaprio J, Rissanen A, et al. Cross-species analyses implicate Lipin 1 involvement in human glucose metabolism. Hum Mol Genet. 2006;15(3):377–86.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  40. 40.

    Kang ES, Park SE, Han SJ, Kim SH, Nam CM, Ahn CW, Cha BS, Kim KS, Lee HC. LPIN1 genetic variation is associated with rosiglitazone response in type 2 diabetic patients. Mol Genet Metab. 2008;95(1-2):96–100.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  41. 41.

    Ryu D, Oh KJ, Jo HY, Hedrick S, Kim YN, Hwang YJ, Park TS, Han JS, Choi CS, Montminy M, et al. TORC2 regulates hepatic insulin signaling via a mammalian phosphatidic acid phosphatase, LIPIN1. Cell Metab. 2009;9(3):240–51.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  42. 42.

    He XP, Xu XW, Zhao SH, Fan B, Yu M, Zhu MJ, Li CC, Peng ZZ, Liu B. Investigation of Lpin1 as a candidate gene for fat deposition in pigs. Mol Biol Rep. 2009;36(5):1175–80.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  43. 43.

    Lixuan WEIZC, He C, Guo B, Lan X, Chen H, Lei C. Association between polymorphism of Lpin1 gene and Bovine economic traits. Acta Agriculturae Boreali Occidentalis Sinica. 2012;21(1):11–5.

    Google Scholar 

  44. 44.

    Liu Y, Liu XL, He H, Gu YL. Four SNPs of insulin-induced gene 1 associated with growth and carcass traits in Qinchuan cattle in China. Genet Mol Res. 2012;11(2):1209–16.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  45. 45.

    Yang T, Espenshade PJ, Wright ME, Yabe D, Gong Y, Aebersold R, Goldstein JL, Brown MS. Crucial step in cholesterol homeostasis: sterols promote binding of SCAP to INSIG-1, a membrane protein that facilitates retention of SREBPs in ER. Cell. 2002;110(4):489–500.

    CAS  PubMed  Article  Google Scholar 

  46. 46.

    Lee PC, DeBose-Boyd RA. Intramembrane glycine mediates multimerization of Insig-2, a requirement for sterol regulation in Chinese hamster ovary cells. J Lipid Res. 2010;51(1):192–201.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  47. 47.

    Gong Y, Lee JN, Lee PC, Goldstein JL, Brown MS, Ye J. Sterol-regulated ubiquitination and degradation of Insig-1 creates a convergent mechanism for feedback control of cholesterol synthesis and uptake. Cell Metab. 2006;3(1):15–24.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  48. 48.

    Sun J, Gao Y, Liu D, Ma W, Xue J, Zhang C, Lan X, Lei C, Chen H. Haplotype combination of the bovine INSIG1 gene sequence variants and association with growth traits in Nanyang cattle. Genome. 2012;55(6):429–36.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  49. 49.

    Li MN, Guo X, Bao PJ, Wu XY, Ding XZ, Chu M, Liang CN, Yan P. Association of genetic variations in the ACLY gene with growth traits in Chinese beef cattle. Genet Mol Res. 2016;15(2). https://doi.org/10.4238/gmr.15028250.

  50. 50.

    Wakil SJ, Abu-Elheiga LA. Fatty acid metabolism: target for metabolic syndrome. J Lipid Res. 2009;50(Suppl):S138–43.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  51. 51.

    Sharpe LJ, Brown AJ. Controlling cholesterol synthesis beyond 3-hydroxy-3-methylglutaryl-CoA reductase (HMGCR). J Biol Chem. 2013;288(26):18707–15.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  52. 52.

    Hegardt FG. Mitochondrial 3-hydroxy-3-methylglutaryl-CoA synthase: a control enzyme in ketogenesis. Biochem J. 1999;338(Pt 3):569–82.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  53. 53.

    Xu H, Luo J, Ma G, Zhang X, Yao D, Li M, Loor JJ. Acyl-CoA synthetase short-chain family member 2 (ACSS2) is regulated by SREBP-1 and plays a role in fatty acid synthesis in caprine mammary epithelial cells. J Cell Physiol. 2018;233(2):1005–16.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  54. 54.

    Eberle D, Hegarty B, Bossard P, Ferre P, Foufelle F. SREBP transcription factors: master regulators of lipid homeostasis. Biochimie. 2004;86(11):839–48.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  55. 55.

    Harvatine KJ, Bauman DE. SREBP1 and thyroid hormone responsive spot 14 (S14) are involved in the regulation of bovine mammary lipid synthesis during diet-induced milk fat depression and treatment with CLA. J Nutr. 2006;136(10):2468–74.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  56. 56.

    Jerez-Timaure NC, Eisen EJ, Pomp D. Fine mapping of a QTL region with large effects on growth and fatness on mouse chromosome 2. Physiol Genomics. 2005;21(3):411–22.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  57. 57.

    Sanchez-Solana B, Li DQ, Kumar R. Cytosolic functions of MORC2 in lipogenesis and adipogenesis. Biochim Biophys Acta. 2014;1843(2):316–26.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  58. 58.

    Yang X, Ouyang H, Chen F, Pang D, Dong M, Yang S, Liu X, Peng Z, Wang F, Zhang X, et al. HMG-CoA reductase is negatively associated with PCV2 infection and PCV2-induced apoptotic cell death. J Gen Virol. 2014;95(Pt 6):1330–7.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  59. 59.

    Soto-Acosta R, Mosso C, Cervantes-Salazar M, Puerta-Guardo H, Medina F, Favari L, Ludert JE, del Angel RM. The increase in cholesterol levels at early stages after dengue virus infection correlates with an augment in LDL particle uptake and HMG-CoA reductase activity. Virol. 2013;442(2):132–47.

    CAS  Article  Google Scholar 

  60. 60.

    Desplanques AS, Pontes M, De Corte N, Verheyen N, Nauwynck HJ, Vercauteren D, Favoreel HW. Cholesterol depletion affects infectivity and stability of pseudorabies virus. Virus Res. 2010;152(1-2):180–3.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  61. 61.

    Levy PL, Duponchel S, Eischeid H, Molle J, Michelet M, Diserens G, Vermathen M, Vermathen P, Dufour JF, Dienes HP, et al. Hepatitis C virus infection triggers a tumor-like glutamine metabolism. Hepatology. 2017;65(3):789–803.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  62. 62.

    Munger J, Bennett BD, Parikh A, Feng XJ, McArdle J, Rabitz HA, Shenk T, Rabinowitz JD. Systems-level metabolic flux profiling identifies fatty acid synthesis as a target for antiviral therapy. Nat Biotechnol. 2008;26(10):1179–86.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  63. 63.

    Raniga K, Liang C. Interferons: Reprogramming the Metabolic Network against Viral Infection. Viruses. 2018;10(1).

  64. 64.

    Sawai H, Hirano A, Mori H, Shinozuka K, Dong B, Silverman RH. Synthesis, characterization, and biological properties of 8-azido- and 8-amino-substituted 2′,5′-oligoadenylates. J Med Chem. 2003;46(23):4926–32.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  65. 65.

    Pari M, Kuusksalu A, Lopp A, Reintamm T, Justesen J, Kelve M. Expression and characterization of recombinant 2′,5′-oligoadenylate synthetase from the marine sponge Geodia cydonium. FEBS J. 2007;274(13):3462–74.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  66. 66.

    Dao CT, Zhang DE. ISG15: a ubiquitin-like enigma. Front Biosci. 2005;10:2701–22.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  67. 67.

    Judith V, Eef P, Bert S, Walter F, Xavier S. Interferon-inducible protein Mx1 inhibits influenza virus by interfering with functional viral ribonucleoprotein complex assembly. J Virol. 2012;86(24):13445–55.

    Article  CAS  Google Scholar 

  68. 68.

    Pavlovic J, Haller O, Staeheli P. Human and mouse mx proteins inhibit different steps of the influenza virus multiplication cycle. J Virol. 1992;66(4):2564–9.

    CAS  PubMed  PubMed Central  Google Scholar 

  69. 69.

    Georg K, Christian J, Heinz H, Otto H. Antivirally active MxA protein sequesters La Crosse virus nucleocapsid protein into perinuclear complexes. Proc Natl Acad Sci U S A. 2002;99(5):3153–8.

    Article  CAS  Google Scholar 

  70. 70.

    Reichelt M, Stertz S, Krijnse-Locker J, Haller O, Kochs G. Missorting of LaCrosse virus nucleocapsid protein by the interferon-induced MxA GTPase involves smooth ER membranes. Traffic. 2010;5(10):772–84.

    Article  CAS  Google Scholar 

  71. 71.

    Datta PK, Rappaport J. HIV and complement: hijacking an immune defense. Biomed Pharmacother. 2006;60(9):561–8.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  72. 72.

    Hovhannisyan LP, Mkrtchyan GM, Sukiasian SH, Boyajyan AS. Alterations in the complement cascade in post-traumatic stress disorder. Allergy Asthma Clin Immunol. 2010;6(1):3.

    PubMed  PubMed Central  Article  Google Scholar 

  73. 73.

    Carroll MC. Complement and humoral immunity. Vaccine. 2008;26(Suppl 8):I28–33.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  74. 74.

    Ricklin D, Hajishengallis G, Yang K, Lambris JD. Complement: a key system for immune surveillance and homeostasis. Nat Immunol. 2010;11(9):785–97.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  75. 75.

    Turner MW. The role of mannose-binding lectin in health and disease. Mol Immunol. 2003;40(7):423–9.

    CAS  PubMed  Article  Google Scholar 

  76. 76.

    Flyvbjerg A. The role of the complement system in diabetic nephropathy. Nat Rev Nephrol. 2017;13(5):311–8.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  77. 77.

    Assogba FG, Couchoud C, Hannedouche T, Villar E, Frimat L, Fagot-Campagna A, Jacquelinet C, Stengel B. Trends in the epidemiology and care of diabetes mellitus-related end-stage renal disease in France, 2007-2011. Diabetologia. 2014;57(4):718–28.

    PubMed  Article  PubMed Central  Google Scholar 

  78. 78.

    Mawatari S, Uto H, Ido A, Nakashima K, Suzuki T, Kanmura S, Kumagai K, Oda K, Tabu K, Tamai T, et al. Hepatitis C virus NS3/4A protease inhibits complement activation by cleaving complement component 4. PLoS One. 2013;8(12):e82094.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  79. 79.

    Sichun Yin HY, Li F, Luo B, Liang P, Hua L, Ye G. Dynamic changes of serum immunoglobulin and complement in children with rotavirus enteritis and their clinical significance. J Chinese Physician. 2006;8(6):834–5.

    Google Scholar 

  80. 80.

    Brown AK, Webb AE. Regulation of FOXO factors in mammalian cells. Curr Top Dev Biol. 2018;127:165–92.

    PubMed  Article  PubMed Central  Google Scholar 

  81. 81.

    Chen W, Wahl SM. TGF-beta: receptors, signaling pathways and autoimmunity. Curr Dir Autoimmun. 2002;5:62–91.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  82. 82.

    Akhurst RJ, Hata A. Targeting the TGFbeta signalling pathway in disease. Nat Rev Drug Discov. 2012;11(10):790–811.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  83. 83.

    Ramezani A, Nikravesh H, Faghihloo E. The roles of FOX proteins in virus-associated cancers. J Cell Physiol. 2019;234(4):3347–61.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  84. 84.

    Grant S, Qiao L, Dent P. Roles of ERBB family receptor tyrosine kinases, and downstream signaling pathways, in the control of cell growth and survival. Front Biosci. 2002;7:d376–89.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  85. 85.

    Sanajou D, Ghorbani Haghjo A, Argani H, Aslani S. AGE-RAGE axis blockade in diabetic nephropathy: current status and future directions. Eur J Pharmacol. 2018;833:158–64.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  86. 86.

    Pan D. The hippo signaling pathway in development and cancer. Dev Cell. 2010;19(4):491–505.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  87. 87.

    Huang J, Wu S, Barrera J, Matthews K, Pan D. The hippo signaling pathway coordinately regulates cell proliferation and apoptosis by inactivating Yorkie, the Drosophila homolog of YAP. Cell. 2005;122(3):421–34.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  88. 88.

    Kim YK, Shin JS, Nahm MH. NOD-like receptors in infection, immunity, and diseases. Yonsei Med J. 2016;57(1):5–14.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  89. 89.

    Altman AM, Mahmud J, Nikolovska-Coleska Z, Chan G. HCMV modulation of cellular PI3K/AKT/mTOR signaling: new opportunities for therapeutic intervention? Antivir Res. 2019;163:82–90.

    CAS  PubMed  Article  Google Scholar 

  90. 90.

    Charpentier T, Hammami A, Stager S. Hypoxia inducible factor 1alpha: a critical factor for the immune response to pathogens and Leishmania. Cell Immunol. 2016;309:42–9.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  91. 91.

    Sun Y, Liu WZ, Liu T, Feng X, Yang N, Zhou HF. Signaling pathway of MAPK/ERK in cell proliferation, differentiation, migration, senescence and apoptosis. J Recept Signal Transduct Res. 2015;35(6):600–4.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  92. 92.

    Chen G, Goeddel DV. TNF-R1 signaling: a beautiful pathway. Science. 2002;296(5573):1634–5.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  93. 93.

    Liu C, Cui SJ. Establishment and application of a SYBR green real-time PCR assay for detection of bovine viral diarrhea virus. J Mod Husbandry Vet Anim. 2018;9:9–13.

    Google Scholar 

Download references

Acknowledgements

We thank CACTUS for language revision of this manuscript and Shanghai Majorbio Bio-pharm Technology Co., Ltd. for their technical support.

Funding

This research was supported by the Agricultural Science and Technology Innovation Program (ASTIP-IAS15).

Author information

Affiliations

Authors

Contributions

YMZ and SJC designed and supervised the study; CL wrote the manuscript. CL, and LL performed the experimental and analyzed the data. YHL analyzed the data and revised the manuscript. All authors have reviewed and approved the final manuscript.

Corresponding authors

Correspondence to Shangjin Cui or Yanming Zhang.

Ethics declarations

Ethics approval and consent to participate

No applicable.

Consent for publication

No applicable.

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary information

Additional file 1: Table S1.

Length distribution of transcripts. Table S2. The annotation of novel transcripts/genes compared with NR, Swiss-Prot, and Pfam databases. Table S3. Nodes centrality analysis of PPI network of group Mock vs. MBV2h. Table S4. Nodes centrality analysis of PPI network of group Mock vs. MBV24h. Table S5. Nodes centrality analysis of PPI network of group MBV2h vs. MBV6h.

Additional file 2: Figure S1.

Venn analysis of functional annotation and classification within GO, COG and KEGG database. Figure S2. Functional annotation of the genes obtained from RNA-Seq. The different colours of the pie chart represent different GO terms, and the area represents the relative proportion of genes/transcripts in the GO Term. The number of genes were shown behind the subsets. Go terms of level 2 were annotated within the three Go categories, biological progress (BP), cellular progress (CC), and molecular function (MF). The functional annotation was performed using Blast2go. Figure S3. Function classification in Clusters of Orthologous Groups of Proteins (COG) of genes from RNA-Seq. Capital letters on x axis indicated the COG categories as listed on the right of histogram; y axis indicated the number of gene. Figure S4. Pathway annotation of genes from RNA-Seq. KEGG Pathway of level 2 within the seven categories of KEGG pathway as listed on the right of histogram were annotated. x axis indicated the number of genes annotated to the pathway; y axis indicated the name of the KEGG pathway. Figure S5. Venn analysis of differential expression genes from group Mock vs. MBV2h, Mock vs. MBV6h, Mock vs. MBV12h, Mock vs. MBV24h. Figure S6. Dynamic changes of BVDV RNA in MDBK cells infected with 10 MOI BVDV. a: BVDV positive strand RNA; b: BVDV negative strand RNA. Copies of BVDV positive or negative strand RNA was quantified by strand specific SYBR Green real-time PCR and shown as copies of positive or negative strand RNA of BVDV per microgram total RNA of cells. Figure S7. GO enrichment analysis of DEGs in comparison groups MBV2h vs. MBV6h(a), MBV6h vs. MBV 12 h(b), MBV12h vs. MBV24h(c). GO terms are on the x axis. Enrichment ratio of genes shown as GO terms for BP, CC, and MF. * means GO terms with significant enrichment. Figure S8. KEGG Pathway enrichment analysis of DEGs in comparison groups MBV2h vs. MBV6h (a), MBV6h vs. MBV12h (b), MBV12h vs. MBV24h (c). The name of KEGG pathway are on the x axis. Enrichment ratio of genes shown as the name of KEGG pathway for seven categories. * means KEGG Pathway with significant enrichment. Figure S9. Protein–Protein Interaction Network Analysis of DEGs with FDR ≤ 0.05 and |Log2FC| ≥ 2 in comparison group Mock vs. MBV2h. OAS1Y, FOS and EGR1 were three genes with the most node degree in the analysis, labelled in red blot, green blot and yellow blot, respectively. Figure S10. Protein–Protein Interaction Network Analysis of DEGs with FDR ≤ 0.05 and |Log2FC| ≥ 2 in comparison group Mock vs. MBV24h. OAS1Y and FOSB were the two genes with the most node degree, labeled in red blot and yellow blot, respectively.

Additional file 3.

DEGs identified in comparison group Mock vs. MBV2h.

Additional file 4.

DEGs identified in comparison group Mock vs. MBV6h.

Additional file 5.

DEGs identified in comparison group Mock vs. MBV12h.

Additional file 6.

DEGs identified in comparison group Mock vs. MBV24h.

Additional file 7.

DEGs identified in comparison group MBV2h vs. MBV6h.

Additional file 8.

DEGs identified in comparison group MBV6h vs. MBV12h.

Additional file 9.

DEGs identified in comparison group MBV12h vs. MBV24h.

Additional file 10.

Detail information of GO enrichment of DEGs in mock vs. MBV12h.

Additional file 11.

Detail information of GO enrichment of DEGs in mock vs. MBV24h.

Additional file 12.

Detail information of KEGG pathway enrichment of DEGs in mock vs. MBV12h.

Additional file 13.

Detail information of KEGG pathway enrichment of DEGs in mock vs. MBV24h.

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Cite this article

Liu, C., Liu, Y., Liang, L. et al. RNA-Seq based transcriptome analysis during bovine viral diarrhoea virus (BVDV) infection. BMC Genomics 20, 774 (2019). https://doi.org/10.1186/s12864-019-6120-4

Download citation

Keywords

  • Bovine viral diarrhoea virus
  • RNA-Seq
  • Transcriptome analysis
  • Pathogen-host interactions