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).
Table1. Overview of the RNA sequencing data
|
Sample
|
Raw Data Reads No.
|
Bases (bp)
|
Q30 (bp)
|
Clean Reads No
|
Clean Bases (bp)
|
Clean Reads %
|
Genome Mapping Rate
|
Sequencing Mode
|
C-1
|
150897246
|
22634586900
|
21331618679(94.24%)
|
137417676
|
20612651400
|
91.06
|
125822533 (91.56%)
|
Paired-end,2×150bp
|
C-2
|
137054082
|
20558112300
|
19407891463(94.4%)
|
126612888
|
18991933200
|
92.38
|
116736730 (92.20%)
|
Paired-end,2×150bp
|
C-3
|
138645590
|
20796838500
|
19632609578(94.4%)
|
124766544
|
18714981600
|
89.98
|
114421545 (91.71%)
|
Paired-end,2×150bp
|
E-1
|
129518036
|
19427705400
|
18289744390(94.14%)
|
117058784
|
17558817600
|
90.38
|
107868528 (92.15%)
|
Paired-end,2×150bp
|
E-2
|
154994906
|
23249235900
|
21875735326(94.09%)
|
139996898
|
20999534700
|
90.32
|
129086227 (92.21%)
|
Paired-end,2×150bp
|
E-3
|
130792272
|
19618840800
|
18520223188(94.4%)
|
119803500
|
17970525000
|
91.59
|
110508187 (92.24%)
|
Paired-end,2×150bp
|
Average
|
140317022
|
21047553300
|
19842970437(94.27%)
|
127609381
|
19141407250
|
90.95
|
117407291(92.01%)
|
|
Q30: nucleotides with a quality value above 30 in reads.
Genome mapping rate: the percentage of reads mapped to the reference genome.
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 4,044 transcripts. Moreover, transcripts with cover > 3 in at least one sample were screened to obtain 4,008 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, 2,640 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 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 down-regulated (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. For example, the most enriched GO–BP terms were regulation of mRNA binding (GO: 1902415), positive regulation of mRNA binding (GO: 1902416), and positive regulation of RNA binding (GO: 1905216). In addition, several terms directly related to virus infection were found, for example, negative regulation of viral transcription (GO: 0032897 top 6), viral transcription (GO: 0019083 top 11), and regulation of viral transcription (GO: 0046782 top 12). The translation initiation factor binding (GO: 0031369), UDP-glucose 4-epimerase activity (GO: 0003978), glutamate–cysteine ligase activity (GO: 0004357), polysome (GO: 0005844), and glutamate–cysteine ligase complex (GO: 0017109), and viral replication complex (GO: 0019034) were the top three enriched GO-MP and GO-CC terms, respectively (Fig. 6C). The KEGG pathway analysis indicated that lncRNA target 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.
Immune-related mRNA and lncRNA analysis
Innate immunity is a crucial defense mechanism of cells against viral infections. We screened the DEGs related to innate immunity in HD11 cells infected with IBV for 36 h. These genes included CSF2, IFIT5, IL15, IL1RAPL1, IL22, IL8, MX1, NR1H4, S100A9, SYK, TRAF5, TRIM67, and ZFPM2. Among these, IL15, IL1RAPL1, and SYK were significantly down-regulated, whereas some anti-viral genes, such as IFIT5 and MX1, and some inflammatory factor genes, such as IL8 and IL22, were significantly up-regulated (Additional File 1).
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.
LncRNA–miRNA–mRNA regulation network analysis
LncRNAs can affect the gene expression through a variety of strategies [13]. We have previously reported that miR-146a-5p and gga-miR-30d had significant regulatory roles in IBV infection in HD11 cells [16, 17]. We screened for lncRNAs that interacted with miR-146a-5p and gga-miR-30d. The results revealed 1,563 lncRNAs to interact with gga-miR-30d. Among these, 30 lncRNAs were differentially expressed after 36 h of IBV infection. A total of 1,563 lncRNAs were found to interact with miR-146a-5p, and the expression of 32 lncRNAs changed significantly after 36 h of IBV infection. We constructed a miRNA–mRNA–lncRNA interaction network on the basis of potential interactions between them (Fig. 8 and Additional File 6). It is believed that these lncRNAs either function alone or compete for miR-146a-5p and gga-miR-30d with three genes (USP47, IRAK2, and TNFRSF18) to regulate IBV infection. In addition, eight lncRNAs (MSTRG.8180.7, MSTRG.4755.14, MSTRG.22271.3, MSTRG.21445.2, MSTRG.15550.10, ENSGALT00000104335, ENSGALT00000095670, and ENSGALT00000094718) were found to interact with both miR-146a-5p and gga-miR-30d (Fig. 8).
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 RNA-sequencing.
Table 2. Primers for qPCR
mRNA/lncRNA
|
Forword primer
|
Reverse primer
|
MSTRG.7587.19
|
GACCGTCGTGAGACAGGTTAGT
|
CTCCTCAGCCAAGCACATACAC
|
MSTRG.22271.3
|
GCAACTTCCAGAGACCACAGAA
|
CTCCAGCCACCAAGCACAAC
|
MSTRG.9359.1
|
GCTACCAAGCAATGTGTTCCAC
|
AGTGAGGCAAGTGAGGAGAAGG
|
MSTRG.6458.14
|
GGTGTGGCTGGTGGACTGTA
|
AGCCGCACCTGTAGTGAGAC
|
MSTRG.26120.58
|
ACGACATTAGGCGGTACGGAAT
|
GAGGCTGGAGTGGCACAAGA
|
IL22
S100A9
OASL
CYP3A5
GRID2
chGAPDH
|
CAGGCTCTCAGATCAGGACACT
TGGTGAAGTGATGCTCCTGAT
GGTCAAGCACTGGTACAAGGAG
ACAGCAGAAGGTGGTGAATGAA
ACTCACATCACCACAACAACCT
ACTGTCAAGGCTGAGAACGG
|
TTCATCATGTAGCAGCGGTTGT
GATGTTGGTGTTGGTGTTGGT
CAGGCGTAGATGGTCAGCAG
GCCATGTCGAGGTATTCCAACT
TGCCGAAGGAGAAGCCAGAA
GCTGAGGGAGCTGAGATGA
|