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

Solexa sequencing analysis

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.

Libraries	Known-miRNAs	novel-miRNAs	Total
PRV-Olf	502	345	847
PRV-Lung	500	440	940
PRV-Spleen	459	402	861
Olf	518	376	894
Lung	479	379	858
Spleen	452	349	801
Total	595	688	1283

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 differentially-expressed 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 differentially-expressed miRNAs were identified with enriched expression in olfactory bulbs, including six down-regulated 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 differentially-expressed miRNAs were identified with enriched expression in the lung, of which 10 were down-regulated miRNAs, 13 were up-regulated miRNAs and 19were non-conserved miRNAs.

In total, 199 differentially-expressed miRNAs were identified in the spleens of infected compared with uninfected rats. Moreover, 171differentially-expressed miRNAs 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 differentially-expressed miRNAs had the same or similar expression patterns (Fig. 6).

Gene annotations of differentially-expressed miRNAs

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.

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 miRNAs 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 miRNAs occurred only in the olfactory bulb. Therefore, we speculated that these miRNAs may be involved in the neurotropism of PRV. MiR-205, which was down-regulated, 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.

Analysis of differentially-expressed miRNAs specific to lungs of PRV-infected rats

PRV can cause respiratory symptoms, such as breathing difficulties and coughing. In the lungs of PRV-infected rats, 69 differentially-expressed miRNAs were identified. Forty-two differentially-expressed miRNAs were identified with enriched expression in the lung, including nine down-regulated miRNAs (miR-9a-3p, miR-653-5p, miR-488-3p, miR-361-5p, miR-329-5p, miR-129-5p, miR-128-1-5p, miR-124-3p, let-7c-1-3p), 10 up-regulated miRNAs (miR-134-5p, miR-223-3p, miR-31a-5p, miR-382-3p, miR-382-5p, miR-451-5p, miR-495, miR-505-5p, miR-6216, miR-672-3p) and 23 non-conserved miRNAs. 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 down-regulation 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 PRV-infected 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 β₃), 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 cell-mediated cellular immunity and humoral immunity. In the spleens of PRV-infected rats, 199 differentially-expressed miRNAs were identified. Moreover, 171 differentially-expressed 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 cell-surface glycoprotein involved in T-cell death, T-cell adhesion and leukocyte migration [46]. MiR-511-3p was down-regulated 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 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 K63-polyubiquitin 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.

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 (CertificateNo. 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.

Materials

PRV-XJ strain, a wild-type PRV strain isolated in Xinjin, China, was propagated in BHK-21 cells cultured in Dulbecco’s Modified Eagle’s Medium (DMEM) supplemented with 10% (v/v) fetal bovine serum. Specific-pathogen-free male rats (180 g bodyweight) were purchased from Chengdu Dashuo Experimental Animal Co. Ltd. All animal experiments were performed in accordance with the guidelines set out by the European Convention for the Protection of Vertebrate Animals used for Experimental and other Scientific Purposes and were approved by the Medical Faculty Experimentation Ethics Committee of Sichuan Agricultural University.

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. 1.5 μg RNA per sample was used as input material to construct small RNA libraries using the small RNA Sample Prep Kit. The 3′ and 5′ ends of small RNAs were connected to adaptors using T4 RNA Ligase 1 and T4 RNA Ligase 2 (truncated), and then cDNA was synthesized using reverse transcriptase. Finally, small RNA libraries were produced using PCR, agarose gel electrophoresis, and gel extraction. Small RNA sequencing was performed using an Illumina HiSeq 2500 Genome Analyzer.

PolyA RT-qPCR assay

Quantitation of mature miRNAs was performed as previously described. Briefly, 2 μg of total RNA was reverse transcribed into cDNA using the miRNA First Strand cDNA Synthesis Kit (Sangon Biotech) and miRNA-specific poly(A) primers. The reverse-transcription reaction was incubated at 37 °C for 60 min and at 85 °C for 5 min. To amplify miRNAs, PCR was performed a using universal reverse primer and a miRNA-specific forward primer. Real-time quantitative PCR (RT-qPCR) was performed using SYBR PremixEx Taq™ (Takara) and a Roche LightCycler® 96 (Swiss) with three replicates for testing expression differences. Triplicate reactions were performed in a total volume of 25 μL containing 12.5 μL of SYBR® Premix Ex Taq™ II (2×), 2 μL of cDNA template, 1 μL of each primer (0.5–2.5 μL), and 8.5 μL of ddH₂O. The thermal cycling conditions were as follows: 95 °C for 1 min, 40 cycles at 95 °C for 15 s and 60 °C for 30 s. The reference rat U6 spliceosomal RNA gene, which showed only slight expression variation in RNA-seq analyses, was used for expression data normalization. Relative miRNA expression was calculated using the comparative ΔΔCt method, and values were expressed as 2^-ΔΔCt.

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 high-quality 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.