Exploring the tissue tropism of pseudorabies virus based on miRNA level analysis

Background Pseudorabies virus (PRV, or suid herpesvirus, SuHV-1), a member of the herpesvirus family, has an extremely broad host range and threatens the pig industry in China. PRV can evade host innate immunity and infect the kidney, lung, brain and other tissues. At the same time, many studies have reported that microRNA (miRNA) can affect the replication of viruses by regulating gene expression levels. Results Here, to identify changes in miRNA expression and post-transcriptional regulation associated with PRV infection in the lung, spleen, and olfactory bulb, we sequenced small RNAs in tissues of rats infected or uninfected with PRV strain XJ (PRV-XJ). Sixty-one, 199 and 29 differentially-expressed miRNAs were identified in the lung, spleen, and olfactory bulb, respectively, of infected compared with uninfected rats. Among the miRNAs differentially-expressed in PRV-infected rats, 36, 171, and 15 miRNAs showed tissue-selective expression in the olfactory bulb, lung and spleen, respectively. All differentially-expressed miRNAs were analyzed for their GO functional annotations and KEGG pathway associations . Conclusions In PRV-XJ-infected rats, miRNAs were differentially expressed in the lung, spleen and olfactory bulb. These miRNAs were involved in regulating various pathways of the nervous, respiratory and immune systems, and may affect the tissue tropism of the virus and play pivotal roles in viral infection and proliferation.


Background
Pseudorabies virus (PRV, or suid herpesvirus, SuHV-1; herpesvirus family, alphaherpesvirus subfamily) shows the typical virion structure of the alphaherpesviruses. Pigs are the natural hosts of PRV and the source of most human infections [1]. The virus is spread by horizontal propagation, through food, water and feces. PRV can infect primary sensory neurons. After proliferating in primary sensory neurons, transmission can occur transsynaptically to second-order neurons [2,3].
PRV infection can result in respiratory diseases, such as hemorrhagic rhinitis and pharyngitis. PRV-infected pigs can develop interstitial pneumonia, hemorrhagic pneumonia, and significantly increasing widths of the alveolar septum and the intralobular interstitium. PRVinfected fetuses sometimes show hemorrhagic necrosis. In piglets, PRV infection can cause pulmonary edema, bleeding, and lobular pneumonia, and occasionally results in gray necrosis on the lung surface. PRV can evade host anti-viral immunity in two major ways. First, PRV escapes humoral immunity by interfering with complement regulatory pathways: for example, PRV glycoprotein C can bind porcine C3, and other viral proteins can block the binding of chemokines to their receptors. Second, PRV escapes cellular immunity: for example, PRV glycoprotein D can regulate CD112 [4], thus inhibiting NK cell-mediated degranulation and lysis, and other viral proteins can inhibit antigen presentation and expression of interferon-regulatory genes [5].
MicroRNAs (miRNAs) are recently-discovered small endogenous non-coding RNA molecules, and are only 21-25 nucleotides in length [6,7]. They are widely found in eukaryotes and serve to repress the expression of target genes through degradation of mRNAs at the transcriptional level after miRNA binding to 3'UTRs. MiRNAs, which account for about 2% of genomic sequences, are widespread in multicellular organisms [8]. They have been highly conserved over biological evolution and regulate a variety of physiological processes through involvement in development, apoptosis, fat metabolism, neuronal development, cell differentiation, and hormone secretion [9,10]. MicroRNAs are essential regulators of central nervous system development and also play an important role in stabilizing cell differentiation [11,12] .Previous research has reported alterations of miRNA expression during PRV infection.
In this study, we infected rats with PRV-XJ strain via intramuscular injection in the thigh. The miRNA content of various tissues was determined by high-throughput sequencing to identify differentially-expressed miRNAs. Our overall goals were to study the regulatory effects of miRNA on PRV infection and proliferation, to understand the tissue tropism of PRV, and to clarify the roles of miR-NAs in regulating host antiviral responses.
Previously-described and novel miRNAs were compared with reference genome sequences and identified using miRDeep2 software (Table 3). Due to the specificity of the Dicer and DCL enzymes, the lengths of mature miRNAs are primarily in the range of 20-24 nt. The length distributions of identified miRNAs and novel unidentified miRNAs are shown in Fig. 1.

MiRNA RT-qPCR assay
MiRNA RT-qPCR was used to verify the expression levels of host rat-encoded candidate miRNAs to validate the deep sequencing results. Ten differentially-expressed miRNAs in infected rats were randomly selected and quantified using this approach; non-conserved miRNAs were excluded. Each group has three repetitions (Fig. 2). Overall, the relative miRNA expression levels determined using deep sequencing were consistent with the levels of these miRNAs quantified using RT-qPCR.

Identification of differentially-expressed miRNAs
Levels of miRNA in various tissues were compared in infected and uninfected rats by analyzing the Solexa sequencing results. Differentially-expressed miRNAs between infected and uninfected samples were identified using IDEG6 [13] software using | log2 (FC) | > = 1 and FDR < = 0.01 as screening criteria. Fold change (FC) represents the ratio of miRNA expression between two samples. P-values can be interpreted as the probability of there being no true difference in miRNA expression given the experimental data. Finally, false discovery rate (FDR) was used as a reference indicator of differential miRNA expression. The number of differentiallyexpressed miRNAs in infected vs. uninfected rats is shown in Fig. 3. Tissue-specific differentially-expressed miRNAs were identified by comparing the sets of differentially-expressed miRNAs in each tissue (Fig. 4). Twenty-nine differentially-expressed miRNAs were identified by comparing miRNA sequences from the olfactory bulbs of infected and uninfected rats, of which six were non-conserved miRNAs (These miRNAs have different functions in different species) and 23 were conserved miRNAs. Excluding the six non-conserved miRNAs, four conserved miRNAs were up-regulated and 19 were down-regulated in infected rats. Fifteen differentiallyexpressed miRNAs were identified with enriched expression in olfactory bulbs, including six downregulated miRNAs, three up-regulated miRNAs and six non-conserved miRNAs.
Sixty-one differentially-expressed miRNAs were identified in the lungs of infected rats. Forty-two differentiallyexpressed miRNAs were identified with enriched expression in the lung, of which 10 were down-regulated miR-NAs, 13 were up-regulated miRNAs and 19were nonconserved miRNAs.
In total, 199 differentially-expressed miRNAs were identified in the spleens of infected compared with uninfected rats. Moreover, 171differentially-expressed miR-NAs were identified with enriched expression in the spleen, of which 58 were down-regulated miRNAs, 31 were up-regulated miRNAs and 82 were non-conserved miRNAs.
Differences in miRNA expression levels between two samples, as well as the statistical significance of these differences, can be seen in the volcano plot in Fig. 5. A cluster analysis was used to assess whether differentiallyexpressed miRNAs had the same or similar expression patterns (Fig. 6).

GO annotations and KEGG pathway analysis of target genes of differentially-expressed miRNAs
Gene Ontology (GO) annotations were determined for miRNA target genes to investigate their biological functions. The results revealed that miRNA target genes belonged to 56 functional subclasses according to the GO database. Secondary functional classes with significant differences included cell aggregation, cell killing, biological phase, biological growth, channel regulator and chemoattractants, all of which were in agreement with the results of a previous study (Fig. 7). We annotated these genes according to their KEGG pathways for statistical analysis to understand metabolic pathway involvement of the targets genes. KEGG annotation results were sorted by type of pathway (Fig. 8). An Enrichment Factor (EF) was calculated to analyze pathway enrichment, and Fisher's method was used to assess enrichment significance (Fig. 9). Larger EF values indicated that the target genes of differentially-expressed miRNAs were significantly enriched in the pertinent pathway.

Discussion
MicroRNAs have been linked to many biological processes, cellular components and molecular functions. Moreover, miRNAs can clearly be involved in viral infection: for instance, high expression of miR-146a can increase viral replication [23]. Many diseases can be clinically diagnosed and treated by reference to a variety of miRNA biomarkers. PRV shows differential infectivity to various host tissues. In this study, using Illumina Solexa deep sequencing technology, six small RNA libraries from the olfactory bulbs, lungs and spleens of PRV-infected and uninfected rats were constructed. Differentially-expressed miRNAs in each tissue were  obtained by comparing samples from infected and uninfected animals. Differentially-expressed miRNAs with enriched expression in specific tissues were identified by comparing levels of each differentially-expressed miRNA.
Analysis of differentially-expressed miRNAs specific to olfactory bulbs of PRV-infected rats In mature neurons, miRNAs can affect viral evolution and regulate viral tissue tropism. Therefore, studying changes in miRNA levels during virus replication is beneficial to understanding host defense mechanisms, as well as the mechanisms that viruses use to resist them. Twenty-nine differentially-expressed miRNAs were identified by comparing olfactory bulbs of infected and uninfected rats, which included six non-conserved miRNAs and 23 conserved miRNAs. Fifteen differentially-expressed miR-NAs with enriched expression in olfactory bulbs were identified, including down-regulated miRNAs (miR-1224, miR-181c-5p, miR-205, miR-298-3p, miR-412-5p, miR-873-5p) and up-regulated miRNAs (miR-26b-3p, miR-881-3p, miR-1-3p). PRV-associated differential expression of these miR-NAs occurred only in the olfactory bulb. Therefore, we speculated that these miRNAs may be involved in the neurotropism of PRV. MiR-205, which was downregulated, was predicted to be involved in regulation of neurotransmitter secretion according to its GO annotation. In Parkinson's disease (PD), Marques et al. observed lower expression of miR-205 in brain regions of 15 PD patients, but not in brains of multiple system atrophy (MSA) patients. Thus, miR-205 level was used as a biomarker to distinguish between PD and MSA [24]. Furthermore, a low level of serum miR-205 expression was found in glioma patients [25]. In our study, miR-26b-3p was up-regulated in olfactory bulbs. According to its GO annotation, miR-26b-3p may be involved in nervous system development, spinal cord motor neuron differentiation and neural tube formation. A high level of miR-26b-3p expression also occurs in other brain diseases, such as Alzheimer's disease. Satoh et al. showed that miR-26b-3p was expressed differentially between patient subgroups in AD using omiRas [26]. In our study, miR-873-5p was determined to be involved in the regulation of many neuronal cell biological processes according to its GO annotation, such as motor neuron axon guidance, chemorepulsion necessary for embryonic olfactory bulb interneuron precursor migration, and retinal ganglion cell axon guidance. Rats infected with PRV developed acute neurological symptoms such as trembling, twitching and muscle spasm. In addition, low levels of miR-873-5p expression were found in the hippocampus of rats in a model of temporal lobe epilepsy [27].  Thus, low levels of miR-873-5p expression may be due to neurological symptoms, which are typically associated with PRV infection.
Interestingly, miR-223-3p, miR-451-5p, miR-451-5p, miR-495, miR-505-5p, miR-361-5p, miR-129-5p, and let-7c-1-3p also showed significant differential expression in lung cancer patients. MiR-495 was up-regulated in lungs of PRV-infected rats. Researchers have found that miR-495 targeted the 3′-UTR of GRP78, leading to significant up-regulation of GRP78, especially in lung cancer. GRP78 (Gene ID: 3309) is a major regulator of the unfolded protein response (UPR) pathway. The UPR plays an important role in maintaining homeostasis within different lung cell types [28]. Studies have shown that miR-495 can target metastasis-associated protein 3 (MTA3), which is overexpressed in non-small cell lung cancer. Thus, miR-495 targeting of MTA3 may be involved in the regulation of lung cancer growth and migration [29]. In our study, miR-361-5p was minimally expressed in the lungs. As a tumor suppressor in lung cancer, miR-361-5p could directly target FOXM1 leading to downregulation of FOXM1 expression in lung cancer cells [30]. Low expression of miR-129-5p, which could directly target DLK1, was found in the lungs of PRVinfected rats. DLK1 belongs to the epidermal growth factor-like family. Apart from its roles in regulating cell differentiation, such as osteogenesis and adipogenesis, DLK1 had been shown to contribute to tumor cell invasion and is a negative regulator of Notch signaling. Furthermore, it was shown to regulate cell stemness [31]. The let-7 miRNA family has been shown to be closely involved in drug resistance [32], carcinogenesis [33], immunocyte differentiation [34,35] and proliferation [36]. Our study found that let-7c-1-3p was poorly expressed in lungs of PRV-infected rats. Zhao confirmed that both ITGB3 and MAP4K3 were target genes responsible for the functions of let-7c [37]. Upon let-7c overexpression, expression of ITGB3 and MAP4K3 were significantly repressed both at the protein and mRNA level. ITGB3 (also known as CD61 or integrin β 3 ), a member of the integrin family, had been shown to be involved in "outside-in" signaling transduction and cell adhesion [38]. MAP4K3, a member of the MAP4K family, is an important component in various signaling pathways [39]. Thus, let-7c-1-3p might play a pivotal role in PRV infection and proliferation.

Analysis of differentially-expressed miRNAs specific to spleens of PRV-infected rats
The spleen has a variety of important physiological functions in both the immune and endocrine systems. Immune function is the most important role of the spleen, both non-specifically and through its effects on T cells, B cellmediated cellular immunity and humoral immunity. In the spleens of PRV-infected rats, 199 differentially-expressed miRNAs were identified. Moreover, 171 differentiallyexpressed miRNAs were identified with enriched expression in the spleen, including 102 down-regulated miRNAs and 69 up-regulated miRNAs. Many of these miRNAs were involved in regulation of immune system development, lymphocyte activation, and the immune response. The miRNA miR-10b-5p could target MAP3K7. MiR-10b overexpression inhibited production of IL-17A by both total CD4 T cells and Th17 cells [40]. MiR-125a-3p was shown to positively regulate apoptosis and differentiation in THP-1 cells. Type I interferon (IFN1) is instrumental to the cellular response against viral infections. The STAT3 signaling pathway was negatively regulated by IFN1 binding to type-I interferon receptor 1 (IFNAR1) [41,42], and matrix metallopeptidase-9 (MMP9) is an important effector molecule in the STAT3 pathway. MiR-93-5p was identified to target IFNAR1, and MMP9 expression was increased by the miR-93-5p-IFNAR1 axis via the STAT3 pathway [43]. MiR-92b-3p directly targeted and reduced Gabra3 expression via the AKT/mTOR and JNK pathways [44]. MiR-664 directly bound to the 3′ UTR of phospholipid protein 2 (PLP2) as shown by luciferase assays, and suppressed the expression of PLP2 at the protein level as shown by western blotting. PLP2 encodes an integral membrane protein that localizes to the endoplasmic reticulum, which can multimerize and may function as an ion channel. PLP2 is necessary for down-regulation of human CD99 in the plasma membrane [45]. This result showed that CD99 expression was upregulated by reduced PLP2 levels. CD99 is a cellsurface glycoprotein involved in T-cell death, T-cell adhesion and leukocyte migration [46]. MiR-511-3p was downregulated in the spleens of PRV-infected rats, and a recent study showed that miR-511-3p may be a potential tumor suppressor in cancer [12,47]. Serine/threonine-specific protein kinase 3 (AKT3), which belongs to the AKT Fig. 5 Volcano plot of differentially-expressed miRNAs. Each point represents a miRNA; the x-axis represents the logarithm of the fold difference in miRNA expression between the two tissues; and the yaxis represents the negative logarithm of the false discovery rate. The blue points represent non-differentially expressed miRNAs, the red points represent up-regulated miRNAs, and the green points represent down-regulated miRNAs superfamily, plays an important role in cell transformation. MiR-511-3p could target the human AKT3 mRNA transcript to regulate cell transformation. MiR-486 can target KLF4 and Cezanne to jointly activate the NF-κB signaling pathway, as Cezanne and KLF4 are important negative regulators of NF-κB signaling: through different mechanisms, they can repress phospho-NF-κB-p65 (the active form of NF-κB). Cezanne, a newly identified member of the A20 family of deubiquitinases, can act as a potential tumor suppressor via NF-κB activity by deconjugating K63polyubiquitin chains from TRAF6 and RIP. However, KLF4, a zinc-finger transcription factor, can interact with phospho-NF-κB-p65 and promote COX-2 activity [48]. All of the differentially-expressed miRNAs in the spleen were reported to be involved in immune regulation through various pathways. Thus, these differentially-expressed miRNAs may be directly involved in the antiviral immune responses of the host.

Conclusions
We conducted a comprehensive analysis of differentiallyexpressed miRNAs in the olfactory bulbs, lungs and spleens of PRV-XJ infected rats and computationally predicted their putative targets. As reflected by our data, differentially-expressed miRNAs induced by PRV infection were involved in regulation of the nervous system, respiratory system and immune system, which may affect the tissue tropism of the virus as well as viral infection and proliferation. This dataset may be useful to understand the tissue tropism and mechanisms of immune escape of other herpesviruses, such as human herpes virus, at the miRNA level.

Animal ethics statement
All animal procedures used in the present study were conducted in accordance with good animal practices as defined by the laboratory animal use license (Certifica-teNo. SYXK (CHUAN) 2014-187). All work with rat was approved and supervised by the Committee on the Care and Use of Laboratory Animals of Sichuan and was conducted in the Animal Biotechnology Center of Sichuan Province at Sichuan Agricultural University. The animal interventions were performed in strict accordance with animal ethical standards. In addition, to limit any suffering, carbon dioxide (CO2) inhalation was used to anesthetize experimental animals in this study. All rats used for tissue collection were euthanized by rapid  decapitation at the endpoint, which could also reduce the hemocytes in the tissues used for sequencing analysis.

Sample collection and RNA extraction
Six rats were infected with PRV-XJ at a multiplicity of infection of 0.01 via intramuscular injection. A control group of six rats received an intramuscular injection with an equal volume of aseptic physiological saline. The infected animals were euthanized by rapid decapitation. The lungs, spleens and olfactory bulbs of rats were harvested 102 h post-infection to allow time for PRV to spread to organs and tissues, and stored in liquid nitrogen until RNA extraction. Total RNA was extracted from tissue samples using the miRNeasy Mini Kit. The miRNAs extracted from the olfactory bulbs, lungs and spleens of the six infected rats were pooled and labeled as samples PRV-Olf, PRV-Lung and PRV-Spleen, respectively. The miRNAs extracted from the olfactory bulbs, lungs and spleens of the six uninfected rats were pooled and labeled as samples Olf, Lung and Spleen, respectively.

Library construction and sequencing
The integrity of all RNA samples was verified via Nanodrop, Qubit 2.0, and Agilent 2100 Bioanalyzer.

Data analysis
The raw sequence data contained adaptor sequences as well as low quality sequences. To ensure miRNA sequence accuracy, we restricted the analysis to highquality sequences by excluding sequences with low quality, length < 18 nt, length > 30 nt as well as sequences containing ambiguous nucleotides ('N'). Using Bowtie(v1.0.0) [49] software, the filtered reads were compared with the Silva, GtRNAdb, Rfam and Repbase databases to filter out rRNA, tRNA, snRNA, RNA (snoRNA) and other ncRNA and repeat sequences. The unannotated reads were compared with the rat reference genome (rattus_norvegicus: ftp://ftp.ensembl.org/pub/release-85/fasta/rattus_norvegicus/) to obtain positional information (mapped reads) using miRDeep2 (v2.0.5) [50] software.

Analysis of differentially-expressed miRNAs
Previously-described and novel miRNAs were compared with sequences of reference genomes and identified using the miRDeep2 software as well as using knowledge of the biological characteristics of miRNAs(Mature sequence, Ring shaped structure, star sequence). The miRNA expression in each sample was statistically analyzed and expression levels were normalized using the TPM algorithm. Differentially-expressed miRNAs in the lungs, spleens and olfactory bulbs of PRV-infected rats were identified by comparing miRNA sequence sets derived from infected and uninfected rats using DESeq2 software. Differentially-expressed miRNAs associated with particular tissues were identified by comparing the miRNA sequence sets in lungs, spleens and olfactory bulbs of rat regardless of infection status. MiRnada(v3.3a) and RNAhybrid(v2.1.1) software were used to predict the target genes of differentially-expressed miRNAs. The Gene Ontology (GO) program Blast2GO (https://www.blas-t2go.com/) was used to annotate potential host target genes and to create GO annotation histograms, including histograms for cell component, biological process and molecular function. Target genes were analyzed for their KEGG Pathway signal transduction pathway using the KEGG (http://www.kegg.jp/kegg/pathway.html) public database and compared with the NR protein database to determine the protein annotations of target genes. An Enrichment Factor (EF) was calculated to analyze pathway enrichment, and Fisher's method was used to assess enrichment significance. Larger EF values indicated that the target genes of differentially-expressed miRNAs were significantly enriched in the pertinent pathway.