Identification and analysis of long non-coding RNAs and mRNAs in chicken macrophages infected with avian infectious bronchitis coronavirus

Background Avian infectious bronchitis virus (IBV) is a gamma coronavirus that severely affects the poultry industry worldwide. Long non-coding RNAs (lncRNAs), a subset of non-coding RNAs with a length of more than 200 nucleotides, have been recently recognized as pivotal factors in the pathogenesis of viral infections. However, little is known about the function of lncRNAs in host cultured cells in response to IBV infection. Results We used next-generation high throughput sequencing to reveal the expression profiles of mRNAs and lncRNAs in IBV-infected HD11 cells. Compared with the uninfected cells, we identified 153 differentially expressed (DE) mRNAs (106 up-regulated mRNAs, 47 down-regulated mRNAs) and 181 DE lncRNAs (59 up-regulated lncRNAs, 122 down-regulated lncRNAs) in IBV-infected HD11 cells. Moreover, gene ontology (GO) and pathway enrichment analyses indicated that DE mRNAs and lncRNAs were mainly involved in cellular innate immunity, amino acid metabolism, and nucleic acid metabolism. In addition, 2640 novel chicken lncRNAs were identified, and a competing endogenous RNA (ceRNAs) network centered on gga-miR-30d and miR-146a-5p was established. Conclusions We identified expression profiles of mRNAs and lncRNAs during IBV infection that provided new insights into the pathogenesis of IBV. Supplementary Information The online version contains supplementary material available at 10.1186/s12864-020-07359-3.

recombination of the IBV genome between viruses with large genetic differences. Therefore, commercial vaccines often fail or only provide partial protection against IBVs [5], posing a major challenge to the poultry industry.
Long non-coding RNAs (lncRNAs) are transcripts longer than 200 nucleotides but without proteincoding capacity [6]. LncRNAs are critical regulators of a wide range of biological processes, including cell proliferation, differentiation, apoptosis, autophagy, tissue repair, and remodeling [7]. Several studies have demonstrated that lncRNAs function at the hostpathogen interface to regulate viral infections either by innate immune responses at several levels including activation of pathogen-recognition receptors or by epigenetic, transcriptional, and posttranscriptional effects [8]. For example, the latest research reported that lncRNA Malat1, a negative regulator of antiviral type I IFN (IFN-I) production, suppressed antiviral innate responses by targeting TDP43 activation via the RNA-RBP interactive network [9]. Although the nature of lncRNAs is well characterized in mammals, little is known about their functions in birds, especially in the field of antivirals for poultry [10]. Furthermore, except for loc107051710, lncRNA L11530, and lncRNA L09863 [11,12], little is known about lncRNA-mediated innate immune response in chicken. The functions of lncRNAs in anti-IBV immune response in chickens remain unclear.
LncRNAs regulate viral infections in multiple ways, such as epigenetic regulation and promotion of viral latency, protein scaffolding and nuclear localization, alternative splicing, and transcriptional regulation of mRNA via miRNA "sponges" [13]. One of the primary mechanisms of functioning of lncRNAs is by competing for shared microRNAs with mRNAs-known as the competing endogenous RNA (ceRNA) hypothesis [14]. We have previously analyzed the miRNAs of IBV-infected chicken kidney tissue and obtained two differentially expressed miRNAs, namely gga-miR-30d and miR-146a-5p, which are encoded by chicken chromosomes 1 and 13, respectively [15]. Next, an in vitro study demonstrated that gga-miR-30d inhibited IBV replication in HD11 cells by targeting USP47 [16], whereas miR-146a-5p promoted IBV replication in HD11 cells by targeting IRAK2 and TNFRSF18 [17].
We further analyzed the expression patterns of lncRNAs and mRNAs in IBV-infected HD11 cells using next-generation high throughput sequencing techniques. Differential expression and co-expression network analysis were conducted to identify interactions between mRNAs and lncRNAs to understand their possible functions in IBV infection. Moreover, a ceRNA network based on gga-miR-30d and miR-146a-5p was established. These results reveal a new data platform to conduct functional studies of chicken lncRNAs and provide valuable information on new therapeutic approaches to control IBV.

Replication status of IBV in HD11 cells
The expression of non-coding RNAs (ncRNAs) in cells is closely related to the stage of virus infection [13]. To confirm the status of IBV in HD11 cells, indirect fluorescent immunoassay (IFA) was performed after IBV infection for 0, 36, or 48 h (Fig. 1). The infection rate of the virus was calculated using the relative ratio of red fluorescence (IBV N protein) and blue fluorescence (cell nucleus) using ImageJ [18]. The results showed that the virus replication started a little after 36 h of infection; however, it replicated vigorously after 48 h of infection (manifested by rupture and collapse of cells, cell aggregation under bright-field microscopy, and significant red fluorescence). No red fluorescence was observed in mock-infected cells. Cytopathies affect the stability of nucleic acids in the cells. Therefore, cells after 36 h of infection were selected to extract total RNA and build the libraries.

RNA libraries establishment and lncRNA identification
Total RNA was extracted from 36 h post-infected HD11 cells (Exp 1, 2, and 3) and mock-infected cells (CK 1, 2, and 3). After high-throughput sequencing, six libraries with an average of 140,317,022 raw reads and 21,047, 553,300 bases were obtained. Nucleotides with a quality value above 30 (Q30) in reads were ranged from 94.09 to 94.4%. After data filtering and quality control, an average of 127,609,381 (90.95%) clean reads with 19,141, 407,250 high-quality bases were retained. Following the removal of rRNAs, clean reads were mapped to the chicken reference genome. The percentage mapping rates of six libraries ranged from 91.565 to 92.24% (Table 1).
The lncRNAs were identified as follows (Fig. 2). We used StringTie (version 1.2.4) software to assemble the transcripts based on the comparison results of HISAT2 (version 2.1.0). Transcripts with uncertain strand orientation were removed. The remaining assembled transcripts for lncRNAs were screened for transcripts with length ≥ 200 nucleotides and exon number ≥ 2 to obtain 59,930 transcripts. Transcripts whose class-code was x/ u/i were screened to obtain 4044 transcripts. Moreover, transcripts with cover > 3 in at least one sample were screened to obtain 4008 transcripts. We used PLEK (version 1.2), Coding-Non-Coding Index (CNCI; version 2.0), and PfamScan (version 1.6) to analyze the coding potential of candidate lncRNAs. All three software revealed that new transcripts without coding potential were high confidence lncRNAs. Ultimately, 2640 lncRNAs were identified.
To further identify the characteristics of lncRNAs, we compared predicted lncRNAs with mRNAs for transcript number, length, and exon number (Fig. 3). The result showed that the number of mRNA transcripts was higher than that of lncRNAs. With respect to the transcript length, the predicted lncRNAs were primarily concentrated between 200 and 3000 bp, whereas mRNAs were mainly of a length between 1400 and 5000 bp. In addition, the majority of lncRNAs contained two to three exons and very few  Q30: nucleotides with a quality value above 30 in reads Genome mapping rate: the percentage of reads mapped to the reference genome had more than 10 exons. However, the majority of mRNAs contained more than 10 exons. In summary, compared with mRNAs, lncRNAs had fewer and shorter transcripts, and fewer exons.

Expression profiles of lncRNAs and mRNAs in IBV-infected HD11 cells
In total, we obtained 15,358 mRNAs and 11,510 lncRNAs. Moreover, 153 mRNAs were differentially expressed, with 106 mRNAs significantly up-regulated and 47 mRNAs significantly down-regulated (Additional File 1). In addition, the expression of 181 lncRNAs changed significantly. Among these, 59 lncRNAs were up-regulated and 122 were downregulated (Additional File 2). Heat map and M-A map enrichment analyses (Fig. 4) revealed that compared with the control group, the expression profiles of lncRNAs and mRNAs changed significantly after 36 h of IBV infection.

LncRNA target gene prediction
LncRNAs regulate genes in a cis or trans manner. We enumerated the top 10 DE lncRNAs and their possible target genes by searching the gene-encoding protein within 100 kb upstream and downstream of lncRNAs. These genes were considered potential cis-regulated target genes corresponding to the lncRNAs. For example, SHISA6, located 19,659 bp downstream of MSTRG8180, was considered a potential target gene of MSTRG8180 (Additional File 3). To predict the trans-regulated target genes of lncRNAs, the top 10 DE lncRNAs and 50 most relevant mRNAs were selected to construct the co-expression network of lncRNA-mRNA pairs based on Pearson's correlation coefficient by Cytoscape (Additional File 4). As shown in Fig. 5, the network contained 174 edges. The majority of lncRNAs had multiple target genes, which were related to other lncRNAs, thus forming a large and complex co-expression network.

Pathway analysis of regulated lncRNAs and mRNAs after IBV infection in HD11
To further explore the functions of these DEGs and DE-lncRNAs following IBV infection, GO categorization and pathway analyses were performed. Significantly enriched GO terms (top 10 biological processes [BP], top 5 molecular functions [MFs], and top 5 cellular components [CCs]) and KEGG terms for mRNAs and lncRNAs are listed in Fig. 6. As shown in Fig. 6a, the GO categorization indicated that DEGs were mainly enriched in biological processes of cellular immunity in response to external stimuli. The top three enriched GO terms were related to defense response to other organisms (GO: 0098542), antimicrobial humoral response (GO: 0019730), and flavonoid metabolism (GO: 0009812). Further, GO-MF and GO-CC analyses identified cytosol and transcription factor binding, respectively. Furthermore, the KEGG pathway analysis revealed that the identified mRNAs mainly participated in protein processing in the AGE-RAGE signaling pathway in diabetic complications (ID: gga04933), cytokine-cytokine receptor interaction (ID: gga04060), and arachidonic acid metabolism (ID: gga00590). The top 20 pathways are shown in Fig. 6b.
The GO enrichment analysis of DE lncRNAs showed them to be mainly involved in the regulation of mRNA and RNA binding during IBV replication.  (Fig. 6c). The KEGG pathway analysis indicated that lncRNA target Fig. 2 LncRNA identification Process, "u" (unknown intergenic transcript), "i" (a transfrag falling entirely within a reference intron), and "x" (exonic overlap with reference on the opposite strand). The number on the right side of the picture represents the number of transcripts filtered out from each step genes mainly participated in alanine, aspartate, and glutamate metabolism (gga00250), amino sugar and nucleotide sugar metabolism (gga00520), and sphingolipid metabolism (gga00600). The enrichment pathways are listed in Fig. 6d.
We next analyzed the DE-lncRNAs and screened their target genes related to immunity. We constructed a network diagram using correlation coefficients (Fig. 7) that revealed a complex regulatory network between lncRNAs and immune genes. One lncRNA participated in the regulation of multiple genes in different ways, and one gene was regulated by multiple lncRNAs. For example, lncRNA MSTRG.14220.1 and MSTRG.21445.2 (Additional File 5) were related to at least 10 or 11 immune genes, and they were speculated to function in immune regulation in HD11 cells.

RT-qPCR validation
To validate the high-throughput sequencing results, we performed qPCR to detect the expression of lncRNAs and mRNAs in HD11 cells. Five lncRNAs and five mRNAs were selected randomly for qPCR to determine their relative expression ( Table 2). The results are shown in Fig. 9. The qPCR results indicated that the expression patterns of these lncRNAs and mRNAs were consistent with those by RNAsequencing.

Discussion
The fact that almost all isolated wild-type IBV strains cannot adapt to cell lines restricts their in-depth research [19]. In 2017, Han et al. reported that the IBV Beaudette strain could be serially passaged in HD11 cells    and triggered cytopathic effects [20]. HD11 is an avian myelocytomatosis virus (MC29) transformed chicken macrophage-like cell line [21]. Because macrophages function in both innate and adaptive immunity [22], it is imperative to study the response of ncRNAs in macrophages to viral stimulation. In the present study, we studied the mRNA-lncRNA regulatory network following IBV infection of HD11 cells.
The development of next-generation, high-throughput sequencing techniques has enabled us to obtain global RNA expression profiling. Recently, several studies have implicated genes and ncRNAs in IBV infection. Cong et al. [23] performed a transcriptome analysis in chicken kidney tissues infected with IBV strain ck/CH/LDL/ 091022. Their results showed that IL6, STAT1, MYD88, IRF1, and NFKB2 as key immune and inflammatory responses genes. We conducted in vitro studies to assess the changes in the expression of genes related to innate immunity in IBV-infected HD11 cells using highthroughput sequencing. The results revealed that certain inflammatory factors (IL15, IL1RAPL1, IL22, and IL8) and some antiviral genes such as IFIT5 and MX1 were differentially expressed following IBV infection. However, in addition to the above genes, the typical signaling pathways for transcriptional activation of chicken IFNs (such as Toll-like receptors), their downstream signaling pathways (such as the JAK/STAT pathway), and interferon-stimulated genes (ISGs) did not change significantly. IBV can use multiple strategies to counteract the IFN response. Certain studies report that IBV delays the activation of the IFN response and inhibits IFN-mediated phosphorylation and translocation of STAT1 [24][25][26]. We did not detect changes in typical innate immune signaling pathway factors at mRNA levels during the early replication stage of IBV, which could be one of the ways by which IBV escapes the host's innate immunity.
Screening of lncRNAs and their comparison between virus-and mock-infected cells revealed viral infectionassociated host lncRNAs. In this study, we analyzed the expression patterns of lncRNAs in IBV-infected and mock-infected HD11 cells using RNA sequencing. The   Fig. 8 The competitive endogenous RNA (ceRNAs) network centered on gga-miR-30d and miR-146a-5p. The relationship between miRNAs and target genes has been shown [16,17]. The expression of the lncRNAs in the picture had changed significantly in IBV-infected HD11 cells, and they were predicted to interact with miRNAs at the nucleic acid sequence level  3, whose target genes (TARDBP, FMR1, and TRIM8) were mainly enriched in terms related to nucleic acid and protein metabolism, and not to terms related to immunity (although certain virus-related terms were also been enriched). The most likely reason could be that nucleic acids and proteins required for virus assembly in the cell are synthesized in large quantities during infection. Protein synthesis in cells is regulated by a variety of factors, including host lncRNAs. We speculated that IBV uses the host's lncRNAs to produce viral nucleic acids and proteins required for replication. Moreover, early studies have demonstrated that the virus hijacks host-encoded lncRNAs to establish persistent infections [29]. Vesicular stomatitis virus, vaccinia virus, and herpes simplex virus 1 could hijack the lncRNA-ACOD1 complex to promote viral replication by activating the metabolic enzyme glutamicoxaloacetic transaminase 2 [30,31]. Furthermore, lncRNA NeST, Lethe, lncRNA-CMPK2, VIN, and NRON are used by the virus to facilitate infection or susceptibility [32].
The ceRNA hypothesis proposes that transcripts with shared microRNA (miRNA)-binding sites compete for post-transcriptional control. Recently, lncRNA HAND2-AS1 was found to work as the ceRNA of miR-3118 to suppress the proliferation and migration in breast cancer by upregulating phlpp2 [33]. LncRNA-Six1 is a target of miR-1611 that functions as a ceRNA to regulate the expression of Six1 protein and fiber type switching in chicken myogenesis [34]. LncRNA MARL functioned as a ceRNA for miR-122 to control protein abundance of MAVS, thereby inhibiting the SCRV replication and promoting antiviral responses [35]. Previously, we demonstrated that gga-miR-30d inhibited IBV replication in HD11 cells by targeting USP47 [16]. MiR-146a-5p promoted IBV replication in HD11 cells by targeting IRAK2 and TNFRSF18 [17]. We further analyzed the changes in lncRNAs in IBV-infected HD11 cells. Based on gga-miR-30d and miR-146a-5p, we constructed the lncRNA-microRNA-mRNA ceRNA network that provided insights into the interplays in lncRNAs, microRNAs, and mRNAs. For example, lncRNA MSTRG.21445.2, identified as a lncRNA with a length of about 4000 bp in our study, was significantly up-regulated following IBV infection. It is located downstream of the USP47 and overlaps with the pellino E3 ubiquitin protein ligase family member 2 (PELI2) gene on chicken chromosome 5. The function of PELI2 in innate immunity has already been revealed [36,37]. LncRNA MSTRG.21445.2 regulates IBV infection by competing with mRNAs of USP47, IRAK2, and TNFRSF18 for gga-miR-30d and miR-146a-5p. Moreover, it may participate in natural immunity against IBV infection by acting on PELI2 alone. Although these results provided new insights for studying the mechanism of IBV infection, bioinformaticsbased predictions should be confirmed using wetlaboratory studies.

Conclusion
We performed a comprehensive analysis of lncRNA and mRNA expression profiles in HD11 cells following IBV infection using RNA sequencing. We identified more than 2000 novel chicken lncRNAs. In addition, we identified several DE genes that may play pivotal roles in IBV infection using co-expression network analysis. We screened several DE lncRNAs related to IBV replication and cellular immunity, e.g., Further, we established a ceRNA network centered on gga-miR-30d and miR-146a-5p. We believe our results provide valuable insights into the mechanisms of IBV pathology.

Cell culture and virus infection
The HD11 cell line was kindly provided by Prof. Xin-An Jiao, Yangzhou University. The cells were cultured in Dulbecco's modified Eagle's medium (DMEM, Gemini, USA) supplemented with 10% fetal bovine serum (FBS, Gemini, USA), 100 IU/mL penicillin, and 100 μg/mL streptomycin sulfate (Solabio, USA). The IBV Beaudette strain (GenBank: DQ001339) was kindly gifted by Prof. Ding-Xiang Liu, Nanyang Technological University. The cells were infected with the IBV Beaudette strain (multiplicity of infection [MOI] = 2) and incubated in 5% CO 2 at 37°C for 1 h. Following the incubation, the cells were rinsed with phosphate-buffered saline (PBS, Gemini) and cultured at 37°C in DMEM supplemented with 2% FBS (5% CO 2 ). Controls included cells that were mock infected.

RNA isolation and library construction
HD11 cells were divided into two groups with three biological replicates in each group: cells infected with IBV represented the experimental group (Exp 1, 2, and 3); mock-infected HD11 cells served as control (CK 1, 2, 3). Following the infection for 36 h, total RNA was prepared from these treatments by the TRIzol method (Invitrogen, USA) according to the manufacturer's protocol. Total RNA quality was assessed on an Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA). Sequencing libraries were generated using the rRNAdepleted RNA with NEBNext® Ultra™ Directional RNA Library Prep Kit for Illumina® (NEB, USA) following the manufacturer's recommendations.

Transcript sequencing and genome mapping
The libraries were paired-end sequenced (PE150, sequencing reads were 150 bp) at Personal Biotechnology Co., Ltd. (Shanghai, China) using the Illumina HiSeq 2000. The information on the quality of raw data in FASTQ format was calculated, following which the raw data was filtered using Cutadapt (v2.7) software. The clean data were obtained by removing the reads containing the adapter, reads containing poly-N, and lowquality reads. All subsequent analyses were based on high-quality clean data. The chicken genome GRCg6a (GenBank assembly accession: GCA_000002315.5) was used as the reference genome. The reference genome index was constructed by Bowtie2 (v2.4.1) and the filtered reads were mapped onto the reference genome using HISAT2 (2.1.0), the default mismatch was no more than two. The alignment region distribution of mapped reads was calculated.
Coding potential and conserved analysis of lncRNAs First, the raw data were filtered to obtain high-quality clean data. The HISAT2 (2.1.0) software was used to compare the filtered reads with the chicken genome as the reference. Based on the results of HISAT2 (2.1.0), StringTie (v1.2.4) software was used to assemble the transcripts. Transcripts with uncertain directions were removed and lncRNAs were screened in the remaining assembled transcript set using the following steps: (1) Transcripts with length ≥ 200 bp and the number of exons ≥2 were selected. (2) Only transcripts class-coded as "u" (unknown intergenic transcript), "i" (a transfrag falling entirely within a reference intron), and "x" (exonic overlap with reference on the opposite strand) were selected as candidate lncRNAs. (3) LncRNAs with cover > 3 in at least one sample, that is, those appearing at least thrice in the transcript of one sample, were screened. (4) To further obtain high-confidence lncRNAs, the coding potential of the candidate lncRNA was analyzed using PLEK (1.2), CNCI (v2) and PfamScan (1.6). (5) The mRNAs were compared with predicted lncRNAs in terms of transcript numbers, length, and exon numbers to describe the characteristics of lncRNAs.

LncRNA target gene prediction
The target gene prediction of lncRNAs was divided into cis-target gene prediction and trans-target gene prediction. We screened for genes with protein-coding capabilities within 100 kb upstream and downstream of lncRNAs for cis-regulated target genes corresponding to lncRNAs. The target genes with trans-regulation were screened by calculating the expression correlation (Pearson's correlation test) or co-expression analysis method between lncRNAs and mRNAs. We used the correlation coefficient > 0.95 and p-value < 0.05 to screen the transacting relationship between lncRNAs and mRNAs. Next, the lncRNA-mRNA co-expression network was constructed using Cytoscape.

Identification and bioinformatics analyses of differentially expressed mRNAs and lncRNAs
StringTie (version 1.2.4) software was used to perform counts at the transcript level to obtain the original expression. Next, FPKM (fragments per kilobase of exon model per million fragments mapped) was used to normalize the expression. In general, FPKM > 1 demonstrated expression. DESeq was used to analyze the DEGs and DE lncRNAs in IBV-infected and control cells.