Identification of Differentially Expressed Non-coding RNA in Porcine Alveolar Macrophages from Tongcheng and Large White Pigs Responded to PRRSV

Porcine reproductive and respiratory syndrome (PRRS) is one of the most ruinous diseases in pig production. Our previous work showed that Tongcheng pigs (TC) were less susceptible to PRRS virus (PRRSV) than Large White (LW) pigs. To elucidate the difference in PRRSV resistance between the two breeds, small RNA-seq and ribo-zero RNA-seq were used to identify differentially expressed non-coding RNAs (including miRNAs and lincRNAs) responded to PRRSV in porcine alveolar macrophages (PAMs) from TC and LW pigs. Totally, 250 known mature miRNAs were detected. For LW pigs, there were 44 down-regulated and 67 up-regulated miRNAs in infection group; while for TC pigs, 12 down-regulated and 23 up-regulated miRNAs in TC infection group were identified. The target genes of the common differentially expressed miRNAs (DEmiRNAs) in these two breeds were enriched in immune-related processes, including apoptosis process, inflammatory response, T cell receptor signaling pathway and so on. In addition, 5 shared DEmiRNAs (miR-181, miR-1343, miR-296-3p, miR-199a-3p and miR-34c) were predicted to target PRRSV receptors, of which miR-199a-3p was validated to inhibit the expression of CD151. Interestingly, miR-378 and miR-10a-5p, which could inhibit PRRSV replication, displayed higher expression level in TC control group than that in LW control group. Contrarily, miR-145-5p and miR-328, which were specifically down-regulated in LW pigs, could target inhibitory immunoreceptors and may involve in immunosuppression caused by PRRSV. This indicates that DEmiRNAs are involved in the regulation of the immunosuppression and immune escape of the two breeds. Furthermore, we identified 616 lincRNA transcripts, of which 48 and 30 lincRNAs were differentially expressed in LW and TC pigs, respectively. LincRNA TCONS_00125566 may play an important role in the entire regulatory network, and was predicted to regulate the expression of immune-related genes through binding with miR-1343 competitively. In conclusion, this study provides an important resource for further revealing the interaction between host and virus, which will specify a new direction for anti-PRRSV research.

SCIEnTIfIC RepoRts | (2018) 8:15621 | DOI: 10.1038/s41598-018-33891-0 virus, PRRSV genome mutates with a high rate. Due to the poor cross-protection of the traditional vaccine for PRRSV variants, clinical prevention of PRRS is quite difficult, thus host genetic improvement of PRRSV resistance would be a better choice. As early in 1998, Harbul et al. reported that genetic differences of host could affect the susceptibility to PRRSV and clinical symptoms under PRRSV infection were different among breeds 8 . In the last decade, more artificial infection experiments were conducted within different pig breeds or populations with different backgrounds, which provides a strong evidence and support for the genetic contributions to PRRSV resistance.
In 2006, a highly pathogenic PRRSV (HP-PRRSV) broke out in China, and the epidemic persisted for a long time 9 , which made pig industry in a serious deficit state. While during the outbreak of HP-PRRSV, the Tongcheng (TC) pigs, a fine local variety in central of China, displayed extremely strong resistance to PRRSV 10 . Our previous study with artificial infection showed that TC pigs were less susceptible to PRRSV than LW pigs, manifesting as less tissue lesions, less virus load in serum, lower level of IL-10 but higher level of anti-viral cytokine interferon-gamma (IFN-γ) 11 . With RNA-sequencing, we compared the transcriptome difference of PAMs between TC and LW pigs, which revealed that TC pigs may promote the extravasation and migration of leukocytes to defend against PRRSV infection 12 .
With the development and widespread application of high-throughput sequencing technology, non-coding RNAs, including microRNA (miRNA) and long intergenic non-coding RNA (lincRNA), have been gradually recognized and found to participate in numerous biological processes. Some miRNAs were reported to regulate PRRSV proliferation. For example, miR-181 could strongly inhibit PRRSV replication through binding with ORF4 and PRRSV receptor CD163 13,14 ; miR-23, miR-378 and miR-505 were verified to directly target PRRSV genomic and subgenomic RNA 15 ; miR-26a could up-regulate the expression level of IFN-I and ISGs 16,17 ; miR-125b suppressed PRRSV replication through down-regulate NF-κB pathways 18 ; while miR-373 facilitated the replication of PRRSV by negative regulation of IFN-β 19 . Except for miRNAs, growing evidence suggested that lincRNAs could serve as ceRNAs through competing with mRNAs by binding to miRNAs and an increasing number of lincRNAs were confirmed to function as regulators of immune system [20][21][22] . Despite of the important roles of non-coding RNAs in immune response regulation, however, our previous work as well as other reported PRRSV-related transcriptome studies was focused on the function of protein-coding genes in PRRSV infection, less is known about the expression pattern of miRNAs and lincRNAs in the process of PRRSV infection. In this study, RNA sequencing with small RNA, and total RNA of ribosomal RNA removal, from PAMs of both TC and LW pigs were performed to obtain profiles of miRNAs and lincRNAs, which we hope and believe will provide another perspective for revealing the PRRSV resistance mechanism.

Materials and Methods
Sample preparation and RNA isolation. Pigs used in this study were selected from our previous research 11 . The HP-PRRSV strain was PRRSV WUH3 with a dose of 10 −5 TCID 50 for intramuscular challenge at 3 mL per 15 kg for experimental pigs. Twelve 5-week-old piglets of Tongcheng pigs (TC, n = 6) and Large White pigs (LW, n = 6) were slaughtered on 7-days post infection and then sampled, three for control and the other three for HP-PRRSV infection, respectively. PAMs were collected by bronchioalveolar lavage from lungs, lysed with TRIzol reagent (Thermol Fisher Scientific, Waltham, MA, USA), frozen in liquid nitrogen and stored at −80 °C until RNA extraction. Tissues including heart, liver, spleen, lung, kidney, brain, testis, mesenteric lymph nodes (MLN) and inguinal lymph nodes (ILN) were also collected for RNA isolation and analysis. All animal procedures were approved by the Ethical Committee for Animal Experiments at Huazhong Agricultural University, Wuhan, China (Animal experiment approval No. HZAUSW-2013-005, 08/27/2013). The animal experiments were performed at the Laboratory Animal Center of Huazhong Agricultural University.
Total RNA was extracted according to the manufacturer's instructions. RNA degradation and contamination were initially monitored on 1% agarose gels. RNA purity was determined using the NanoPhotometer ® spectrophotometer (IMPLEN, CA, USA). Then, the quantity of RNA was measured through Qubit ® RNA Assay Kit in Qubit ® 2.0 Flurometer (Life Technologies, CA, USA). RNA integrity was assessed using the RNA Nano 6000 Assay Kit of the Agilent Bioanalyzer 2100 system (Agilent Technologies, CA, USA).
Library construction for small RNA sequencing and data analysis. Library construction and sequencing were conducted in the Beijing Novogene Technology Co. Ltd. Company (Beijing, China). A total of 3 μg RNA per sample was used to generate a small RNA library using NEBNext® Multiplex Small RNA Library Prep Set for Illumina® (NEB, USA.) following manufacturer's recommendations and index codes were added to attribute sequences to each sample. Library quality was assessed on the Agilent Bioanalyzer 2100 system. Totally, twelve libraries were constructed. After the clustering of the index-coded samples on a cBot Cluster Generation System using TruSeq SR Cluster Kit v3-cBot-HS (Illumina), the libraries were sequenced on an Illumina Hiseq. 2000 platform and 50 bp single-end reads were generated.
To guarantee the quality of subsequent analysis, raw reads were processed to obtain clean reads as described: (1) Discard reads with poly-N > 10%, or poly-A, T, C, G and low-quality reads; (2) Remove reads containing 5′adapter contaminants; (3) Filter out reads without 3′adapters or insert tags; (4) Trim 3′adaptor sequences. Length-filtered clean reads were mapped to the reference pig genome (Sscrofa 10.2) by Bowtie. Then the aligned reads were compared with pig miRNA precursors in miRBase 21.0 to find known miRNA. Custom scripts were used to obtain miRNA counts, which were later normalized by TPM (transcript per million): mapped read count × 10 6 /total read count 23 . MiRNAs with at least five reads coverage in either library of the twelve were considered to be expressed and used for the following analysis. Differential expression analysis was performed by SARTools package 24 , setting p value < 0.05 and Fold Change (abbreviated as FC) ≥ 2 as threshold.
Target genes of miRNAs were predicted by software miRanda 25 . Genome sequence of PRRSV-WUH3 were downloaded from NCBI (accession number: HM853673.2) and used for miRNA binding sites prediction. Meanwhile, 3′UTRs of pig mRNA from Ensembl database were extracted to find underlying host target genes of expressed miRNAs. Also, an integrated analysis between differentially expressed (DE) miRNAs and DEmRNAs 12 was performed.
Library construction for ribo-zero RNA sequencing and data analysis. The detailed procedure of library construction for ribo-zero RNA sequencing is described in previous publication 12 . In brief, after removal of rRNA, twelve strand-specific cDNA libraries were constructed and 100 bp paired-end raw reads were generated. Clean data were obtained by removing reads containing adapter, reads containing ploy-N and low quality reads. Clean reads were mapped to the reference pig genome (Sscrofa 10.2) by TopHat v2.0.9 26 . The aligned reads were assembled into transcripts with Cufflinks and merged together with Cuffmerge 27 . To identify lincRNA, four steps were taken: (1) transcripts with single exon or shorter than 200nt were removed; (2) coding potential score, generated with CPC software, lower than −1 was supposed to be non-coding 28 ; (3) transcripts with FPKM value in all samples less than 0.01 were discarded; (4) transcripts overlapping with known gene or falling into 1 kb of any protein-coding genes were also removed 29 . Cuffdiff was used to identify DElincRNAs between groups, setting p value < 0.05 and |FC| ≥ 2 as threshold. DEmRNAs 12 within 500 kb of lincRNA were regarded as the cis target genes, while DEmRNAs co-expressed with DElincRNAs were considered as trans target genes. For the predicted target genes of DEmiRNAs and DElincRNAs, gene ontology (GO) and KEGG pathway analysis were performed using DAVID (https://david.ncifcrf.gov/), considering p value < 0.05 as the condition of significantly enrichment.
Construction of lincRNA-miRNA-mRNA regulatory network. Based on the sequence complementary and negative correlation between miRNA and target mRNA or target lincRNA, as well as the co-expression relationship between lincRNA and mRNA, the regulatory network of lincRNA-miRNA-mRNA was analyzed.
qRT-PCR validation of differentially expressed miRNA and lincRNA. First strand cDNA synthesis kit (Takara, Dalian, China) was used for reverse transcription, and common qRT-PCR and stem-loop qRT-PCR were conducted respectively to determine the relative expression of lincRNAs and miRNAs. Pig U6 snRNA and GAPDH were used as internal control for miRNA and lincRNA, respectively. Primers used for qRT-PCR are listed in Table S1. The qRT-PCR was performed in a total volume of 10 μL containing 5 μL 2 × SYBR Green Real-time PCR Master Mix, 0.2 μM of each primer, 1 μL cDNA and 3.6 μL ddH 2 O. The qRT-PCR reaction was conducted at 95 °C for 3 min, followed by 40 cycles of 95 °C for 15 s, 60 °C for 15 s, 72 °C for 20 s. All of the reactions were run in triplicate. Relative gene expression level was calculated using the 2 −ΔΔCT method 30 . Student's t-test was performed for data analysis, p value < 0.05 and p value < 0.01 was considered as the statistical threshold of significant difference and highly significant difference, respectively.

Results
Overview of miRNA sequencing. To reveal the expression difference of non-coding RNAs in PAMs from PRRSV infected TC and LW pigs, twelve small RNA libraries and twelve cDNA libraries were constructed and sequenced. The experiment was divided into four groups: TC_CON, TC_INF, LW_CON and LW_INF. For small RNA sequencing, the numbers of raw reads which covered at least 51.6 million bases ranged from 10M to 13M (Table 1). Q20 was above 97.89% and Q30 was above 94.70%, which indicated a high accuracy of sequencing data. Clean reads accounted for 96.92% to 98.53% of the raw reads. The length of clean reads was mainly from 20 nt to 24 nt, which is in accord with the general length range of miRNA. Of the 411 annotated mature miRNA, 208 miR-NAs were transcribed in all sequenced individuals. We analyzed the first nucleotide bias of the identified miRNA, which revealed that it tended to be uracil at most cases (Fig. S1).

Identification of DEmiRNAs and qRT-PCR validation.
Setting p value < 0.05 and |FC| ≥ 2 as the criteria for screening DEmiRNAs, 67 up-regulated and 44 down-regulated miRNAs were identified to be differentially expressed after PRRSV infection in LW pigs. Also, compared with TC_Control group, 23 up-regulated and 12 down-regulated miRNAs were identified in TC_Infection group (Table S2). Additionally, there were 31 and 30 miRNAs differentially expressed in the control and infection groups between these two breeds, respectively (Table S2; Fig. 1a). It is worth noting that miR-451, miR-486, miR-199b-3p, miR-199a-3p, miR-199b-5p and miR-31 were differentially expressed in all the four paired-comparisons (Table S2): CON_TC vs CON_LW, INF_TC vs INF_LW, LW_CON vs LW_INF and TC_CON vs TC_INF. Heatmap for all DEmiRNAs in the four groups revealed that TC_CON and LW_CON displayed similar expression pattern and clustered together, and the two infection groups clustered together. However, within one cluster, no matter control or infection group, some miR-NAs displayed extreme expression difference between the two breeds (Fig. 1b). miR-378, miR-378b-3p, miR-582 and miR-7135-3p were specifically down-regulated in PAMs of TC pigs, while there were 80 specific DEmiRNAs in PAMs of LW pigs. To validate the analysis results, six DEmiRNAs (miR-146b, miR-335, miR-378, miR-451, miR-532-5p and miR-9-1) were randomly selected for stem-loop qRT-PCR assay. The results indicated that trends of relative expression of qRT-PCR were consistent with small RNA-seq (Fig. 2). Functional analysis of the target genes of DEmiRNAs. miRanda was used to predict target genes of DEmiRNAs against PRRSV genome, PRRSV receptor and other host genes. Of the shared DEmiRNAs before and after infection in the two breeds, 23 miRNAs could bind to PRRSV genome (Table S3). Of which, miR-1343, miR-296-3p, miR-199a-3p and miR-34c could target PRRSV receptor CD151; miR-181a and miR-181b could target CD163 and miR-34c could target CD209. They were differentially expressed between control and infection group in both breeds (Table S2), which indicated that these common DEmiRNAs might play a key role in the regulation of PRRSV infection in pigs. In addition, four specific DEmiRNAs in TC pigs (miR-378, miR-378b-3p, miR-7135-3p and miR-582) were predicted to have ten binding sites on PRRSV genome. Then miR-199a-3p was randomly selected to validate the result of bioinformatics analysis. As shown in Fig. 3, over-expression of miR-199a-3p by mimics downregulated the expression of CD151.

Identification and qRT-PCR validation of differentially expressed lincRNAs.
Overview of ribo-zero RNA sequencing data were described in another publication 12 . After series of rigorous screening, 616 lincRNA transcripts were identified eventually, which were distributed on each chromosome except chromosome Y. The average length of these lincRNA transcripts is 2709 nt and more than 60% of the lincRNAs contain only two exons. Compared with protein-coding genes, the identified lincRNA has fewer exon number, lower coding potential and lower expression level (Fig. 5). Compared with TC control group, there were 11 significantly down-regulated and 19 up-regulated lincRNAs in TC infection group. Also, 29 lincRNAs and 19 lincRNAs were found to be down-regulated and up-regulated in LW infection group, respectively (Table S5). 11 DElincRNAs were randomly selected to validate through qRT-PCR. The results showed that the relative expression of lincRNAs was consistent with RNA-seq analysis result (Fig. 6a). Additionally, the expression level of 5 DElincRNAs was detected in ten tissues by RT-PCR, which revealed that the selected DElincRNAs expressed highly in immune system except TCONS_00146873 (Fig. 6b).

Functional annotation of the target genes of DElincRNAs.
To explore potential cis-regulatory lincRNAs, DEmRNAs within 500 kb upstream and downstream of DElincRNAs were analysed 31 . Totally, 150 DEmRNAs were found and they were enriched in regulation of signal transduction and reproduction-associated biological processes (Table S6). Based on trans-acting, Pearson Correlation Coefficient of DElincRNAs and DEmRNAs was calculated. Result showed that 915 DEmRNAs co-expressed with lincRNAs were found, which were significantly enriched in multiple immune-related biological processes (Table 3). Target genes of DElincRNAs between control and infection group in both breeds were enriched in the similar biological processes and KEGG pathways.
Construction of lincRNA-miRNA-mRNA regulatory network. As down-regulated miRNAs closely related to immunity, they were used to do the following analysis. Based on the sequence complementarity and negative correlation between miRNA and its targets, as well as the co-expression relationship between lincRNA and mRNA, the regulatory network of lincRNA-miRNA-mRNA was analysed (Fig. 7). The results showed that lincRNA TCONS_00125566 (TC: log 2

Discussion
PRRS has been pandemic for over 30 years, but it's still one of the main enemies of large-scale pig farms. Clinical prevention is difficult due to its high variability and the pressure of immune selection. Therefore, it is particularly important to enhance the genetic resistance to PRRSV and further improve the genetic structure. In this study, small RNA-seq and ribo-zero RNA-seq were performed to study the regulation of miRNA and lincRNA in the interaction between virus and host, which not only provides a new facet to investigate the differences of disease resistance of TC and LW pigs, but also lays foundation for further studying the antiviral function of miRNA and lincRNA. As a result, we identified some miRNAs and lincRNAs which may play important roles in PRRSV defending in two breeds, including miR-181, miR-296-3p, miR-744, miR-185, let-7c, miR-145-5p, miR-328, etc., as well as TCONS_00074289, TCONS_00037786, TCONS_00125566 and so on. As an important regulatory factor, miRNAs participate in diverse biological process, such as development, metabolism, cell proliferation and apoptosis, also in the intricate networks of host-pathogen interactions and innate immunity. In accordance with previous research, miR-181, validated to inhibit PRRSV replication 13,14 , was up-regulated post infection compared with control groups in both breeds. miR-296-3p and miR-744 was  decreased more than 5-fold in both breeds, and they could target LDOC1 (up-regulated post-infection), which served as a negative regulator of NF-κB 32 and was involved in anti-inflammatory response in both two breeds. Interestingly, miR-451 was up-regulated in both breeds, while it was expressed higher in TC pigs, even 17 times to LW infection group. Previous research indicated that miR-451 could inhibit the expression of pro-inflammatory factors 33 . Thus, its higher expression level probably contributed to relatively weakened tissue inflammation in TC pigs. It is also well known that PRRSV is an immunosuppressive disease, and IL-10 is a vital immunosuppressive factor during PRRSV infection 34 . In the current study, miR-185 and let-7c, which were specifically down-regulated in LW pigs, targeted SIGLEC5 that was up-regulated during PRRSV infection. Previous study showed that mouse SIGLEC5 enhanced IL-10 production while inhibiting TNF-α production in macrophages 35 . As a member of Siglec family, SIGLEC5 could facilitate the escape of pathogenic organisms from the control of the natural immune system 36 . Thus, LW pigs with higher expression of SIGLEC5 and IL10 manifested a state of immunosuppression. In addition, the specifically down-regulated miR-145-5p in LW pigs could target Cytotoxic T-lymphocyte associated protein 4 (CTLA-4), which had a significant lower expression in TC infection group than that in LW infection group. CTLA-4, as a negative regulator of T-cell activation, could reduce response to antigen 37 . Besides, miR-328, specifically up-regulated in LW pigs, could target programmed death ligand-1 Graft-versus-host disease 6 Allograft rejection 6 Autoimmune thyroid disease 6 Primary immunodeficiency 5 (PD-L1), which was up-regulated in LW infection group compared with TC infection group. Research showed that increased PD-L1 expression on antigen presenting cells (APCs) could aid in virus survival and decrease T-cell activity 38 . TIM-3, which could inhibit Th1 cells' immune activity by mediating apoptosis and thus inducing immune suppression, might be targeted by miR-362 and miR-365-3p. To sum up, the above mentioned DEmiRNAs were specifically down-regulated in LW pigs, while all of their target genes were up-regulated in LW pigs. Interestingly, these target genes together with IL-10 participated in immunosuppression, which indicated that some DEmiRNAs were involved in the regulation of immunosuppression in LW pigs. As a common DEmiRNAs of two breeds, miR-199 was effective in targeting and regulating HBV (Hepatitis B virus), thus curing the disease caused by HBV 39 . In this study, we found miR-199a-3p could inhibit CD151 expression at protein level. Since CD151 serves as the receptor of PRRSV, this would be an appropriate way to prevent PRRSV by miR-NAs through modulating their targeted genes. Apart from miRNAs, lincRNA is also an indispensable part of transcriptome. However, pig lincRNAs of immune system have been rarely reported. In the current study, we identified 616 differentially expressed lin-cRNA transcripts in PAMs, between control and infection group within the each pig breed or between control groups of the two pig breeds.
Co-expression analysis based on Pearson correlation coefficient identified many DEmRNAs, which were involved in immune-related process and pathway, including inflammatory response, apoptosis, cytokine-cytokine receptor interaction and so on. The regulation of these DEmRNAs by lincRNA was mostly mediated by trans-acting, which was consistent with a previous research 40 . Interestingly, one DElincRNA-TCONS_00074289-was expressed only in PAMs and lung, which indicated its function in the development of lung and PAMs. Recent reports have suggested that lincRNAs can potentially interact with other classes of non-coding RNAs including miRNAs and then regulate the expression of mRNA 41 .
In this study, we systematically analyzed the complex effects of the interactions between miRNAs and their target genes and provided lncRNA-miRNA-mRNA networks. TCONS_00037786 and TCONS_00125566 were recognised as important competing endogenous RNA. Of the target genes, GADD45B was an anti-apoptosis factor 42 ; SLAMF7 could inhibit the production of proinflammatory cytokines 43 ; Belong to chemokine receptor family, CXCR3 and CCR5 were involved in inflammatory response; DUSP4 is closely related with cell proliferation, differentiation and apoptosis through negatively regulating MAPK pathway 44 . Therefore, we speculated that these two lincRNAs together with miR-296-3p and miR-1343 took part in the regulation of inflammination and apoptosis in pigs, which of course need further experimental validation.
In summary, miRNAs were involved in PRRSV defending, and due to the difference of genetic background, some of them displayed specific expression pattern in one breed, and also conbined with lincRNAs to modulate physiological process such as anti-virus and regulation for host immune, which probably provides a new evidence for the genetic contributions to PRRSV resistance.

Conclusions
In this study, the expression and regulation of miRNA and lincRNA were analyzed in TC and LW pigs responded to PRRSV. Down-regulated miRNAs were involved in the regulation of PRRSV proliferation by participating in immune-related biological processes and pathway. Also, miRNA could target immunosuppressive receptor and ligand genes, leading to a stronger degree of immunosuppression in Large White pigs. Moreover, network interaction analysis was performed and some functional non-coding RNAs were found. This study lays the foundation for exploring the interaction mechanism between host and PRRSV and further revealing the disease resistance mechanisms of Tongcheng pigs.