Determination of antiviral action of long non-coding RNA loc107051710 during infectious bursal disease virus infection due to enhancement of interferon production

ABSTRACT The functions and profiles of lncRNAs during infectious bursal disease virus (IBDV) infection have not been determined, yet. The objectives of this study were to determine the antiviral action of loc107051710 lncRNA during IBDV infection by investigating the relationship between loc107051710 and IRF8, Type I IFN, STATs, and ISGs. DF-1 cells were either left untreated as non-infected controls (n = 1) or infected with IBDV (n = 3). RNA sequencing was applied for analysis of mRNAs and lncRNAs expression. Differentially expressed genes were verified by RT-qPCR. Then identification, of 230 significantly different expressed genes (182 mRNAs and 48 lncRNA) by pairwise comparison of the infected and control groups, was carried out. The functions of differentially expressed lncRNAs were investigated by selection of lncRNAs and mRNAs significantly enriched in the aforementioned biological processes and signaling pathways for construction of lncRNA-mRNA co-expression networks. The techniques of gene ontology and Kyoto Encyclopedia of Genes and Genomes pathways were applied. It was suggested that these differentially expressed genes were involved in the interaction between the host and IBDV. Loc107051710 was found to have potential antiviral effects. RT-qPCR and western blot were applied and revealed that loc107051710 was required for induction of IRF8, type I IFN, STAT, and ISG expression, and its knockdown promoted IBDV replication. By fluorescence in situ hybridization, it was found that loc107051710 was translocated from the nucleus to the cytoplasm after infection with IBDV. Overall, loc107051710 promoted the production of IFN-α and IFN-β by regulating IRF8, thereby promoting the antiviral activity of ISGs.


Introduction
Infectious bursa disease (IBD) is a highly contagious and immunosuppressive disease that affects young chickens. It causes high mortality rates and large economic losses to the poultry industry [1,2]. Infectious bursal disease virus (IBDV) mainly causes the destruction of B-lymphocytes in the bursa of Fabricius, which can lead to severe immunosuppression and secondary infections in infected chickens [3][4][5]. The virus (IBDV) is non-enveloped with an icosahedral capsid. It belongs to the family Birnaviridae. The genome of IBDV is bi-segmented, double-stranded RNA [6].
While the importance of protein encoding genes has long been known, the various types of noncoding RNAs are increasingly receiving attention [7]. Long noncoding RNAs (lncRNAs), which are composed of more than 200 nucleotides, in particular, are involved in apoptosis during organ development and tissue differentiation [8]. They also play regulatory roles in many diseases, particularly, human cancers [9][10][11], and in the adaptive and innate immune responses [8,12]. Moreover, the lncRNAs have antiviral functions by changing the expression levels of IFNs or ISGs against Theiler's virus, encephalomyocarditis virus, hepatitis C virus, influenza A virus, and HIV [13][14][15][16][17].
Although DF-1 cells, a continuous line of chicken embryo fibroblasts, have been used to study the transcriptional changes and differentially expressed proteins of host cells in response to IBDV infections [18][19][20][21][22][23][24], the role of host lncRNAs in response to IBDV infection is still unknown. Application, of some sophisticated techniques such as techniques of in depth bioinformatics, RNA sequencing, RT-qPCR, Fluorescence in situ hybridization, and western blotting, is expected to provide a theoretical basis for novel preventive and therapeutic strategies against CONTACT Guangpeng Ma maguangpeng1977@163.com; Lijie Tang tanglijie@163.com *Xuewei Huang, Yigang Xu and Qingyu Lin contributed equally to this work.
Supplemental data for this article can be accessed here.
IBDV. The objective of this study was to determine the antiviral action of loc107051710 during IBDV infection in cultured DF-1 cells after investigation of the relationship between the loc107051710 and IRF8, type I IFNs, STATs and ISGs. The antiviral action of loc107051710 was determined by application of RNA-seq for analysis of mRNAs and lncRNAs expression. Next, identification of differentially expressed genes by pairwise was carried out. After that, the functions of differentially expressed lncRNAs were investigated by selection of lncRNAs and mRNAs significantly enriched in the aforementioned biological processes and signaling pathways for construction of lncRNA-mRNA co-expression networks. The techniques of GO and KEGG pathways were applied. Subsequently, RT-qPCR and western blot were applied to investigate the relationship between loc107051710 and IRF8, type I IFNs, STATs, and ISGs expression. Eventually, fluorescence in situ hybridization was applied to determine the location of loc107051710.

Experimental design and RNA extraction for sequencing
The DF-1 cell cultures were divided into two groups. The first group was infected with IBDV (n = 3). The second group was kept without infection as a control (n = 1). The total RNA of the control (n = 1) and IBDV-infected (n = 3) groups was extracted using TRIzol reagent (Invitrogen, Carlsbad, CA, USA) according to the manufacturer's instructions. The purity, concentration, and integrity of the total RNA were analyzed using an Agilent Bioanalyzer 2100 (Agilent Technologies, Santa Clara, CA, USA). Samples with an RNA integrity number (RIN) value > 9.7 and an optical density 260:280 ratio > 2.0 were used for library construction and deep sequencing. RNA-seq was performed by Novel Bioinformatics Co., Ltd.

Quality control of RNA readings
The raw readings were preprocessed with custom Perl scripts to ensure the high quality of RNA for subsequent analyzes. Adapter-polluted readings, low-quality readings (in which more than 15% of bases had a Phred score ≤ 19, or more than 5% were undetermined, N), and readings matching rRNA were removed.
Alignment, mapping of RNA-seq to reference genome and transcriptome assembly Using TopHat version 2.0.12 [25] the RNA-seq readings were aligned and mapped to the reference genome obtained from NCBI Gallus_gallus-5.0 (ftp://ftp.ncbi. nlm.nih.gov/genomes). Using the Cufflinks 2.2.1 program [26], the transcriptome of each sample was assembled independently with the help of the reference annotation-based transcript assembly technique.

Selection of lncRNAs
Based on the assembly results, some transcripts were removed such as transcripts with readings of 0. For selection of lncRNAs, the following transcripts were removed: Transcripts which were shorter than 200 nucleotides and having less than two exons, transcripts which were only present in one sample, and transcripts which were encoding a protein family or were a known mRNA transcript. The Coding-Non-Coding Index (CNCI) [27], Coding Potential Calculator (CPC) [28], and Coding-Potential Assessment Tool (CPAT) [29] were used to evaluate the coding potential of the transcripts. Default parameters were used with all the software. After filtering, transcripts without coding potential were considered candidate lncRNAs.

Quantification of gene expression
The read count for each gene in each sample was determined by HTSeq v0.6.0 [30], and the number of fragments per kilobase of transcript sequence per million base pairs sequenced [31] was calculated to measure the expression levels of both mRNAs and lncRNAs in each sample. DESeq v1. 16.0 [32] was used for the differential expression analysis of two groups. Genes, with a p value <0.05 and a fold change >1.5, were considered differentially expressed in the two groups.
Determination of significantly enriched biological functions and pathways in mRNAs differentially expressed in IBDV infected DF-1 cells and control DF-1 cells Using fisher.test and p.adjust routines, 2 analyses were performed with custom R scripts. These 2 analyses were gene ontology (GO) and Kyoto encyclopedia of genes and genomes (KEGG) pathway enrichment. These analyses were used to determine biological functions and pathways significantly enriched in mRNAs differentially expressed in IBDV infected DF-1 cells and control DF-1 cells. GO categories and KEGG pathways with p values <0.05 were considered significantly enriched.
Protein-protein and lncRNA-mRNA co-expression network analyses We constructed gene co-expression networks to identify the interactions among differentially expressed genes. The gene co-expression networks were built using the normalized signal intensity of specific expression genes. For each pair of genes, the Pearson correlation coefficient was calculated, and significantly correlated pairs were selected to construct the networks. In network analysis, degree centrality* was the simplest and most important measure of the relative importance of a gene within a network. Moreover, to analyze certain properties of the networks, k-cores**, from graph theory, were introduced for simplifying graph topology analysis. In the present study, the purpose of the network structure analysis was to locate genes in one network. When analyzing the different networks, genes with the largest degree of difference between the two classes were selected. (Note: *Degree centrality is defined as the number of links connecting one node to other nodes; **A k-core of a network is a subnetwork in which all nodes are connected to at least k other genes in the subnetwork. Accordingly, a k-core of a protein-protein interaction network contains cohesive groups of proteins).

Gene silencing of loc107051710
A specific siRNA for loc107051710 was designed by Genepharma Co., Ltd. (Shanghai, China). The sequences of the specific and control siRNA are listed in Table 1. In total, 5 × 10 5 DF-1 cells were seeded in 6-well plates, and when they reached 50-60% confluence, they were transfected with 100 nM negative siCont or siloc107051710 using Lipofectamine 2000 (Invitrogen, Carlsbad, CA, USA) according to the manufacturer's instructions. After transfection for 24 h, the cells were infected for 24 h with IBDV (MOI = 1).

Quantitative reverse transcription PCR (Rt-qPCR)
Total RNA was extracted from control and IBDVinfected groups using TRIzol according to procedures previously described and reverse-transcribed to cDNA using M-MLV Reverse Transcriptase (Promega, Madison, WI, USA) according to the manufacturer's procedures. Twenty-four differentially expressed genes were selected to confirm the accuracy of RNA-seq by RT-qPCR, including 12 lncRNAs and 12 mRNAs. Moreover, RT-qPCR was applied for detection of the expression of type I IFNs, STATs and ISGs. Primers were designed using oligo6 software (Tables 2-4). For RT-qPCR, the ABI 7500 Real-Time PCR system (Applied Biosystems, Carlsbad, CA, USA) and SYBR Green Master (TaKaRa, Dalian, China) were used; the reactions were performed in a final volume of 20 μL.
The GAPDH gene was used as an internal control. The relative fold change was calculated by the 2 −ΔΔCt method [33]. Experiments were repeated three times.

Western blotting
Total proteins were extracted, subjected to 15% SDS-PAGE and transferred to PVDF membranes. After blocking with  Table 2. RT-qPCR primers used for verification of mRNA results.

Fluorescence in situ hybridization (FISH)
Using FISH technique, loc107051710 was detected on DF-1 cells. The biotin-labeled antisense and sense probes of loc107051710 were synthesized by RiboBio (Guangzhou, China) and the sequences are listed in Table 1. The FISH assay was performed using a Fluorescent In Situ Hybridization Kit (RiboBio, Guangzhou, China) according to the manufacturer's procedures. In brief, DF-1 cells were fixed with 4% polyoxymethylene. After washing with PBS, the cells were treated with 0.5% Triton X-100 at 4°C for 5 min. Next, the cells were treated with pre-hybridization lncRNA FISH probe mix at 37°C for 30 min. Hybridization was performed by the adding antisense and sense loc107051710 FISH probe mix (10 mM ATP, 10 mM CTP, 10 mM GTP, 6.5 mM UTP and 3.5 mM Biotin-16-UTP), followed by incubation at 37°C overnight. After washing with 4×, 2×, and 1× SSC, the cells were stained with DAPI, washed three times with PBS, and visualized by laser scanning confocal microscopy.

Statistical analyses
The data were analyzed using GraphPad Prism5 software. One-way ANOVA was used to determine significance. P values <0.05 were considered statistically significant.

Overview of lncRNA and mRNA expression profiles
In total, 361 100 380 raw reads were obtained from the infected and control DF-1 cells (NCBI SRA Run Selector, accession number SRP145165). After quality control, 308 057 830 clean reads (readings after removal of Adapter-polluted readings, low-quality readings and readings matching rRNA), were obtained, and the percentage of clean reads ranged from 83.74% to 87.27%. The percentage of clean reads with a Phred quality value of more than 30 ranged from 95.03% to 95.54% (Table 5). Overall, 86-87% of the clean reads aligned with the chicken reference genome ( Table 6).

Identification of differentially expressed lncRNAs and mRNAs
To analyze the differential expression of mRNAs and lncRNAs between the control and IBDV-infected groups, heat maps (Figure 1(a)) and volcano plots (Figure 1(b)) were generated. Given the criteria of p value <0.05 and fold change >1.5, 91 mRNAs were upregulated and 91 were downregulated (Table 7, supplementary materials 1) . Furthermore, 13 lncRNAs were upregulated, and 35 lncRNAs were downregulated in IBDV-infected DF-1 cells (Table 7, supplementary materials 2). Table 3. RT-qPCR primers used for verification of lncRNA results.

Functional annotation of differentially expressed mRNAs
In total, 182 differentially expressed mRNAs were subjected to GO and KEGG pathway enrichment analyses (Figures 2 and 3) to obtain a better understanding of the potential roles of host factors in IBDV infection. We found that the GO terms were significantly enriched in differentially expressed mRNAs, as well as for just upregulated mRNAs, were mainly involved in biological processes and pathways that may be related to the replication of IBDV, including the type I IFN signaling pathway, defense response to virus, IFN-γmediated signaling pathway, and positive regulation of the type I IFN production signaling pathway (Figure 2  and supplementary materials 3). Moreover, KEGG analysis showed that differentially expressed genes were involved in the RIG-I-like receptor, NF-κB, Toll-like receptor, influenza A, and hepatitis C signaling pathways ( Figure 3, supplementary materials 4).

Protein-protein interaction network analysis
mRNA-mRNA functional interaction network analysis demonstrated that IFIT5, OASL, USP18, TLR3, IRF1, IFI6, SAMD9L, EPSTI1, and several other mRNAs were at the core of the network (Figure 4). Most of these mRNAs were upregulated and significantly enriched in the type I IFN signaling pathway and the RIG-I, NF-κB, and Toll-like receptor signaling pathways.

lncRNA-mRNA co-expression network analysis
To further examine the functions of the differentially expressed lncRNAs, lncRNA-mRNA co-expression networks were constructed. We first identified and selected the differentially expressed mRNAs and lncRNAs significantly enriched in GO terms and pathways related to the replication of IBDV. We then constructed and analyzed 10 lncRNA-mRNA co-expression networks ( Figure 5). The lncRNA-mRNA gene pair loc107051710-IRF8 had the best correlation.

Validation of RNA-seq data by Rt-qPCR
To confirm the accuracy of RNA-seq, several differentially expressed mRNAs and lncRNAs were selected and analyzed by RT-qPCR. The RT-qPCR results were consistent with those from RNA-seq (Figure 6(a,b)), indicating the validity of the RNA-seq data.

Effect of silencing lncRNA loc107051710 on IBDV replication
To further illustrate the association of the loc107051710 with IRF8, the lncRNA was knocked-down in DF-1 cells using a specific siRNA. The siRNA for loc107051710 was effective in reducing the expression of the lncRNA (Figure 7(a)). Moreover, the decrease in loc107051710 expression was accompanied by a decrease in IRF8 expression (Figure 7(b)). To determine whether the reduction in loc107051710 expression resulted in the reduction of type I IFNs (IFN-α and IFN-β), STATs (STAT1 and STAT2), and ISGs (OAS, Mx1, and PKR) expression, their expression levels were measured by RT-qPCR. Surprisingly, silencing loc107051710 decreased the levels of type I IFNs, STATs, and ISGs (Figures 8 and 9). To verify the antiviral effect of loc107051710, we infected siControl or siloc107051710-transduced DF-1 cells and measured the expression of IBDV. The results revealed that IBDV replication was increased by loc107051710 knockdown (Figure 10).

Determination of location of loc107051710 in DF-1 cells by application of fluorescence in situ hybridization
Using FISH, loc107051710 was detected in both the cytoplasm and nucleus of uninfected DF-1 cells (Figure 11(ac)). However, it was mainly detected in the cytoplasm after infection with IBDV ( Figure 11 d-f), and the results showed that Loc107051710 did not bind to the sense probe in IBDV-infected DF-1 cells (g-i).

Discussion
IRFs are a family of transcription factors that were originally isolated as positive and negative regulators of IFNs and IFN-responsive genes [34]. More recently, extensive studies of the IRF family have confirmed that IRFs have gained great attention due to their diverse roles such as in the initiation of antiviral responses, regulation of inflammatory cytokine expression, and control of cell cycle and apoptosis [35,36]. Numerous pattern recognition receptor that can detect cytosolic nucleic acids have been identified and most of them can evoke type I IFN responses via activation of IRFs [37]. The biological activities of type I IFN are initiated by binding to the type I IFN receptor. This leads to the activation of the receptor associated tyrosine kinases JAK1 and Tyk2, which phosphorylate STAT1 and STAT2. Phosphorylated STAT1 and STAT2 interact strongly with DNA-binding protein IRF-9. The IRF-9-STAT1-STAT2 heterotrimer is called ISGF3, which participates in the induction of the expression of ISGs, including OAS, Mx1 and PKP [38]. These ISGs can directly inhibit the replication, assembly and proliferation of the virus. To date, eight IRF homologues have been identified in chickens, though their functions are not yet well defined [39]. IRF8 (also known as IFN consensus sequence-binding protein [ICSBP]), a member of the IRF family, is expressed by B cells, dendritic cells (DCs) and macrophages, and has a role in the regulation of innate immune responses [40,41].
As an important transcription mediator, IRF8 is also required for type I IFN induction in virus-stimulated DCs [42]. It appears that IRF8 is involved in the transcriptional regulation of type I IFN genes; IRF8 can bind to promoters of IFN-α/β genes to promote the transcription of IFN-α/β [43]. Accordingly, it plays a critical role in the initiation of host innate immune responses to viral infection. It is well known that host cells can regulate the expression levels of various mRNAs during viral infection. In the current study, differentially expressed  Previous studies have demonstrated that IBDV infections induce the expression of different cytokine genes [18]. Furthermore, the genes of cytokines that initiate inflammatory responses, including the IL8, nitric oxide synthase, and cyclooxygenase-2 genes, were reported to be upregulated [44,45]. In the present study, we similarly found that the expression of STAT1, IL8L1, and IL8L2 was increased in DF-1 cells after infection with IBDV. Regulatory factors that could modulate the expression levels of these cytokine genes were also over-expressed, including IFN-induced GTP-binding Figure 5. LncRNA-mRNA co-expression networks. A dot represents a differentially expressed mRNA, and a triangle represents a differentially expressed lncRNA. Genes in the same k-core are indicated with the same color. The progression from gray to red represents increasing k for the k-core. The color scale bar is shown.  protein Mx1 and IFN-induced protein with tetratricopeptide (IFIT). IFIT5 is a member of the IFIT family, which can be triggered by viral infection [46] and enhances the antiviral response by promoting IRF3and NF-κB-mediated gene expression [47]. Our results showed that IFIT5 expression was 23.8-fold higher and that expression of the NF-κB-mediated genes IL8L1 and IL8L2 were 1.58-fold and 1.52-fold higher, respectively, in the infected group than in the control groups.
Radical S-adenosyl methionine (RSAD2) is an endoplasmic reticulum-associated virus inhibitory protein that can be induced by double-stranded RNA viral infection [48], as well as in an IFN-independent pathway by IRF1 [49]. In our study, IRF1 and RSAD2 expression were 3.01-fold and 58.4-fold higher, respectively. Although the specific antiviral mechanism of RSAD2 has not yet been determined, it is thought to indirectly inhibit viral replication by regulating cell survival [50]. Moreover, the upregulation of the tripartite motif-containing 25 (TRIM25) protein at the beginning of the infection may be associated with RIG-I in the detection of viral RNA intermediates [51]. A previous study has demonstrated that RIG-I can upregulate type I IFNs and reduce viral gene expression [52]. Our results showed that TRIM25 expression was increased 1.6-fold, and the RIG-I-like receptor signaling pathway was significantly enriched during infection. The differential expression patterns of these   mRNAs were consistent with those reported in previous studies. We believe that these over-expressed mRNAs may enhance the antiviral ability of DF-1 cells during the early stage of IBDV infection.
In this study, the type I IFN signaling pathway showed the highest enrichment. Type I IFNs can up-regulate the expression of hundreds of ISGs by activating the wellcharacterized JAK-STAT pathway [53]. Several studies have shown that ISGs induced by type I IFNs contribute to the antiviral response [54,55]. To determine whether differentially expressed lncRNAs play an antiviral role in IBDV-infected DF-1 cells by regulating interferon expression, we first screened out differentially expressed mRNAs significantly enriched in the GO terms and KEGG signaling pathways with antiviral activity. Next, we selected all the differentially expressed lncRNAs and mRNAs enriched in the GO terms and KEGG pathways to construct 10 lncRNA-mRNA co-expression networks. While there were many regulatory relationships in the networks, loc107051710, loc107052218, loc107052243, and loc107052781 had a close relationship with the IRF8 gene, of which loc107051710 and IRF8 gene pair had the closest relationship.
We confirmed the cascade relationship between loc107051710 and IRF8 and found that loc107051710 positively regulates the expression of IRF8. In addition, we found that loc107051710 acts as a positive transcriptional regulator of the antiviral-related IFN-α, IFN-β, STAT1, STAT2, OAS, Mx1 and PKR. Therefore, we speculate that loc107051710 promotes the production of IFN-α and IFN-β by regulating IRF8, thereby promoting ISGs antiviral activity. However, whether loc107051710 directly regulates IRF8 expression remains unknown.
The loc107051710 seems to possess antiviral function. However, it was unclear where loc107051710 was localized in uninfected cells and where it exerted antiviral effects after DF-1 cells were infected with IBDV. FISH showed that the amount of loc107051710 increased and shifted from the nucleus to cytoplasm during infection, indicating that it plays an antiviral role not only at the transcriptional level, but also at the post-transcriptional level. Therefore, the antiviral potential of the loc107051710 may be greater than the other lncRNAs analyzed in this study.
In conclusion, we profiled lncRNAs and mRNAs from IBDV-infected DF-1 cells, and based on GO and KEGG pathway enrichment analyses, we identified the differentially expressed mRNAs with antiviral functions. Notably, lncRNA-mRNA co-expression analysis revealed the potential antiviral function of the loc107051710, which involved the regulation of type I IFN, STAT, and ISG production to prevent IBDV infection. In brief, it was concluded that downregulation of the long non-coding RNA loc107051710 enhanced the replication of infectious bursal disease virus by reducing interferon production. Vice versa up-regulation of the long noncoding RNA loc107051710 suppressed the replication of infectious bursal disease virus by elevated interferon production. Our findings provide a valuable foundation for future studies of the molecular mechanisms of innate anti-IBDV responses in DF-1 cells.