RESEARCH ARTICLE

Transcriptome sequencing reveals non-coding RNAs respond to porcine reproductive and respiratory syndrome virus and Haemophilus parasuis co-infection in Kele piglets

Jing Zhang1,#https://orcid.org/0000-0002-7281-0219, Chunping Zhao1,#https://orcid.org/0000-0002-1042-1138, Min Yao2https://orcid.org/0000-0003-4880-9397, Jing Qi1https://orcid.org/0000-0001-8191-5314, Ya Tan1https://orcid.org/0000-0003-2128-2396, Kaizhi Shi1,*https://orcid.org/0000-0001-7255-7180, Jing Wang1https://orcid.org/0000-0002-8884-3277, Sixuan Zhou1https://orcid.org/0000-0002-2589-136X, Zhixin Li3https://orcid.org/0000-0003-4484-0952
Author Information & Copyright
1Institute of Animal Husbandry and Veterinary Science, Guizhou Academy of Agricultural Sciences, Guiyang 550002, China
2Inspection and Testing Department, Guizhou Testing Center for Livestock and Poultry Germplasm, Guiyang 550002, China
3College of Animal Science, Guizhou University, Guiyang 550002, China

#These authors contributed equally to this work.

*Corresponding author: Kaizhi Shi, Institute of Animal Husbandry and Veterinary Science, Guizhou Academy of Agricultural Sciences, Guiyang 550002, China. Tel: +86-851-85400451, E-mail: nkyxms6462@163.com

© Copyright 2024 Korean Society of Animal Science and Technology. This is an Open-Access article distributed under the terms of the Creative Commons Attribution Non-Commercial License (http://creativecommons.org/licenses/by-nc/4.0/) which permits unrestricted non-commercial use, distribution, and reproduction in any medium, provided the original work is properly cited.

Received: Mar 02, 2024; Revised: May 11, 2024; Accepted: May 13, 2024

Published Online: Jul 31, 2024

Abstract

Co-infection with porcine reproductive and respiratory syndrome virus (PRRSV) and Haemophilus parasuis (HPS) has severely restricted the healthy development of pig breeding. Exploring disease resistance of non-coding RNAs in pigs co-infected with PRRSV and HPS is therefore critical to complement and elucidate the molecular mechanisms of disease resistance in Kele piglets and to innovate the use of local pig germplasm resources in China. RNA-seq of lungs from Kele piglets with single-infection of PRRSV or HPS and co-infection of both pathogens was performed. Two hundred and twenty-five differentially expressed long non-coding RNAs (DElncRNAs) and 30 DEmicroRNAs (DEmiRNAs) were identified and characterized in the PRRSV and HPS co-infection (PRRSV–HPS) group. Compared with the single-infection groups, 146 unique DElncRNAs, 17 unique DEmiRNAs, and 206 target differentially expressed genes (DEGs) were identified in the PRRSV–HPS group. The expression patterns of 20 DEmiRNAs and DElncRNAs confirmed by real-time quantitative polymerase chain reaction (RT-qPCR) were consistent with those determined by high-throughput sequencing. In the PRRSV–HPS group, the target DEGs were enriched in eight immune Gene Ontology terms relating to two unique DEmiRNAs and 16 DElncRNAs, and the unique target DEGs participated the host immune response to pathogens infection by affecting 15 immune-related Kyoto Encyclopedia of Genes and Genomes enrichment pathways. Notably, competitive endogenous RNA (ceRNA) networks of different groups were constructed, and the ssc-miR-671-5p miRNA was validated as a potential regulatory factor to regulate DTX4 and AEBP1 genes to achieve innate antiviral effects and inhibit pulmonary fibrosis by dual-luciferase reporter assays. These results provided insight into further study on the molecular mechanisms of resistance to PRRSV and HPS co-infection in Kele piglets.

Keywords: Kele piglets; Porcine reproductive and respiratory syndrome virus (PRRSV); Haemophilus parasuis; Co-infection; Non-coding RNAs; ceRNAs

INTRODUCTION

The Kele pig is a precious Chinese indigenous pig breed with disease resistance and economic value found in Guizhou Province located in the Karst mountainous area. Porcine respiratory disease complex (PRDC) is the most important and economically significant disease in the Chinese pig industry, and the economic loss caused by this disease is also the most serious in pig farms. As far as the incidence of pig farms worldwide is concerned, PRDC is generally a mixed infection of multiple pathogens [1]. The dual infection of porcine reproductive and respiratory syndrome virus (PRRSV) and Haemophilus parasuis (HPS) is epidemic and severe in China [2], resulting in the excessive use of antimicrobials and antivirals and causing reduced pig survival, feed conversion, and pork quality. Previous studies that focused on single pathogen invasion [3,4] and host defense mechanisms [5,6] were mainly concerned with protein-coding genes. Most studies on PRRSV and HPS co-infection only focus on molecular biological detection [7,8], comparisons of clinical symptoms and pathological changes of target organs [9], and immune responses in vitro [10,11]. According to our previous clinical observation and research, the results preliminarily suggested that Kele pigs have some ability to resist PRRSV and HPS concurrent invasion [12]. More information is still needed to better understand the host immune responses to PRRSV and HPS co-infection in Kele pigs.

Non-coding RNAs, including microRNAs (miRNAs) and long non-coding RNAs (lncRNAs), have been gradually recognized as key regulatory factors involved in the pathogenesis of PRDC. Some miRNAs were reported in pigs infected with PRRSV [13], porcine circovirus type 2 (PCV2) [14], Mycoplasma hyopneumoniae [15], and HPS [16]. The above studies have shown that miRNAs play critical roles in pathogen invasion, gene expression regulation, cell growth, differentiation, and apoptosis, immune response pathways. An increasing number of lncRNAs in the studies of HPS [17], PCV2 [18], and PRRSV [19,20] were revealed to function as crucial regulators of an immune response, and dysregulation of lncRNAs may also lead to disease. Furthermore, based on the competitive endogenous RNA (ceRNA) hypothesis, lncRNAs harboring miRNA response elements (MREs) act as miRNA “sponges”, which can compete with miRNA target genes for shared miRNAs. These MRE-sharing elements form the post-transcriptional ceRNA network to regulate mRNA expression [21]. Therefore, an efficient way to infer the potential function of miRNAs and lncRNAs related to pathogen infection and the host immune response is by exploring their relationships with annotated mRNAs.

To elucidate the mechanisms of resistance to PRDC pathogen invasion in Kele pigs for disease resistance breeding, it is necessary to research the Kele pig immune response to these two pathogens and the potential effect of differentially expressed transcripts on the development of a protective immune response against PRRSV and HPS co-infection. In this study, with an established PRRSV and HPS infection model, transcriptome sequencing was performed to explore the differences in miRNA and lncRNA expression patterns and functions between PRRSV and HPS infection alone and co-infection, and particularly the changes after co-infection. The functions of the target genes of the differentially expressed miRNAs (DEmiRNAs) and differentially expressed lncRNAs (DElncRNAs) were analyzed and a ceRNA regulatory network was constructed in the PRRSV and HPS co-infection group. Furthermore, the selected miRNAs and lncRNAs were validated by real-time quantitative polymerase chain reaction (RT-qPCR) and dual-luciferase reporter assays. This revealed a novel posttranscriptional regulation perspective that could help us understand the crucial roles of non-coding RNAs for regulating gene expression in immune resistance to PRRSV and HPS co-infection.

MATERIALS AND METHODS

Ethics statement

All animal procedures were conducted in accordance with the National Research Council Guide for the Care and Use of Laboratory Animals and approved by the Animal Ethics and Research Committee of Guizhou University (EAE-GZU-2020-P008).

Virus and bacterium

The viruses and bacteria used in this study were the same as in a previously published study [12]. The titer of the PRRSV genotype 2 GZBJ12 strain was 105.675 Tissue Culture Infective Dose (TCID50)/mL. The concentration of the HPS serotype 4 GZ strain inoculum was 1.66 × 109 CFU/mL.

Animal artificial challenge and sample collection

Thirteen 5-week-old Kele piglets were selected and divided into four groups (There were four piglets in the PRRSV–HPS group and three piglets in the HPS, PRRSV and control groups, respectively.). The piglets were inoculated with pathogens and managed for feeding as described in a previous study [12]. All animals were euthanized after 10 d post-infection. Pentobarbital sodium salt (Sigma, MO, USA) was used to relieve the pain of experimental animals. Lung tissues were collected, homogenized, and mixed with Trizol (TaKaRa, Kusatsu, Japan) for total RNA and stored at −80°C.

RNA and small RNA sequencing and analysis

Total RNA was extracted from 13 frozen lung tissue samples using Trizol (Invitrogen, Carlsbad, CA, USA) according to the manufacturer’s instructions. The RNA was used for miRNA and RNA sequencing on an Illumina platform at Novogene Technology (Beijing, China). After quality validation using an Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA) and gel electrophoresis, the small RNA library was constructed using the NEBNext® Multiplex Small RNA Library Prep set for Illumina® (New England Biolabs, New Ipswich, MA, USA), and its quality was assessed on an Agilent 2100 Bioanalyzer system (Agilent Technologies). Then, clustering of the index-coded samples was performed on a cBot Cluster Generation system using a TruSeq SR Cluster kit v3-cBot-HS (Illumina, San Diego, CA, USA) and the preparation libraries were subjected to deep sequencing. The detailed procedure of the library construction and quality control for RNA sequencing were carried out according to our related report [12].

The small RNA sequencing clean data were obtained by removing reads containing ploy-N, with 5’ adapter contaminants, without 3’ adapter or the insert tag, containing ploy A or T or G or C and low quality reads from raw data by using custom scripts. And then unique sequences with length in 18–26 nucleotide were mapped to the reference genome (Suscrofa 10.2.72) by Bowtie (-v 0 -k 1) [22] for further annotation analysis. Next, miRBase20.0 [23] was used as reference, the known miRNAs were analyzed by miRDeep2 (quantifier.pl -p -m -r -y -g 0 -T 10) [24]. After removing protein-coding genes, repeat sequences, rRNA, tRNA, snRNA, and snoRNA, small RNA tags by RepeatMasker (-species -nolow -no_is -norna -pa 8) and Rfam databases, the novel miRNAs were predicted using miREvo (-i -r -M -m -k -p 10 -g 50000) [25].

The raw RNA-sequencing data were filtered and quality control was carried out according to our previous study [12]. Novel potential lncRNA transcripts were identified following below steps: (1) Filter the transcripts with shorter than 200 base pair (bp); (2) Filter the transcripts with single exon; (3) Filter the known coding protein transcript; (4) Filter the transcripts with Fragments Per Kilobase of per Million fragments mapped (FPKM) < 0.5 and p < 0.05; (5) Coding-Non-Coding-Index (CNCI, score < 0) [26], Coding Potential Calculator (CPC, score < 0) [27], Pfam-scan (-E 0.001 --domE 0.001) [28], and phylogenetic codon substitution frequency (phyloCSF, --orf = ATGStop -frames = 3) [29] were used to distinguish mRNAs from lncRNAs, transcripts with coding potential predicted by any one of these four tools were filtered. Finally, those passing through all steps were filtered as the candidate set of lncRNAs.

Identification of DEmiRNAs and DElncRNAs

The miRNA and lncRNA expression levels were estimated by Transcripts Per Million reads (TPM) and FPKM, respectively. Differential expression analysis between challenged and normal groups was performed using DESeq2 [30]. MiRNAs with a log2|fold change| > 1 and p < 0.05 between the two groups (PRRSV–HPS vs. Control, PRRSV vs. Control, and HPS vs. Control) were identified as differentially expressed. DElncRNAs were identified between the two groups (PRRSV–HPS vs. Control, PRRSV vs. Control, and HPS vs. Control) by log2|fold change| > 2 and q < 0.05 according to the methods of our related transcriptome sequencing study, and differentially expressed genes (DEGs) were used in this study from our published data [12].

Prediction and functional annotation of the DEmiRNA and DElncRNA target genes

Predicting the target genes of the DEmiRNAs was performed by miRanda [31] with parameters: -sc 140 -en 10 -scale 4 -strict; and RNAhybrid [32] with parameters: -e 10 -p 0.05 -m 50,000; and the common targets were the final results. Target gene prediction of DElncRNAs was performed according to two roles: (1) Cis-role of the target gene prediction: a cis-role is the lncRNA acting on neighboring target genes. We searched coding genes 100,000 bp upstream and downstream of the lncRNA. (2) Trans-role of the target gene prediction: a trans-role is the lncRNA identifying non-neighboring genes by the expression level of correlation with Pearson correlation coefficient (PCC) |R| > 0.95 and p < 0.05. The expressed correlation between lncRNAs and coding genes was calculated by “cor.test” in R.

The Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses were performed using KOBAS software [33] with the hypergeometric test and p < 0.05 considered significantly enriched. Furthermore, Cytoscape3.7.2 [34] ClueGO plug-in was used to build the GO-target gene networks.

Construction of the ceRNA regulatory network

Clearly establishing the potential regulatory relationship between lncRNA-miRNA-mRNA gene pairs was based on the sequence complementarity and PCC (R > 0.5 or R < −0.5, and p < 0.05) between the shared miRNA and target mRNA or target lncRNA, as well as the co-expression relationship between the lncRNA and mRNA. The regulatory network of the ceRNAs was analyzed, and the network was constructed by Cytoscape3.7.2.

RT-qPCR validation of DEmiRNAs and DElncRNAs

The RT-qPCR validation was performed with the total RNA from lungs used in the transcriptome sequencing analysis with a PrimeScriptTM RT reagent kit (Perfect Real Time) (Takara). cDNA was used for RT-qPCR with TB Green® Premix Ex Taq™ II (Tli RNaseH Plus, Takara), following the instruction manual. Reactions were performed on an Eppendorf realplex Sequence Detection system (Eppendorf, Hamburg, Germany). Triplicate wells of reactions (25 µL) contained 12.5 µL of TB Green Premix Ex Taq II (2×), 2 µL of 50 ng/µL cDNA, 1 µL of 10 µM of each primer, and 8.5 µL of ddH2O. The RT-qPCR conditions were 95°C for 30 s, followed by 40 cycles at 95°C for 10 s, 60°C for 30 s, and 72°C for 20 s, followed by a melt curve. Ten DEmiRNAs and ten DElncRNAs were selected randomly from the three infected groups after being filtered and identified by sequencing. U6 [16] and GAPDH [20] were chosen as the internal reference genes for miRNAs and lncRNAs, and the 2−∆∆Ct method was used to calculate the fold change for gene expression. All primers of the selected genes were synthesized by Sangon Biotech (Shanghai, China) and are shown in Supplementary Table S1.

Dual-luciferase reporter assay

To verify the effects of ssc-miR-671-5p on target genes, the full length sequence of DTX4 and AEBP1 3’-UTRs were cloned into pmirGLO vector (Promega, Madison, WI, USA) downstream of the firefly luciferase gene, respectively. The Renilla luciferase gene was expressed as a reference reporter in the pmirGLO vector. Luciferase activity was measured using a Dual-Luciferase Reporter Assay System (Promega) according to the manufacturer’s instructions. Briefly, PmirGLO-DTX4-3’-UTR wild type (wt), PmirGLO-DTX4-3’-UTR mutated type (mt), PmirGLO-AEBP1-3’-UTR wt, or PmirGLO-AEBP1-3’-UTR mt was transfected into HEK293 cells along with ssc-miR-671-5p mimics or negative control miRNA mimics (mimics NC) (Genecrete, Wuhan, China) in 24-well plates using Lipofectamine 2000 (Invitrogen) according to the manufacturer’s instructions. The Dual-Luciferase Reporter Assay kit (Beyotime, Shanghai, China) and GloMax 96 Microplate Luminometer (Biosino, Beijing, China) were used to detect firefly and Renilla luciferase activities. Luciferase activity was normalized to Renilla luciferase activity at 48 hours after transfection.

Statistical analysis

Three replicates of RT-qPCR and dual luciferase reporter assays were performed for each sample. The results of dual-luciferase reporter assays were presented as means ± SD, p values were calculated using Student’s t test between two groups, and values of p < 0.05 were considered to indicate a statistically significant difference.

RESULTS

Overview of transcriptomics data

To reveal the expression difference of non-coding RNAs in lungs from PRRSV- and HPS-infected pigs, 13 small RNA libraries and 13 cDNA libraries were constructed and sequenced. For small RNA sequencing, the number of raw reads ranged from 0.786 G to 1.113 G (Supplementary Table S2). Q20 was above 99.19 % and Q30 was above 97.34 %, which indicated high accuracy of the sequencing data. Clean reads accounted for 98.85%–99.32% of the raw reads. The length of the clean reads was mainly from 20 nt to 24 nt, which is in accordance with the general length range of miRNAs. Of the 360 annotated mature miRNAs and 135 novel miRNAs, 326 miRNAs were transcribed in all sequenced individuals.

The lncRNA sequencing showed that 12,029 annotated lncRNAs and 17,338 novel lncRNAs were identified in all sequenced individuals. These identified lncRNAs were characterized as shorter in exon length, with fewer exons, shorter open reading frames (ORFs), and lower expression levels than protein-coding genes (Supplementary Figs. S1A, S1B, S1C, and S1D). Among the lncRNAs, 51.1% were intronic lncRNAs, 34.4% were lncRNAs, and 14.5% accounted for antisense lncRNAs (Supplementary Fig. S2).

Identification of DElncRNAs and DEmiRNAs in the different groups

DEmiRNAs and DElncRNAs were identified to understand the non-coding RNA changes in different groups. Compared with controls, 13 significant DEmiRNAs (10 upregulated miRNAs and 3 downregulated miRNAs) and 159 significant DElncRNA (66 upregulated lncRNAs and 93 downregulated lncRNAs) were identified in the HPS group (Figs. 1A and 1D). There were 20 significant DEmiRNAs (11 upregulated miRNAs and 9 downregulated miRNAs) and 122 significant DElncRNAs (79 upregulated lncRNAs and 43 downregulated lncRNAs) in the PRRSV group (Figs. 1B and 1E), and 36 significant DEmiRNAs (20 upregulated miRNAs and 10 downregulated miRNAs) and 225 significant DElncRNAs (144 upregulated lncRNAs and 81 downregulated lncRNAs) in the PRRSV–HPS group (Figs. 1C and 1F).

jast-66-4-663-g1
Fig. 1. Volcano plots of differentially expressed miRNAs and lncRNAs from the HPS (A, D), PRRSV (B, E), and PRRSV–HPS (C, F) groups. down: downregulated, blue dot; non: no differential expression, grey dot; up: upregulated, red dot. miRNAs, microRNAs; lncRNAs, long non-coding RNAs; HPS, Haemophilus parasuis; PRRSV, porcine reproductive and respiratory syndrome virus.
Download Original Figure

As shown in the Venn diagrams of the differentially expressed miRNAs and lncRNAs, there were 15 shared DEmiRNAs in the three infected groups, and 17 unique DEmiRNAs, including two novel miRNAs (novel_209, novel_494), in the PRRSV–HPS group (Fig. 2A), 92 shared DElncRNAs in the three infected groups, and 146 unique DElncRNAs in the PRRSV–HPS group (Fig. 2B). This showed that unique miRNAs and lncRNAs were regulated in the PRRSV–HPS group, which may indicate a potential function with PRRSV and HPS dual infection.

jast-66-4-663-g2
Fig. 2. Venn diagrams of differentially expressed miRNAs (A) and lncRNAs (B) from HPS, PRRSV, and PRRSV–HPS groups. PRRSV, porcine reproductive and respiratory syndrome virus; HPS, Haemophilus parasuis; miRNAs, microRNAs.
Download Original Figure
Target genes of DEmiRNA and DElncRNA prediction in different groups

Non-coding RNAs play a role by affecting gene expression and by regulating biological processes and signaling pathways. Therefore, the primary method to explore DEmiRNAs and DElncRNAs is to predict their target genes and functions. The target genes of miRNAs were identified by miRanda and RNAhybrid, with 16, 58, and 115 target DEGs of DEmiRNAs in the HPS (SupplementaryTable S3A), PRRSV (Supplementary Table S3B), and PRRSV–HPS groups (Supplementary Table S3C), respectively. Meanwhile, potential cis-regulated and trans-regulated target genes of DElncRNAs were explored. In total, 30 cis-regulated and 90 trans-regulated DEGs were found in the HPS group (Supplementary Table S4A), 27 cis-regulated and 37 trans-regulated DEGs were found in the PRRSV group (Supplementary Table S4B), and 71 cis-regulated and 135 trans-regulated DEGs were found in the PRRSV–HPS group (Supplementary Table S4C). The Circos plots of DEmiRNAs, DElncRNAs, and differentially expressed target genes in the HPS (Supplementary Fig. S3A), PRRSV (Supplementary Fig. S3B), and PRRSV–HPS (Fig. 3A) groups were visualized, which showed a more complicated non-coding RNA and target gene regulation pattern in the PRRSV–HPS group.

jast-66-4-663-g3
Fig. 3. The positions of the DEmiRNAs, DElncRNAs, and target DEGs in the chromsomes and the relationship of the DEmiRNAs, DElncRNAs, and target DEGs in the PRRSV–HPS group. (A) The Circos plots of DEmiRNAs, DElncRNAs, and target differentially expressed genes in the PRRSV–HPS group. From outside in, the first layer is the pig genome chromosomes. The second layer shows the gene labels. The third layer represents the heatmap of DElncRNAs and DEGs by red and blue bars, respectively. The respective fold changes of the DElncRNAs and DEGs are shown in the fourth and fifth layers. The fifth circle represents the fold changes of DEmiRNAs. The network in the center of the Circos plot represents the relationship of DEmiRNA, DElncRNA and DEmRNA location. The red lines indicate the linked RNAs in the same chromosome, while blue is from different chromosomes. (B) The regulation network of DEmiRNAs, DElncRNAs, and target DEGs. The green rectangles indicate DEGs, the pink hexagons represent DElncRNAs, and the blue circles indicate DEmiRNAs. DEmiRNAs, differentially expressed micro RNAs; DElncRNAs, differentially expressed long non-coding RNAs; DEGs, differentially expressed genes; PRRSV, porcine reproductive and respiratory syndrome virus; HPS, Haemophilus parasuis.
Download Original Figure

Furthermore, compared with the HPS and PRRSV single-infection groups, 206 target DEGs were predicted by the unique 146 DElncRNAs and 17 DEmiRNAs in the PRRSV–HPS group (Fig. 3B), and some DEGs related to pathogen infection and the immune system were regulated. For example, matrix metalloproteinases (such as MMP8; log2 (fold change) = 22.825, q = 5.33e-06) and their inhibitors (such as TIMP1; log2 (fold change) = 3.981, q = 0.0017), associated with pulmonary fibrosis, were regulated by LNC_002890, LNC_005631, LNC_008554, LNC_009017, LNC_010161, LNC_011347, LNC_015493, LNC_015814, and LNC_017117; CD163 molecule (CD163; log2 (fold change) = 2.154, q = 0.0017), the membrane protein receptor of PRRSV, was regulated by LNC_008554 and LNC_017117; T lymphocyte surface glycoprotein beta chain (CD8B; log2 (fold change) = 2.178, q = 0.044), calcium/calmodulin-dependent protein kinase II alpha (CAMK2A; log2 (fold change) = -4.913, q = 0.039), and complement c1q A chain (C1QA; log2 (fold change) = 8.760, q = 0.00012) were regulated by ssc-miR-671-5p; and interleukin 21 receptor (IL21R; log2 (fold change) = 2.260, q = 0.033) was regulated by LNC_000282, LNC_015493, LNC_017117, and ssc-miR-7137-3p. The results indicated that the regulated non-coding RNAs might potentially regulate the expression of genes to influence pathogen invasion and the host immune response against PRRSV and HPS infection.

Differentially expressed lncRNAs and miRNAs are involved in the immune response pathways in the concurrent infection group

Next, we investigated whether the target DEGs of DElncRNAs and DEmiRNAs from PRRSV–HPS group influenced the immune response. The possible enrichment pathways and biological processes of the target DEGs were predicted and analyzed. The GO enrichments demonstrated that more immune-related biological processes were present in the PRRSV–HPS group compared with the HPS and PRRSV single-infection groups. The top 20 enrichments are shown in Supplementary Table S5. Eight immune biological processes were significantly enriched in the unique DElncRNAs and DEmiRNAs of the PRRSV–HPS group, including complement activation (GO: 0006956), interaction with host (GO: 0051701), and B cell-mediated immunity (GO: 0019724) (Fig. 4). These immune processes were activated by two special miRNAs (novel_209 and ssc-miR-671-5p) and 16 lncRNAs (including LNC_012285, LNC_014095, and LNC_015466), which regulated 13 DEGs, including C1S (log2 (fold change) = 2.071, q = 0.00054), CD46 (log2 (fold change) = 2.071, q = 0.00054), C1QA (log2 (fold change) = 8.760, q = 0.00012), CXCL8 (log2 (fold change) = 2.833, q = 0.0013), and CD163 (log2 (fold change) = 2.154, q = 0.0017).

jast-66-4-663-g4
Fig. 4. The enriched biological process items of targets of unique DElncRNAs and DEmiRNAs determined by the Cytoscape ClueGO plug-in for the PRRSV–HPS group. The darker the red color, the smaller the p value. Genes in blue font are related to immune system biological processes. DElncRNAs, differentially expressed long non-coding RNAs; DEmiRNAs, differentially expressed micro RNAs; GO, Gene Ontology; PRRSV, porcine reproductive and respiratory syndrome virus; HPS, Haemophilus parasuis.
Download Original Figure

KEGG pathway enrichment showed that important immune pathways were regulated by the target genes of regulated non-coding RNAs in the three infection groups (Supplementary Table S6). For example, in the networks of the miRNA/lncRNA-KEGG pathway, the HPS group was enriched in adherens junction (ssc-miR-21-3pPTPRF [log2 (fold change) = 3.121, q = 2.88e-08]), and viral protein interaction with cytokine and cytokine receptor and Toll-like receptor signaling pathway (LNC_001726IL6 [log2 (fold change) = 4.929, q = 0.0045]; LNC_015081, LNC_015082CXCL10 [log2 (fold change) = 2.57, q = 0.0011], CXCL11 [log2 (fold change) = 3.104, q = 0.040]) (Fig. 5A). Cell adhesion molecules (LNC_001035CD274 [log2 (fold change) = 5.544, q = 0.0091] and LNC_014466SLA-5 [log2 (fold change) = −20.992, q = 0.00017]), antigen processing and presentation (LNC_002022CTSL [log2 (fold change) = 3.882, q = 0.00023]), and salmonella infection (ssc-miR-4331-3pCCL3L1 [log2 (fold change) = 2.114, q = 1.46e-06]) were significantly enriched in the PRRSV group (Fig. 5B). The KEGG pathways of the HPS and PRRSV single-infection groups indicated that the HPS or PRRSV infection influenced the immune system through regulated miRNAs and lncRNAs.

jast-66-4-663-g5
Fig. 5. The KEGG pathways for regulated non-coding RNAs in different groups. The circles indicate KEGG pathways, the diamond shapes represent DElncRNAs, and the hexagons show DEmiRNAs. A gradient from blue to red means downregulated to upregulated, respectively. (A) HPS group. (B) PRRSV group. (C) PRRSV–HPS group. KEGG, Kyoto Encyclopedia of Genes and Genomes; DElncRNAs, differentially expressed long non-coding RNAs; DEmiRNAs, differentially expressed micro RNAs; HPS, Haemophilus parasuis; PRRSV, porcine reproductive and respiratory syndrome virus.
Download Original Figure

Compared with the HPS and PRRSV single-infection groups, more DEmiRNAs and DElncRNAs regulated immune-related KEGG pathways in the PRRSV–HPS group (Fig. 5C), such as ECM–receptor interaction, bacterial invasion of epithelial cells, and complement and coagulation cascades, which were enriched by novel_209, ssc-miR-671-5p, ssc-miR-7137-3p, LNC_015466, and LNC_017117. Detailed descriptions of the immune-related pathways in the PRRSV–HPS group by targets of unique DElncRNAs and DEmiRNAs are listed in Table 1, suggesting that they may affect complement and coagulation cascades and the HIF-1 signaling pathway to activate the host immune defense mechanisms. Across these results are regulated non-coding RNAs that influence or take part in the immune biological processes and signaling pathways to defend against PRRSV and HPS invasion, and more immune response pathways were enriched in the PRRSV–HPS group compared with the single-infection groups.

Table 1. The immune-related KEGG pathways by target genes of unique DEmiRNAs and DElncRNAs in the PRRSV–HPS group
KEGG pathway term Count p-value Gene name
Complement and coagulation cascades 6 0.000004240 C4BPA, CD46, C1QA, F13A1, C1S, C1R
HIF-1 signaling pathway 6 0.000202924 EIF4E2, TIMP1, HK3, CAMK2A, INSR, PGK1
ECM-receptor interaction 4 0.004097015 AGRN, ITGA10, ITGB5, DAG1
Staphylococcus aureus infection 3 0.008168818 C1S, C1R, C1QA
Endocytosis 6 0.009587744 GRK6, ZFYVE16, ARPC4, WIPF2, SH3GL1, SH3GL2
Chemokine signaling pathway 5 0.011332147 GRK6, VAV1, CXCL13, CCL8, CXCL8
PI3K-Akt signaling pathway 7 0.013666818 EIF4E2, ITGB5, ERBB4, RXRA, INSR, ITGA10, TSC2
Yersinia infection 3 0.015345267 VAV1, WIPF2, CXCL8
Focal adhesion 4 0.016119951 VAV1, ITGB5, ITGA10, PARVG
Adherens junction 2 0.016137438 PTPRF, INSR
PPAR signaling pathway 3 0.019938900 CYP7A1, FADS2, RXRA
Salmonella infection 2 0.023451319 CXCL8, ARPC4
Phagosome 3 0.023528555 CTSL, ITGB5, C1R
Viral protein interaction with cytokine and cytokine receptor 3 0.026487035 CXCL8, CXCL13, CCL8
mTOR signaling pathway 4 0.027305425 FZD9, INSR, EIF4E2, TSC2

KEGG, Kyoto Encyclopedia of Genes and Genomes.

Download Excel Table
Construction of lncRNA–miRNA–mRNA network

Next, we investigated the regulatory network between differentially expressed non-coding RNAs and target genes based on the “ceRNA hypothesis.” The lncRNA–miRNA–mRNA networks with the sequence complementarity and high PCC of expression levels between miRNAs and their targets (lncRNA or mRNA) were constructed for the HPS, PRRSV, and PRRSV–HPS groups (Supplementary Fig. S4), which revealed the more complicated lncRNA–miRNA–mRNA network in the PRRSV–HPS group.

Furthermore, based on the sequence complementarity and negative correlation between DEmiRNAs and their targets, as well as the co-expression relationship between lncRNA and mRNA, the ceRNA network was constructed and analyzed in the PRRSV-HPS group (Fig. 6). The result showed that ssc-miR-671-5p (log2 [fold change] = 1.183, p = 0.0018) was the hub miRNA of the regulatory network, which regulated four downregulated lncRNAs (LNC_001016, LNC_002194, LNC_003637, and LNC_009394) and five downregulated mRNAs, including DTX4 (log2 [fold change] = −26.266, q = 6.09e-08), CAMK2A (log2 [fold change] = −4.913, q = 0.039), FAM160A1 (log2 [fold change] = −25.821, q = 1.05e-07), SLC2A9 (log2 [fold change] = −6.597, q = 0.027), and AEBP1 (log2 [fold change] = −27.228, q = 1.76e-08). The novel_207 (log2 [fold change] = −2.936, p = 0.018) was the only miRNA with downregulated expression that regulated ADM5 (log2 [fold change] = 2.370, q = 0.000026). In addition, the effect of ssc-miR-106a, ssc-miR-20b, ssc-miR-7139-5p, and ssc-miR-2411 on the expression of ten lncRNAs (such as LNC_009394, LNC_003226, and LNC_003637) formed a closed regulatory network.

jast-66-4-663-g6
Fig. 6. Regulatory network of lncRNA–miRNA–mRNA of the PRRSV–HPS group. The rectangles indicate miRNAs, the hexagons represent lncRNAs, and the circles indicate mRNAs. A gradient from green to red means downregulated to upregulated, respectively. lncRNA, long non-coding RNA; miRNA, micro RNA; PRRSV, porcine reproductive and respiratory syndrome virus; HPS, Haemophilus parasuis.
Download Original Figure
RT-qPCR validation of differentially expressed lncRNAs and miRNAs

To validate the expression pattern of DEmiRNAs and DElncRNAs identified from RNA-seq, ten DEmiRNAs and ten DElncRNAs were selected for RT-qPCR assays with the U6 and GAPDH genes as the internal controls. Among the selected DEmiRNAs and DElncRNAs, seven DEmiRNAs (including ssc_miR-21-3p, ssc_miR-371-5p, and novel_209) and eight lncRNAs (including LNC_015980, LNC_002561, and LNC_004879) showed a consistent upregulated trend in the results of both RT-qPCR and RNA-seq. Three DEmiRNAs (ssc_miR-202-5p, novel_611, and novel_207) and two DElncRNAs (LNC_009367 and LNC_002788) were downregulated in both RT-qPCR and RNA-seq (Fig. 7A). Although the extent of the fold change varied between RT-qPCR and RNA-seq, the RT-qPCR assay confirmed the results of RNA-seq, and the methods displayed a strong correlation (miRNAs: R2 = 0.853, lncRNAs: R2 = 0.936) (Figs. 7B and 7C), indicating the high reliability of RNA-seq data.

jast-66-4-663-g7
Fig. 7. Identification of selected miRNAs and lncRNAs from RNA-seq by RT-qPCR. (A) Log2 (fold change) obtained from RT-qPCR and RNA-seq data. The x-axis shows the name of the selected differentially expressed non-coding RNAs, and the Y-axis shows the value of the log2 (fold change). The correlation of selected miRNAs (B) and lncRNAs (C) between RT-qPCR and RNA-seq data. FC, fold change. RT-qPCR, real-time quantitative polymerase chain reaction; lncRNA, long non-coding RNA; miRNA, micro RNA.
Download Original Figure
Ssc-miR-671-5p targets the 3’-UTR of DTX4 and AEBP1 genes

To validate the predicted ssc-miR-671-5p target genes (DTX4 and AEBP1) (Fig. 8A), luciferase reporter assays were performed in HEK293 cells. Ssc-miR-671-5p mimics cotransfection with pmirGLO- DTX4 3’-UTR wt resulted in a significant decrease in luciferase activity, whereas the opposite effect was observed when mimics NC cotransfection with pmirGLO- DTX4 3’-UTR wt (Fig. 8B). Moreover, mutations in the predicted DTX4 3’-UTR binding sites abolished this effect. The suppression function of ssc-miR-671-5p on AEBP1 gene was consistent with DTX4 gene (Fig. 8C). Therefore, DTX4 and AEBP1 are target genes of ssc-miR-671-5p.

jast-66-4-663-g8
Fig. 8. Ssc_miR_671_5p reduces the levels of DTX4 and AEBP1. (A) Schematic representation of wild-type and mutant pmir GLO-DTX4/AEBP1-3’-UTR miRNA expression vectors used in luciferase reporter assays. The nucleotides altered in the mutant binding site are colored red. Relative luciferase activities in HEK293 cells corresponding to DTX4-3’UTR (B) and AEBP1-3’UTR (C) after cotransfection with the pmirGLO-DTX4/AEBP1-3’-UTR reporter and the mimic NC, ssc_miR_671_5p mimic for 48 hours. Data are presented as the mean ± SD of three independent experiments (**p < 0.01 and ***p < 0.001, Student’s t test).
Download Original Figure

DISCUSSION

As a local pig breed in China, the Kele pig is famous for its resistance to diseases and for its high-quality pork, making it excellent breeding material. PRDC is a severe ailment in pig breeding in China. Our previous study reported that Kele piglets challenged with PRRSV and HPS showed changes in clinical symptoms, body temperature, histopathology, and transcriptome levels in the lung, inferring that reactive oxygen species (ROS) production and pulmonary fibrosis were affected in Kele piglets in response to PRRSV and HPS co-infection [12]. This study aimed to further explore the effects of PRRSV and HPS infection on the miRNA and lncRNA levels in Kele piglets, which improved our understanding of the immune response to these pathogens in the host. Based on the comparative analysis, 122 DElncRNAs and 20 DEmiRNAs were found in the PRRSV group. The number of DElncRNAs was greater than that reported by a previously published study, but the number of DEmiRNAs was less [20]. We inferred that the different pig species and pathogens are key factors influencing the amount of non-coding RNA identified in studies, which suggest that experiments with different breeds infected with the same pathogens could be taken in further studies to compare the immune responses of different breeds. We also found that more DElncRNAs and DEmiRNAs were identified in the co-infection group than these in the single infection groups, implying that multi-pathogens co-infection tends to stimulate more changes in non-coding RNAs, which may lead to more complex post-transcriptional immune regulation. Moreover, compared with protein-coding genes, lncRNAs in lung tissue had the same characteristics in different tissues and cells [20,35].

As an important regulatory factor of post-transcriptional regulation, miRNAs play an important role in many biological processes. There were eight DEmiRNAs shared by the three infected groups, among which ssc-miR-21-5p, ssc-miR-20b, ssc-miR-106a, ssc-miR-7-5p, and ssc-miR-21-3p were upregulated. The expression of ssc-miR-202-5p, novel_611, and novel_190 was downregulated. High expression of the shared DEmiRNA miR-21-5p can inhibit H2O2 (PTEN/AKT)-induced apoptosis of alveolar epithelial cells, thus improving acute lung injury induced by hyperoxia in rats [36], and it is also related to the epithelial–mesenchymal transition (EMT) in humans [37]. These studies showed that ssc-miR-21-5p was highly expressed in the lungs of pigs infected with pathogens, and the highest expression level in the co-infection group (log2 [fold change] = 1.7005) may be related to pulmonary fibrosis after acute lung injury. MiR-20b can induce polarization of M2 alveolar macrophages in mice. Macrophages present in alveoli and the mesenchyme are activated to form M2 alveolar macrophages, participating in the process of pulmonary fibrosis after lung tissue is continuously stimulated by antigen [38]. Ssc-miR-20b also had the highest expression in the co-infection group, indicating that both pathogens can lead to pulmonary fibrosis formation after continuous infection, and PRRSV and HPS co-infection promoted the process. Both miR-106a and miR-7-5p were found to be associated with human lung cancer [39,40]. Ssc-miR-21-3p, shown to promote PRRSV replication in MARC-145 cells in a previous study [41], was upregulated post-infection in the three treatment groups compared with the control. Ssc-miR-202-5p was downregulated about 5-fold in the PRRSV and PRRSV–HPS groups compared with the control group, and is related to EMT [42], apoptosis [43] and innate immunity [44]. In the PRRSV–HPS group, among 21 unique DEmiRNAs, novel_209, novel_494, ssc_miR-371-5p, and novel_207 were differentially expressed more than 2-fold compared with the single-infection groups. The overexpression of ssc_miR-371-5p may inhibit proliferation and promote apoptosis of mesangial cells by directly targeting HIF-1α in humans [45], suggesting that upregulation of ssc-miR-371-5p after infection might also induce apoptosis.

Increasingly more evidence shows that miRNAs can regulate expression of their targets through lncRNAs, thus affecting gene function [20,35,46,47]. Although many lncRNAs have been found, the potential role of lncRNAs in the pig respiratory system is far from understood. A large number of target genes of DEmiRNAs and DElncRNAs participating in host immune regulation were found in this research. In the PRRSV–HPS group, the CD163 molecule, having long been confirmed to be the receptor protein of PRRSV [48], was regulated by two DElncRNAs (LNC_008554 and LNC_017117), so we speculated that these two DElncRNAs may participate in the invasion process of PRRSV. MMPs and TIMPs were regulated by nine DElncRNAs and are related to pulmonary fibrosis [49], indicating that these nine DElncRNAs may be involved in the regulation of pulmonary fibrosis. Some studies have shown that C1QA [50], CD8B [51], and IL21R [52] are related to the classical complement pathway and immune response of the host. In our study, C1QA and CD8B were regulated by ssc-miR-671-5p, IL21R were regulated by three DElncRNAs (LNC_000282, LNC_015493, LNC_017117) and one DEmiRNA (ssc-miR-7137-3p), which showed that these DEmiRNAs and DElncRNAs may regulate the expression of genes.

We also found that there were many differentially expressed non-coding RNAs regulating immune-related biological processes and immune-related pathways. Particularly in the co-infection group, two special DEmiRNAs and 16 DElncRNAs were involved in regulating eight immune-related biological processes, such as complement activation, interaction with host, and B cell-mediated immunity. Compared with PRRSV or HPS infection alone, KEGG pathways in the PRRSV–HPS group were more complex and were mainly involved in the regulation of complement and coagulation cascades and the HIF-1 signaling pathway. Many lncRNAs and miRNAs reported to be related to the porcine respiratory system were also involved in host immune response [20], consistent with our results. Therefore, we suggested that these specific differentially expressed non-coding RNAs play critical roles in the innate immune and adaptive immune processes.

In the ceRNA network, we observed that ssc-miR-671-5p interacted with AEBP1, DTX4, SLC2A, FAM160A1, CAMK2A, and four DElncRNAs in the PRRSV–HPS group. Moreover, the results of the dual-luciferase reporter assays showed that the relative luciferase activities were significantly lower in cells transfected with the ssc-miR-671-5p mimics than in cells transfected with mimics NC in DTX4-wt and AEBP1-wt groups. This validated the role of ssc-miR-671-5p in regulating the expression of the DTX4 and AEBP1 genes. It was reported that silencing AEBP1 markedly suppressed the proliferation, migration, invasion, metastasis, and EMT of guanine cytosine (GC) cells [53,54]. However, overexpression of ssc-miR-671-5p led to the downregulation of AEBP1 expression in the present study, which suggests that ssc-miR-671-5p inhibits pulmonary fibrosis by regulating the expression of the AEBP1 gene. The published research showed that NLRP4-DTX4 mediated TBK1 degradation to promote virus infection, while the DTX4 gene knock-down eliminated TBK1 ubiquitination and degradation, and enhanced TBK1 and transcription factor IRF3 phosphorylation to achieve innate antiviral effects [55]. It is speculated that the expression of DTX4 is downregulated under the influence of ssc-miR-671-5p after persistent co-infection, and the body produces an antiviral response. Therefore, ssc-miR-671-5p may be a key promoter of the host’s inhibition of lung diseases and its antiviral mechanism, which deserves further study. Remarkably, novel_207 was the only downregulated miRNA in the ceRNA network, targeting the ADM5 gene. ADM5 may induce cell migration and invasion [56]. Therefore, the downregulation of novel_207 may involve cell migration and invasion by interacting with ADM5 after infection with PRRSV and HPS.

CONCLUSION

Differentially expressed lncRNAs and miRNAs have been identified between PRRSV and HPS single-infection and co-infection groups compared with the control group, the complex relationships between DEmiRNAs, DElncRNAs, and their target genes were analyzed in three infected groups. The ceRNA regulatory networks were constructed for the different groups. Two miRNAs (ssc-miR-671-5p and novel_207) were found that play important roles after PRRSV and HPS concurrent infection. Ssc-miR-671-5p is a regulator for DTX4 and AEBP1 genes, which was verified by dual-luciferase reporter assays. This explained the potential relationship between the two pathogenic infections of Chinese local pigs at the transcriptome level.

SUPPLEMENTARY MATERIALS

Competing interests

No potential conflict of interest relevant to this article was reported.

Funding sources

This research was funded by the funds of Guizhou Science and Technology Department (grant number [2021] 5616 and [2022] KEY032) and the youth fund of Guizhou Academy of Agricultural Sciences (grant number [2022]23).

Acknowledgements

Not applicable.

Availability of data and material

All RNA and small RNA sequencing data generated in this study are available in the NIH short read archive under accession numbers PRJNA661122 and PRJNA930985, respectively.

Authors’ contributions

Conceptualization: Shi K.

Data curation: Yao M, Tan Y.

Formal analysis: Zhang J, Qi J.

Methodology: Shi K.

Software: Zhao C, Yao M.

Validation: Zhang J, Zhao C, Qi J.

Investigation: Wang J, Zhou S, Li Z.

Writing - original draft: Zhang J, Zhao C.

Writing - review & editing: Zhang J, Zhao C, Yao M, Qi J, Tan Y, Shi K, Wang J, Zhou S, Li Z.

Ethics approval and consent to participate

All animal procedures were conducted in accordance with the National Research Council Guide for the Care and Use of Laboratory Animals and approved by the Animal Ethics and Re-search Committee of Guizhou University (EAE-GZU-2020-P008).

REFERENCES

1.

Cohen LM, Grøntvedt CA, Klem TB, Gulliksen SM, Ranheim B, Nielsen JP, et al. A descriptive study of acute outbreaks of respiratory disease in Norwegian fattening pig herds. Acta Vet Scand. 2020; 62:35

2.

Lin WH, Shih HC, Lin CF, Yang CY, Chang YF, Lin CN, et al. Molecular serotyping of Haemophilus parasuis isolated from diseased pigs and the relationship between serovars and pathological patterns in Taiwan. PeerJ. 2018; 6e6017

3.

Ke H, Han M, Kim J, Gustin KE, Yoo D. Porcine reproductive and respiratory syndrome virus nonstructural protein 1 beta interacts with nucleoporin 62 to promote viral replication and immune evasion. J Virol. 2019; 93:e00469-19

4.

Li X, Chi SQ, Wu LY, Liu C, Sun T, Hong J, et al. PK/PD modeling of Ceftiofur Sodium against Haemophilus parasuis infection in pigs. BMC Vet Res. 2019; 15:272

5.

Wang J, Liu JY, Shao KY, Han YQ, Li GL, Ming SL, et al. Porcine reproductive and respiratory syndrome virus activates lipophagy to facilitate viral replication through downregulation of NDRG1 expression. J Virol. 2019; 93:e00526-19

6.

Yue C, Li J, Jin H, Hua K, Zhou W, Wang Y, et al. Autophagy is a defense mechanism inhibiting invasion and inflammation during high-virulent Haemophilus parasuis infection in PK-15 cells. Front Cell Infect Microbiol. 2019; 9:93

7.

Cheong Y, Oh C, Lee K, Cho K. Survey of porcine respiratory disease complex-associated pathogens among commercial pig farms in Korea via oral fluid method. J Vet Sci. 2017; 18:283-9

8.

Sunaga F, Tsuchiaka S, Kishimoto M, Aoki H, Kakinoki M, Kure K, et al. Development of a one-run real-time PCR detection system for pathogens associated with porcine respiratory diseases. J Vet Med Sci. 2020; 82:217-23

9.

Li J, Wang S, Li C, Wang C, Liu Y, Wang G, et al. Secondary Haemophilus parasuis infection enhances highly pathogenic porcine reproductive and respiratory syndrome virus (HP-PRRSV) infection-mediated inflammatory responses. Vet Microbiol. 2017; 204:35-42

10.

Kavanová L, Prodělalová J, Nedbalcová K, Matiašovic J, Volf J, Faldyna M, et al. Immune response of porcine alveolar macrophages to a concurrent infection with porcine reproductive and respiratory syndrome virus and Haemophilus parasuis in vitro. Vet Microbiol. 2015; 180:28-35

11.

Kavanová L, Matiašková K, Levá L, Nedbalcová K, Matiašovic J, Faldyna M, et al. Concurrent infection of monocyte-derived macrophages with porcine reproductive and respiratory syndrome virus and Haemophilus parasuis: a role of IFNα in pathogenesis of co-infections. Vet Microbiol. 2018; 225:64-71

12.

Zhang J, Wang J. Transcriptome profiling identifies immune response genes against porcine reproductive and respiratory syndrome virus and Haemophilus parasuis co-infection in the lungs of piglets. J Vet Sci. 2022; 23e2

13.

Do DN, Dudemaine PL, Mathur M, Suravajhala P, Zhao X, Ibeagha-Awemu EM. miRNA regulatory functions in farm animal diseases, and biomarker potentials for effective therapies. Int J Mol Sci. 2021; 22:3080

14.

Li C, Sun Y, Li J, Jiang C, Zeng W, Zhang H, et al. PCV2 regulates cellular inflammatory responses through dysregulating cellular miRNA-mRNA networks. Viruses. 2019; 11:1055

15.

Mucha SG, Ferrarini MG, Moraga C, Di Genova A, Guyon L, Tardy , et al. Mycoplasma hyopneumoniae J elicits an antioxidant response and decreases the expression of ciliary genes in infected swine epithelial cells. Sci Rep. 2020; 10:13707-29

16.

Fu S, Liu J, Xu J, Zuo S, Zhang Y, Guo L, et al. The effect of baicalin on microRNA expression profiles fo long in porcine aortic vascular endothelial cells infected by Haemophilus parasuis. Mol Cell Biochem. 2020; 472:45-56

17.

Guo L, Liu J, Zhang Y, Fu S, Qiu Y, Ye C, et al. The effect of baicalin on the expression profiles of long non-coding RNAs and mRNAs in porcine aortic vascular endothelial cells infected with Haemophilus parasuis. DNA Cell Biol. 2020; 39:801-15

18.

He J, Leng C, Pan J, Li A, Zhang H, Cong F, et al. Identification of lncRNAs involved in PCV2 infection of PK-15 cells. Pathogens. 2020; 9:479-90

19.

Wu J, Peng X, Qiao M, Zhao H, Li M, Liu G, et al. Genome-wide analysis of long noncoding RNA and mRNA profiles in PRRSV-infected porcine alveolar macrophages. Genomics. 2020; 112:1879-88

20.

Zhen Y, Wang F, Liang W, Liu J, Gao G, Wang Y, et al. Identification of differentially expressed non-coding RNA in porcine alveolar macrophages from tongcheng and large white pigs responded to PRRSV. Sci Rep. 2018; 8:15621

21.

Sen R, Ghosal S, Das S, Balti S, Chakrabarti J. Competing endogenous RNA: the key to posttranscriptional regulation. Sci World J. 2014; 2014:896206

22.

Langmead B. Aligning short sequencing reads with bowtie. Curr Protoc Bioinformatics. 2010; 32:11.7.1-14

23.

Kozomara A, Griffiths-Jones S. miRBase: annotating high confidence microRNAs using deep sequencing data. Nucleic Acids Res. 2014; 42:D68-73

24.

Friedländer MR, Mackowiak SD, Li N, Chen W, Rajewsky N. miRDeep2 accurately identifies known and hundreds of novel microRNA genes in seven animal clades. Nucleic Acids Res. 2012; 40:37-52

25.

Wen M, Shen Y, Shi S, Tang T. miREvo: an integrative microRNA evolutionary analysis platform for next-generation sequencing experiments. BMC Bioinform. 2012; 13:140

26.

Sun L, Luo H, Bu D, Zhao G, Yu K, Zhang C, et al. Utilizing sequence intrinsic composition to classify protein-coding and long non-coding transcripts. Nucleic Acids Res. 2013; 41e166

27.

Kang YJ, Yang DC, Kong L, Hou M, Meng YQ, Wei L, et al. CPC2: a fast and accurate coding potential calculator based on sequence intrinsic features. Nucleic Acids Res. 2017; 45:W12-6

28.

El-Gebali S, Mistry J, Bateman A, Eddy SR, Luciani A, Potter SC, et al. The Pfam protein families database in 2019. Nucleic Acids Res. 2019; 47:D427-32

29.

Lin MF, Jungreis I, Kellis M. PhyloCSF: a comparative genomics method to distinguish protein coding and non-coding regions. Bioinformatics. 2011; 27:i275-82

30.

Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014; 15:550

31.

Enright AJ, John B, Gaul U, Tuschl T, Sander C, Marks DS. MicroRNA targets in Drosophila. Genome Biol. 2003; 5:R1

32.

Krüger J, Rehmsmeier M. RNAhybrid: microRNA target prediction easy, fast and flexible. Nucleic Acids Res. 2006; 34:W451-4

33.

Wu J, Mao X, Cai T, Luo J, Wei L. KOBAS server: a web-based platform for automated annotation and pathway identification. Nucleic Acids Res. 2006; 34:W720-4

34.

Bindea G, Mlecnik B, Hackl H, Charoentong P, Tosolini M, Kirilovsky A, et al. ClueGO: a Cytoscape plug-in to decipher functionally grouped gene ontology and pathway annotation networks. Bioinformatics. 2009; 25:1091-3

35.

Zhang K, Ge L, Dong S, Liu Y, Wang D, Zhou C, et al. Global miRNA, lncRNA, and mRNA transcriptome profiling of endometrial epithelial cells reveals genes related to porcine reproductive failure caused by porcine reproductive and respiratory syndrome virus. Front Immunol. 2019; 10:1221

36.

Qin S, Wang H, Liu G, Mei H, Chen M. miR-21-5p ameliorates hyperoxic acute lung injury and decreases apoptosis of AEC II cells via PTEN/AKT signaling in rats. Mol Med Rep. 2019; 20:4953-62

37.

Li Q, Li B, Li Q, Wei S, He Z, Huang X, et al. Exosomal miR-21-5p derived from gastric cancer promotes peritoneal metastasis via mesothelial-to-mesenchymal transition. Cell Death Dis. 2018; 9:854

38.

Lou J, Wang Y, Zhang Z, Qiu W. MiR-20b inhibits mycobacterium tuberculosis induced inflammation in the lung of mice through targeting NLRP3. Exp Cell Res. 2017; 358:120-8

39.

Tian Y, Sun C, Zhang L, Pan Y. Clinical significance of miRNA-106a in non-small cell lung cancer patients who received cisplatin combined with gemcitabine chemotherapy. Cancer Biol Med. 2018; 15:157-64

40.

Li Q, Wu X, Guo L, Shi J, Li J. MicroRNA-7-5p induces cell growth inhibition, cell cycle arrest and apoptosis by targeting PAK2 in non-small cell lung cancer. FEBS Open Bio. 2019; 9:1983-93

41.

Liang Z, Wang L, Wu H, Singh D, Zhang X. Integrative analysis of microRNA and mRNA expression profiles in MARC-145 cells infected with PRRSV. Virus Genes. 2020; 56:610-20

42.

Mody HR, Hung SW, Pathak RK, Griffin J, Cruz-Monserrate Z, Govindarajan R. miR-202 diminishes TGFβ receptors and attenuates TGFβ1-induced emt in pancreatic cancer. Mol Cancer Res. 2017; 15:1029-39

43.

Chen J, Sun Q, Liu GZ, Zhang F, Liu CY, Yuan QM, et al. Effect of miR-202-5p-mediated ATG7 on autophagy and apoptosis of degenerative nucleus pulposus cells. Eur Rev Med Pharmacol Sci. 2020; 24:517-25

44.

Liu W, Jin Y, Zhang W, Xiang Y, Jia P, Yi M, et al. miR-202-5p inhibits RIG-I-dependent innate immune responses to RGNNV infection by targeting TRIM25 to mediate RIG-I ubiquitination. Viruses. 2020; 12:261

45.

Yao F, Sun L, Fang W, Wang H, Yao D, Cui R, et al. Hsa-miR-371-5p inhibits human mesangial cell proliferation and promotes apoptosis in lupus nephritis by directly targeting hypoxia-inducible factor 1α. Mol Med Rep. 2016; 14:5693-8

46.

Li N, Guo X, Liu L, Wang L, Cheng R. Molecular mechanism of miR-204 regulates proliferation, apoptosis and autophagy of cervical cancer cells by targeting ATF2. Artif Cells Nanomed Biotechnol. 2019; 47:2529-35

47.

Yan Y, Yu J, Liu H, Guo S, Zhang Y, Ye Y, et al. Construction of a long non-coding RNA-associated ceRNA network reveals potential prognostic lncRNA biomarkers in hepatocellular carcinoma. Pathol Res Pract. 2018; 214:2031-8

48.

Yang H, Zhang J, Zhang X, Shi J, Pan Y, Zhou R, et al. CD163 knockout pigs are fully resistant to highly pathogenic porcine reproductive and respiratory syndrome virus. Antiviral Res. 2018; 151:63-70

49.

Odaka C, Tanioka M, Itoh T. Matrix metalloproteinase-9 in macrophages induces thymic neovascularization following thymocyte apoptosis. J Immunol. 2005; 174:846-53

50.

Zhang J, Li G, Liu X, Wang Z, Liu W, Ye X. Influenza A virus M1 blocks the classical complement pathway through interacting with C1qA. J Gen Virol. 2009; 90:2751-8

51.

Luo J, Carrillo JA, Menendez KR, Tablante NL, Song J. Transcriptome analysis reveals an activation of major histocompatibility complex 1 and 2 pathways in chicken trachea immunized with infectious laryngotracheitis virus vaccine. Poult Sci. 2014; 93:848-55

52.

Li N, Zhu Q, Li Z, Han Q, Chen J, Lv Y, et al. IL21 and IL21R polymorphisms and their interactive effects on serum IL-21 and IgE levels in patients with chronic hepatitis B virus infection. Hum Immunol. 2013; 74:567-73

53.

Liu J, Jiang L, Liu JJ, He T, Cui Y, Qian F, et al. AEBP1 promotes epithelial-mesenchymal transition of gastric cancer cells by activating the NF-κB pathway and predicts poor outcome of the patients. Sci Rep. 2018; 8:11955

54.

Gerhard GS, Hanson A, Wilhelmsen D, Piras IS, Still CD, Chu X, et al. AEBP1 expression increases with severity of fibrosis in NASH and is regulated by glucose, palmitate, and miR-372-3p. PLOS ONE. 2019; 14e0219764

55.

Cui J, Li Y, Zhu L, Liu D, Songyang Z, Wang HY, et al. NLRP4 negatively regulates type I interferon signaling by targeting the kinase TBK1 for degradation via the ubiquitin ligase DTX4. Nat Immunol. 2012; 13:387-95

56.

Oulidi A, Bokhobza A, Gkika D, Vanden Abeele F, Lehen’kyi V, Ouafik L, et al. TRPV2 mediates adrenomedullin stimulation of prostate and urothelial cancer cell adhesion, migration and invasion. PLOS ONE. 2013; 8e64885