Differential expression profile and in-silico functional analysis of long noncoding RNA and mRNA in duck embryo fibroblasts infected with duck plague virus
BMC Genomics volume 23, Article number: 509 (2022)
Duck plague virus (DPV), belonging to herpesviruses, is a linear double-stranded DNA virus. There are many reports about the outbreak of the duck plague in a variety of countries, which caused huge economic losses. Recently, increasing reports revealed that multiple long non-coding RNAs (lncRNAs) can possess great potential in the regulation of host antiviral immune response. Furthermore, it remains to be determined which specific molecular mechanisms are responsible for the DPV-host interaction in host immunity. Here, lncRNAs and mRNAs in DPV infected duck embryonic fibroblast (DEF) cells were identified by high-throughput RNA-sequencing (RNA-seq). And we predicted target genes of differentially expressed genes (DEGs) and formed a complex regulatory network depending on in-silico analysis and prediction.
RNA-seq analysis results showed that 2921 lncRNAs were found at 30 h post-infection (hpi). In our study, 218 DE lncRNAs and 2840 DE mRNAs were obtained in DEF after DPV infection. Among these DEGs and target genes, some have been authenticated as immune-related molecules, such as a Macrophage mannose receptor (MR), Anas platyrhynchos toll-like receptor 2 (TLR2), leukocyte differentiation antigen, interleukin family, and their related regulatory factors. Furthermore, according to the Kyoto Encyclopedia of Genes and Genomes (KEGG) and Gene Ontology (GO) enrichment analysis, we found that the target genes may have important effects on biological development, biosynthesis, signal transduction, cell biological regulation, and cell process. Also, we obtained, the potential targeting relationship existing in DEF cells between host lncRNAs and DPV-encoded miRNAs by software.
This study revealed not only expression changes, but also the possible biological regulatory relationship of lncRNAs and mRNAs in DPV infected DEF cells. Together, these data and analyses provide additional insight into the role of lncRNAs and mRNAs in the host's immune response to DPV infection.
Duck plague (DP) is caused by duck herpesvirus type I . The causative agent of the disease belongs to the species Anatid herpesvirus I, genus Mardivirus, subfamily Alphaherpesvirinae, family Herpesviridae . As waterfowl migrate, the disease can spread widely inside and between continents around the world, causing huge economic losses [3,4,5,6]. After DPV infection, the main injured organs in the disease seem to be the brain, liver, kidney, and lung, and the main pathological features are bleeding, inflammation, and necrosis of tissues, and organs [1, 7]. In China, the duck plague was first reported in 1959 and verified in many regions of China. And now the duck plague has not been wholly conquered in China, and there are still reports once in a while .
The definition of Long non-coding RNAs (lncRNAs) is non-protein-coding transcripts with more than 200 nucleotides, mainly classified into antisense lncRNAs, sense RNAs, bidirectional lncRNAs, large intergenic non-coding RNAs (lincRNAs), and intronic transcript [9, 10]. It has been found that lncRNAs have a variety of functions, including epigenetic regulation, transcriptional regulation, and post-transcriptional regulation, and unexpectedly, some transcripts annotated as lncRNAs can encode small proteins [10,11,12]. In addition, with the deepening of research, more researchers have started to pay more attention to that lncRNAs can also participate in the interaction between virus and host. The lnczc3h7a existed in mouse macrophage RAW264.7, can be used as a scaffold to stabilize retinoid acid-inducible gene I (RIG-I) -Tripartite Motif Containing 25 (TRIM25) complex and promote TRIM25-mediated ubiquitination of RIG-I, finally to participate in type I interferon (IFN) signal transduction and regulate RNA virus replication both in vivo and in vitro . Similarly, MALAT1 can bind to TAR DNA-binding protein 43 (TDP43) in resting macrophages and inhibit IFNs expression induced by interferon regulatory Factor 3 (IRF3) . In addition, there is a unique mechanism by which the virus absorbs host lncRNA to assist virus replication. After hepatitis C virus (HCV) infection, lncRNA EGOT was found to be significantly up-regulated. And Previous studies showed that EGOT had a highly significant negative correlation with innate immune-related genes that may block the antiviral response to promote virus replication .
Currently, there are increasing studies on lncRNAs related to avian diseases, such as the H5N1 influenza virus , duck enteritis salmonella , and duck Tembusu virus (DTMUV) . However, there is a lack of reports on host lncRNAs expression induced by DPV infection. To solve this problem, the transcriptome of DPV infected duck embryonic fibroblasts (DEF) was sequenced by RNA-seq. Here we identified the differentially expressed (DE) lncRNAs and mRNAs compared with the infected and uninfected DEF cells, and predicted the target genes of DE lncRNAs using bioinformatics analysis. Based on the enrichment analysis of Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG), we identified that differentially expressed genes (DEGs) in DEF cells were significantly associated with the biological processes and signal pathways.
Characterization of DPV infection in DEFs
At a multiplicity of infection (MOI) of 1.0, duck embryonic fibroblasts (DEF) cells were infected with DPV. In our study, CPE was observed in DPV-infected DEFs after 24 hpi, and DPV propagation in DEF cell culture began to slow down at 30 hpi (Fig. 1a, b), similar to previous studies [19, 20]. As a result, 30 hpi and an MOI of 1.0 were chosen as the condition for preparing cell samples for future high-throughput RNA-sequencing (RNA-seq).
After a series of data sorting and filtering processes, the number of high-quality reads per each sample is about 66,715,946–77,362,586 (Additional file 1: Table S1 and Additional file 2: Table S2). The filtered data were aligned to the reference genome. The percentage of the total aligned reads in the control group was more than 90%, while in the infection group it was more than 70% (Additional file 3: Table S3). We screened lncRNAs in the assembled transcript set [10, 21, 22]: (1) Screening the transcripts with length > = 200 bp and exon number > = 2; (2) Screening the transcripts with Class x/u/i (Anti-sense lncRNA/Intergenic lncRNA/Intronic lncRNA); (3) Screening for transcripts with coverage > 3 in at least one sample. That is, it appeared at least three times in one sample’s transcripts.
A total of 61,391 long transcripts were screened out in the first step, 5994 long transcripts were screened out in the second step, and finally, 5941 long transcripts were screened out for follow-up analysis (Fig. 2a). Further analysis of the coding potential of the candidate lncRNAs was performed to get high-credibility lncRNAs. A total of 2921 novel lncRNAs with high-reliability were found by the following three methods: plek, CNCI, and pfamscan (Fig. 2a, b).
Comparative analysis of lncRNAs and mRNAs
The structures of known and unknown lncRNAs and mRNAs were compared. On the one hand, the result showed the molecular differences between lncRNAs and mRNAs, on the other hand, it verified whether the predicted lncRNAs accord with the general characteristics. As shown in Fig. 2, the results showed that on the whole, the number of mRNA transcripts was widely distributed, ranging from 1 to 10, mainly in the range of 1–2, whereas the number of lncRNAs transcripts was only in the range of 1–6, and mainly a single transcript (Fig. 3a). In terms of length analysis, the major length of mRNAs was in the range of > = 5000 and 800–1600, while lncRNAs were chiefly distributed in the range of 200–300 and showed a stepwise downward trend in the range of 400–4800, but there was a certain quantitative trend rising in the range of > = 5000 (Fig. 3b). Most mRNAs contained a large number of exons, mainly in the range of > 10, while the number of exons in lncRNAs was small, only 2–3 (Fig. 3c). Based on the results of FPKM (fragments per kilobases per million fragments) density distribution, we observed that the FPKM distribution of mRNAs was gathered at 1–2, while the expression level of lncRNA is lower (Fig. 3d, e).
Deep sequencing results
We also performed PCA (Principal Component Analysis) to check the correlation between samples. And it was found that there were greater differences in the total lncRNAs of biological repeat samples in each treatment group, compared with mRNAs (Fig. 4a, b).
The calculation method of P value of DESeq2 is very strict. Therefore, we choose |log2 fold change (infected/mock in expression) |> 1 and p-value < 0.05 as the screening standard, hoping to get more DEGs for analysis. The differential expression analysis results indicated that compared with the control group (Fig. 4c), we identified 2840 DE mRNAs in DEF after DPV infection, of which 1386 were up-regulated and 1454 were down-regulated (Additional file 4: Table S4), and 218 DE lncRNAs, of which 169 were up-regulated and 49 were down-regulated (Additional file 5: Table S5 and Additional file 6: Table S6). Besides, we used volcano plots and heat maps to show DE lncRNAs and DE mRNAs (Fig. 5a-d).
Among them, the mRNA with the most significant difference in expression level are actin and myosin-related factors. Further, some DE mRNAs have been identified as immune-related molecules, which can be classified into the following three categories: pattern recognition receptors (PRRs) (such as Macrophage mannose receptor 1 and Anas platyrhynchos toll-like receptor 2), leukocyte differentiation antigens (like CD81, CD82, CD36, and CD2 associated protein) and immune cell-derived cytokines and their related regulatory factors (such as interleukin 15, interleukin 21 receptor, interleukin 6 family cytokine, interleukin 23 receptor, interleukin-9 receptor-like gene, interleukin 18 receptor accessory protein, interleukin 19, interleukin 6 signal transducer, cytokine-inducible SH2 containing protein, and Cytokine-like 1 gene). Interestingly, most of these genes were up-regulated in the infection group (Additional file 4: Table S4).
Potential interacting genes prediction
Then, we predicted potential interacting genes of lncRNA in cis- and trans-regulation, as shown in figures (Fig. 6 and Fig. 7). The previous study has shown that lncRNA can regulate the expression of adjacent genes by recruiting regulatory factors and competing for TF (Transcription Factor) and Pol ll (RNA polymerase II) required for transcription. Or the transcription of lncRNA can interfere with the cis-regulatory element of mRNA, contained in the lncRNA encoded region . Besides, lncRNAs can also regulate chromatin state and structure, organization of nuclear, and trans-acting factors to regulate distal gene expression . And it was found that in the predicted results, 3762 lncRNAs (of a total 3814) were predicted to have cis regulatory function, and1388 lncRNAs could act in trans-regulation. Interestingly, 1262 lncRNAs have a large number of possible associated genes, whether in cis-regulation or trans-regulation (Additional file 7: Table S7 and Additional file 8: Table S8).
Functional enrichment analysis of DE lncRNAs
The GO and KEGG analyses were conducted to predict the biological regulatory function of DE duck lncRNAs (Fig. 8a-d). Based on the GO enrichment analysis, we found cis-associated genes were significantly referred to the cell biological items (p < 0.05) such as hormone activity, double-stranded DNA binding, protein phosphorylation, microtubule anchoring, and transferase activity, while trans-acting potential target genes were significantly enriched in (p < 0.05) RNA process, filopodium and cell cycle, involving cell biological regulation, signal transduction, biosynthesis, response to external stimulus, and cell process (Additional file 9: Table S9 and Additional file 11: Table S11). Besides, there are some GO terms that may be associated with the process of virus infection, such as Wnt signaling pathway (GO:0,030,111, p < 0.05; GO:0,030,177, p < 0.05; GO:0,016,055) in trans-regulation, and IL-2 terms (GO: 0,004,911, p < 0.05; GO: 0,019,976, p < 0.05), JAK/STAT signaling pathway (GO: 0,007,259, GO: 0,097,696) and Wnt signaling pathway (GO:0,030,111, GO:0,030,178) in cis-regulation.
At the same time, KEGG analysis showed that multiple signal pathways associated with inflammation and host antiviral immune response were enriched in (Additional file 10: Table S10 and Additional file 12: Table S12), for instance, Toll-like receptor signaling pathway, C-type lectin receptor signaling pathway, RIG-I-like receptor signaling pathway, and NOD-like receptor signaling pathway were found in both trans- and cis-regulation. Also, the mitogen-activated protein kinase (MAPK) signaling pathway and Cytokine-cytokine receptor interaction were enriched at cis-regulation.
Functional enrichment results of DE mRNAs
In addition, we performed functional analyses of DE mRNAs (Fig. 9a, b). As shown in the results, DE mRNA was significantly related to (p < 0.05) biological development and biosynthesis (Additional file 13: Table S13). There are also a large number of enriched GO terms related to immunity, such as regulation of defense response (GO:0,031,347, p < 0.05), response to virus (GO:0,009,615, p < 0.05), defense response to virus (GO:0,051,607, p < 0.05), and negative regulation of I-kappaB kinase/NF-kappaB signaling (GO:0,043,124, p < 0.05).
The results of KEGG pathway analyses revealed that the enrichment results were significantly related to (p < 0.05) the signal transduction, metabolism, and biomembrane (Additional file 14: Table S14). And the top enriched KEGG pathway was MAPK signaling pathway (p < 0.05). Similar to the enrichment results of lncRNA predicted target genes, Cytokine-cytokine receptor interaction, Wnt signaling pathway, NOD-like receptor signaling pathway, Toll-like receptor signaling pathway, and RIG-I-like receptor signaling pathway were all present in the KEGG pathway analyses.
Establishment of predicted competitive endogenous RNA (ceRNA) network
LncRNA-miRNA plays a significant role in gene expression regulation. LncRNA with microRNA response elements (MREs) can be used as a competing endogenous RNA (ceRNA) to inactivate miRNA by binding to miRNA . Because most lncRNAs are structurally similar to mRNAs, miRNA may also negatively regulate the expression of lncRNA through a similar mechanism of action on mRNA, and then plays a series of biological roles . At the same time, miRNA can also promote the expression of specific lncRNA . In this study, we obtained the prediction results of the interaction between DE lncRNAs and miRNAs encoded by DPV virus . And we found that each miRNA was predicted to interact with a large number of lncRNAs, and many lncRNAs could also be targeted by multiple miRNAs at the same time (Additional file 15: Table S15 and Additional file 16: Figure S1), for instance, 24 viral miRNAs were all predicted to target ENSAPLT00020019865.
Validation of DE lncRNAs and mRNAs by qRT-PCR
To validate the DEGs in RNA-seq results, we randomly selected 12 DEGs (6 lncRNAs and 6 mRNAs) for qRT-PCR analysis. And it was found that among the 6 lncRNAs, there were 5 up-regulated lncRNAs and 1 down-regulated lncRNA (MSTRG.24727.2) by qRT-PCR analysis, differential expression patterns of which were consistent with RNA-seq results (Additional file 17: Figure S2). Similarly, the expression trend of the chosen mRNAs was in accordance with the results of RNA-seq analysis (Additional file 18: Figure S3).
With an in-depth understanding of lncRNA, it has been found that virus infection can affect the expression of host lncRNAs, and lncRNAs can also play a unique role in the process of virus infection. For example, RIG-I-dependent IAV upregulated noncoding RNA (RDUR), induced by Influenza A virus (IAV) infection, can positively regulate the signal transduction of type I interferon (IFN-1) by affecting IRF3 and further upregulates the expression of several vital antiviral molecules. It is concluded that RDUR is a critical regulator of innate immunity against viral infection . And lncRNA-GM can bind to glutathione S-transferase M1 (GSTM1), block the interaction between GSTM1 and TANK binding kinase 1 (TBK1), and reduce GSTM1-mediated glutathione. Finally, the deficiency of lncRNA-GM in mice increases sensitivity to viral infection and interferes with the production of IFN-I . There are also lncRNA-ACOD1, which are effectively induced by several viruses. This lncRNA can promote viral replication by activating the glutamic acid transaminase 2 (GOT2), a metabolic enzyme. Interestingly, the enhancement of GOT2 activity by lncRNA-ACOD1 is necessary for virus replication . DPV is a pathogen causing a hot, acute, and septic infectious disease in ducks, geese, and swans, with a great number of economic losses . Moreover, the regulatory mechanism between DPV and host remains to be discussed. Here we analyzed the transcriptome of DPV infected DEF and obtained 218 DE lncRNAs and 2840 DE mRNAs at 30 hpi.
In our study, we found that some DE mRNAs were closely related to immune regulation, such as toll-like receptor 2 (TLR2), mannose receptor (MR) family, T cell surface antigen, interleukin family, and their related regulatory factors. There are many reports showed that the above factors are involved in the regulation of the host antiviral immune response. For example, the coordinated activation of TLR2 and TLR9 in trigeminal ganglion exists in the defence response of mice to the herpes simplex virus-1 (HSV-1), while the expression of TLRs in HSV-1 infected mice is significantly increased . And it has been reported that innate immune recognition of vaccinia virus (VV) is mediated by TLR2 . However, it is also reported that the above immune factors can participate in the interaction between the virus and the host. For example, macrophage mannose receptor (MMR) can participate in the infection of mouse macrophages by influenza virus, which is necessary for endocytic uptake of virus . Similarly, CD81 can promote the entry of the HCV into host cells , and CD2 related protein (CD2AP) can promote the assembly of HCV . Meanwhile, the differential expression level of actin-related factors changed significantly. Actin is abundant in eukaryotic cells and previous studies have shown that actin could intimately interact with lentiviruses . Also, cofilin, one of the actin-binding proteins, was found to be involved in the process of HSV infection . The findings suggested that DE mRNAs may similarly regulate the pathogenicity of DPV.
In addition, the functional analysis indicated that the target genes and DEGs significantly related to (p < 0.05) biological development, biosynthesis, metabolism, and signal transduction. At the same time, some targets were enriched in IL-2 terms, Wnt signaling pathway, JAK/STAT pathway terms, and other important GO terms. Interleukin-2 (IL-2) is a growth factor of various cells, which is closely related to the growth and differentiation of T cells, the development of functional B cells, and the activation of NK cells . Studies have proved that duck IL-2 is an important immunomodulatory molecule, which has been fully verified in chickens and mammals . Similarly, JAK/STAT pathway is an important cytokine signal transduction pathway, which is activated by a variety of cytokines, growth factors, and receptors. And it involves many biological processes such as immune regulation, inflammatory response, cell proliferation, differentiation, apoptosis, angiogenesis, oxidative stress, and tumor formation . Also, it is related to the process of a variety of viruses infection, such as porcine deltacoronavirus (PDCoV) , hepatitis B virus (HBV) . And the Wnt pathway is an important signal pathway, which participates in the regulation of a variety of cell processes, which is closely related to the development of the body, the maintenance of homeostasis, and the occurrence and development of a variety of diseases, including viral diseases [42,43,44]. Therefore, the above results indicated that there may be IL-2, Wnt pathway or JAK/STAT pathway-dependent manners in DEF after DPV infection to deal with virus infection. In addition, after the KEGG pathway enrichment analysis, we also found target genes were significantly enriched in MAPK signaling pathway (p < 0.05), RIG-I-like receptor signaling pathway, Toll-like receptor signaling pathway, and other immune-related pathways. It has been reported that MAPK pathway activation is important for the replication of many viruses, such as foot-and-mouth disease virus (FMDV) , bovine herpesvirus type 1 (BoHV-1) [46, 47], and KS associated herpesvirus (KSHV) . Perhaps there is a similar mechanism in DPV infected DEF cells. Further, both RIG-I and TLR can participate in NF-κB signal transduction after detecting viruses to regulate host antiviral immunity [49, 50]. This is also consistent with the enrichment of NF-κB signaling pathway-related GO terms (GO:0,043,124, p < 0.05) in the functional analysis of DEGs.
In addition, we also analyzed the potential viral miRNAs that may interact with DE lncRNAs based on the DPV-encoded miRNA sequence spectrum of Wu et al. . Previous studies have also shown that viral miRNAs can also play an important role in virus invasion and replication process [51,52,53], and that lncRNA could post-transcriptionally regulate target genes by microRNA response elements (MREs) . Here, we would like to further understand whether DE lncRNAs may be involved in a similar process during DPV infection. However, these analysis results need further experimental exploration to prove.
In our study, we screened out numerous DE lncRNAs and mRNAs by RNA-seq, formed a potential target interaction network in silico analysis, and predicted potential biological functions. LncRNA continues to be hot in the field of scientific research in recent years because of its various gene regulation mechanisms and a wide range of biological functions. We obtained a total of 218 DE lncRNAs and 2840 DE mRNAs, which were beneficial to understand the cellular mechanism of DPV infection with DEF.
Cell and virus
Primary DEF cells were obtained from 9-day-old duck embryos after trypsin digestion for subsequent experiments and cultured in Dulbecco's modified eagle medium (DMEM) supplemented with 10% fetal bovine serum (Gibco, USA) at 37 °C with 5% CO2. The DPV CHv (GenBank accession No. JQ647509) used in this study was kept in our laboratory. At a multiplicity of infection (MOI) of 1.0, duck embryonic fibroblasts (DEF) cells were infected with DPV or mock-infected with DMEM as control.
RNA isolation and RNA-Seq
All experiments were carried out in triplicate. At 30 h post-infection (hpi) [19, 20], total cellular RNAs were extracted using the RNAiso plus reagent (TaKaRa, Japan). The purity, integrity, and concentration of total RNA samples were evaluated and tested by NanoDrop NC2000 Spectrophotometer (Thermo Fisher Scientific, USA) and Agilent RNA 6000 nano kit (Agilent Technologies, USA). Rin value (RNA integrity number) was automatically generated by Agilent 2100 Bioanalyzer (Agilent Technologies, USA) which represents the quantitative value of RNA integrity. The quality test results showed that the RIN values of the samples were 10, which met the requirements of database construction and sequencing. All the infected and uninfected samples were sequenced by Shanghai Personal Biotechnology (Shanghai, China).
Reads filtering and mapping
We calculate the data of each sample (Raw Data) separately, including sample names, Q30, percentage of ambiguous reads, Q20 (%), and Q30 (%). Remove low-quality reads using Cutadapt software (http://cutadapt.readthedocs.io/en/stable/). Then the filtered reads were aligned to the reference genome (Anas_platyrhynchos.ASM874695v1.dna.toplevel.fa, http://asia.ensembl.org/Anas_platyrhynchos/Info/Annotation) using the HISAT2 software (http://ccb.jhu.edu/software/hisat2/index.shtml) . Based on the preliminary gene mapping results, we assembled transcripts by StringTie software (http://ccb.jhu.edu/software/stringtie/) . We used the following three methods for coding potential analysis: plek (https://sourceforge.net/projects/plek/files/) , CNCI (https:// github.com/www-bioinfo-org/CNCI) , pfamscan (https://www.ebi.ac.uk/seqdb/confluence/display/JDSAT/PfamScan+Help+and+Documentation) , and considered that the new transcripts without coding potential determined by the three software were high-reliability lncRNAs.
Identification of differentially expressed RNAs
The correlation of gene expression level between biological repeat samples is an important index to test the reliability of experiments and the rationality of sample selection. We used DESeq Software package(http://www.bioconductor.org/packages/release/bioc/html/DESeq.html) to carry out PCA (principal component analysis) for each sample according to the expression amount [59, 60]. Then, we used HTSeq statistics (https://htseq.readthedocs.io/en/master/) to obtain the reads count on each gene as the original expression level . In order to facilitate comparison of gene expression of different genes in different samples, we used FPKM (fragments per kilobases per million fragments) as the standard of measurement by the Cufflinks program (http://cole-trapnell-lab.github.io/cufflinks/) . We also used Stringtie software to obtain the reads count of lncRNAs and then homogenized the expression level with FPKM. The differential expression analyses of RNAs were performed by DESeq (http://www.bioconductor.org/packages/release/bioc/html/DESeq.html) [59, 60]. The conditions for screening DEGs were as follows: |log2 fold change (infected/mock in expression) |> 1 and p-value < 0.05. We used the phatmap software package (https://CRAN.R-project.org/package=pheatmap) to conduct a two-way cluster analysis on the Union and samples of different genes in all comparison groups. According to the expression level of the same genes in different samples and the expression patterns of different genes in the same sample, the Euclidean method was used to calculate the distance, and the hierarchical clustering longest distance method (complete linkage) was used for clustering.
Target genes prediction of lncRNAs
Increasing researches showed that the cis-regulatory function of lncRNAs is associated with the protein-coding genes near their genomic locations [10, 19, 63]. Therefore, to explore the regulatory function of lncRNAs, the protein-coding genes within 100 kb upstream and downstream of the lncRNA genes were searched and thought to be the cis-regulatory target genes of their corresponding lncRNAs.
In contrast, the basic principle of trans-regulatory target genes prediction holds that the function of lncRNA is dependent on its co-expressed protein-coding genes [10, 19, 63]. The conditions of correlation coefficient | correlation |> 0.95 and p-value < 0.05 were used to screen the trans-regulatory relationship between lncRNAs and mRNAs.
In-silico gene functional prediction analysis
We used the TopGO (http://www.bioconductor.org/packages/release/bioc/html/RamiGO.html) for GO enrichment analysis of DEGs and target genes . Based on the DEGs and target genes annotated by the GO database (http://geneontology.org/), then we calculated the p-value of DEGs and target genes by hypergeometric distribution method (p-value < 0.05 is regarded as the standard of significant enrichment), and took the whole genome as the background to screen the GO terms with significant enrichment of DEGs and target genes, to determine the main biological functions performed by DEGs and target genes . Furthermore, we predicted the associated pathways of DEGs and target genes by the KEGG database (http://www.kegg.jp/) [66, 67]. Similarly, the criterion for significant enrichment is p-value < 0.05.
Competitive endogenous RNAs (ceRNAs) regulatory network construction
We used miRanda software (https://www.cs.kent.ac.uk/people/staff/dat/miranda/) and the DPV-encoded miRNA sequence spectrum to predict target miRNAs for lncRNAs respectively, and finally obtained the predicted miRNA-lncRNA pairs to have a targeting relationship [24, 68]. The final lncRNA-miRNA regulatory networks were visualized with Cytoscape software (https://cytoscape.org/) .
QRT -PCR Analysis
We randomly selected several DEGs and determined the expression changes by Quantitative Real-time PCR (qRT-PCR) to verify the validity of sequencing data. The extracted total RNA was used to reverse transcribed to cDNA by PrimeScript™ RT reagent Kit with gDNA Eraser (Perfect Real Time) (TaKaRa, Japan). The expression of DEGs was detected by the CFX Connect real-Time PCR Detection System (Bio-Rad, USA) and TB Green® Premix Ex Taq™ II(Tli RNaseH Plus) (TaKaRa, Japan). The β-actin gene was used to standardize the expression of DEGs, and calculate relative expression levels by the 2-△△Ct method . The primers used in this study were synthesized by Huada Company (Guangdong, China) and shown in Table 1.
Availability of data and materials
All data generated or analysed during this study are included in this article and its additional files. All RNA-seq data analysed during this study has been deposited in in SRA database (https://www.ncbi.nlm.nih.gov/sra) under accession number: SRR16573633, SRR16573634, SRR16573635, SRR16573636, SRR16573637, SRR16573638.
Duck plague virus
Long non-coding RNA
Duck embryonic fibroblast
Differentially expressed genes
Toll-like receptor 2
Kyoto Encyclopedia of Genes and Genomes
Duck viral enteritis
Large intergenic non-coding RNA
Retinoid acid-inducible gene I
Tripartite Motif Containing 25
TAR DNA-binding protein 43
Interferon regulatory Factor 3
Hepatitis C virus
Duck Tembusu virus
Dulbecco's modified eagle medium
Fetal bovine serum
Phosphate buffered saline
Multiplicity of infection
Quantitative Real-time PCR
Pattern recognition receptors
Mitogen-activated protein kinase
Competing endogenous RNA
MicroRNA response elements
Short internal repeat
Short terminal repeat
Herpes simplex virus-1
CD2 related protein
Porcine deltacorona virus
Hepatitis B virus
Foot-and-mouth disease virus
Bovine herpesvirus type 1
KS associated herpesvirus
RIG-I-dependent IAV upregulated noncoding RNA
Influenza A virus
Glutathione S-transferase M1
TANK binding kinase 1
Glutamic acid transaminase 2
Principal component analysis
- Pol ll:
RNA polymerase II
Dhama K, Kumar N, Saminathan M, Tiwari R, Karthik K, Kumar MA, et al. Duck virus enteritis (duck plague) - a comprehensive update. Vet Q. 2017;37(1):57–80.
Matthews REF. Virus Taxonomy: Classification and Nomenclature of Viruses: Ninth Report of the International Committee on Taxonomy of Viruses. Intervirology. 1979;12:129–296.
Kaleta EF, Kuczka A, Kühnhold A, Bunzenthal C, Bönner BM, Hanka K, et al. Outbreak of duck plague (duck herpesvirus enteritis) in numerous species of captive ducks and geese in temporal conjunction with enforced biosecurity (in-house keeping) due to the threat of avian influenza A virus of the subtype Asia H5N1. Dtsch Tierarztl Wochenschr. 2007;114(1):3–11.
Campagnolo ER, Banerjee M, Panigrahy B, Jones RL. An outbreak of duck viral enteritis (duck plague) in domestic Muscovy ducks (Cairina moschata domesticus) in Illinois. Avian Dis. 2001;45(2):522–8.
Davison S, Converse KA, Hamir AN, Eckroade RJ. Duck viral enteritis in domestic muscovy ducks in Pennsylvania. Avian Dis. 1993;37(4):1142–6.
Burgess EC, Ossa J, Yuill TM. Duck plague: a carrier state in waterfowl. Avian Dis. 1979;23(4):940–9.
Metwally SA. Duck Virus Enteritis (Duck Plague). Diseases of Poultry. 2013:431–40.
Wang G, Qu Y, Wang F, Hu D, Liu L, Li N, et al. The comprehensive diagnosis and prevention of duck plague in northwest Shandong province of China. Poult Sci. 2013;92(11):2892–8.
Meng XY, Luo Y, Anwar MN, Sun Y, Gao Y, Zhang H, et al. Long Non-Coding RNAs: Emerging and Versatile Regulators in Host-Virus Interactions. Front Immunol. 2017;8:1663.
Kopp F, Mendell JT. Functional Classification and Experimental Dissection of Long Noncoding RNAs. Cell. 2018;172(3):393–407.
Fortes P, Morris KV. Long noncoding RNAs in viral infections. Virus Res. 2016;212:1–11.
Ali T, Grote P. Beyond the RNA-dependent function of LncRNA genes. Elife. 2020;9: e60583.
Lin H, Jiang M, Liu L, Yang Z, Ma Z, Liu S, et al. The long noncoding RNA Lnczc3h7a promotes a TRIM25-mediated RIG-I antiviral innate immune response. Nat Immunol. 2019;20(7):812–23.
Liu W, Wang Z, Liu L, Yang Z, Liu S, Ma Z, et al. LncRNA Malat1 inhibition of TDP43 cleavage suppresses IRF3-initiated antiviral innate immunity. Proc Natl Acad Sci U S A. 2020;117(38):23695–706.
Carnero E, Barriocanal M, Prior C, Pablo Unfried J, Segura V, Guruceaga E, et al. Long noncoding RNA EGOT negatively affects the antiviral response and favors HCV replication. EMBO Rep. 2016;17(7):1013–28.
Lu C, Xing Y, Cai H, Shi Y, Liu J, Huang Y. Identification and analysis of long non-coding RNAs in response to H5N1 influenza viruses in duck (Anas platyrhynchos). BMC Genomics. 2019;20(1):36.
Zhang Y, Dong X, Hou L, Cao Z, Zhu G, Vongsangnak W, et al. Identification of Differentially Expressed Non-coding RNA Networks With Potential Immunoregulatory Roles During Salmonella Enteritidis Infection in Ducks. Front Vet Sci. 2021;8: 692501.
Lin Y, Yang J, He D, Li X, Li J, Tang Y, et al. Differently Expression Analysis and Function Prediction of Long Non-coding RNAs in Duck Embryo Fibroblast Cells Infected by Duck Tembusu Virus. Front Immunol. 2020;11:1729.
Kocan RM. Duck plague virus replication in muscovy duck fibroblast cells. Avian Dis. 1976;20(3):574–80.
Breese SS Jr, Dardiri AH. Electron microscopic characterization of duck plague virus. Virology. 1968;34(1):160–9.
Qian X, Zhao J, Yeung PY, Zhang QC, Kwok CK. Revealing lncRNA Structures and Interactions by Sequencing-Based Approaches. Trends Biochem Sci. 2019;44(1):33–52.
Cabili MN, Trapnell C, Goff L, Koziol M, Tazon-Vega B, Regev A, et al. Integrative annotation of human large intergenic noncoding RNAs reveals global properties and specific subclasses. Genes Dev. 2011;25(18):1915–27.
Thomson DW, Dinger ME. Endogenous microRNA sponges: evidence and controversy. Nat Rev Genet. 2016;17(5):272–83.
Chiyomaru T, Fukuhara S, Saini S, Majid S, Deng G, Shahryari V, et al. Long non-coding RNA HOTAIR is targeted and regulated by miR-141 in human cancer cells. J Biol Chem. 2014;289(18):12550–65.
Braconi C, Kogure T, Valeri N, Huang N, Nuovo G, Costinean S, et al. microRNA-29 can regulate expression of the long non-coding RNA gene MEG3 in hepatocellular cancer. Oncogene. 2011;30(47):4750–6.
Wu X, Jia R, Zhou J, Wang M, Chen S, Liu M, et al. Virulent duck enteritis virus infected DEF cells generate a unique pattern of viral microRNAs and a novel set of host microRNAs. BMC Vet Res. 2018;14(1):144.
Chen Y, Hu J, Liu S, Chen B, Xiao M, Li Y, et al. RDUR, a lncRNA, Promotes Innate Antiviral Responses and Provides Feedback Control of NF-κB Activation. Front Immunol. 2021;12: 672165.
Wang Y, Wang P, Zhang Y, Xu J, Li Z, Li Z, et al. Decreased Expression of the Host Long-Noncoding RNA-GM Facilitates Viral Escape by Inhibiting the Kinase activity TBK1 via S-glutathionylation. Immunity. 2020;53(6):1168-81.e7.
Wang P, Xu J, Wang Y, Cao X. An interferon-independent lncRNA promotes viral replication by modulating cellular metabolism. Science. 2017;358:1051–5.
Zolini GP, Lima GK, Lucinda N, Silva MA, Dias MF, Pessoa NL, et al. Defense against HSV-1 in a murine model is mediated by iNOS and orchestrated by the activation of TLR2 and TLR9 in trigeminal ganglia. J Neuroinflammation. 2014;11:20.
Zhu J, Martinez J, Huang X, Yang Y. Innate immunity against vaccinia virus is mediated by TLR2 and requires TLR-independent production of IFN-beta. Blood. 2007;109(2):619–25.
Upham JP, Pickett D, Irimura T, Anders EM, Reading PC. Macrophage receptors for influenza A virus: role of the macrophage galactose-type lectin and mannose receptor in viral entry. J Virol. 2010;84(8):3730–7.
Zona L, Lupberger J, Sidahmed-Adrar N, Thumann C, Harris HJ, Barnes A, et al. HRas signal transduction promotes hepatitis C virus cell entry by triggering assembly of the host tetraspanin receptor complex. Cell Host Microbe. 2013;13(3):302–13.
Zhang H, Zhang C, Tang H, Gao S, Sun F, Yang Y, et al. CD2-Associated Protein Contributes to Hepatitis C, Virus Propagation and Steatosis by Disrupting Insulin Signaling. Hepatology. 2018;68(5):1710–25.
Cantin R, Méthot S, Tremblay MJ. Plunder and stowaways: incorporation of cellular proteins by enveloped viruses. J Virol. 2005;79(11):6577–87.
Zheng K, Xiang Y, Wang X, Wang Q, Zhong M, Wang S, et al. Epidermal growth factor receptor-PI3K signaling controls cofilin activity to facilitate herpes simplex virus 1 entry into neuronal cells. mBio. 2014;5(1):e00958-13.
Smith KA. Interleukin-2: inception, impact, and implications. Science. 1988;240(4856):1169–76.
Zhou JY, Wang JY, Chen JG, Wu JX, Gong H, Teng QY, et al. Cloning, in vitro expression and bioactivity of duck interleukin-2. Mol Immunol. 2005;42(5):589–98.
Gurzov EN, Stanley WJ, Pappas EG, Thomas HE, Gough DJ. The JAK/STAT pathway in obesity and diabetes. FEBS J. 2016;283(16):3002–15.
Zhu X, Wang D, Zhou J, Pan T, Chen J, Yang Y, et al. Porcine Deltacoronavirus nsp5 Antagonizes Type I Interferon Signaling by Cleaving STAT2. J Virol. 2017;91(10):e00003-17.
Zhu X, Xie C, Li YM, Huang ZL, Zhao QY, Hu ZX, et al. TMEM2 inhibits hepatitis B virus infection in HepG2 and HepG2.2.15 cells by activating the JAK-STAT signaling pathway. Cell Death Dis. 2016;7(6):e2239.
Croce JC, McClay DR. Evolution of the Wnt pathways. Methods in molecular biology (Clifton, NJ). 2008;469:3–18.
Zhu X, Wen L, Wang W, Xiao Q, Li B, He K. PCV2 inhibits the Wnt signalling pathway in vivo and in vitro. Vet Microbiol. 2020;247: 108787.
Zwezdaryk KJ, Combs JA, Morris CA, Sullivan DE. Regulation of Wnt/β-catenin signaling by herpesviruses. World journal of virology. 2016;5(4):144–54.
Zhu Z, Li W, Zhang X, Wang C, Gao L, Yang F, et al. Foot-and-Mouth Disease Virus Capsid Protein VP1 Interacts with Host Ribosomal Protein SA To Maintain Activation of the MAPK Signal Pathway and Promote Virus Replication. J Virol. 2020;94(3):e01350-e1419.
Zhu L, Yuan C, Ding X, Jones C, Zhu G. The role of phospholipase C signaling in bovine herpesvirus 1 infection. Vet Res. 2017;48(1):45.
Zhu L, Ding X, Zhu X, Meng S, Wang J, Zhou H, et al. Biphasic activation of PI3K/Akt and MAPK/Erk1/2 signaling pathways in bovine herpesvirus type 1 infection of MDBK cells. Vet Res. 2011;42(1):57.
Zhu X, Zhou F, Qin D, Zeng Y, Lv Z, Yao S, et al. Human immunodeficiency virus type 1 induces lytic cycle replication of Kaposi’s-sarcoma-associated herpesvirus: role of Ras/c-Raf/MEK1/2, PI3K/AKT, and NF-κB signaling pathways. J Mol Biol. 2011;410(5):1035–51.
Zhang HX, Liu ZX, Sun YP, Zhu J, Lu SY, Liu XS, et al. Rig-I regulates NF-κB activity through binding to Nf-κb1 3’-UTR mRNA. Proc Natl Acad Sci U S A. 2013;110(16):6459–64.
Kawai T, Akira S. Signaling to NF-kappaB by Toll-like receptors. Trends Mol Med. 2007;13(11):460–9.
Židovec Lepej S, Matulić M, Gršković P, Pavlica M, Radmanić L, Korać P. miRNAs: EBV Mechanism for Escaping Host’s Immune Response and Supporting Tumorigenesis. Pathogens (Basel, Switzerland). 2020;9(5):353.
Zhu ZJ, Teng M, Li HZ, Zheng LP, Liu JL, Chai SJ, et al. Marek’s Disease Virus (Gallid alphaherpesvirus 2)-Encoded miR-M2-5p Simultaneously Promotes Cell Proliferation and Suppresses Apoptosis Through RBM24 and MYOD1-Mediated Signaling Pathways. Front Microbiol. 2020;11: 596422.
Yang X, Li H, Sun H, Fan H, Hu Y, Liu M, et al. Hepatitis B Virus-Encoded MicroRNA Controls Viral Replication. J Virol. 2017;91(10):e01919-e2016.
Zhang Y, Park C, Bennett C, Thornton M, Kim D. Rapid and accurate alignment of nucleotide conversion sequencing reads with HISAT-3N. Genome Res. 2021;31(7):1290–5.
Pertea M, Pertea GM, Antonescu CM, Chang TC, Mendell JT, Salzberg SL. StringTie enables improved reconstruction of a transcriptome from RNA-seq reads. Nat Biotechnol. 2015;33(3):290–5.
Mistry J, Bateman A, Finn RD. Predicting active site residue annotations in the Pfam database. BMC Bioinformatics. 2007;8:298.
Sun L, Luo H, Bu D, Zhao G, Yu K, Zhang C, et al. Utilizing sequence intrinsic composition to classify protein-coding and long non-coding transcripts. Nucleic Acids Res. 2013;41(17): e166.
Li A, Zhang J, Zhou Z. PLEK: a tool for predicting long non-coding RNAs and messenger RNAs based on an improved k-mer scheme. BMC Bioinformatics. 2014;15(1):311.
Wang L, Feng Z, Wang X, Wang X, Zhang X. DEGseq: an R package for identifying differentially expressed genes from RNA-seq data. Bioinformatics. 2010;26(1):136–8.
Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15(12):550.
Anders S, Pyl PT, Huber W. HTSeq–a Python framework to work with high-throughput sequencing data. Bioinformatics. 2015;31(2):166–9.
Trapnell C, Hendrickson DG, Sauvageau M, Goff L, Rinn JL, Pachter L. Differential analysis of gene regulation at transcript resolution with RNA-seq. Nat Biotechnol. 2013;31(1):46–53.
Chen LL. Linking Long Noncoding RNA Localization and Function. Trends Biochem Sci. 2016;41(9):761–72.
Ashburner M, Ball CA, Blake JA, Botstein D, Butler H, Cherry JM, et al. Gene Ontology: tool for the unification of biology. Nat Genet. 2000;25(1):25–9.
Götz S, García-Gómez JM, Terol J, Williams TD, Nagaraj SH, Nueda MJ, et al. High-throughput functional annotation and data mining with the Blast2GO suite. Nucleic Acids Res. 2008;36(10):3420–35.
Kanehisa M, Araki M, Goto S, Hattori M, Hirakawa M, Itoh M. KEGG for linking genomes to life and the environment. Nucleic Acids Res. 2008;36:D480–4.
Moriya Y, Itoh M, Okuda S, Yoshizawa AC, Kanehisa M. KAAS: an automatic genome annotation and pathway reconstruction server. Nucleic Acids Res. 2007;35(Web Server issue):W182–5.
Enright AJ, John B, Gaul U, Tuschl T, Sander C, Marks DS. MicroRNA targets in Drosophila. Genome Biol. 2003;5(1):R1.
Shannon P, Markiel A, Ozier O, Baliga NS, Wang JT, Ramage D, et al. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. 2003;13(11):2498–504.
Livak KJ, Schmittgen TD. Analysis of relative gene expression data using real-time quantitative PCR and the 2(-Delta Delta C(T)) Method. Methods (San Diego, Calif). 2001;25(4):402–8.
This work was supported by the National Natural Science Foundation of China (31872475), Sichuan Veterinary Medicine and Drug Innovation Group of China Agricultural Research System (CARS-SVDIP), China Agricultural Research System (CARS-42–17), and Sichuan Province Program (2022NSFSC0078).
Ethics approval and consent to participate
All animal experimental procedures were approved by the Institutional Animal Care and Use Committee of Sichuan Agriculture University in Sichuan, China (Protocol Permit Number: SYXK(川)2019–187). Experiments were conducted in accordance with guidelines and regulations of the Institutional Animal Care and Use Committee of Sichuan Agriculture University. This study was carried out in compliance with the ARRIVE guidelines.
Consent for publication
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Table S1. Summary of sequencing data.
Table S2. Summary of Filtered data.
Table S3. Summary of RNASeq Mapping.
Table S4. Summary of differentially expressed mRNA in DPV infected DEF cells.
: Table S5. Summary of differentially expressed lncRNA in DPV infected DEF cells.
Table S6. genomic position of lncRNAs.
Table S7. Prediction of cis-regulatory targets of lncRNA.
: Table S8. Prediction of trans-regulatory targets of lncRNA.
Table S9. GO enrichment of cis-regulatory targets of DE lncRNA.
KEGG enrichment of cis-regulatory targets of DE lncRNA.
Table S11. GO enrichment of trans-regulatory targets of DE lncRNA.
Table S12. KEGG enrichment of trans-regulatory targets of DE lncRNA.
Table S13. GO enrichment of DE mRNA.
Table S14. KEGG enrichment of targets of DE mRNA.
Table S15. Prediction of potential targeting viral miRNAs.
Figure S1. The predicted targeting regulatory networks of lncRNA-viral miRNA in DPV infected DEFs.
: Figure S2. Validation of the differential expression of 6 DE lncRNAs by qRT-PCR.
Figure S3. Validation of the differential expression of 6 DE mRNAs by qRT-PCR.
About this article
Cite this article
Wu, Z., Zeng, Y., Cheng, A. et al. Differential expression profile and in-silico functional analysis of long noncoding RNA and mRNA in duck embryo fibroblasts infected with duck plague virus. BMC Genomics 23, 509 (2022). https://doi.org/10.1186/s12864-022-08739-7
- Duck plague virus
- Functional analysis