Skip to main content
Download PDF

Genome-wide analysis of differentially expressed mRNAs, lncRNAs, and circRNAs in chicken bursae of Fabricius during infection with very virulent infectious bursal disease virus

BMC Genomics volume 21, Article number: 724 (2020) Cite this article

Abstract

Background

Infectious bursal disease virus (IBDV) causes acute, highly contagious, immunosuppressive, and lethal infectious disease in young chickens and mainly infects the bursa of Fabricius (BF). To investigate interactions between IBDV and its host, RNA sequencing was applied to analyze the responses of the differentially expressed transcriptional profiles of BF infected by very virulent IBDV (vvIBDV).

Results

In total, 317 upregulated and 94 downregulated mRNAs were found to be significantly differentially expressed in infected chickens, compared to controls. Long non-coding RNA (lncRNA) and circular RNA (circRNA) alterations were identified in IBDV-infected chickens, and significantly different expression was observed in 272 lncRNAs and 143 circRNAs. Gene Ontology and Kyoto Encyclopedia of Genes and Genomes pathway enrichment analyses were performed to assess the functions of significantly dysregulated genes, which showed that the JAK-STAT signaling pathway, the NOD-like receptor signaling pathway, and apoptosis may be activated by IBDV infection. We predicted interactions between differentially expressed genes and produced lncRNA-mRNA and circRNA-miRNA-mRNA regulator network.

Conclusions

The present study identified the expression profiles of mRNAs, lncRNAs, and circRNAs during vvIBDV infection and provides new insights into the pathogenesis of IBDV and antiviral immunity of the host.

Background

Infectious bursal disease virus (IBDV), a non-enveloped double-stranded RNA virus, is a member of the family Birnaviridae; it can cause acute, highly contagious, and immunosuppressive disease in chickens aged 3–6 weeks, leading to high mortality and considerable economic losses [1, 2]. Serotype-I and -II strains of IBDV are recognized, and serotype I, which includes attenuated, classically virulent, and very virulent (vv) variants, can cause different degrees of pathogenicity and mortality in chickens [3]. IBDV predominantly targets the precursors of B lymphocytes, particularly those in the bursa of Fabricius (BF), an important immune organ, the infection of which may lead to B lymphocyte depletion and BF deterioration [4, 5]. The severe immunosuppression caused by IBDV increases the susceptibility of chickens to other infectious diseases and reduces immune responses to vaccination [6]. Therefore, elucidating the interactions between IBDV and its host is a matter of urgency to determine effective strategies for preventing and controlling IBDV infections.

High-throughput sequencing technology and proteomic approaches have been used for this purpose, and various types of non-coding RNAs (ncRNAs) are also gaining increasing attention in this regard [7], including long non-coding RNAs (lncRNAs) [8] and circular RNAs (circRNAs) [9]. lncRNAs are longer than 200 nucleotides in length and can regulate gene expression at various levels, including epigenetics and transcriptional and post-transcriptional regulation. lncRNAs typically occur at low abundance and are frequently not conserved between species [10]. Moreover, lncRNAs have been displayed to affect viral replication by regulating the expression of antiviral-related key genes [11,12,13], indicating that lncRNA might play a crucial role in antiviral responses. CircRNAs form covalently closed continuous loops, which have been observed to be widely expressed in plants [14] and animals [15, 16]. More importantly, studies have demonstrated that circRNAs influence viral infection by inhibiting viral factors [17, 18]. Therefore, it is important to identify lncRNAs and circRNAs as well as their targets to understand the dynamics of gene regulation and effectively control the occurrence and development of disease.

Although many studies have been conducted to assess the effects of host lncRNAs on IBDV infection and its underlying molecular mechanisms [11, 19,20,21], the genome-wide expression effects and functional roles of lncRNAs and circRNAs during vvIBDV infection have not been examined so far. We investigated the expression profiles of mRNAs, lncRNAs, and circRNAs associated with vvIBDV infection of chicken’ BF and constructed lncRNA-mRNA and circRNA-miRNA-mRNA co-expression networks, which may provide valuable information for new therapeutic approaches to control this disease.

Results

Identification of lncRNAs and circRNAs in chicken bursa tissue

In total, 589,776,342 raw reads were obtained from the control and IBDV-infected bursa tissues. After data filtering and quality control, 584,284,990 clean reads of high quality were retained, and the proportion of clean reads ranged from 98.99 to 99.13% (Additional file 1: Table S1). Subsequently, reads mapping to ribosomal RNA (rRNA) were removed, and the proportion of retained reads ranged from 86.69 to 99.54% (Additional file 2: Table S2). After removing the rRNA, clean reads were mapped to the chicken reference genome (Table 1). The transcripts were screened using Cufflinks software (v. 2.1.1) [22] based on the location of the transcripts on the reference genome, a transcript length of ≥200 bp, and an exon count of ≥2. In total, 4324 known lncRNAs transcripts (Additional file 3) were obtained from chicken bursa tissue, including 1957 (45.3%) intergenic lncRNAs, 706 (16.3%) bidirectional lncRNAs, 1442 (33.3%) sense lncRNAs, and (5.1%) 219 lncRNAs of other types. Nevertheless, intronic lncRNA and antisense lncRNA were not detected in our study. Moreover, 1086 novel lncRNAs were screened based on three protein-coding potential software analyses (Fig. 1a; Additional file 4). The 1086 novel lncRNAs comprised 610 (56.2%) intergenic lncRNAs, 212 (19.5%) sense lncRNAs, 109 (10.0%) bidirectional lncRNAs, 76 (7.0%) antisense lncRNAs, and 79 (7.3%) lncRNAs of other types (Fig. 1b). No intronic lncRNAs were detected in the current study. Besides, anchor reads were mapped to the chicken reference genome, 7808 novel circRNAs were identified in the study (Additional file 5; Table 2).

Table 1 High quality clean reads compared with the reference genome
Full size table
Fig. 1
figure 1

Screening and classification of lncRNAs in chicken BF. a The novel lncRNAs were mainly classified as intergenic lncRNAs, bidirectional lncRNAs, and sense lncRNAs. b CNCI, CPC, and the SwissProt database were used to analyze the coding potential of the novel lncRNAs. RNAs identified by all three analytical tools were chosen as candidate lncRNAs. c The novel lncRNAs were mainly classified as intergenic lncRNAs, sense lncRNAs, bidirectional lncRNAs, and antisense lncRNAs

Full size image
Table 2 Anchor reads compared with the reference genome
Full size table

Differentially expressed profiles of mRNAs, lncRNAs, and circRNAs

The number of mRNAs, lncRNAs, and circRNAs shared between the infection group and the control group is shown in Fig. 2. Based on the criteria p < 0.05 and fold change > 1.5 or 2, 411 mRNAs (Additional file 6), 272 lncRNAs (Additional file 7), and 143 circRNAs (Additional file 8) were considered differentially expressed (Table 3). Moreover, 317 mRNAs were upregulated and 94 were downregulated in the treatment group (Figs. 3a and 4a). Moreover, 172 upregulated and 100 upregulated lncRNAs were identified (Figs. 3b and 4b), and 63 and 80 circRNAs were upregulated and downregulated, respectively (Figs. 3c and 4c). Differentially expressed mRNAs, lncRNAs, and circRNAs were used for cluster analysis. A heat map indicated that the control and IBDV-infected individuals produced two separate clusters (Fig. 5). The expression patterns of the mRNAs, lncRNAs, and circRNAs thus differed between the two groups, suggesting that the differentially expressed genes (DEGs) in the chicken bursa tissue infected with IBDV were significantly different from those in the non-infected tissue.

Fig. 2
figure 2

Venn diagram showing the number of overlapping genes in the IBDV-infected group and the control group. a mRNAs, b lncRNAs, and c circRNAs

Full size image
Table 3 The number of differentially expressed mRNA and lncRNA in IBDV treated group
Full size table
Fig. 3
figure 3

Histogram of the differentially expressed mRNAs (a), lncRNAs (b) and circRNAs (c) in the two groups. The red and green columns indicate the significantly upregulated and downregulated genes (p < 0.05), respectively

Full size image
Fig. 4
figure 4

Volcano plots of the differentially expressed mRNAs (a), lncRNAs (b), and circRNAs (c). The red and green dots indicate the significantly upregulated and downregulated genes (p < 0.05), respectively. The black dots indicate the genes that were not significantly differentially expressed (p > 0.05)

Full size image
Fig. 5
figure 5

Heatmap of differentially expressed mRNAs (a), lncRNAs (b), and circRNAs (c). The coloration gradient from blue to red indicates low to high expression level. All biological replicates were pooled to identify the DEGs based on a threshold fold change > 2 (mRNAs and circRNAs) or fold change > 1.5 (lncRNAs) at p < 0.05

Full size image

Comparison of mRNA and lncRNA characteristics

In total, 44188 mRNAs and 5410 lncRNAs were identified in all samples. To examine the differences in the mRNAs and lncRNAs, genetic structure and sequence conservation was compared, and the distribution of the transcript length of the lncRNAs differed slightly from that of the mRNAs (Fig. 6a). However, there were fewer exons in the lncRNAs than in the mRNAs (Fig. 6b). In addition, the open reading frames of the lncRNAs were shorter than those of the mRNAs (Fig. 6c).

Fig. 6
figure 6

Comparison of mRNA and lncRNA characteristics. a Number of exons in the mRNAs and lncRNAs. b Distribution of transcript lengths in the mRNAs and lncRNAs. The horizontal axis indicates the length of the transcripts, and the vertical axis indicates the abundance. c Number of open reading frames (ORFs) in the mRNAs and lncRNAs

Full size image

Functional annotation of DEGs

Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses were performed to analyze the 411 differentially expressed mRNAs and the target genes of the differentially expressed lncRNAs and circRNAs to examine the functions of gene dysregulation during IBDV infection. GO enrichment of the biological processes (BPs), cellular components (CCs), and molecular functions (MFs) of the differentially expressed mRNAs, lncRNAs, and circRNAs is shown in Figs. 7 and 8. BPs associated with the enriched mRNAs predominantly included biological regulations, cellular processes, and single-organism processes; enriched mRNA CCs included cells, cell parts, and organelles; and significantly enriched MFs included binding and catalytic activities (Fig. 7; Additional file 9). Based on the GO analyses of the differentially expressed lncRNA target genes, the most enriched BPs were cellular processes, biological regulation, and single-organism processes; the most enriched CCs were cells, cell parts, and organelles; and the most enriched MFs were binding and catalytic activity (Fig. 8a). GO enrichment analysis was performed for the antisense, cis, and trans roles of the target genes of the lncRNA, showing that target genes were also mainly enriched in cellular processes, cells, and binding (Additional file 10: Figure S1). The circRNA results were consistent with those of the lncRNAs (Fig. 8b).

Fig. 7
figure 7

Gene ontology (GO) analysis of the differentially expressed mRNAs in IBDV-infected chicken BF. a-c Directed acyclic graph showing the significantly enriched biological processes, cellular components, and molecular functions of the differentially expressed mRNAs. d Number of genes in GO terms are shown in the histogram

Full size image
Fig. 8
figure 8

Gene ontology analysis of the differentially expressed lncRNAs (a) and circRNAs (b) in the two groups. The green, red, and blue column indicate biological processes, cellular components, and molecular functions, respectively

Full size image

KEGG is the main public pathway-related database, and it has been used previously to determine the significantly enriched pathways of dysregulated gene products during IBDV infection [23, 24]. The top 20 pathways associated with the mRNAs, lncRNAs, and circRNAs are shown in Fig. 9. The results show that mRNAs differentially expressed owing to IBDV infection were associated with the JAK-STAT signaling pathway, the NOD-like receptor signaling pathway, and the cytokine-cytokine receptor interaction signaling pathway, among others (Fig. 9a; Additional file 11). According to the KEGG analyses of the target genes of the differentially expressed lncRNAs, spliceosome, JAK-STAT signaling pathway, ribosome, and Toll-like receptor signaling pathway were enriched (Fig. 9b). The target genes of lncRNA antisense, cis, and trans-regulation were mainly enriched in the notch signaling pathway, the glycosphingolipid biosynthesis-lacto complement, and spliceosome, respectively (Additional file 12: Figure S2). In the circRNAs, the MAPK signaling pathway, protein processing in the endoplasmic reticulum, and ubiquitin-mediated proteolysis were identified as predominantly enriched KEGG pathways (Fig. 9c).

Fig. 9
figure 9

Kyoto Encyclopedia of Genes and Genomes pathway enrichment of the differentially expressed mRNAs (a), lncRNAs (b), and circRNAs (c) in the two groups. The vertical axis indicates the pathways, and the horizontal axis represents the Rich factor. The dot size indicates the number of DEGs in the pathway, and the coloration corresponds to the Q-value range

Full size image

Target predictions and lncRNA-mRNA and circRNA-miRNA-mRNA regulatory network analysis

Generally, lncRNAs may exert their effects by regulating their target genes. In the study, we predicted the potential target genes of lncRNAs and constructed lncRNA-mRNA regulatory networks (Fig. 10). A total of 2101 pairs of lncRNA-target genes containing 51 lncRNAs and 342 mRNAs (Additional file 13) were detected. In the lncRNA-mRNA network, LOC112532624 (XR_003075708.1) with the most significant difference (fold change = 703.3). Moreover, LOC107053928 (XR_001467739.2), LOC107054815 (XR_001469507.2), LOC107053352 (XR_001466515.2), and LOC107053557 (XR_001466920.2 97) was predicted to regulate the expression of 110, 100, 100 and 97 target genes, respectively, involved in antiviral responses, including STAT1, STAT3, STAT4, TRIM25, and IFIH1 (Additional file 14: Table S3). More importantly, a total of 44 differentially expressed lncRNAs were found to be co-expressed with STAT1, a key antiviral marker molecule. The target prediction indicated that a lncRNA had multiple target genes, and a target gene was also targeted by multiple lncRNAs.

Fig. 10
figure 10

LncRNA-mRNA regulatory networks. The red ellipses indicate the differentially expressed lncRNAs, and the green ellipses indicate the differentially expressed mRNAs

Full size image

Recent studies evidenced that circRNAs commonly play an important regulatory role as miRNA sponges in several diseases [16, 25,26,27,28]. Here, the potential miRNA targets of 30 differentially expressed circRNAs were predicted based on complementary base pairing, with 36 miRNAs identified (Fig. 11; Additional file 15). In the circRNA-miRNA-mRNA network, circRNAs would indirectly regulate 25 chicken mRNAs, such as STAT1/4 and IRF1/7, indicating that these circRNAs might play a critical role in regulating vvIBDV-infection.

Fig. 11
figure 11

CircRNA-miRNA-mRNA regulatory networks. The red triangles, blue diamonds, and purple ellipses indicate the differentially expressed circRNA, miRNAs, and mRNAs, respectively

Full size image

Validation of differential gene expression by quantitative PCR

To validate the accuracy of the RNA-sequencing results, 10 differentially expressed mRNAs, lncRNAs, and circRNAs were selected and quantified by reverse-transcription quantitative PCR (RT-qPCR; Fig. 12); these RT-qPCR results showed trends similar to those of RNA sequencing.

Fig. 12
figure 12

Validation of the differentially expressed mRNAs (a), lncRNAs (b), and circRNAs (c) by RT-qPCR. RT-qPCR experiments were performed in triplicate

Full size image

Discussion

lncRNAs and circRNAs has emerged as regulator of gene expression and play important roles in various diseases [29, 30]. Many lncRNAs and circRNAs of chicken have been identified and analyzed by high-throughput sequencing technology, and suggesting that aberrantly expressed lncRNAs and circRNAs in cells or tissues play a crucial role in virus-host interactions [31,32,33,34]. Interesting, classical IBDV infection affects lncRNAs expression in DF-1 cell has been testified [11]. Infectious bursal disease caused by IBDV is one of the most important infectious diseases that severely affect the poultry industry. Infection with vvIBDV strains results in high mortality in chickens, and vvIBDV appears to have emerged as the predominant clinical disease type in nearly all poultry-producing regions of the world. Therefore, the BF of experimentally vvIBDV-infected chickens and control individuals were collected to study the differential expression profiles of mRNAs, lncRNAs, and circRNAs. In this study, 411 mRNAs, 272 lncRNAs, and 143 circRNAs were considered differentially expressed. The results indicated that these DEGs may play a significant role in the vvIBDV infection process, suggesting that they may include potential therapeutic targets for treating IBDV infections.

Host cells counteract the invasion of viral particles through a cascade of mechanisms, and IBDV infection elicits the expression of various genes. Previous study reported that the expression of interferon regulatory factor 8 (IRF8), TRIM25, IFIT1, Mx1; of STAT1, STAT4, and of Toll-like receptors, including TLR3 and TLR7, is increased in IBDV-infected chickens or DF-1 or DT40 cells [4, 7, 11, 35,36,37,38]. In our study, the JAK-STAT signaling pathway, the NOD-like receptor signaling pathway, the cytokine-cytokine receptor interaction signaling pathway, apoptosis, the chemokine signaling pathway, and the Toll-like receptor signaling pathway were significantly enriched according to the KEGG enrichment analyses of the respective differentially expressed mRNAs. IBDV exploits these pathways to induce the expression of STAT1, STAT3, STAT4, TRIM25, IFIT1, and Mx1 in the bursal tissue, and our results were in line with those of previous studies. In STAT1, a member of the STAT family, phosphorylation induces the expression of interferon-stimulating genes through a series of signal transduction steps to elicit antiviral mechanisms [39, 40]. Interestingly, STAT1 can be regulated by many differentially expressed lncRNAs, suggesting that STAT1 may be an important regulator during IBDV infection in chickens. Additionally, STAT3, STAT4, TRIM25, IFIT1, and Mx1 may be involved in interactions between vvIBDV and the host.

The potential functions of lncRNAs are commonly predicted according to their target genes. Increasing evidence suggests that lncRNAs have important roles in adaptive or innate immune responses [41,42,43]. In a previous study, we found that loc107051710 suppressed the replication of IBDV by upregulating type I interferon production [11]. In the current study, it is worth noting that LOC107053928, LOC107054815, LOC107053352, and LOC107053557 were identified as regulated various target genes associated with immunomodulation, including STAT1, STAT3, STAT4, TRIM25, and IFIH1. TRIM25 is a member of the tripartite motif family of E3 ubiquitin ligases and has been demonstrated to play an important role in RIG-I antiviral pathways. Reportedly, TRIM25 can promote transcriptional upregulation of type I interferons (IFNs) by binding viral RNA-activated RIG-I [44, 45]. It has been well established that type I IFNs plays a crucial role in the antiviral processes, where type I IFNs can activate the STAT1 phosphorylation after binding its ligand to induce interferon-stimulating genes (ISGs) expression, such as Mx1, OAS, and IFIH1, through the JAK-STAT signaling pathway [40, 46]. Therefore, this implied that these lncRNAs and their target genes, STAT1, STAT3, STAT4, IFIH1 and TRIM25, might play a vital role in the IBD anti-viral response. Intriguingly, gene targeting studies have shown that STAT1 target genes can promote antagonizing proliferation and inflammation; however, STAT3 was the opposite [47, 48]. Therefore, the expression levels of STAT1 and STAT3 may reflect the balanced response of the organism. We believe that further insight into the roles of LOC107053928, LOC107054815, LOC107053352, and LOC107053557 is very important for understanding the complex regulatory mechanism of gene expression in response to vvIBDV infection in chicken BF.

CircRNAs, a newly discovered class of ncRNAs, can affect the prognosis of diseases, especially tumors [49, 50] by regulating the activity of target genes and by participating in the regulation of gene transcription and protein production [51, 52]. In the current study, 63 upregulated and 80 downregulated circRNAs were identified, and their expression levels were generally lower than those of the mRNAs and lncRNAs. Most circRNAs in normal and cancer tissues in humans also occur at low abundance and may thus be by-products of pre-mRNA splicing [53, 54]. Additionally, circRNA can act as a competing endogenous RNA (ceRNAs) that impair miRNA activity through sequestration, thereby upregulating miRNA target gene expression [55]. In the study, circRNA-miRNA-mRNA network was constructed and showed that the regulatory relationships between circRNAs, miRNAs and target mRNAs were complex. In the network, circRNAs novel_circ_000574 and novel_circ_001469 were expressed 131,072-fold and 286,862-fold higher levels, respectively, during IBDV infection. Their potentially affected target genes involved in immune-related included STAT1 and IRF7 by binding to miR-1587-x and miR-4507-y respectively, suggesting that these circRNAs also play an important role in IBDV pathogenicity.

Conclusion

Many signaling pathways were found to be involved in IBDV infection, particularly the JAK-STAT signaling pathways; however, further research is needed to assess their effects on the pathogenesis of IBDV infection. Additionally, we constructed lncRNA-mRNA and circRNA-miRNA-mRNA co-expression networks and predicted that LOC107053928, LOC107054815, LOC107053352, and LOC107053557 may affect viral replication by regulating STAT1, STAT3, STAT4, IFIH1 and TRIM25 expression. Our results provide new insights into the pathways and mechanisms that mediate host immune responses to vvIBDV. Further studies are needed to explore the biological functions of LOC107053928, LOC107054815, LOC107053352, and LOC107053557 during vvIBDV infection.

Methods

Study animals and viruses

Three-week-old specific-pathogen-free (SPF) White Leghorn chickens were obtained from the Harbin Institute of Veterinary Medicine (Harbin, China). A vvIBDV strain isolated from chicken BF and maintained in our laboratory was used in this study (strain LJ-5).

Sample collection and preparation

Eighteen chickens were randomly assigned to two groups (control group and vvIBDV infection group) with nine individuals per group, and the two groups were housed independently. Chickens in the control group were injected with phosphate-buffered saline (PBS; PH = 7.4), and the vvIBDV infection group chickens were inoculated with the LJ-5 strain virus at a dose of 103 × 50% embryo lethal dose (ELD) per 0.2 mL via eye-nose drops [56, 57]. The chickens were housed at an animal facility under pathogen-free conditions and were provided a standard diet and water. Moreover, these chickens were observed 2–3 times daily. On the third day after inoculation, three chickens from each group were randomly selected and euthanized via T-61 intravenously (0.4 ml/kg) to isolate the BF. The animal remains were treated innocuously according to the requirements of school animal welfare management. The procedures used in this experiment were approved by the Institutional Animal Care and Use Committee of Northeast Agricultural University (2016NEFU-315, 13 April 2017). All sections of this report adhere to the ARRIVE Guidelines for reporting animal research (Additional file 16). The National Centre for the Replacement, Refinement, and Reduction of Animals in Research (NC3Rs) is an independent scientific organization that minimize the use of animals in research and improve animal welfare and help overcome the challenges and limitations of the use of animals in research to the benefit of the whole bioscience community [58]. Samples were placed in liquid nitrogen and stored at − 80 °C. All instruments were treated with DEPC before use to remove RNases.

RNA isolation, library preparation, and sequencing

Total RNA was isolated using Trizol reagent (Life Technologies, Carlsbad, USA), and the purity and integrity of the RNA were measured using a Nanodrop micro-spectrophotometer (Thermo Fisher Scientific, Waltham, USA) and an Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, USA), respectively. RNA samples qualified for analyses under the following conditions: concentration ≥ 100 ng/μL, RNA integrity number ≥ 7.0, and RNA 260/280 ratio of 1.8–2.0. For RNA sequencing, 3 μg RNA per sample was used, and ribosomal RNA was removed to generate sequencing libraries using a NEBNext® Ultra™ RNA Library Prep kit (NEB E7490L; New England Biolabs, Ipswich, USA). First-strand complementary DNA was synthesized using random hexamer primers and ProtoScript II Reverse Transcriptase, and the second strand was generated using a Second Strand Synthesis Enzyme Mix. Uracil N Glycosylase was used to digest complementary DNA, and Agencourt AMPure XP Beads (Beckman Coulter, Brea, USA) were used to purify library fragments to retain DNA fragments of 150–200 bp. The quality of the libraries was evaluated using a High Sensitivity DNA assay Kit (Agilent Technologies), and the libraries were sequenced on an Illumina HiSeq 4000 platform by a commercial service provider (Gene Denovo Biotechnology Co.; Guangzhou, China). For circRNA sequencing, RNase R (EPICENTRE Biotechnologies, Madison, USA) was used on rRNA-depleted RNAs to remove linear RNAs before sequencing library preparation.

Filtration of raw sequencing reads

Raw reads containing adapters or bases of low quality affect assembly and analysis; therefore, to obtain clean high-quality sequences, raw reads were filtered to remove adapters, low-quality reads, and poly-N reads. Clean high-quality data was used for further analyses.

RNA sequence mapping and transcriptome assembly

The short reads alignment tool Bowtie 2 (v. 2.2.8) [59] was used to map the sequences using an rRNA database. After removing the rRNA reads, the sequences were mapped to the Gallus gallus GRCg6a reference genome using TopHat2 (v. 2.1.1) [60]; Cufflinks (v. 2.1.1) [22] was used to assemble mapped reads independently with a reference annotation-based transcript assembly technique. Raw sequencing data were made available in the NCBI short reads archive.

lncRNA and circRNA predictions

CNCI (version: 2.0) [61] and CPC [62] were used to assess the potential protein-coding functions of novel lncRNAs. Transcripts were mapped using the SwissProt database to assess protein annotation. The intersection of the respective results was chosen as lncRNAs. Phast (v. 1.3) [63], phyloFit [64], and PhastCons [65] were used to evaluate sequence conservation in the transcripts, calculate phylogenetic models among species, and compute the conservation scores of coding genes and lncRNAs, respectively. For circRNA, unmapped reads were extracted from the above results, and the ends of the unmapped reads were intercepted (default 20 bp) to identify the anchor reads, which were then mapped to the reference genome. The results were processed using find_circ software to predict circRNAs.

Prediction of lncRNA and circRNA target genes

lncRNAs are involved in many post-transcriptional regulation processes, as are some other small RNAs such as microRNA (miRNA), the regulation of which depends on complementary base pairing. Some antisense lncRNAs may regulate gene silencing, transcription, and mRNA stability. To assess the interactions between lncRNAs and mRNAs, RNAplex software [66] was used to predict the complementary correlation of antisense lncRNAs and mRNAs. This software program includes the Vienna RNA package, and best base pairing predictions were based on a calculation of minimum free energy in the thermodynamic structure. Moreover, the cis and trans target genes of the differentially expressed lncRNAs were predicted. For the cis target genes, the mRNAs and the genomic location of the lncRNAs were mapped. We searched coding genes less than 100 kb up- or downstream of each lncRNA and then analyzed their functions. For the trans targets, correlations between the lncRNA and protein-coding gene expression was analyzed using Pearson’s correlation coefficient, and protein-coding genes with absolute correlation values > 0.94 were screened. The differentially expressed lncRNAs (fold change > 1.5 and p < 0.05) and predicted target genes (fold change > 2 and p < 0.05) were chosen to construct a co-expression regulatory network, which was visualized using Cytoscape software (v3.6.0) [67]. The binding sites of the miRNAs on the circRNAs and mRNAs were predicted using mirTarBase software [68], and the target relationships of miRNAs-mRNAs and miRNAs-circRNAs were assessed accordingly. Subsequently, we synthesized targeted relationships between the miRNAs and mRNAs and between the miRNAs and circRNAs to identify the mRNAs associated with the circRNAs.

Functional annotation of DEGs

To understand the functions of the differentially expressed transcripts, including the mRNA and the target genes of the lncRNAs and circRNAs, the GO (http://www.geneontology.org/) and KEGG (http://www.kegg.jp/kegg/) databases were used to perform GO and KEGG enrichment analyses, respectively. GO terms and KEGG pathways with p values < 0.05 were considered significantly enriched.

RT-qPCR validation of the RNA-sequencing results

To assess the accuracy of the sequencing results, five upregulated and five downregulated mRNAs, lncRNAs, and circRNAs were selected and quantified by RT-qPCR. Primers were designed using oligo6 software (Tables 4, 5 and 6). Total RNA was extracted from vvIBDV-infected and uninfected bursal tissue. RT-qPCR was performed using SYBR Green Master (Roche, Mannheim, Germany) and an ABI 7500 Real-Time PCR system (Applied Biosystems, Waltham, USA). β-Actin was used as an internal control of the mRNA. The 2−ΔΔCt method was used to calculate the relative expression levels of the target genes. Experiments were conducted using three replicates.

Table 4 RT-qPCR primers used for verification of mRNA results
Full size table
Table 5 RT-qPCR primers used for verification of lncRNA results
Full size table
Table 6 RT-qPCR primers used for verification of circRNA results
Full size table

Statistical analyses

GraphPad Prism5 software was used to analyze the RT-qPCR results with one-way ANOVA. Data are shown as the means ± standard error of the mean. Statistical significance was considered at p < 0.05.

Availability of data and materials

The National Center for Biotechnology Information (NCBI) Nucleotide database accession number for IBDV LJ-5 strain is MT133301 (A fragments: https://www.ncbi.nlm.nih.gov/nuccore/MT133301.1/) and MT446362 (B fragments: https://www.ncbi.nlm.nih.gov/nuccore/MT446362).

The raw data sets supporting the results of this article are available in the NCBI short reads archive and accession number is PRJNA635782. For information linking and citing, please refer to: https://www.ncbi.nlm.nih.gov/search/all/?term=PRJNA635782.

Abbreviations

IBDV:

Infectious bursal disease virus

BF:

Bursa of Fabricius

lncRNAs:

Long non-coding RNAs

circRNAs:

Circular RNAs

NC3Rs:

The National Centre for the Replacement, Refinement, and Reduction of Animals in Research

DEGs:

Differentially expressed genes

GO:

Gene Ontology

KEGG:

Kyoto Encyclopedia of Genes and Genomes

BP:

Biological processe

CC:

Cellular component

MF:

Molecular function

IRF8:

Regulatory factor 8

IFNs:

Interferons

ISGs:

Interferon-stimulating genes

ceRNAs:

Competing endogenous RNA

miRNA:

microRNA

References

  1. Li Z, Wang Y, Li X, Li X, Cao H, Zheng SJ. Critical roles of glucocorticoid-induced leucine zipper in infectious bursal disease virus (IBDV)-induced suppression of type I interferon expression and enhancement of IBDV growth in host cells via interaction with VP4. J Virol. 2013;87(2):1221–31.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  2. Vukea PR, Willows-Munro S, Horner RF, Coetzer TH. Phylogenetic analysis of the polyprotein coding region of an infectious South African bursal disease virus (IBDV) strain. Infect Genet Evol. 2014;21:279–86.

    Article  CAS  PubMed  Google Scholar 

  3. McFerran JB, McNulty MS, McKillop ER, Connor TJ, McCracken RM, Collins DS, Allan GM. Isolation and serological studies with infectious bursal disease viruses from fowl, turkeys and ducks: demonstration of a second serotype. Avian Pathol. 1980;9(3):395–404.

    Article  CAS  PubMed  Google Scholar 

  4. Ruby T, Whittaker C, Withers DR, Chelbi-Alix MK, Morin V, Oudin A, Young JR, Zoorob R. Transcriptional profiling reveals a possible role for the timing of the inflammatory response in determining susceptibility to a viral infection. J Virol. 2006;80(18):9207–16.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  5. Becht H, Muller H. Infectious bursal disease--B cell dependent immunodeficiency syndrome in chickens. Behring Inst Mitt. 1991;89:217–25.

    Google Scholar 

  6. Berg TP. Acute infectious bursal disease in poultry: a review. Avian Pathol. 2000;29(3):175–94.

    Article  CAS  PubMed  Google Scholar 

  7. Li YP, Handberg KJ, Juul-Madsen HR, Zhang MF, Jorgensen PH. Transcriptional profiles of chicken embryo cell cultures following infection with infectious bursal disease virus. Arch Virol. 2007;152(3):463–78.

    Article  CAS  PubMed  Google Scholar 

  8. Bhan A, Mandal SS. Long noncoding RNAs: emerging stars in gene regulation, epigenetics and human disease. ChemMedChem. 2014;9(9):1932–56.

    Article  CAS  PubMed  Google Scholar 

  9. Qu S, Yang X, Li X, Wang J, Gao Y, Shang R, Sun W, Dou K, Li H. Circular RNA: a new star of noncoding RNAs. Cancer Lett. 2015;365(2):141–8.

    Article  CAS  PubMed  Google Scholar 

  10. Derrien T, Johnson R, Bussotti G, Tanzer A, Djebali S, Tilgner H, Guernec G, Martin D, Merkel A, Knowles DG, et al. The GENCODE v7 catalog of human long noncoding RNAs: analysis of their gene structure, evolution, and expression. Genome Res. 2012;22(9):1775–89.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  11. Huang X, Xu Y, Lin Q, Guo W, Zhao D, Wang C, Wang L, Zhou H, Jiang Y, Cui W, et al. Determination of antiviral action of long non-coding RNA loc107051710 during infectious bursal disease virus infection due to enhancement of interferon production. Virulence. 2020;11(1):68–79.

    Article  CAS  PubMed  Google Scholar 

  12. Nishitsuji H, Ujino S, Yoshio S, Sugiyama M, Mizokami M, Kanto T, Shimotohno K. Long noncoding RNA #32 contributes to antiviral responses by controlling interferon-stimulated gene expression. Proc Natl Acad Sci U S A. 2016;113(37):10388–93.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  13. Carnero E, Barriocanal M, Segura V, Guruceaga E, Prior C, Börner K, Grimm D, Fortes P. Type I interferon regulates the expression of long non-coding RNAs. Front Immunol. 2014;5:548.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  14. Sanger HL, Klotz G, Riesner D, Gross HJ, Kleinschmidt AK. Viroids are single-stranded covalently closed circular RNA molecules existing as highly base-paired rod-like structures. Proc Natl Acad Sci U S A. 1976;73(11):3852–6.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  15. Capel B, Swain A, Nicolis S, Hacker A, Walter M, Koopman P, Goodfellow P, Lovell-Badge R. Circular transcripts of the testis-determining gene Sry in adult mouse testis. Cell. 1993;73(5):1019–30.

    Article  CAS  PubMed  Google Scholar 

  16. Memczak S, Jens M, Elefsinioti A, Torti F, Krueger J, Rybak A, Maier L, Mackowiak SD, Gregersen LH, Munschauer M, et al. Circular RNAs are a large class of animal RNAs with regulatory potency. Nature. 2013;495(7441):333–8.

    Article  CAS  PubMed  Google Scholar 

  17. Tagawa T, Gao S, Koparde VN, Gonzalez M, Spouge JL, Serquiña AP, Lurain K, Ramaswami R, Uldrick TS, Yarchoan R, et al. Discovery of Kaposi’s sarcoma herpesvirus-encoded circular RNAs and a human antiviral circular RNA. Proc Natl Acad Sci U S A. 2018;115(50):12805–10.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  18. Wang M, Yu F, Wu W, Zhang Y, Chang W, Ponnusamy M, Wang K, Li P. Circular RNAs: a novel type of non-coding RNA and their potential implications in antiviral immunity. Int J Biol Sci. 2017;13(12):1497–506.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  19. Ouyang W, Wang YS, Du XN, Liu HJ, Zhang HB. gga-miR-9* inhibits IFN production in antiviral innate immunity by targeting interferon regulatory factor 2 to promote IBDV replication. Vet Microbiol. 2015;178(1–2):41–9.

    Article  CAS  PubMed  Google Scholar 

  20. Wei O, Jing Q, Qun-xing P, Jing-yu W, Xing-xia X, Xiao-li W, Yu-mei Z, Yong-shan W. gga-miR-142-5p attenuates IRF7 signaling and promotes replication of IBDV by directly targeting the chMDA5′s 3′ untranslated region. Vet Microbiol. 2018;221:74–80.

    Article  CAS  Google Scholar 

  21. Wang Y, Sun H, Shen P, Zhang X, Xia X, Xia B. Effective inhibition of replication of infectious bursal disease virus by miRNAs delivered by vectors and targeting the VP2 gene. J Virol Methods. 2010;165(2):127–32.

    Article  CAS  PubMed  Google Scholar 

  22. Trapnell C, Roberts A, Goff L, Pertea G, Kim D, Kelley DR, Pimentel H, Salzberg SL, Rinn JL, Pachter L. Differential gene and transcript expression analysis of RNA-seq experiments with TopHat and Cufflinks. Nat Protoc. 2012;7(3):562–78.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  23. Mao X, Cai T, Olyarchuk JG, Wei L. Automated genome annotation and pathway identification using the KEGG Orthology (KO) as a controlled vocabulary. Bioinformatics. 2005;21(19):3787–93.

    Article  CAS  PubMed  Google Scholar 

  24. Kanehisa M, Araki M, Goto S, Hattori M, Hirakawa M, Itoh M, Katayama T, Kawashima S, Okuda S, Tokimatsu T, et al. KEGG for linking genomes to life and the environment. Nucleic Acids Res. 2008;36(Database issue):D480–4.

    CAS  PubMed  Google Scholar 

  25. Hansen TB, Jensen TI, Clausen BH, Bramsen JB, Finsen B, Damgaard CK, Kjems J. Natural RNA circles function as efficient microRNA sponges. Nature. 2013;495(7441):384–8.

    Article  CAS  PubMed  Google Scholar 

  26. Lukiw WJ. Circular RNA (circRNA) in Alzheimer’s disease (AD). Front Genet. 2013;4:307.

    PubMed  PubMed Central  Google Scholar 

  27. Xu S, Zhou L, Ponnusamy M, Zhang L, Dong Y, Zhang Y, Wang Q, Liu J, Wang K. A comprehensive review of circRNA: from purification and identification to disease marker potential. PeerJ. 2018;6:e5503.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  28. Zhang HD, Jiang LH, Sun DW, Hou JC, Ji ZL. CircRNA: a novel type of biomarker for cancer. Breast Cancer. 2018;25(1):1–7.

    Article  PubMed  Google Scholar 

  29. Elling R, Chan J, Fitzgerald KA. Emerging role of long noncoding RNAs as regulators of innate immune cell development and inflammatory gene expression. Eur J Immunol. 2016;46(3):504–12.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  30. Seeger C, Mason WS. Molecular biology of hepatitis B virus infection. Virology. 2015;479-480:672–86.

    Article  CAS  PubMed  Google Scholar 

  31. Imam H, Bano AS, Patel P, Holla P, Jameel S. The lncRNA NRON modulates HIV-1 replication in a NFAT-dependent manner and is differentially regulated by early and late viral proteins. Sci Rep. 2015;5:8639.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  32. Li Z, Chao TC, Chang KY, Lin N, Patil VS, Shimizu C, Head SR, Burns JC, Rana TM. The long noncoding RNA THRIL regulates TNFα expression through its interaction with hnRNPL. Proc Natl Acad Sci U S A. 2014;111(3):1002–7.

    Article  CAS  PubMed  Google Scholar 

  33. Wang L, You Z, Wang M, Yuan Y, Liu C, Yang N, Zhang H, Lian L. Genome-wide analysis of circular RNAs involved in Marek’s disease tumourigenesis in chickens. RNA Biol. 2020;17(4):517–27.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  34. Zhang Y, Wang L, Qiu L, Pan R, Bai H, Jiang Y, Wang Z, Bi Y, Chen G, Chang G. Expression patterns of novel circular RNAs in chicken cells after avian leukosis virus subgroup J infection. Gene. 2019;701:72–81.

    Article  CAS  PubMed  Google Scholar 

  35. Eldaghayes I, Rothwell L, Williams A, Withers D, Balu S, Davison F, Kaiser P. Infectious bursal disease virus: strains that differ in virulence differentially modulate the innate immune response to infection in the chicken bursa. Viral Immunol. 2006;19(1):83–91.

    Article  CAS  PubMed  Google Scholar 

  36. Hui RK, Leung FC. Differential expression profile of chicken embryo fibroblast DF-1 cells infected with cell-adapted infectious bursal disease virus. PLoS One. 2015;10(6):e0111771.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  37. Ou C, Wang Q, Zhang Y, Kong W, Zhang S, Yu Y, Ma J, Liu X, Kong X. Transcription profiles of the responses of chicken bursae of Fabricius to IBDV in different timing phases. Virol J. 2017;14(1):93.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  38. Quan R, Zhu S, Wei L, Wang J, Yan X, Li Z, Liu J. Transcriptional profiles in bursal B-lymphoid DT40 cells infected with very virulent infectious bursal disease virus. Virol J. 2017;14(1):7.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  39. Horvath CM, Stark GR, Kerr IM, Darnell JE Jr. Interactions between STAT and non-STAT proteins in the interferon-stimulated gene factor 3 transcription complex. Mol Cell Biol. 1996;16(12):6957–64.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  40. Schoggins JW, Wilson SJ, Panis M, Murphy MY, Jones CT, Bieniasz P, Rice CM. A diverse range of gene products are effectors of the type I interferon antiviral response. Nature. 2011;472(7344):481–5.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  41. Heward JA, Lindsay MA. Long non-coding RNAs in the regulation of the immune response. Trends Immunol. 2014;35(9):408–19.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  42. Wang P, Xue Y, Han Y, Lin L, Wu C, Xu S, Jiang Z, Xu J, Liu Q, Cao X. The STAT3-binding long noncoding RNA lnc-DC controls human dendritic cell differentiation. Science. 2014;344(6181):310–3.

    Article  CAS  PubMed  Google Scholar 

  43. Liu W, Ding C. Roles of LncRNAs in viral infections. Front Cell Infect Microbiol. 2017;7:205.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  44. Gack MU, Shin YC, Joo CH, Urano T, Liang C, Sun L, Takeuchi O, Akira S, Chen Z, Inoue S, et al. TRIM25 RING-finger E3 ubiquitin ligase is essential for RIG-I-mediated antiviral activity. Nature. 2007;446(7138):916–20.

    Article  CAS  PubMed  Google Scholar 

  45. Gack MU, Kirchhofer A, Shin YC, Inn KS, Liang C, Cui S, Myong S, Ha T, Hopfner KP, Jung JU. Roles of RIG-I N-terminal tandem CARD and splice variant in TRIM25-mediated antiviral signal transduction. Proc Natl Acad Sci U S A. 2008;105(43):16743–8.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  46. Randall RE, Goodbourn S. Interferons and viruses: an interplay between induction, signalling, antiviral responses and virus countermeasures. J Gen Virol. 2008;89(Pt 1):1–47.

    Article  CAS  PubMed  Google Scholar 

  47. Meraz MA, White JM, Sheehan KC, Bach EA, Rodig SJ, Dighe AS, Kaplan DH, Riley JK, Greenlund AC, Campbell D, et al. Targeted disruption of the Stat1 gene in mice reveals unexpected physiologic specificity in the JAK-STAT signaling pathway. Cell. 1996;84(3):431–42.

    Article  CAS  PubMed  Google Scholar 

  48. Durbin JE, Hackenmiller R, Simon MC, Levy DE. Targeted disruption of the mouse Stat1 gene results in compromised innate immunity to viral disease. Cell. 1996;84(3):443–50.

    Article  CAS  PubMed  Google Scholar 

  49. Kong Z, Wan X, Zhang Y, Zhang P, Zhang Y, Zhang X, Qi X, Wu H, Huang J, Li Y. Androgen-responsive circular RNA circSMARCA5 is up-regulated and promotes cell proliferation in prostate cancer. Biochem Biophys Res Commun. 2017;493(3):1217–23.

    Article  CAS  PubMed  Google Scholar 

  50. Cui X, Niu W, Kong L, He M, Jiang K, Chen S, Zhong A, Li W, Lu J, Zhang L. hsa_circRNA_103636: potential novel diagnostic and therapeutic biomarker in major depressive disorder. Biomark Med. 2016;10(9):943–52.

    Article  CAS  PubMed  Google Scholar 

  51. Chen LL, Yang L. Regulation of circRNA biogenesis. RNA Biol. 2015;12(4):381–8.

    Article  PubMed  PubMed Central  Google Scholar 

  52. Li J, Yang J, Zhou P, Le Y, Zhou C, Wang S, Xu D, Lin HK, Gong Z. Circular RNAs in cancer: novel insights into origins, properties, functions and implications. Am J Cancer Res. 2015;5(2):472–80.

    CAS  PubMed  PubMed Central  Google Scholar 

  53. Guo JU, Agarwal V, Guo H, Bartel DP. Expanded identification and characterization of mammalian circular RNAs. Genome Biol. 2014;15(7):409.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  54. Salzman J, Chen RE, Olsen MN, Wang PL, Brown PO. Cell-type specific features of circular RNA expression. PLoS Genet. 2013;9(9):e1003777.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  55. Thomson DW, Dinger ME. Endogenous microRNA sponges: evidence and controversy. Nat Rev Genet. 2016;17(5):272–83.

    Article  CAS  PubMed  Google Scholar 

  56. Liu J, Zhou J, Kwang J. Antigenic and molecular characterization of recent infectious bursal disease virus isolates in China. Virus Genes. 2002;24(2):135–47.

    Article  PubMed  Google Scholar 

  57. Ou CB, Pan Q, Chen X, Hou N, He C. Protocatechuic acid, a new active substance against the challenge of avian infectious bursal disease virus. Poult Sci. 2012;91(7):1604–9.

    Article  CAS  PubMed  Google Scholar 

  58. Burden N, Chapman K, Sewell F, Robinson V. Pioneering better science through the 3Rs: an introduction to the national centre for the replacement, refinement, and reduction of animals in research (NC3Rs). J Am Assoc Lab Anim Sci. 2015;54(2):198–208.

    PubMed  PubMed Central  Google Scholar 

  59. Langmead B, Salzberg SL. Fast gapped-read alignment with bowtie 2. Nat Methods. 2012;9(4):357–9.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  60. Kim D, Pertea G, Trapnell C, Pimentel H, Kelley R, Salzberg SL. TopHat2: accurate alignment of transcriptomes in the presence of insertions, deletions and gene fusions. Genome Biol. 2013;14(4):R36.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  61. Sun L, Luo H, Bu D, Zhao G, Yu K, Zhang C, Liu Y, Chen R, Zhao Y. Utilizing sequence intrinsic composition to classify protein-coding and long non-coding transcripts. Nucleic Acids Res. 2013;41(17):e166.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  62. Kong L, Zhang Y, Ye ZQ, Liu XQ, Zhao SQ, Wei L, Gao G. CPC: assess the protein-coding potential of transcripts using sequence features and support vector machine. Nucleic Acids Res. 2007;35(Web Server issue):W345–9.

    Article  PubMed  PubMed Central  Google Scholar 

  63. Siepel A, Bejerano G, Pedersen JS, Hinrichs AS, Hou M, Rosenbloom K, Clawson H, Spieth J, Hillier LW, Richards S, et al. Evolutionarily conserved elements in vertebrate, insect, worm, and yeast genomes. Genome Res. 2005;15(8):1034–50.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  64. Darling AE, Jospin G, Lowe E, Matsen FA, Bik HM, Eisen JA. PhyloSift: phylogenetic analysis of genomes and metagenomes. PeerJ. 2014;2:e243.

    Article  PubMed  PubMed Central  Google Scholar 

  65. Ramani R, Krumholz K, Huang YF, Siepel A. PhastWeb: a web interface for evolutionary conservation scoring of multiple sequence alignments using phastCons and phyloP. Bioinformatics. 2019;35(13):2320–2.

    Article  CAS  PubMed  Google Scholar 

  66. Tafer H, Hofacker IL. RNAplex: a fast tool for RNA-RNA interaction search. Bioinformatics. 2008;24(22):2657–63.

    Article  CAS  PubMed  Google Scholar 

  67. Shannon P, Markiel A, Ozier O, Baliga NS, Wang JT, Ramage D, Amin N, Schwikowski B, Ideker T. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. 2003;13(11):2498–504.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  68. Chou CH, Shrestha S, Yang CD, Chang NW, Lin YL, Liao KW, Huang WC, Sun TH, Tu SJ, Lee WH, et al. miRTarBase update 2018: a resource for experimentally validated microRNA-target interactions. Nucleic Acids Res. 2018;46(D1):D296–d302.

    Article  CAS  PubMed  Google Scholar 

Download references

Acknowledgments

We thank Novel Bioinformatics Ltd., Co for the support of bioinformatics analysis and we also thank GENE DENOVO for the support of uploading the sequencing data.

Funding

This work is supported by the National Science and Technology Support Program in Rural Areas of the12th Five-Year Plan [Grant 2015BAD12B01]. These funding agencies provided great support in funding, and played no role in the design of the study and collection, analysis, and interpretation of data and in writing the manuscript.

Author information

Authors and Affiliations

  1. College of Veterinary Medicine, Northeast Agricultural University, Changjiang Road No. 600, Xiang Fang District, Harbin, People’s Republic of China

    Xuewei Huang, Junyan Zhang, Zengsu Liu, Meng Wang, Xiaolong Fan, Li Wang, Han Zhou, Yanping Jiang, Wen Cui, Xinyuan Qiao, Yigang Xu, Yijing Li & Lijie Tang

  2. Heilongjiang Key Laboratory for Animal Disease Control and Pharmaceutical Development, Harbin, People’s Republic of China

    Li Wang, Han Zhou, Yanping Jiang, Wen Cui, Xinyuan Qiao, Yigang Xu, Yijing Li & Lijie Tang

Authors
  1. Xuewei Huang
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. Junyan Zhang
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. Zengsu Liu
    View author publications

    You can also search for this author in PubMed Google Scholar

  4. Meng Wang
    View author publications

    You can also search for this author in PubMed Google Scholar

  5. Xiaolong Fan
    View author publications

    You can also search for this author in PubMed Google Scholar

  6. Li Wang
    View author publications

    You can also search for this author in PubMed Google Scholar

  7. Han Zhou
    View author publications

    You can also search for this author in PubMed Google Scholar

  8. Yanping Jiang
    View author publications

    You can also search for this author in PubMed Google Scholar

  9. Wen Cui
    View author publications

    You can also search for this author in PubMed Google Scholar

  10. Xinyuan Qiao
    View author publications

    You can also search for this author in PubMed Google Scholar

  11. Yigang Xu
    View author publications

    You can also search for this author in PubMed Google Scholar

  12. Yijing Li
    View author publications

    You can also search for this author in PubMed Google Scholar

  13. Lijie Tang
    View author publications

    You can also search for this author in PubMed Google Scholar

Contributions

XWH and JYZ designed and conducted research experiments. XWH, ZSL, MW and XLF prepared the bursa of Fabricius for RNA sequencing. XYQ, LW, and YPJ sorted out the sequencing results. WC and HZ performed qRT-PCR, LJT, YJL, XWH and YGX analyzed the data and wrote the paper. All authors read and approved the final manuscript.

Corresponding author

Correspondence to Lijie Tang.

Ethics declarations

Ethics approval and consent to participate

Animal experiments were carried out in accordance with the recommendations in the institutional and national guidelines for animal care and use. The protocol was approved by the Committee on the Ethics of Animal Experiments of Northeast Agricultural University, Harbin, China (2016NEFU-315, 13 April 2017).

Consent for publication

Not applicable.

Competing interests

The authors have no conflicts of interest to declare.

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.

Data quality of lncRNA and mRNA profiles.

Additional file 2: Table S2.

High quality clean reads compared with the ribosomal RNA.

Additional file 3.

The information of all known lncRNAs identified in this study.

Additional file 4.

The information of all novel lncRNAs identified in this study.

Additional file 5.

The information of all novel circRNAs identified in this study.

Additional file 6.

All the differentially expressed mRNAs in this study.

Additional file 7.

All the differentially expressed lncRNAs in this study.

Additional file 8.

All the differentially expressed circRNAs in this study.

Additional file 9.

GO enrichment analysis of mRNAs.

Additional file 10: Figure S1.

Gene ontology enrichment analysis for the antisense, cis, and trans roles of the differentially expressed lncRNAs in chicken BF between the two groups; a antisense; b cis; and c trans. The green, red, and blue columns indicate biological processes (BPs), cellular components (CCs), and molecular functions (MFs), respectively.

Additional file 11.

KEGG pathway analysis of mRNAs.

Additional file 12: Figure S2.

Kyoto Encyclopedia of Genes and Genomes pathway enrichment for the antisense, cis, and trans roles of the differentially expressed lncRNAs in chicken BF between the two groups; a antisense; b cis; and c trans. The vertical axis shows the pathways, and the horizontal axis indicates the Rich factor. The dot size indicates the number of differentially expressed genes in the pathway, and the coloration corresponds to the Q-value range.

Additional file 13.

LncRNA-mRNA co-expression relationships.

Additional file 14: Table S3.

The target genes of lncRNAs.

Additional file 15.

CircRNA/mRNA-miRNA co-expression relationships.

Additional file 16.

ARRIVE guidelines checklist.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. 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 in a credit line to the data.

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Cite this article

Huang, X., Zhang, J., Liu, Z. et al. Genome-wide analysis of differentially expressed mRNAs, lncRNAs, and circRNAs in chicken bursae of Fabricius during infection with very virulent infectious bursal disease virus. BMC Genomics 21, 724 (2020). https://doi.org/10.1186/s12864-020-07129-1

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12864-020-07129-1

Keywords

  • Infectious bursal disease virus
  • Bursa of Fabricius
  • mRNAs
  • lncRNAs
  • circRNAs

BMC Genomics

ISSN: 1471-2164

Contact us