Construction and Analysis of ceRNA Networks Reveal the Key Genes Associated with Bovine Herpesvirus Type 1 Infection

Introduction

BoHV-1, Alphaherpesvirinae subfamily of the family Herpesviridae¹ with double-stranded DNA in the capsule, is a primary virulent pathogen of cattle diseases.² The genome of BoHV-1 is about 135 KB, its diameter of virus particles is about 150 to 220 nm, enabling it to invade mucosal epithelial cells of the upper respiratory tract and those of the reproductive tract of cattle and proliferate in large numbers at the site of infection, resulting in cytopathic effects (CPE).³ BoHV-1 mainly infects cattle and causes respiratory and reproductive diseases and is associated with other clinical symptoms and systemic diseases of newborn calves,⁴ which causes significant economic losses to beef and dairy cattle industry all over the world.

Competing endogenous RNA (ceRNA) is the RNA in the complex transcriptional regulation network of an organism, which mainly includes messenger RNA (mRNA), long-chain non-coding RNA (lncRNA), microRNA (miRNA), circular RNA (circRNA) and pseudogenes.^5,6 The regions of these RNAs that can be bound by miRNAs are called microRNA response elements (MRE).⁷ These RNAs compete through common MREs to bind common miRNAs, affecting and regulating the expression of target gene transcripts. In organisms, these ceRNA interactions play important biological functions, such as lncRNA five prime to Xist (LncRNA FTX) can regulate the progression of colorectal cancer (CRC) through miR-192-5p/EIF5A2 Axis.⁸ LncRNA34977 is a lncRNA that is highly expressed in canine mammary tumors. It can promote progression of canine mammary tumors through regulating the expression of miR-8881/ELAVL4,⁹ LncRNA DANCR functions as a miR-33b sponge to positively regulate proliferation and metastasis in pancreatic cancer.¹⁰ Suman Ghosal et al¹¹ identified a large number of non-coding transcripts as immune response-related genes through the ceRNA network. Rahmani-Kukia N and Abbasi A¹² identified the potential application of circRNAs as an antiviral vaccine such as SARS-CoV-2 through ceRNA network analysis.

In recent years, RNA sequencing (RNA-Seq) based on next-generation sequencing technology has become an effective tool for screening genomic variation associated with different phenotypic traits and identifying molecular mechanisms of viral infection,¹³ but to the best of our knowledge, there is no report on ceRNA networks in BoHV-1 virus infected cells. MDBK cells were infected and sequenced by next-generation sequencing in this study. lncRNA-miRNA-mRNA networks were constructed using these high-throughput sequencing data to analyze the potential specific molecular functions and mechanisms. Then, function enrichment analyses and annotation were used to explore the importance of ncRNA ceRNAs in MDBK cells of BoHV-1 infection. From this, ceRNA function in virus-infected cells was analyzed and this enhanced the understanding of ncRNA regulatory networks.

Materials and Methods

Cells and Viruses

The Madin-Darby bovine kidney (MDBK) cell line (NBL-1) (ATCC CCL-22) was stored in the Cell Biology Laboratory of Heilongjiang Bayi Agricultural University. Bovine herpesvirus 1 (BoHV-1) was presented by Professor Zhu Zhanbo of Heilongjiang Bayi Agricultural University and stored in the Cell Biology Laboratory of Heilongjiang Bayi Agricultural University.

Cell Culture

The frozen MDBK cells (13 passages) were removed from liquid nitrogen and seeded onto T25 breathable cell culture flask (NEST, China) containing Dulbecco’s modified eagle medium with 10% Fetal bovine serum (Gibco, USA) at 37°C in 5% CO₂ atmosphere for 24 to 48 h in the cell incubator. When the cell confluence reached 80 to 90%, the cells were digested with Trypsin-EDTA (0.25%) without phenol red (Solarbio, China), then passaged with the ratio of 1:2 in 10 cm cell culture dishes.

When the concentration of cells reached 1 × 10⁷ cells per dish, 1 MOI of BoHV-1 was inoculated to MDBK cells and cultured in DMEM media (Gibco, USA) without serum at 37°C in 5% CO2 atmosphere. Cytopathic effects (CPE) were observed by light microscopy when virus-infected cells were cultured for 18, 24, and 33 h. The samples were harvested when the CPE reached varying degrees of monolayer at 18, 24 and 33 h. The sample from each dish was washed three times with Hank’s solution (Gibco, USA) and digested with 1 mL TRIzol reagent (Life, USA). The sample was then collected into nuclease-free cell cryopreservation tubes and sealed with All-Purpose Parafilm M Laboratory Film (Parafilm, USA) before being stored in liquid nitrogen for quick freezing. Finally, the samples were sequenced by Guangzhou Genedenovo Biotechnology Co., Ltd (Guangzhou, China).

Next Generation Sequencing for Cell Samples

The extraction of cell total RNA refers to the Experimental methods of Yongfei Bai et al¹⁴ and the quality of the RNA was assessed by an Agilent 2100 Bioanalyzer and agarose gel electrophoresis without RNA enzyme. Illumina HiSeqTM 4000 platform was used for PCR amplification and sequencing by Guangzhou Genedenovo Biotechnology Co., Ltd (Guangzhou, China). Specific methods.

Data Processing

Data processing methods are mainly based on the methods of Chen S and Langmead B et al, including the filtering of read by fastp¹⁵ (version 0.18.0) and the reads were compared to the ribosomal RNA (rRNA) database using the short read comparison tool Bowtie2¹⁶ (version 2.2.8). The alignment with the reference genome was determined by using HISAT2 (version 2.1.0) and detailed steps refer to the methods of Kim D et al.¹⁷ Transcripts were reconstructed using the Stringtie software (version 1.3.4).^18,19 Stringtie software was used to analyze transcript abundances in a reference-based approach. The edgeR package was used to analyze the differentially expressed transcripts between the two groups (http://www.bioconductor.org/packages/release/bioc/html/edgeR.html)²⁰. The mRNA and lncRNA were identified with a fold change greater than or equal to two and a false discovery rate (FDR) of less than 0.05 in a comparison of significantly differentially expressed genes (DEGs), and miRNA with a fold change greater than or equal to two and p < 0.05. The genes that had been rising and falling were selected in the experiment to prepare for later analysis. Mireap (Version 0.20), miRanda (Version 3.3a) and TargetScan (Version 7.0) were used for miRNA target prediction, and miRNA sequences and family information were obtained from the TargetScan website (http://www.targetscan.org/).²¹ The results were visualized by using R (version 3.6.1) software with Pheatmap package.²²

Construction of the lncRNA-miRNA-mRNA Network

All co-expressed competing triples identified based on ceRNA theory were used to construct the mRNA-lncRNA-miRNA network and visualized using Cytoscape software (version 3.6.0).^23,24 The networks were analyzed by using “Network Analyzer” which is a plug-in unit of Cytoscape software.²⁵

Functional Enrichment Analysis

To assess the function of differentially expressed genes, GO biological processes term and KEGG pathway analyses were performed by using mRNAs in the ceRNA network. All ceRNAs were mapped to GO terms in the gene ontology database (http://www.geneontology.org/), gene numbers were calculated according to the method of Zhan S and Yu.^26,27 The results were visualized by R (version 3.6.1) software with ggplot 2 package.²⁷

Validation of Differentially Expressed lncRNAs by RT−qPCR

In order to verify the accuracy of sequencing data, RT−qPCR was used to analyze DEGs. Twelve DEGs, namely, 6up, 7up, 12up, 15up, 16up, 17up, 19up, 20up, Hub1, 3D, 4D and 5D (Table 1) were randomly selected for RT-qPCR validation. The total RNA 20 ng was used to carry out reverse transcription and RT-qPCR reactions. Some differentially expressed lncRNA genes were randomly selected, and their primers (Table 1) were designed using Primer Express software. Ubiquitin C-terminal Hydrolase L5 (UCHL5) gene was used as a reference gene.

Table 1 The Primer Sequences of RT-qPCR

Results

States of MDBK Cells by BoHV-1 Infection at Different Points in Time

To understand the morphology of MDBK cells which infected by BoHV-1 virus, the morphology of MDBK cells was observed as uninfected and then, post-infection at 18, 24, and 33 h (Figure 1). In the uninfected group, the edges of MDBK cells were evident, and there was no cell disruption and necrosis present (Figure 1a). At 24 h post-infection, the cell boundary became blurred, appeared necrotic and plaques appeared (Figure 1c) and at 33 h post-infection, the plaques became bigger and cell necrosis was severe (Figure 1d).

Figure 1 States of MDBK cells by BoHV-1 infection at different points in time. (a) MDBK cells uninfected with BoHV-1 were used as a negative control; (b) MDBK cells post-infection 18 h, characterized by focal rounding; (c) MDBK cells post-infection 24 h, detachment of cells from the surface of the cell-culture dish; (d) MDBK cells post-infection 33 h, aggregation of the cells into vacuoles. All cell morphology under the light microscope in multiple rate 100-fold. Arrows indicate the state of cytopathic effects (CPE).

Identification of Differentially Expressed mRNAs, miRNAs, and lncRNAs

We conducted differential expression analysis for each two groups, aiming to observe the differential expression of lncRNA, micRNA and mRNA at different periods. Analysis of differential expression transcripts showed that 984 lncRNAs, 4584 mRNAs, and 374 miRNAs were differentially expressed in this experiment (Figure 2). Moreover, the number of differentially expressed RNAs in the 33-hour group was higher than that in any other group, indicating that the number of differentially expressed RNAs increased with the time of virus infection. These changed RNAs may be involved in virus–host interactions. lncRNAs, mRNAs, and miRNAs with significantly different expression levels in this study were selected from their expression plot on heat maps (Figure 3). From the plot, 68 lncRNAs (46 up-regulated and 22 down-regulated), 162 mRNAs (89 up-regulated and 73 down-regulated) and 43 miRNAs (26 up-regulated and 17 down-regulated) were found. These lncRNAs, mRNAs and miRNAs were continuously up-regulated or down-regulated and conducted for follow-up analysis. The horizontal coordinate indicates that samples clustered into one class are more similar, and the vertical coordinate indicates clusters of different RNAs that have the same or similar functions.

Figure 2 Differential expression analysis of lncRNA, mRNA and miRNA in every two groups The horizontal coordinate is the number of differentially expressed genes, The ordinates represent different groups. Red represents up-regulated genes and blue represents down-regulated genes. The differential expression analysis in two groups was taken by method as follow: lncRNA (a) and mRNA (b) with a fold change ≥ 2 and a false discovery rate (FDR) < 0.05 were chosen, and miRNA (c) with a fold change ≥ 2 and p < 0.05, and visualized by R software with ggplot2 package. B: MDBK cells uninfected with BoHV-1 were used as a negative control; Bo18h:MDBK cells post-infection 18 h; Bo24h:MDBK cells post-infection 24 h; Bo33h:MDBK cells post-infection 33 h.

Figure 3 Heat map analysis of differential expression of lncRNA, mRNA and miRNA for each sample The differential expression heat map in every sample was visualized by R software with a pheatmap package. (a) Heat map of differentially expressed lncRNAs; (b) Heat map of differentially expressed mRNAs; (c) Heat map of differentially expressed miRNAs. Each column represents one sample, and Each row represents a gene. The color legend is on the top-right of the figure, and compound relative abundances were standardized prior to unsupervised hierarchical clustering of samples (rows). Red indicates the expression of genes is up-regulated. Blue indicates the expression of genes is down-regulated. The darker the color, the more up-regulated the gene expression. The color legend named group indicates the sample group.

Construction and Visualization of mRNA-lncRNA-miRNA Networks

All co-expressed competing triples identified based on ceRNA theory were used to construct the mRNA-lncRNA-miRNA network and visualized using Cytoscape software (version 3.6.0). In this study, up-regulated lncRNAs and down-regulated miRNAs and up-regulated mRNA were integrated into the ceRNA network, down-regulated lncRNAs and up-regulated miRNAs and down-regulated mRNA were integrated into the ceRNA network. In up-regulated lncRNAs network (Figure 4a), 46 lncRNA nodes, 89 mRNA nodes, 26 miRNA nodes, and 500 edges were present in this mRNA-lncRNA-miRNA network. In down-regulated lncRNAs network (Figure 5a), 22 lncRNAs, 73 mRNAs, 17 miRNAs, and 255 edges were present in this mRNA-lncRNA-miRNA network. After the network analyzer calculation, the centrality indices of “betweenness centrality” (BC) and “closeness centrality” (CC) were used to determine the importance of individual node with basic functions in complex networks. Through the network analysis, significant differences (p < 0.05) were found in BC and CC among lncRNAs, mRNAs and miRNAs (Figures 4b, c and 5b, c), which showed that lncRNA played a very important roles in mRNA-lncRNA-miRNA networks.

Figure 4 Network of up-regulated lncRNAs with down-regulated miRNAs and the up-regulated mRNAs. (a) The mRNA-lncRNA-miRNA network in upregulated lncRNAs. Every node represents one gene, and Each line represents an interaction between two genes. The red triangle is used to represent lncRNAs, the green circle is used to represent miRNAs and the purple square is used to represent mRNAs. The higher the correlation, the larger the size. (b) The difference of betweenness centrality among lncRNAs, miRNAs, and mRNAs. In this network, betweenness centrality of the lncRNA nodes was higher than that of miRNA and mRNA nodes. (c) The difference of closeness centrality among lncRNAs, miRNAs, and mRNAs. In this network, the closeness centrality of lncRNA nodes was higher than that of miRNA and mRNA nodes.

Figure 5 Network of down-regulated lncRNAs with up-regulated miRNAs and down-regulated mRNA. (a) The mRNA-lncRNA-miRNA network in downregulated lncRNAs. Every node represents one gene, and Each line represents an interaction between two genes. The red triangle is used to represent lncRNAs, the green circle is used to represent miRNAs and the purple square is used to represent mRNAs. The higher the correlation, the larger the size. (b) The difference of betweenness centrality among lncRNAs, miRNAs, and mRNAs. In this network, betweenness centrality of the lncRNA nodes was higher than that of miRNA and mRNA nodes. (c) The difference of closeness centrality among lncRNAs, miRNAs, and mRNAs. In this network, the closeness centrality of lncRNA nodes was higher than that of miRNA and mRNA nodes.

Functional Enrichment Analysis of lncRNA-mRNA-miRNA Networks

To assess functional enrichment, GO and KEGG pathway analyses were performed by using DEGs in the ceRNA network. The top 20 significant GO terms of every group were taken and shown in Figure 6A and b. On the molecular function level, these genes were mostly enriched in enzyme binding, zinc ion binding, transporter activity, transmembrane transporter activity, substrate-specific channel activity, passive transmembrane transporter activity, channel activity and inorganic molecular entity in up-regulated group (Figure 6a). In down-regulated group, these genes were mostly enriched in DNA binding, nuclease activity, catalytic activity, acting on DNA and hydrolase activity, and acting on ester bonds (Figure 6b). The top 20 significant KEGG pathways of DEGs were taken to visualize by ggplot 2 package (Figure 7a and b). These DEGs were mostly enriched in proteoglycans in cancer, MAPK signaling pathway, cytokine-cytokine receptor interaction and viral carcinogenesis in up-regulated group (Figure 7a). In down-regulated group, the DEGs were mostly enriched in Herpes simplex virus 1 infection, microRNAs in cancer, TGF-beta signaling pathway and Thyroid hormone signaling pathway (Figure 7b).

Figure 6 Gene ontology (GO) enrichment analysis of lncRNA-mRNA-miRNA networks. (a) Top 20 significant GO terms on the molecular function level in up-regulated group; (b) top 20 significant GO terms on the molecular function level in down-regulated group; The count and q value were used to measure the size and color of bubbles, indicating the amount and significance of DEGs enriched in these terms.

Figure 7 Kyoto Encyclopedia of Genes and Genome (KEGG) enrichment analysis of lncRNA-mRNA-miRNA networks. lncRNA-mRNA-miRNA networks. (a) Top 20 significant KEGG pathways in up-regulated group. (b) Top 20 significant KEGG pathways in down-regulated group. The count and q value were used to measure the size and color of bubbles, indicating the amount and significance of DEGs enriched in these terms.

RT-qPCR Analysis

To verify the accuracy and reproducibility of transcriptome sequencing data, twelve differential expression genes (6up, 7up, 12up, 15up, 16up, 17up, 19up, 20up, Hub1, 3D, 4D and 5D) were randomly selected for RT-qPCR analysis (Figure 8a and b). The results were consistent with the transcriptome sequencing data, suggesting transcriptome sequencing were reliable. The correlation of RNA-seq and RT-qPCR was measured. As shown in Figure 9, the results of RT-qPCR showed a high correlation with RNA-seq, and the correlation coefficient (R2) was as high as 0.8274. These results indicated that most of the identified lncRNAs was truly expressed in BoHV-1 infected MDBK cells.

Figure 8 Expression level of randomly selected gene validated by RT-qPCR. The relative expression of 12 randomly selected lncRNA genes which were up-regulated (a) and down-regulated in RNA-seq (b) was performed by RT-qPCR. UCHL5 gene was used as internal control, and the comparative 2-ΔΔCt method was used to calculate the relative expression of each gene (fold change). The relative expression values shown are the same as the mean SD values. ⁎p < 0.05, ⁎⁎p < 0.01, ⁎⁎⁎p < 0.001 and ⁎⁎⁎⁎p < 0.0001 using by Student’s t-test to evaluate differences between every two groups. Treat: BHV-1 infected MDBK cells for 33 hours; Blank: MDBK cells not infected by the virus.

Figure 9 Correlation analysis of qRT-PCR vs RNA-seq. The correlation of log2 (fold change) between RT-qPCR and RNA-Seq data was analyzed by GraphPad software 6.0 (San Diego, CA), The abscissa represents the RT-qPCR data, the ordinate represents RNA-Seq data. The correlation coefficient (R2) was as high as 0.8274.

Discussion

Bovine herpesvirus type 1 (BoHV-1) is a member of the subfamily α-herpesvirinae and a virulent pathogen causing significant economic losses to the worldwide cattle industry,^28,29 but knowledge about the molecular mechanism of RNA in BoHV-1 infection such as the expression level of mRNAs in immune pathway, non-coding RNAs expression stats are still unclear. The use of RNA-Seq analysis is becoming one of the main techniques for transcript analyses. Recently, RNA-Seq has been mostly used in some virus infection studies,^30–32 such as bluetongue virus (BTV), peste des petits ruminants virus (PPRV) and bovine viral diarrhea virus (BVDV), but there were fewer reports about the changes of RNA level after the infection of BoHV-1. More in-depth studies of the interaction between BoHV-1 and the host will be of great significance for elucidating the pathogenesis and the progression of persistent BoHV-1 infection. In this study, MDBK cells were infected with BoHV-1 for 18, 24, 33 hours and sequenced by next-generation sequencing. The reason for choosing these time points was to look at the expression of lncRNAs in the early, middle, and late stages of viral infection. After the virus infects MDBK cells, it generally shows complete cytopathic effect within 48 hours, accompanied by mass cell death and shedding. In order to prevent RNA degradation due to excessive cell death, we chose 33 hours when cytopathic effects had already occurred but were not severe. The ceRNA hypothesis has been recently suggested as a novel layer of gene regulation. The mRNA-miRNA-lncRNA network was constructed by using raw-data, which had been sequenced and the lncRNAs, miRNAs and mRNAs in this network were analyzed to reveal the key genes related to BoHV-1 infection. By constructing mRNA-miRNA-lncRNA networks, this study found that lncRNA XR_003029469.1 had a higher degree in up-regulated lncRNAs network and bound four miRNAs including bta-miR-23b-5p, bta-miR-30b-3p, miR-12032-x, and miR-25-x. Sequencing results showed that MSTRG.8067.1 was a novel lncRNA whose function was unknown. Through mRNA-miRNA-lncRNA network construction, it was found that it formed a network with bta-miR-10174-5p, XM_005222015.4 and NM_001037628.2. This indicated that the function of MSTRG.8067.1 was related to bta-miR-10174-5p, XM_005222015.4 and NM_001037628.2. In down-regulated lncRNAs network, lncRNA XR_003036761.1 had a higher degree and bound three miRNAs, including miR-2370-y, miR-30-x and bta-miR-12059. In turn, miR-30-x bound to NM_001076487. This indicated that the function of lncRNA XR_003036761.1 might be exerted by the miR-30-x targeting mRNA NM_001076487.

In ceRNA networks, nodes with high connectivity often have important biological significance, and these genes are called hub genes.^33,34 In most studies,^35–37 hub genes were focused on some cancers and diseases such as prostate cancer, large B-Cell lymphoma and non-small cell lung cancer. However, hub genes predicted based on bioinformatics technology may be caused by errors in algorithms or association networks and the function of these genes is not clear and requires further research.

To assess functional enrichment, GO and KEGG pathway analyses were performed by using DEGs in the ceRNA networks. This study’s results indicated that the DEGs were mainly enriched in enzyme binding, substrate-specific channel activity, DNA binding and nuclease activity on the molecular function level. The studies³⁸ have shown that the BoHV-1 genome is GC-rich, with SP1 consensus binding sites and other GC-rich Moines, which may promote binding to specific KLF transcription factors and stimulate viral transcription early in latency reactivation.⁴ In the KEGG pathway, these DEGs were mostly enriched in the proteoglycans in cancer, cytokine-cytokine receptor interaction and viral carcinogenesis, mitogen-activated protein kinase and calcium signaling pathway. MAPK signaling pathways play an important role in pathogen–host interaction, some viruses such as Epstein-Barr virus (EBV) and herpes simplex virus (HSV) affect host cell survival or cell cycle or viral replication or alter immune responses by MAPK-ERK cascade.³⁹ Some studies^40–42 reported that virus infection with hepatitis C virus (HCV), human immunodeficiency virus (HIV) and pseudorabies virus (PRV) can activate the calcium signaling and MAPK signaling pathways. Especially for Pseudorabies virus (PRV), an alpha herpesvirus which can enter the mammalian nervous system in a similar way to BoHV-1, but whether the molecular mechanism of the BoHV-1 infection is the same as that of Pseudorabies virus (PRV) needs further study.

In summary, this study revealed that BoHV-1 infection can affect the expression of RNA in MDBK cells and the regulation of ceRNA network to carry out corresponding biological functions in the host, but further experimental studies are still necessary to prove the hub genes function in ceRNA network and the molecular mechanism in BoHV-1 infection.