Analysis of expression profiles of long noncoding RNAs and mRNAs in brains of mice infected by rabies virus by RNA sequencing

Rabies, caused by rabies virus (RABV), is still the deadliest infectious disease. Mechanism of host immune response upon RABV infection is not yet fully understood. Accumulating evidences suggest that long noncoding RNAs (lncRNAs) plays key roles in host antiviral responses. However, expression profile and function of lncRNAs in RABV infection remain unclear. In the present study, expression profile of lncRNAs and mRNAs profiles were investigated in RABV-infected brain tissues of mice by RNA sequencing. A total of 140 lncRNAs and 3,807 mRNAs were differentially expressed in RABV-infected animals. The functional annotation and enrichment analysis using Gene Oncology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) revealed that differentially expressed transcripts were predominantly involved in signaling pathways related to host immune response. The expression profiles of the selected lncRNAs in brains of mice during RABV infections were verified by quantitative real time polymerase chain reaction (qRT-PCR). To our knowledge, this is the first report to profile the lncRNA expression in RABV infected mice. Our findings provide insights into understanding the role of lncRNAs in host immune response against RABV infection.

Rabies is one of the deadliest zoonosis disease caused by rabies virus (RABV) 1 . It is nearly 100% fatal once clinical symptoms develop 2 . Rabies claims more than 60,000 human deaths annually, which is more than any other single zoonotic disease in the world. More than 80% of the deaths occurred in countries in Asia. China is the second most burden countries in the world. It showed that 40% of the deaths are children and 99% of the cases are resulted from bites of infected dogs 3 . Meanwhile, in developed countries like USA and Canada, bat RABV poses a serious threat to public health 4 .
RABV is a negative-stranded RNA virus that belongs to the family Rhabdoviridae, genus Lyssavirus, and species Rabies lyssavirus. Genome of RABV is approximately 12 kb and encodes five structural proteins, i.e. nucleoprotein (N), phosphoprotein (P), matrix protein (M), glycoprotein (G) and RNA-dependent RNA polymerase (L) 5 . Most RABV infections start from a dermal or muscular wound. RABV replicates locally in muscle tissue and then enters a neuron and spreads to motor neurons through synapses between muscles and motor neurons. It transports to central neural system (CNS) by retrograde axonal transport. Displaying of clinical symptoms means RABV reached the CNS 6 , where RABV elicit neuronal dysfunction and ultimately lead to death 7 .
Interferon (IFN)-mediated immune response is essential for protection against RABV infection 8 . Studies have shown that IFN-stimulated genes (ISGs), which were the effector of type I IFN response, exerted diverse antiviral effects 9,10 . Previous studies have demonstrated that deficiency in IFN production increased susceptibility to RABV in mouse model 11 . Although much advances have been achieved in prevention of RABV, the mechanism by which RABV causes fatal disease remains unclear.
Long noncoding RNAs (lncRNAs) are transcripts longer than 200 nucleotides and incapable of coding functional proteins. Most lncRNAs are capped at the 5′-end and polyadenylated at the 3′-end 12 . According to their genomic position, lncRNAs are generally classified as intergenic, intronic, bidirectional, antisense and pseudogene 13 . In the recent years, increasing evidences suggested that lncRNAs regulated numerous physiological processes, such as differentiation 14 , apoptosis 15 , development 16 , and immune responses 17 . In 2006, Rangarajan et al. 18 first reported a virus-induced lncRNA (VINC) in the CNS of mouse after Japanese encephalitis infection. Since then, many viral infections such as influenza (IAV) 19 , HIV 20 , hepatitis B 21 were reported to induce specific lncRNAs. LncRNA NRAV is downregulated during IAV infection and negatively regulates the transcription of ISGs 22 . Meanwhile, NRAV is the first lncRNA that is involved in inhibiting HIV-1 replication and facilitates the expression of antiviral genes during influenza virus and herpes simplex virus infection 23 . However, little is known about lncRNA expression profile and their regulating roles in immune responses during RABV infection.
To explore the role of lncRNAs during RABV infection, we analyzed the lncRNA expression profile in brain tissues of mice infected by RABV strain CVS-11 utilizing RNA sequencing (RNA-Seq). Our results indicated that RABV induced significant changes in lncRNA expression. Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis revealed that differentially expressed lncRNAs regulated immune response against RABV infection. To our knowledge, this is the first study to report profile the lncRNA expression in RABV infected mice. Our findings provide insights into understanding the role of lncRNAs in host immune response against RABV infection.

RNA-seq and identification of differentially expressed lncRNA.
To investigate lncRNA expression profile in mice infected with RABV, high-throughput RNA sequencing was performed on CVS11 infected brain tissues of mice. We sequenced 15 rRNA-deprived total RNA samples, including 5 brain tissues of mock-infected mice and 10 brain tissues of CVS-11 infected of mice. Each assay was duplicates. Average 80 million raw reads were produced for each sample using Illumina HiSeq platform by two-pair end sequencing. After removing the low-quality and adaptor sequences, clean reads were further analyzed.
Hierarchical clustering was used to analyze the lncRNA expression profiles in mock-or RABV-infected mice. As it was observed, the lncRNA expression profiles were significantly modified after RABV infection ( Fig. 2A). A total of 140 lncRNAs were differentially expressed in mice at days post infection (dpi) 8, with 38 lncRNAs up-regulated and 102 lncRNAs down-regulated (Fig. 2B). Of the dysregulated lncRNAs, 20 lncRNAs were changed with a fold change (FC) of more than 5.0, compared with mock infected group ( Table 1). The most up-regulated lncRNA was AW112010, with a FC of more than 140, and the most down-regulated transcript was a novel lncRNA, termed LNC_000415 with a FC of more than 9 ( Table 1).
The differentially expressed lncRNA in RABV-infected mice were widely scattered in all chromosomes, while the numbers were various in different chromosomes. Chromosome 7, 12 and 16 had the largest number of altered lncRNAs, while 18,19 and x had the least altered lncRNAs (Fig. 1C). Differential expression of mRNAs in brain tissues of mice between mock-and RABV-infected groups. We also examined the change of mRNA expression in brain tissues of mice post CVS-11 infection.
Hierarchical cluster analysis showed that mRNA expression profile was significantly changed in mice after RABV infection compared with mock-infected controls (Fig. 3A). A total of 3,807 lncRNAs were differentially expressed in the CVS-11 infected mice (FC ≥ 2 and P < 0.05), including 2,187 up-regulated and 1,620 down-regulated (Fig. 3B). To our surprise, 67 genes were upregulated with a fold change of more than 100 after infection. The most up-regulated gene was Cyba (FC = 4.75E + 30). The most down-regulated genes in our study are Alb, with a fold change of 31.16. The top 20 differentially expressed genes were listed in Table 2.  Similar to the distribution pattern of lncRNAs, the differentially expressed mRNAs in RABV infected mice were not equally scattered among chromosomes. The chromosome 2, 7 and 11 had the most differentially expressed mRNAs and chromosome 16 and 18 have the least numbers, while Y chromosome was absent of differentially expressed mRNAs (Fig. 3C).

Genomic features of lncRNAs and mRNAs in mice.
Then we systematically analyzed the basic features of the lncRNAs and compared them with protein-coding genes. As shown in Fig. 4A, the average expression levels of lncRNAs were lower than those of mRNAs. The length of transcripts of lncRNAs was shorter than   those of mRNAs (Fig. 4B). Moreover, the exon number of lncRNA was also less than those of mRNAs (Fig. 4C). Furthermore, most of the mRNAs had a longer Open Reading Frames (ORFs) than those of lncRNAs (Fig. 4D).

Functional prediction of RABV-induced lncRNAs.
To better understand the functions of differentially expressed lncRNAs in RABV infected mice, GO term and KEGG pathway analysis was performed to predict the functions of cis-and trans-target genes of differentially expressed lncRNAs. We found that the target genes of differentially expressed lncRNAs were highly enriched in biological processes like Intracellular signal transduction, Immune response and Synaptic transmission. The top 20 significant GO biological terms were presented in Fig. 5A. The targets of differentially expressed lncRNAs were involved in important signaling pathways, such as TNF signaling pathway, Toll-like receptor signaling pathway, NF-κB signaling pathway and MAPK signaling pathway. The top 20 significant enriched pathways were presented in Fig. 5B. These findings suggested that lncR-NAs regulate the immune responses during RABV infection.

Validation of differentially expressed lncRNAs by quantitative PCR. RNA-seq analysis indicated
that 140 lncRNA were differentially expressed post RABV infection. To validate the RNA-seq data, we investigated the expression levels of the eight most up-regulated lncRNAs at four time points after RABV infection, i.e. dpi 0, 3, 6 and 8, using quantitative real-time PCR (qRT-PCR). The results showed that expression patterns of these eight selected lncRNAs were consistent with RNA-seq data (Fig. 6). Moreover, the levels of these lncRNAs continuously increased from dpi 3 to dpi 8, which may reflect their correlations with progression of clinical symptom.

Discussion
RABV is still one of the deadliest zoonoses and remains as an important threat to public health in the world. Currently, although rabies is prevented by giving post-exposure prophylaxis (PEP) promptly, it lacks curable treatment. Effective protection of exposed subjects of rabies correlates with the induction of rabies-specific virus-neutralizing antibodies (VNAs). However, current vaccine not only requires multiple injections but also time-consuming and expensive, thus prevent many rabies exposed subjects away from timely vaccinated. As a result, rabies still cause around 70,000 deaths annually around the world despite efficacious vaccines are available 24 . Therefore, it is an emergent need to develop a cost-effective vaccine which elicits long-lasting immunity by a single vaccination and could ideally clear virus infection from the CNS. During infection, virus is detected by pattern-recognition receptors (PPRs), either canonical or non-canonical, which activate nuclear factor κB (NF-κB) and interferon regulatory transcription factors (IRFs) and induce the expression of type I interferons 25,26 . When RABV infections occur, the innate immune responses are promptly induced. PRRs are activated in the periphery and RABV is recognized in the CNS by retinoic-acid-inducible gene I (RIG-I), which sequentially activate NF-κB and type I IFN-regulated responses 27,28 . However, although much advance has been achieved on RABV biology and anti-RABV immune response, the mechanism underlies how RABV causes fatal disease is not fully understood. Previously, we found that protein-coding gene profile of host cell was significantly changed after RABV infection. We have identified some genes that function against viral replication, i.e. interferon-stimulated genes 15 (ISG15) and ubiquitin-like modifier-activating enzyme 7(UBA7) 29,30 . Recently, many studies have suggested that lncRNAs played key roles in the host immune response against viral infections 31,32 . However, the role of lncRNAs in RABV infection remained unclear. In the present study, we examined the expression profile of lncRNAs and mRNAs in brain tissues of mice after RABV infection. We identified 140 lncRNAs that were significantly differentially expressed between mock-or RABV-infected mice. To be noted, several lncRNAs, i.e. AW112010, AU020206, AI662270 and Ifi30 were up-regulated with a fold change of more than 30. The expression of lncRNAs has been confirmed by qRT-PCR. The dynamic change of lncRNA expression in brain tissues of mice further suggested that lncRNAs might play significant biological roles in RABV infection. Meanwhile, our results showed that 3,807 mRNAs were differentially expressed after infected with CVS-11, including 2,187 up-regulated and 1,620 down-regulated. We also characterized the genomic feature of lncRNAs in brain tissues of mice. Compared with mRNAs, lncRNAs are less enriched in expression, shorter in length, have fewer exons [33][34][35] . Previous studies have demonstrated that lncRNAs have poor primary sequence conservation compared to protein-coding genes 36 . It has been reported that less than 6% of zebrafish lncRNAs exhibited sequence conservation with lncRNAs of human or mouse and the sequence conservation of lncRNAs between human and other species were only about 12% 37,38 . We also evaluated the sequence conservation of the differentially expressed lncRNAs identified from RABV-infected mouse and found that only about 15% of the lncRNAs appeared to be conserved in human.
Unlike protein-coding genes or microRNAs, the sequences or structures of lncRNAs were currently uninformative for predicting its function 39 . In the present study, the function of lncRNAs was predicted according to their cis-or trans-target genes. GO terms were significantly enriched in biological processes like Intracellular signal transduction, Regulation of molecular function, Immune system process, Synaptic transmission. It suggested that lncRNAs induced by RABV infection may regulate the immune responses against RABV. KEGG pathway analysis showed that target genes of differentially expressed lncRNA were enriched in the pathways like NF-κB signaling pathway, Toll-like receptor signaling pathway, T cell receptor signaling pathway and TNF signaling pathway, which suggested that lncRNAs take part in host immune response against virus infection through various pathways.
In conclusion, the present study is for the first time to report the expression profile of lncRNAs upon RABV infection in mice. The results suggested that lncRNAs might have key roles in regulating immune responses post RABV infection and exert important biological effects.

RABV infection.
Six-to-eight-week-old male BALB/c mice were purchased from the Guangdong Medical Laboratory Animal Center. Mice were kept in an animal room with stable temperature and light, freely fed and drink. Mice were randomly assigned to two groups: ten for CVS-11 infected group and another ten for mock infected group. For virus infected group, mice were inoculated intracranially with 200 plaque-forming units (PFU) of CVS-11 in 50 µl DMEM, whereas the mock infected group was injected with equal volume of DMEM. On day 8 post infection, mice infected with CVS-11 showed clinical signs, i.e. disordered movement, hunched back, trembling and shaking. Mice from both groups were euthanized, and brain samples were collected. All animal experiments were performed following the National Institute of Health Guide for the Care and Use of Laboratory Animals, and the experimental protocols were approved by the Ethical Committee of Meizhou People's Hospital (Huangtang Hospital), Meizhou Hospital Affiliated to Sun Yat-sen University, Guangdong, China. All virus experiments were performed at Biosafety Level 2 laboratory.
Total RNA extraction. Total RNA was extracted from brain tissues of mice using RNeasy Kit (TianGene, Beijing, China) according to manufacturer's protocol. The quantity and purity of total RNA were evaluated by Nanodrop 2000. The ratio of A260/A280 should be from 1.8 to 2.0. RNA integrity was analyzed by the Bioanalyzer 2100 system (Agilent Technologies, CA, USA). Bioinformatic analysis. Raw data were filtered by removing the adaptors, low-quality reads and poly-N reads to obtain clean data using the SOAPnuke. The Q20, Q30, and GC information were calculated to evaluate the clean data. Then the filtered reads were mapped to the mice reference genome (version: mus_musculus. GRCm38) by Tophat 2. The transcripts were assembled with the mapped reads by reference annotation based transcripts (BRAT) method using Cufflink 40 .
Expressed profile of lncRNAs and mRNAs in brains of mice upon RABV infection were shown in Tables S1 and S2.

GO and KEGG pathway analyses.
To predict the target genes of differentially expressed lncRNAs, cis-and trans-analyses were performed. The genes located within a 10 kb window upstream or downstream of lncRNAs were classified as the cis target genes. The trans target genes were predicted on the expression levels of coding genes.
GO enrichment analyses were performed to identify biological processes associated with cis-or trans-target genes of lncRNA. KEGG was used to analyze the associated pathways of cis-or trans-target genes of the lncRNAs. A false discovery rate (FDR) was used to correct the P values. A corrected P value (Q values) <0.05 were considered significant.
Real-time RT-PCR assay. Real-time RT-PCR was performed to detect express of the selected lncRNAs using Luna ® Universal One-Step RT-qPCR Kit (New England Biolabs, USA) according to manufacturer's instructions. Primers used for validation of lncRNA expression were shown in Table 3. The amplify program is as follow: 95 °C for 30 s, 40 cycles (95 °C for 5 s, 60 °C for 30 s, and 72 °C for 30 s). The specificity of the amplified products was evaluated using dissociation curves. Relative expression of lncRNA were normalized to Glyceraldehyde 3-phosphate dehydrogenase (GAPDH) using the 2 −ΔΔCt method. The tests were triplicated. Statistical analysis. Data were analyzed using Student's t-test or one-way Analysis of variance (ANOVA) followed by Dunnett's multiple comparison test (compare all groups to the control group). All data are demonstrated as the means ± S.D. (*P < 0.05, **P < 0.01, ***P < 0.001). For correlation studies, a two-tailed non-parametric Spearman analysis was used. P ≤ 0.05 were considered as significant.