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)
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).
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.
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.
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 . 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 . 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 , including long non-coding RNAs (lncRNAs)  and circular RNAs (circRNAs) . 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 . 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  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.
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)  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).
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.
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).
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).
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).
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.
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.
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.
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 . 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 . 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 . 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.
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.
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 . 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)  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) ; Cufflinks (v. 2.1.1)  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)  and CPC  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) , phyloFit , and PhastCons  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  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) . The binding sites of the miRNAs on the circRNAs and mRNAs were predicted using mirTarBase software , 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.
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.
Infectious bursal disease virus
Bursa of Fabricius
Long non-coding RNAs
The National Centre for the Replacement, Refinement, and Reduction of Animals in Research
Differentially expressed genes
Kyoto Encyclopedia of Genes and Genomes
Regulatory factor 8
Competing endogenous RNA
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.
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.
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.
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.
Becht H, Muller H. Infectious bursal disease--B cell dependent immunodeficiency syndrome in chickens. Behring Inst Mitt. 1991;89:217–25.
Berg TP. Acute infectious bursal disease in poultry: a review. Avian Pathol. 2000;29(3):175–94.
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.
Bhan A, Mandal SS. Long noncoding RNAs: emerging stars in gene regulation, epigenetics and human disease. ChemMedChem. 2014;9(9):1932–56.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
Lukiw WJ. Circular RNA (circRNA) in Alzheimer’s disease (AD). Front Genet. 2013;4:307.
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.
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.
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.
Seeger C, Mason WS. Molecular biology of hepatitis B virus infection. Virology. 2015;479-480:672–86.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
Heward JA, Lindsay MA. Long non-coding RNAs in the regulation of the immune response. Trends Immunol. 2014;35(9):408–19.
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.
Liu W, Ding C. Roles of LncRNAs in viral infections. Front Cell Infect Microbiol. 2017;7:205.
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.
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.
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.
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.
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.
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.
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.
Chen LL, Yang L. Regulation of circRNA biogenesis. RNA Biol. 2015;12(4):381–8.
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.
Guo JU, Agarwal V, Guo H, Bartel DP. Expanded identification and characterization of mammalian circular RNAs. Genome Biol. 2014;15(7):409.
Salzman J, Chen RE, Olsen MN, Wang PL, Brown PO. Cell-type specific features of circular RNA expression. PLoS Genet. 2013;9(9):e1003777.
Thomson DW, Dinger ME. Endogenous microRNA sponges: evidence and controversy. Nat Rev Genet. 2016;17(5):272–83.
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.
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.
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.
Langmead B, Salzberg SL. Fast gapped-read alignment with bowtie 2. Nat Methods. 2012;9(4):357–9.
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.
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.
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.
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.
Darling AE, Jospin G, Lowe E, Matsen FA, Bik HM, Eisen JA. PhyloSift: phylogenetic analysis of genomes and metagenomes. PeerJ. 2014;2:e243.
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.
Tafer H, Hofacker IL. RNAplex: a fast tool for RNA-RNA interaction search. Bioinformatics. 2008;24(22):2657–63.
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.
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.
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.
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.
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
The authors have no conflicts of interest to declare.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Data quality of lncRNA and mRNA profiles.
High quality clean reads compared with the ribosomal RNA.
The information of all known lncRNAs identified in this study.
The information of all novel lncRNAs identified in this study.
The information of all novel circRNAs identified in this study.
All the differentially expressed mRNAs in this study.
All the differentially expressed lncRNAs in this study.
All the differentially expressed circRNAs in this study.
GO enrichment analysis of mRNAs.
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.
KEGG pathway analysis of mRNAs.
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.
LncRNA-mRNA co-expression relationships.
The target genes of lncRNAs.
CircRNA/mRNA-miRNA co-expression relationships.
ARRIVE guidelines checklist.
About this article
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
- Infectious bursal disease virus
- Bursa of Fabricius