RNA-sequence analysis of primary alveolar macrophages after in vitro infection with porcine reproductive and respiratory syndrome virus strains of differing virulence.

Porcine reproductive and respiratory syndrome virus (PRRSV) mainly infects porcine alveolar macrophages (PAMs), resulting in porcine reproductive and respiratory syndrome (PRRS) in pigs. Most of the transcriptomic studies on PAMs infected with PRRSV conducted thus far have made use of microarray technology. Here, we investigated the transcriptome of PAMs in vitro at 12 h post-infection with two European PRRSV strains characterized by low (Lelystad, LV) and high (Lena) virulence through RNA-Seq. The expression levels of genes, isoforms, alternative transcription start sites (TSS) and differential promoter usage revealed a complex pattern of transcriptional and post-transcriptional gene regulation upon infection with the two strains. Gene ontology analysis confirmed that infection of PAMs with both the Lena and LV strains affected signaling pathways directly linked to the innate immune response, including interferon regulatory factors (IRF), RIG1-like receptors, TLRs and PKR pathways. The results confirmed that interferon signaling is crucial for transcriptional regulation during PAM infection. IFN-β1 and IFN-αω, but not IFN-α, were up-regulated following infection with either the LV or Lena strain. The down-regulation of canonical pathways, such as the interplay between the innate and adaptive immune responses, cell death and TLR3/TLR7 signaling, was observed for both strains, but Lena triggered a stronger down-regulation than LV. This analysis contributes to a better understanding of the interactions between PRRSV and PAMs and outlines the differences in the responses of PAMs to strains with different levels of virulence, which may lead to the development of new PRRSV control strategies.

Introduction

Porcine reproductive and respiratory syndrome (PRRS) is a viral disease that causes substantial economic losses to the swine industry [1], [2]. The PRRS virus (PRRSV) consists of a positive-sense single-stranded RNA molecule that is14.5 kb in length and has nine open reading frames. PRRSV strains can be divided into two groups, the European and North American genotypes, based on genetic and antigenic properties [3]. The European PRRSV genotype is divided into three subtypes: the mild pan-European subtype 1 (prototype virus: Lelystad virus, LV) and the more virulent eastern European subtypes 2 and 3 (prototype virus: Lena) [4]. The Lena virus has been found to be more virulent than the Lelystad virus [5]. Infection studies showed that the expression levels of IFN-α, IL-10, IL-12 and TNF-α were higher in Lena-infected pigs than in LV-infected pigs [6]–[8]. Differences in virulence between the PRRSV strains are related to their immunomodulatory properties and replication capacity [9], [10]. Nevertheless, Morgan et al. [11] found that the “SU1-bel, subtype 3”strain from Belarus exhibited high virulence as a consequence of an enhanced inflammatory response, rather than an increased replication capacity.

PRRSV induces a strong humoral response; however, this response is not effective at controlling the virus, and a persistent infection frequently develops [12], [13]. The genetic diversity of PRRSV strains has been extensively characterized, and the correlation between their genetic diversity and geographical/temporal distances has been investigated to reveal the mechanisms of viral spreading [14]–[16]. Many strategies for controlling PRRSV transmission have been proposed but have generally shown little success, which has stimulated the search for new ways to control PRRSV transmission, including the possibility of using genetic breeding to achieve resistance to PRRSV [17].

Pulmonary alveolar macrophages (PAMs) are the main target cells for PRRSV infection, and many gene expression studies have explored the immune response of PAMs to PRRSV. Such studies have shown that the expression levels of MX1, USP, IFN-β, IL-10 and TNF-α are affected by PRRSV infection [18], [19]. Overall, these analyses suggest that PRRSV subverts host defenses by inhibiting the expression of pro-inflammatory cytokines [20]–[22] and stimulating weak production of IFN-α [reviewed in 23].

RNA-Seq is increasingly being used to study gene expression, as it provides unbiased profiles, compared to microarrays [24], [25], and can be extremely accurate if a sufficient level of coverage is obtained [26]. Validation techniques, such as qPCR [27] and spike-in RNA [28], have corroborated the accuracy of RNA-Seq. One advantage of RNA-Seq over microarrays is that, in addition to gene expression, RNA-Seq data can be used to identify transcription start sites (TSS), splicing variants and differential promoter usage [28]. Annotations of the human genome suggest that a high proportion of protein-coding genes exhibit alternative promoters, and aberrant promoter usage has been found to be associated with various diseases [29].

In this study, RNA-Seq was used to characterize gene expression, splice variants, TSSs and differential promoter usage in PAMs infected with either LV (subtype 1, low virulence) or Lena (subtype 3, high virulence). The PAM signaling pathways that were affected by PRRSV infection were investigated and were found to differ following infection with the two PRRSV strains. The obtained data confirmed that many signaling pathways are potentially involved in the host immune response to PRRSV infection and that the regulation of interferon signaling is strongly modulated during PAM infection. In contrast to infection with LV, infection with Lena down-regulated the pathways involved in the interplay between the innate and adaptive immune response, cell death and TLR3/TLR7 signaling. Additionally, we also highlight many differences in splicing isoforms, TSSs and differential promoter usage between LV and Lena infections.

Materials and Methods

Animals, cells and library preparation

Pulmonary alveolar macrophages were collected from three 3-week-old piglets, which were the offspring of hybrid sows (JSR Genepacker 90 English: Landrace x Large White) and Pietrain boars from a PRRS-negative herd. The PRRSV- and PCV2-negative status of the piglets was confirmed via an immunoperoxidase monolayer assay (IPMA) and PCR, respectively. The piglets were treated daily with 1 ml of enrofloxacin (5% solution) and 1 ml of lincospectin/spectinomycin (5 or 10% solution) for 3 days to eliminate bacterial pathogens and then were sacrificed one week later via intravenous injection with a lethal dose of Na-pentobarbital (20 ml/animal) in the jugular vein. The sacrifice procedure was carried out by a technician certified in Laboratory Animal Science (Felasa Category C). The PAMs were collected by broncho-alveolar lavage and frozen in liquid nitrogen as described by Wensvoort et al. [30].

The animals were housed at the Faculty of Veterinary Medicine of Ghent University (Belgium). All experimental procedures and the sacrifice of the animals were conducted in compliance with the relevant legislation on animal experiments (EU Directive 2010/63) and with the approval of the local ethical committee.

Prior to infection, the PAMs were thawed and cultured for 48 h, as described by Delputte et al. [31]. The primary culture from each of the three animals was split into three fractions: one fraction was infected at a multiplicity of infection (MOI) of 0.5 with the LV strain; the second was inoculated at a MOI of 0.5 with the Lena strain; and the third was maintained as a control (mock inoculated). The cells were collected at 12 h post-infection for transcriptome analysis. All of the infection experiments and controls were performed in duplicate. One of each duplicate was used to assess the percentage of infected cells. After 12 h, the cells were fixed in acetone-100% methanol at −20°C, and the infected cells were detected via immunoperoxidase staining with the monoclonal anti-nucleocapsid protein antibody P3/27 (48), followed by incubation with a horseradish peroxidase-labeled goat anti-mouse secondary antibody and development with a substrate solution containing 3-amino-9-ethylcarbazole. Viral antigen-positive cells and total cells were counted using an Olympus light microscope (Olympus Optical Co., Hamburg, Germany), and the percentage of infected cells was calculated. Three microscopic fields and a minimum of 200 cells per field were counted for each experimental condition.

Total RNA was extracted from the PAMs using TRIzol (Invitrogen Life Technologies, Milan-Italy) and RNeasy columns (Qiagen). RNA quality was assessed via microcapillary electrophoresis on an Agilent 2001 Bioanalyzer (Agilent Technologies) with RNA 6000 Nanochips. RNA was quantified using spectrophotometry (ND-1000; NanoDrop Technologies).

cDNA libraries were constructed using RNAfrom the LV-, Lena- or mock-infected PAMs. A total of nine cDNA libraries were obtained using TruSeq Sample Prep Kits (Illumina Inc., San Diego, CA) and were sequenced via 2×100 paired-end sequencing on an IlluminaHiSeq 2000 instrument. The obtained sequences were submitted to the Sequence Read Archive (SRA) at NCBI (SRX352447).

The Laboratory of Virology has been acknowledged by the Belgian Federal Agency for the Safety of the Food Chain (FASFC) for its work with experimental animals, including pigs (national approval number LA1400076). As euthanasia was the only action performed on the animals, no approval was required by the local ethical committee (Article 1, 5f, EU Directive 2010/63).

Quality check, mapping, assembly and visualization of the sequenced reads

Sequence visualization and statistical analyses were performed with FASTQC (http://www.bioinformatics. bbsrc.ac.uk/projects/fastqc). Sickle software (https://github.com/najoshi/sickle) was employed to retrieve the paired sequences with quality scores above a threshold of 20 and a length greater than 50 bp. The reads passing the quality criteria were processed and aligned to the reference pig genome (Build Sus_scrofa.Sscrofa10.2.71.) using TopHat v2.0.8 [32], and a maximum of 2 mismatches were allowed when mapping the reads to the reference genome. The default parameters for TopHat were used [32]. The mapping results were then used to identify “islands” of expression, which can be interpreted as potential exons. The mapping results were visualized using the UCSC genome browser [33] and the Integrative Genomics Viewer (IGV) software, which is available at http://www.broadinstitute.org/igv/.

Transcript assembly and abundance estimation

Transcript assembly and abundance estimation were performed using two methods: one for the genes and the other for the remaining features (isoforms, TSS and promoters).

The aligned reads were further processed using Cufflinks v2.1.1 [34] (Figure S4-A) to infer the splicing structure of each gene using the optimal reconstruction of the gene model, which was determined using a rigorous statistical model. The reads were assembled into transcripts by forcing Cufflinks to use not just the existing gene annotations during the assembly of transcripts but instead to construct the transcripts based on both the reference gene annotations and the sequence data. This approach allows Cufflinks to identify new potential alternative transcription and splicing events. Confidence intervals for the FPKM estimates were calculated using a Bayesian inference method [35]. Once all of the short read sequences were assembled with Cufflinks, the output files were merged via the Cuffmerge function and processed using Cuffcompare, along with a reference GTF annotation file downloaded from the Ensembl database.

The HTseq python package (http://www-huber.embl.de/users/anders/HTSeq/doc/count.html) was employed to generate an unambiguous table of counts per gene to estimate expression levels. The HTseq-count function was employed, with ‘union’ in overlapping mode and ‘gene’ as a feature.

Differential expression analysis

The cells obtained from each pig were divided into three batches to represent the two treatment conditions (LV-infected orLena-infected) and the control (mock-infected). To remove baseline differences between the pigs, a paired design was employed, in which treatment values were subtracted from control values for each pig prior to the computation of the average and variance. Differential expression analysis was performed using the EdgeR package [36], which employs the table of counts per gene generated by HTSeq as an input.

Using edgeR, a general linear model was generated to test for a treatment (infection with Lena or LV) effect. This method is based on the negative binomial distribution and uses the conditional weighted likelihood to moderate the level of overdispersion across the sequences. Three pairwise comparisons of gene expression levels were performed for the three combinations of treatments: i) LV vs. mock, ii) Lena vs. mock and iii) LV vs. Lena.

The Cuffdiff function was used to assess the presence and frequency of different isoforms, TSSs and differential promoter usage. Briefly, Cuffdiff starts by modeling the variability in the fragment count for each gene across the replicates corresponding to the mock, LV and Lena samples. Second, the counts for each isoform are estimated in each replicate, taking into consideration the uncertainty arising from ambiguously mapped reads. Finally, Cuffdiff estimates the count variances for each transcript in each sample, and this information is used in statistical testing to report the differentially expressed transcripts, including the isoforms, TSSs and differential promoter usage [34]. Cuffdiff was used to perform three pairwise comparisons of the observed isoforms, TSSs and promoter usage: i) LV vs. mock ii) Lena vs. mock and iii) LV vs. Lena. For all of the analyses, a feature was considered significant if the false discovery rate (FDR) was less than 0.05 and the fold change (FC) difference was greater than or equal to 1.5.

Functional analysis of gene lists through Ingenuity Pathways Analysis (IPA)

IPA (http://www.ingenuity.com) was used to identify the most significant molecular networks, biological functions and canonical metabolic pathways. The p-values associated with the biological processes or pathway annotation were calculated according to the right-tailed Fisher's Exact Test. Affected genes were first matched to their corresponding gene names in IPA. The obtained lists were then subjected to IPA to reveal the canonical pathways and biological functions that were significantly associated with the gene lists, using the default “universe” of genes and endogenous chemicals in the IPA library.

Results

1. Expressed and differentially expressed genes, isoforms and TSSs, as well as differential promoter usage,thein LV vs. mock, Lena vs. mock and LV vs. Lena analyses

The percentages of PAMs infected by LV and Lena were 21% and 16% on average, respectively. This result was confirmed using quantitative PCR; the threshold cycle (Ct) values obtained were 15.6 and 20.5 for LV and Lena, respectively.

Nine libraries obtained from PAMs corresponding to LV (3 samples), Lena (3 samples) and mock infections (3 samples) were sequenced. The mean numbers of reads produced per strand were 14,483,298, 17,381,381 and 10,473,546 for the LV, Lena and Mock groups, respectively (Table S1-A). Data filtering using FASTQC and Sickle excluded between 11% and 17% of the reads (Table S1-B). Between 73% and 87% of the total paired reads passing the quality test were mapped to the reference pig genome (Build Sus_scrofa.Sscrofa10.2.71).

In total, 29,900 annotated genes and 58,347 isoforms were expressed in both the infected and non-infected PAMs. Of these genes, 0.46% were expressed at more than 10,000 FPKM, 0.6% between 1000 and 10,000 FPKM, 3.14% between 100 and 1000 FPKM and 40.64% between 1 and 100 FPKM. The genes with the highest expression were XLOC_009444 (14,122,900 FPKM in the LV-infected cells), XLOC_027735 (13,467,000FPKM in the Lena-infected cells) and SCARNA16 (13,133,300 FPKM in the LV-infected cells). The differentially expressed features and the intersection of the differentially expressed genes, isoforms, TSSs and differential promoter usage observed in the LV vs. mock, Lena vs. mock and LV vs. Lena comparisons are shown in Figure 1 (1.I, 1.II, 1.III and 1.IV).

Figure 1

Venn diagram illustrating the significantly affected genes (I), isoforms (II), TSSs (III) and differential promoter use (VI) found in the LV vs. mock (A), Lena vs. mock (B) and LV vs. Lena (C) comparisons.

I. Venn diagram illustrating the significantly affected genes found in the LV vs. mock (A, 446), Lena vs. mock (B, 153) and LV vs. Lena (C, 241) comparisons. II. Venn diagram illustrating the significantly affected isoforms found in the LV vs. mock (A, 187), Lena vs. mock (B, 72) and LV vs. Lena (C, 34) comparisons. III. Venn diagram illustrating the significantly affected TSSs found in the LV vs. mock (A, 240), Lena vs. mock (B, 93) and LV vs. Lena (C, 50) comparisons. VI. Venn diagram illustrating the significantly affected promoters found in the LV vs. mock (A, 14), Lena vs. mock (B, 6) and LV vs. Lena (C, 7) comparisons. The blue circle, “A”, corresponds to the LV vs. mock case (I. significantly affected genes, II. significantly affected isoforms, III. significantly affected TSSs and VI. significantly affected promoters). The green circle, “B”, corresponds to the Lena vs. mock case (I. significantly affected genes, II. significantly affected isoforms, III. significantly affected TSSs and VI. significantly affected promoters). The purple circle, “C”, corresponds to the LV vs. Lena case (I. significantly affected genes, II. significantly affected isoforms, III. significantly affected TSSs and VI. significantly affected promoters).

1.1. Canonical pathways and biological functions differentially expressed in the LV- vs. mock-infected samples

A total of 446 differentially expressed genes were identified between the mock-infected and LV-infected PAMs (Table S2). Among these genes, 107 presented at least two isoforms that were differentially expressed between the mock- and LV-infected PAMs, and one was differentially regulated at the promoter level (Figure 2.I). In total, 413 genes could be mapped to the IPA database to assign gene ontology. The top five canonical pathways identified as being differentially expressed in the LV vs. mock via IPA were “Interferon Signaling”, “Activation of IRF by Cytosolic Pattern Recognition Receptors”, the “Role of Pattern Recognition Receptors in Recognition of Bacteria and Viruses”, the “Role of PKR in Interferon Induction and Antiviral Response” and “Retinoic Acid Mediated Apoptosis Signaling” (Table 1A). All of these pathways are interconnected and are directly involved in the immune response to viral infection. In the infected cells, the “Interferon Signaling” pathway, which is central to the innate defense mechanisms of the host against viral infection, included 12 up-regulated genes: IFIT1, IFIT3, IFITM1, IFNβ1, IRF1, JAK2, MX1, OAS1, PSMB8, SOCS1, STAT1 and STAT2. An additional 14 up-regulated genes were involved in the activation of IRF by the cytosolic pattern recognition receptor pathway (ADAR, DDX58, DHX58, IFIH1, IFIT2, IFNβ1, IL10, IRF7, ISG15, NFKBIA, STAT1, STAT2, TNF and ZBP1). The “Role of Pattern Recognition Receptors in Recognition of Bacteria and Viruses” pathway, which generally leads to the serine phosphorylation of IRF3 and IRF7 [37], involved IRF7, but not IRF3. The “Role of PKR in Interferon Induction and Antiviral Response” and “Retinoic Acid Mediated Apoptosis Signaling” pathways were up-regulated and shared four genes (CASP8, BID, IFNβ1 and IRF1).

Figure 2

Venn diagram illustrating the significantly affected genes (A, chartreuse color), isoforms (B, violet color) and differential promoter use (C, orange color) found in the LV vs. mock (I), Lena vs. mock (II) and LV vs. Lena (III) comparisons.

I. Venn diagram illustrating the significantly affected genes (A, 446), isoforms (B, 187) and differential promoter use (C, 14) found in the LV vs. mock comparison. II. Venn diagram illustrating the significantly affected genes (A, 153), isoforms (B, 72) and promoters (C, 6) found in Lena vs. mock comparison. III. Venn diagram illustrating the significantly affected genes (A, 241), isoforms (B, 34) and promoters (C, 7) found in LV vs. Lena comparison. The chartreuse circle, “A”, corresponds to the differentially expressed genes (I. LV vs. mock, II. Lena vs. mock and III. LV vs. Lena). The violet circle, “B”, represents the differentially expressed isoforms (I. LV vs. mock, II. Lena vs. mock and III. LV vs. Lena). The orange circle, “C”, corresponds to the differentially used promoters(I. LV vs. mock, II. Lena vs. mock and III. LV vs. Lena).

Table 1

Lists of the 5 top canonical pathways and the corresponding affected genes in the LV vs. mock (A), Lena vs. mock (B) and LV vs. Lena (C) comparisons.

A. LV vs. mock infection
Ingenuity Canonical Pathways	−log(p-value)*	Ratio£	Molecules&
Interferon Signaling	1.08E+01	3.87E-01	IFIT3,SOCS1,IFIT1,OAS1,MX1,IFNB1, IFITM1,STAT2,PSMB8,JAK2,STAT1,IRF1
Activation of IRF by Cytosolic Pattern Recognition Receptors	9.97E+00	2.64E-01	DHX58,IL10,ZBP1,IFNB1,ADAR,ISG15, IFIH1,IRF7,NFKBIA,DDX58,STAT2,IFIT2, STAT1,TNF
Role of Pattern Recognition Receptors in Recognition of Bacteria and Viruses	8.55E+00	1.78E-01	OAS1,OAS2,IL10,PIK3R1,IFNB1,RNASEL,IFIH1,IRF7,DDX58,TLR7,CASP1,IL1B,PIK3CD, EIF2AK2,TLR3,TNF
Role of PKR in Interferon Induction and Antiviral Response	7.03E+00	2.5E-01	NFKBIA,IFNB1,BID,EIF2AK2,TLR3,STAT1, CASP8,RNASEL,TNF,IRF1
Retinoic Acid Mediated Apoptosis Signaling	6.96E+00	2.16E-01	BID,CASP8,IFNB1,IRF1,PARP9,PARP11, PARP14,TIPARP,TNFRSF10A,TNFSF10, ZC3HAV1

*probability that a canonical pathway is significantly involved in the analysis.

number of genes in the list involved in a certain canonical pathway divided by the number of genes in the IPA database involved in the same canonical pathway.

genes involved in the corresponding canonical pathway.

The IFN-β gene was common to all of the top canonical pathways, which is consistent with the previously reported fundamental role of this gene in the response of pigs to PRRSV (reviewed in [23]). Interestingly, IFN-β and IFN-αω were both up-regulated following LV and Lena infection, but the expression of IFN-β was higher in the PAMs infected with LV (Figure 3).

Figure 3

Gene expression of IFN-β and IFN-αω genes in FPKM (Y axis).

Top left: IFN-β expression in the mock, LV and Lena groups. IFN-β was differentially expressed between both the LV and Lena groups when compared to the mock infection group. IFN-β was also expressed at a significantly higher level in the LV group than in the Lena group. Lower right: IFN-αω expression in the mock, LV and Lena groups. IFN-αω was differentially expressed between the LV and Lena groups from one side and the mock group from the other side. IFN-β was also expressed at a significantly higher level in the LV group than in the Lena group.

Among the top ten biological functions that differed between the LV-and mock-infected samples (Table 2A), the categories “Cell death of immune cells”, “Leukocyte infection”, other functions related to the hematological system (the numbers of leukocytes, mononuclear leukocytes, lymphocytes and T lymphocytes) and the “Cell death/survival” pathway were all up-regulated. Only two biological functions (apoptosis of phagocytes and cellular homeostasis) were down-regulated.

Table 2

Lists of affected biological functions and the corresponding affected genes in the LV vs. mock (A) and Lena vs. mock (B) comparisons.

A. LV vs. mock
Function Annotation	p-Value*	Predicted Activation State£	Activation z-score&	Molecules§
A. LV vs. mock
cell death of immune cells	1.89E-13	Increased	2.147	BID,CASP1, CASP8, CCL4, CCR5, CD274, CD38, CD5, DHX58, EIF2AK2, HAVCR2, HCAR2, HSH2D, IDO1, IFNB1, IL10, IL15, IL15RA, IL1B, IL27, IL7R, IRF1, IRF2, IRF4, JAK2, LGALS9, MCL1, MX1, NAMPT, NFKBIA, PARP14, PIK3CD, PIK3R1, PIM1, PMAIP1, PPBP, PRDM1, PYCARD, RHOH, RNASEL, SEMA4D, SOCS1, STAT1, TLR3, TNF, TNFAIP3, TNFSF10, TNFSF13B
quantity of leukocytes	2.62E-13	Increased	2.297	ABCA1, ARNTL, B3GNT5, BID, CASP8, CCL2, CCL3L1/CCL3L3, CCR2, CCR5, CD274, CD38, CD5, CXCL10, CYP27B1, DAPP1, DDX58, ELF1, ENTPD7, HAVCR2, HBEGF, HLA-B, IDO1, IFNB1, IL10, IL15, IL15RA, IL1B, IL27, IL7R, IRF1, IRF2, IRF4, ISG15, JAK2, LGALS9, MAEA, MAP3K8, MCL1, MXD1, NBN, NFKBIA, PDE4B, PIK3CD, PIK3R1, PIM1, PIM3, PLAC8, PML, PPM1D, PRDM1, PTPN13, PYCARD, RFX5, RHOH, RICTOR, SEMA4D, SOCS1, SOX4, STAT1, TCF4, TEC, TIMD4, TLR3, TNF, TNFAIP3, TNFRSF10A, TNFSF10, TNFSF13B, VCAM1
Cell death	2.62E-13	Increased	2.297	ABCA1, ARNTL, B3GNT5, BID, CASP8, CCL2, CCL3L1/CCL3L3, CCR2, CCR5, CD274, CD38, CD5, CXCL10, CYP27B1, DAPP1, DDX58, ELF1, ENTPD7, HAVCR2, HBEGF, HLA-B, IDO1, IFNB1, IL10, IL15, IL15RA, IL1B, IL27, IL7R, IRF1, IRF2, IRF4, ISG15, JAK2, LGALS9, MAEA, MAP3K8, MCL1, MXD1, NBN, NFKBIA, PDE4B, PIK3CD, PIK3R1, PIM1, PIM3, PLAC8, PML, PPM1D, PRDM1, PTPN13, PYCARD, RFX5, RHOH, RICTOR, SEMA4D, SOCS1, SOX4, STAT1, TCF4, TEC, TIMD4, TLR3, TNF, TNFAIP3, TNFRSF10A, TNFSF10, TNFSF13B, VCAM1
quantity of lymphocytes	1.61E-12	Increased	2.703	BID, CASP1, CASP8, CCL2, CCL4, CCR5, CD274, CD38, CD5, DHX58, EIF2AK2, HAVCR2, HCAR2, HLA-B, HSH2D, IDO1, IFNB1, IL10, IL15, IL15RA, IL1B, IL27, IL7R, IRF1, IRF2, IRF4, JAK2, LGALS9, MCL1, MX1, NAMPT, NFKBIA, PARP14, PIK3CD, PIK3R1, PIM1, PMAIP1, PPBP, PRDM1, PYCARD, RHOH, RNASEL, SEMA4D, SOCS1, STAT1, TLR3, TNF, TNFAIP3, TNFSF10, TNFSF13B
quantity of mononuclear leukocytes	5.16E-12	Increased	3.212	ARNTL, B3GNT5, BID, CASP8, CCL2, CCR2, CCR5, CD274, CD38, CD5, CXCL10, CYP27B1, DAPP1, DDX58, ELF1, ENTPD7, HLA-B, IFNB1, IL10, IL15, IL15RA, IL27, IL7R, IRF1, IRF2, IRF4, ISG15, JAK2, LGALS9, MAP3K8, MCL1, NBN, NFKBIA, PIK3CD, PIK3R1, PIM1, PIM3, PPM1D, PRDM1, PTPN13, PYCARD, RFX5, RHOH, SEMA4D, SOCS1, SOX4, STAT1, TCF4, TEC, TIMD4, TLR3, TNF, TNFAIP3, TNFRSF10A, TNFSF13B, VCAM1
infection of leukocytes	5.16E-12	Increased	3.212	ARNTL, B3GNT5, BID, CASP8, CCL2, CCR2, CCR5, CD274, CD38, CD5, CXCL10, CYP27B1, DAPP1, DDX58, ELF1, ENTPD7, HLA-B, IFNB1, IL10, IL15, IL15RA, IL27, IL7R, IRF1, IRF2, IRF4, ISG15, JAK2, LGALS9, MAP3K8, MCL1, NBN, NFKBIA, PIK3CD, PIK3R1, PIM1, PIM3, PPM1D, PRDM1, PTPN13, PYCARD, RFX5, RHOH, SEMA4D, SOCS1, SOX4, STAT1, TCF4, TEC, TIMD4, TLR3, TNF, TNFAIP3, TNFRSF10A, TNFSF13B, VCAM1
proliferation of B lymphocytes	7.80E-12	Increased	3.218	ARNTL, B3GNT5, BID, CASP8, CCL2, CCR2, CCR5, CD274, CD38, CD5, CXCL10, CYP27B1, DAPP1, DDX58, ELF1, ENTPD7, HLA-B, IFNB1, IL10, IL15, IL15RA, IL27, IL7R, IRF1, IRF2, IRF4, ISG15, JAK2, LGALS9, MAP3K8, MCL1, NBN, NFKBIA, PIK3CD, PIK3R1, PIM1, PIM3, PML, PPM1D, PRDM1, PTPN13, PYCARD, RFX5, RHOH, SEMA4D, SOCS1, SOX4, STAT1, TCF4, TEC, TIMD4, TLR3, TNF, TNFAIP3, TNFRSF10A, TNFSF13B, VCAM1
quantity of T lymphocytes	7.80E-12	Increased	3.218	ARNTL, B3GNT5, BID, CASP8, CCL2, CCR2, CCR5, CD274, CD38, CD5, CXCL10, CYP27B1, DAPP1, DDX58, ELF1, ENTPD7, HLA-B, IFNB1, IL10, IL15, IL15RA, IL27, IL7R, IRF1, IRF2, IRF4, ISG15, JAK2, LGALS9, MAP3K8, MCL1, NBN, NFKBIA, PIK3CD, PIK3R1, PIM1, PIM3, PML, PPM1D, PRDM1, PTPN13, PYCARD, RFX5, RHOH, SEMA4D, SOCS1, SOX4, STAT1, TCF4, TEC, TIMD4, TLR3, TNF, TNFAIP3, TNFRSF10A, TNFSF13B, VCAM1
apoptosis of phagocytes	1.47E-10	Decreased	−2.201	CASP1, CCL4, CCR5, EIF2AK2, HAVCR2, IL15, IRF1, LGALS9, MGLL, PIK3R1, SAMHD1, STAT1, TNF
Cellular homeostasis	1.47E-10	Decreased	−2.201	CASP1, CCL4, CCR5, EIF2AK2, HAVCR2, IL15, IRF1, LGALS9, MGLL, PIK3R1, SAMHD1, STAT1, TNF
B. Lena vs. mock
cell death of immune cells	4.78E-08	Increased	2.276	CCL4, CCR5, CD274, DHX58, EIF2AK2, HAVCR2, HCAR2, IFNB1, IL15RA, IL27, IL7R, IL8, MX1, PARP14, PIK3CD, PPBP, RNASEL, SOCS1, TNF, TNFSF10
Cell death	8.23E-08	Increased	2.575	CCL2, CCL4, CCR5, CD274, DHX58, EIF2AK2, HAVCR2, HCAR2, IFNB1, IL15RA, IL27, IL7R, IL8, MX1, PARP14, PIK3CD, PPBP, RNASEL, SOCS1, TNF, TNFSF10
Cell death of lymphocytes	6.11E-07	Increased	2.921	CCL2, CCL4, CCL8, CCR5, CCRL2, CXCL10, EDN1, IFNB1, IL8, PDE4B, PIK3CD, PPBP, TNF
apoptosis of lymphocytes	6.11E-07	Increased	2.921	CCL2, CCL4, CCL8, CCR5, CCRL2, CXCL10, EDN1, IFNB1, IL8, PDE4B, PIK3CD, PPBP, TNF
apoptosis of T lymphocytes	6.11E-07	Increased	2.921	CCL2, CCL4, CCL8, CCR5, CCRL2, CXCL10, EDN1, IFNB1, IL8, PDE4B, PIK3CD, PPBP, TNF
invasion of phagocytes	6.11E-07	Increased	2.921	CCL2, CCL4, CCL8, CCR5, CCRL2, CXCL10, EDN1, IFNB1, IL8, PDE4B, PIK3CD, PPBP, TNF
chemotaxis of antigen presenting cells	2.43E-06	Increased	2.987	ADAR, CCL2, CD274, DDX58, DUSP5, EIF2AK2, IFNB1, IL15RA, IL27, IL7R, IL8, IRF7, JAG1, PIK3CD, PML, RSAD2, SEMA4A, SOCS1, TNF
hematopoiesis of myeloid progenitor cells	2.43E-06	Increased	2.987	ADAR, CCL2, CD274, DDX58, DUSP5, EIF2AK2, IFNB1, IL15RA, IL27, IL7R, IL8, IRF7, JAG1, PIK3CD, PML, RSAD2, SEMA4A, SOCS1, TNF
cell movement of PBMCs	2.43E-06	Increased	2.987	ADAR, CCL2, CD274, DDX58, DUSP5, EIF2AK2, IFNB1, IL15RA, IL27, IL7R, IL8, IRF7, JAG1, PIK3CD, PML, RSAD2, SEMA4A, SOCS1, TNF
inflammatory response	2.53E-06	Increased	2.635	CCL2, CCL4, CCL8, CCR5, CXCL10, IFNB1, IL8, PIK3CD

*probability that a biological function is significantly involved in the analysis.

activation state of the biological function; three cases are possible: activation, inhibition or “none”, if the statistical significance is not sufficient for a clear conclusion regarding the activation/inhibition.

prediction of the activation state of the biological function; IPA uses z-scores to evaluate the activation/inhibition of the biological function. A negative score is indicative of inhibition, while a positive score is indicative of activation. If the statistical significance is not sufficient for a clear conclusion regarding activation/inhibition, we discuss the tendency toward activation/inhibition on the basis of the positive/negative activation z-score.

genes involved in the corresponding biological function.

1.2. Canonical pathways and biological functions differentially expressed in the Lena- vs. mock-infected samples

Three of the top five canonical pathways modulated by the infection of PAMs with Lena were shared with those affected by LV infection: “Activation of IRF by Cytosolic Pattern Recognition Receptors”, “Interferon Signaling” and the “Role of Pattern Recognition Receptors in Recognition of Bacteria and Viruses” (Table 1B). The two canonical pathways specific to Lena infection were the “Role of RIG1-like Receptors in Antiviral Innate Immunity”, which is involved in the recognition of viral RNA [38]–[40], and the “Role of Hypercytokinemia/hyperchemokinemia in the Pathogenesis of Influenza”.

Among the top ten biological functions, the “cell death of immune cells and lymphocytes and apoptosis of lymphocytes and T lymphocytes” was up-regulated (Table 2B). Other up-regulated functions related to virus infection included the invasion of phagocytes, chemotaxis of antigen-presenting cells and the inflammatory response. Cell death/survival pathways were less affected by Lena infection than by LV infection and involved different sets of genes (Figure S1-A and B).

1.3. Canonical pathways and biological functions differentially expressed in the LV-and Lena-infected samples

A direct comparison between the LV and Lena infections was performed at the gene, isoform and promoter levels (Table S2). The top five canonical pathways (not taking into account the canonical pathways already found to be differentially expressed in the LV vs. mock and Lena vs. mock) were “IL-15 production”, “TREM1 signaling”, “Communication between innate and adaptive Immune Cells”, “Crosstalk between dendritic cells and natural killer cells” and the “Role of JAK2 in hormone-like cytokine signaling” (Table 1C). Lena infection triggered the up-regulation of RPL31 (FC = 243.87), IL-8 (FC = 144), PA2G4 (FC = 75) and CXCL2 (FC = 9.18), while other genes were up-regulated following LV infection, including P2RY2 (FC = 29.24), CCL2 (FC = 13.18), WNT4 (FC = 11.31), ACADVL (FC = 7.26), CXCL10(FC = 3.97) and IFNβ1 (FC = 3.48).

Interestingly, RPL31 and IL8 were almost exclusively expressed after Lena infection (Figure S2-A and B), while other genes, including CCL2 and ACADVL, were differentially regulated by infection with both strains, although to different extents (Figure S2-C and D). CXCL10 exhibited three isoforms, but only one isoform (TCONS_00046382) was differentially expressed (Figure S3-A) between the two strains. This gene was also regulated at the TSS level; we identified 3 TSSs, of which one (TSS31082) was differentially regulated between the LV and Lena infections (Figure S3-B).

2. Transcriptional and post-transcriptional regulation and differential promoter usage in PAMs after infection with the LV and Lena PRRSV strains

Differentially expressed isoforms with different TSSs are transcriptionally regulated, while differentially expressed isoforms with the same TSSs are regulated at the post-transcriptional level [28]. For the genes involved in the top canonical pathways identified in the previous analyses (LV vs. mock, Lena vs. mock and LV vs. Lena), the regulation of PRRSV infection was investigated at the levels of the transcripts, isoforms and TSSs (Figure S4-B). Three groups of genes were defined:

Genes with one isoform and one TSS; these genes were classified as “un-spliced and transcriptionally regulated” genes and included TNF, BID, CASP8, ZC3HAV1, IL8, CCR5, CCL4, CCL2 and TLR4 (Figure S5-Group 1).

Genes with more than one isoform and one TSS;these genes were classified as “Spliced and post-transcriptionally regulated” genes and included IFNβ1, STAT1, TIPARP, MICB, IRF7, CASP1, IL1B, OAS2 and STAT2 (Figure S5-Group 2).

Genes with more than one isoform and more than one TSS. In this group, some isoforms of a given gene are transcriptionally regulated, while others are post-transcriptionally regulated. This group was classified as “Spliced and both transcriptionally and post-transcriptionally regulated”. This group could be further differentiated into two sub-groups (Figure S5-Group 3): (a) “Spliced genes whose isoforms are all transcriptionally regulated”, which included IFIT3, IFIT1, IRF1, ZBP1, ISG15, OAS1, IFITM1, TLR7, TNFSF10, TNFRSF10A and CXCL10; and (b) “Spliced genes whose isoforms are either transcriptionally or post-transcriptionally regulated”, which included IL10, NFKBIA, EIF2AK2, PARP14, PSMB8, DHX58, ADAR and RNASEL.

The analysis of differential promoter usage offers an additional perspective to the overview of gene regulation triggered by PRRSV infection. Therefore, we considered the different isoforms that showed the same TSS yields and presented transcripts produced from different start sites [Figure S4-B]. This analysis suggested that 14, 6 and 7 genes were differentially regulated at the promoter level in the LV vs. mock, Lena vs. mock and LV vs. Lena comparisons, respectively (Table S2). Two genes (PSMB8 and TIPARP7) expressed two isoforms, of which one was exclusively expressed in one PRRSV strains (Lena or LV) but not in the other (Figure S5-Group 3 and 2). Moreover, many genes exhibited different levels of expression of the various isoforms under infection with one PRRSV strain vs. the other (PARP14, SOCS1, IFIT1, MX1, IFOTM1, RNASEL, PIKSCD and TLR3; Figure S5).

Discussion

This study provides insights into the transcriptome of PAMs infected with the LV and Lena PRRSV strains.

Canonical pathways and biological functions differentially expressed in the LV vs. mock, Lena vs. mock and LV vs. Lena analyses

The biological annotation of the immune response-related genes that were differentially expressed in the PAMs during LV and Lena infection can be interpreted using two sequential steps.

i) The role of the TLRs in mounting an immune response during PRRSV infection has been confirmed [41], [42]. Once activated, IRF1, IRF2 and IRF7 translocate into the nucleus and induce the transcription of IFN-β and IFN-αω, but not IFN-α. IRF1, IRF2 and IRF7 were all affected by LV infection, but only IRF7 was found to be modulated after Lena infection. The fact that the PRRSV strains provoked the expression of IFN-β, but not IFN-α, is in agreement with previous experiments performed in PAMs and dendritic cells in vitro [43]–[45].

(ii) Once secreted, IFN-β and most likely IFN-αω are exported to the extracellular space, where they bind to the IFN receptors, provoking the recruitment of STAT1/STAT2 and inducing the subsequent IFN-stimulated genes: IFIT1, IFIT3, IFITM1, ISG12, ISG15, ISG20, PKR, IRF1, JAK2, MX1, OAS1, PSMB8, RNaseL and SOCS1;all of these genes play a fundamental role in the host immune response. Specifically, PKR dimerizes and phosphorylates EIF2α to inhibit the translational process [45], [46]. OAS1 is an IFN inducer that is capable of inhibiting protein synthesis and viral growth by degrading viral and cellular RNA [47]. ISG12, ISG15 and ISG20 are ubiquitin homologs involved in the regulation of innate immunity [48]. MX1 exhibits antiviral activity [49].

The biological functions altered by Lena and LV infections were focused on cell death/survival, which is considered an innate defense mechanism. We observed considerable interplay between cell death and cell survival, which might explain the discrepancies reported in previous PRRS studies concerning apoptosis [50]–[52].

Comparison between Lena and LV infections

In this study, the replication rate was slightly higher in LV than Lena strain. Similar results were obtained by Morgan et al. [11] who reported no obvious differences in PRRSV replication between the LV and the “SU1-bel, subtype 3” strain in vitro. This differs from results of in vivo studies, in which a higher replication rate was reported for Lena compared to other PRRSV subtype 1 strains [6]. Frydas et al [53] have shown that, while LV can replicate in a restricted subpopulation of monocytic cells, Lena appears to target an ample range of cell subpopulations to disseminate within the mucosa. This finding might explain the high replication rate of the Lena strain in vivo, as airway mucosal surface is a common entry site for the virus [53].

i. Canonical pathways

Two canonical pathways might explain the differences in cell death/survival that were reported for the LV and Lena PRRSV strains: “IL-15 production” and “Role of JAK2 in Hormone-like Cytokine Signaling”.“IL-15 production” involved five genes (IFNβ1, IL15, IRF1, JAK2 and STAT1), all of which were up-regulated following LV infection compared to Lena infection. This pathway is one hallmark of the survival and differentiation of NK and T cells [54]. Furthermore, the “Role of JAK2 in Hormone-like Cytokine Signaling” is involved in cell survival and differentiation [55].

The “TREM1 signaling” pathway involved seven genes (CCL2, IL8, JAK2, TLR3, TLR4, TLR7 and TNF), among which IL-8 and TLR4 were up-regulated during Lena infection. It has been reported that TREM1 activation is associated with increases in TLR4 [56] and interleukin-8 (IL-8) expression [57]. Interestingly, IL-8 was exclusively expressed during Lena and not LV infection, while TLR4 expression in the LV group was close to null. This finding suggests that TREM1 signaling is mostly specific to Lena infection. Based on a meta-analysis, Badaoui et al. [58] reported the possible involvement of TREM1 signaling during PRRSV infection. This finding could partially explain the high inflammatory response reported following Lena infection (Table S2), as TREM1 is directly implicated in the enhancement of inflammation [59]. From a comparative analysis of PRRSV type 1 and PRRSV type 2, Lee & Lee [60] reported that the up-regulation of IL-8 only occurred 12 h after infection with PRRSV type 2 and was not observed following infection with PRRSV type 1. As IL-8 is known to activate lymphocytes, basophils and neutrophils, its modulation might play an important role in controlling the host immune response.

The interplay between the innate and adaptive immune response is a fundamental aspect of host-pathogen interactions [61]. In this context, two related canonical pathways were found to be differentially expressed in LV- and Lena-infected samples: “Communication between Innate and Adaptive Immune Cells” and “Crosstalk between Dendritic Cells and Natural Killer Cells”. Most of the genes involved in these two canonical pathways were up-regulated to a greater extent by LV infection than by Lena infection. This result suggests that a partial down-regulation of the crosstalk between the innate and adaptive immune responses is triggered by Lena infection. The connection between NK cells and the DC pathway is required for optimal immune cell expansion and activation [62], which might explain the relative decrease in the immune response observed in the Lena-infected samples compared to the LV-infected cells.

ii. Biological functions

The cell death/survival biological functions were strongly activated during LV infection compared to Lena infection (Figure S1). Moreover, screening the network interactions among the genes involved in cell death/survival for both the LV and Lena strains (Figure S1) showed that the canonical apoptosis pathway signals (BID, CASP8, MCL1 and NFKB) were activated upon LV infection, but not following Lena infection.

Weesendorp et al. [6] observed a higher apoptosis rate following Lena infection than following LV infection in vivo. In other studies, PRRSV was found to induce the apoptosis of lymphocytes and macrophages; this would lead to a decreased host immune response and allow the virus to persist in the host [63], [64]. In addition to the difficulty in directly comparing in vitro and in vivo approaches, these discrepancies could be explained by ability of Lena to exploit novel cell receptors to target a wider population of cells compared to LV, leading to higher viral replication and higher apoptosis rate in vivo.

iii. Individual genes

TLR7 and TLR3 were both activated upon LV infection, but not upon Lena infection (Figure S6). The lack of induction of TLR3 and/or TLR7 might be a strategy developed by Lena to escape the host immune response. Su et al. [65] showed that inhibiting the TLR3/TLR7signaling pathway hinders an adequate host response to PRRSV infection. However, Zhang et al. [66] found that TLR3 and TLR7 were more highly expressed in cells infected with a highly pathogenic PRRSV isolate compared to a weakly pathogenic isolate. This discrepancy from our results might be due to the differences in the immune responses triggered by European PRRSV genotypes 1 and 3, which were used in this study, compared to PRRSV genotype 2, which was employed in [66]. Furthermore, although TLR4 was expressed in small quantities (average FPKM = 0.4 in Lena), it was shown to be more highly expressed in the Lena group than in the LV group (FC = 5.98) and was expressed at an even lower level in the LV-infected cells compared to mock-infected cells (Figure S6). TLR4 has been shown to intervene in viral infection with respiratory syncytial virus (RSV) and mouse mammary tumor virus [67].

IFN-β and IFN-αω were up-regulated following infection with either LV or Lena, but IFN-β was more highly expressed in the PAMs infected with LV (Figure 3). The difference in IFN-β expression levels between LV and Lena had thus no apparent impact on their replication rate in vitro, which was even slightly lower for Lena. This would point to the previously proposed mechanisms of impairment of the signaling pathway of IFN- β as a main strategy used by PRRSV to block or delay apoptosis until sufficient progeny have been produced [44], [68]. Such impairment could act via individual signaling pathways like RIG-I and Toll like receptor or via a cross talk among several signaling pathways. In addition, some structural proteins of the virus could explain the delayed expression of INF- β [68].

It is worth mentioning that the involvement of IFN-αω in macrophage infection by PRRSV is reported here for the first time, likely because this gene was only recently annotated in the context of the structural and functional annotation of the porcine immunome [69].

IL-10 was observed in the LV group, but not in the Lena group (Table S2). Furthermore, while TNF-α was up-regulated in both the LV and Lena groups, it was significantly up-regulated in the LV group compared to the Lena group (Table S2). The induction of TNF-α and IL-10 is used as a criterion for the classification of PRRSV isolates [70]. The relative inhibition of TNF-α and IL-10 might be another mechanism employed by the Lena strain to increase its virulence, as suggested for other PRRSV strains [71], [72].

In an in vivo transcriptomic analysis of PAMs [73], of the potential antiviral functions reported to be significant, at least two of the biological functions were also identified in this study: the “up-regulating IFN-induced genes” and “increasing intracellular zinc ion concentration” categories. The intersection between the genes found to be differentially expressed in vivo [73] and those that were differentially expressed in this study between the mock and the LV/Lena groups highlighted 22 genes (BRCA1, CCL2, CCR5, DTX3L, GBP1, GBP2, HERC6, IFIH1, IFIT1, IFIT2, IFIT3, IRF7, ISG15, MX1, PAPD5, PARP9, PLAC8, PPBP, RSAD2, SOCS1, USP18 andXAF1). These genes form a highly significant network centered on IFN signaling (Figure S7). This group of genes could therefore be considered a hallmark of infection by type 1 and type 3 European strains, independent of the infection type (in vivo/in vitro) and strain virulence (LV/Lena).

Transcriptional and post-transcriptional regulation and differential promoter use upon infection of the PAMs with the LV and Lena PRRSV strains

Most of the genes that were differentially expressed in the PAMs following PRRSV infection were expressed as several isoforms that were subjected to transcriptional/post-transcriptional regulation and/or differential promoter usage. In the results section, we classified these genes into three main groups (genes with one isoform and one TSS, Figure S5-Group 1; genes with more than one isoform and one TSS, Figure S5-Group 2; and genes with more than one isoform and more than one TSS, Figure S5-Group 3).

The first two groups included genes crucial for the innate immune response, which may be under stronger selection to prevent the emergence of new isoforms and/or post-transcriptional regulation. However, most genes (19) belonged to the third group, suggesting that they could be subjected to positive selection at the transcriptional and post-transcriptional regulation levels.

The analysis of alternative promoter usage highlighted the complexity of transcriptional regulation during the infection of PAMs with PRRSV. SLC22A15, a regulator of transmembrane transport, was the only differentially expressed gene that appeared to be significantly regulated at the promoter level, and this gene was identified only in the LV vs. mock analysis. However, seven other genes (MPDZ, ANAPC13, GPAT, ENSSSCG00000029349, LIN54, ENSSSCG00000025855 and UBXN11) were differentially regulated at the promoter level in the LV vs. Lena comparison. MPDZ is directly involved in host-pathogen interactions [74], whereas ANAPC13 and LIN54 are involved in cell cycle progression [75], [76].UBXN1 is involved in cytoskeleton remodeling [77], and GPAT plays a role in lipid synthesis [78].

Conclusion

In this study, we used RNA-Seq to profile PAMs following infection with two different strains of PRRSV to obtain a broad overview of the response of PAMs to PRRSV infection. We characterized and compared the cell responses to the two strains in terms of canonical pathways and biological functions, as well as in terms of individual genes. Therefore, our results increase our knowledge of PRRSV infection by consolidating and enriching the previous information generated using microarray technologies.

To the best of our knowledge, this is the first study in which genes, isoforms, TSSs and differential promoter usage were estimated based on RNA-Seq data obtained from PAMs infected with PRRSV. However, further targeted experiments are needed to validate and clarify these types of regulation and to determine their potential effects during macrophage infection.

Top network formed by genes that were differentially expressed in (A) LV vs. mock and (B) Lena vs. mock comparisons that are directly involved in cell death/survival. These two biological functions are highly interconnected and are highlighted in the networks as “cell death” and “cell survival”. In the LV vs. mock set, most genes were up-regulated, except forCD5, PYCARD and ARG1, which were down-regulated. In the Lena vs. mock set, except for HSPB1, all of the genes were up-regulated. The networks were constructed by using focus molecules as “seeds” that were connected together to form a network including the genes in the list. If needed, other non-focus molecules from the dataset were then added to complete the network. The resulting networks were scored and then sorted based on the score. The network scores represent the negative log of the p-value of the likelihood that the network molecules were found together by chance. Therefore, a high score represents a san index indicating that the interconnection of the molecules within the network is more likely to be true.

(PDF)

Click here for additional data file.

Example of top genes that were differentially expressed between the LV and Lena groups (among the largest fold changes). (I) Expressed only in the Lena and not in the LV group: (A) RPL31 and (B) IL8. (II) Expressed in high quantities in both the LV and Lena groups: (C) CCL2. (III) Expressed in the LV group, but in a very low quantity in the Lena group: (D)ACADVL.

(TIFF)

Click here for additional data file.

Example of the top (among the largest fold changes) differentially expressed genes (CXCL10) that were regulated transcriptionally and post-transcriptionally. CXCL10was presentin three isoforms, among which “TCONS_00046382” was differentially expressed between the LV and Lena groups. CXCL10 was also regulated at the TSS level, as it displayed three TSSs, among which “TSS31082” was differentially regulated between the LV and Lena groups.

(TIFF)

Click here for additional data file.

Cufflinks approaches for estimating (A) transcripts and their abundances and (B) transcriptional and post-transcriptional regulatory effects on overall transcript output. In Figure S4-A: For the whole analysis, Cufflinks used the paired-end reads that were aligned to the reference genome (Build Sus_scrofa.Sscrofa10.2.71.),using the TopHat software, to perform the spliced alignments (a). Cufflinks starts by connecting the compatible fragments in an overlap graph. Thus, Cufflinks applies Dilworth's Theorem, which yields a minimal set of paths that cover all of the fragments in the overlap graph, by finding the largest set of reads meeting the criterion that no two reads could have originated from the same isoform (b,c). Subsequently, Cufflinks estimates transcript abundance using a statistical model in which the probability of observing each fragment is a linear function of the abundances of the transcripts from which it could have originated (d). The last step consists of maximizing the likelihood function for all possible sets of relative transcript abundance to determine the set that best explains the observed fragments (e). In Figure S4-B: (a) When the abundance of isoforms A, B and C are grouped by TSS, the changes in the relative abundance of the TSS groups indicate transcriptional regulation (A+B vs. C). Post-transcriptional effects are observed as changes in the levels of the isoforms ina single TSS group (A vs. B) (Adapted from Trapnell et al. 2012).

(PDF)

Click here for additional data file.

Transcriptional/post-transcriptional regulation of the genes involved in the top canonical pathways in the LV vs. mock, Lena vs. mock and LV vs. Lena comparisons. (A) Un-spliced and transcriptionally regulated genes,(B) spliced and post-transcriptionally regulated genes and(C) spliced and both transcriptionally and post-transcriptionally regulated genes. For each transcript,the “XLOC”, “TSS” and “TCONS” suffixes correspond to the genes, TSSs and isoforms, respectively. Differentially expressed isoforms with different TSSs are transcriptionally regulated, while isoforms with the same TSS are regulated at the post-transcriptional level (Figure S4).

(PDF)

Click here for additional data file.

Differential expression ofTLR3, TLR4 and TLR7 between the LV and Lena groups. TLR7 and TLR3 were significantly up-regulated in the LV group, while neither TLR7 nor TLR3wasmodulated in the Lena group. Although TLR4 was expressed in small quantities (medium FPKM = 0.4 in Lena), it was more highly expressed in the Lena group than in the LV group (FC = 5.98) and was expressed at an even lower level during LV infection than during mock infection.

(PDF)

Click here for additional data file.

Network formed by common genes that were differentially expressed in the LV vs. mock and Lena vs. mock comparisons from one side and in an in vivo study performed by Zhou et al. (2001) from the other side. The canonical pathways that were significantly affected by this group of genes are highlighted with blue squares and include interferon signaling, the activation of IRF by cytosolic pattern recognition receptors, the role of hypercytokinemia/hyperchemokinemia in the pathogenesis of influenza, the role of RIG1-like receptors in antiviral innate immunity and the role of PI3K/AKT signaling in the pathogenesis of influenza. The networks were constructed using focus molecules as “seeds” that were connected together to form a network using the genes in the list. If needed, other non-focus molecules from the dataset were then added to complete the network. The resulting networks were scored and then sorted based on the score. The network scores represent the negative log of the p-value of the likelihood that the network molecules were found together by chance. Therefore, a high score represents an index indicating that the interconnection of the molecules within the network is more likely to be true.

(PDF)

Click here for additional data file.

Summary of the reads generated per sample (mock, LV and Lena). A. Numbers of RNA-Seq reads generated per sample (mock, LV and Lena). B. RNA-Seq reads eliminated following quality evaluation via FASTQC and Sickle. C. RNA-Seq reads mapped to the pig genome (Build Sus_scrofa.Sscrofa10.2.71.) using the TopHat v2.0.8 algorithm.

(DOC)

Click here for additional data file.

Genes, isoforms, TSSs and promoter usage showing differential expression in the LV vs. mock, Lena vs. mock and LV vs. Lena comparisons. Genes_Mock_vs_LV: differentially expressed genes between the LV and mock infections; Genes_Mock_vs_Lena: differentially expressed genes between the Lena and mock infections; Genes_LV_vs_Lena: differentially expressed genes between the LV and Lena infections. Isoforms_Mock_vs_LV: differentially expressed isoforms between the LV and mock infections; Isoforms_Mock_vs_Lena: differentially expressed isoforms between the Lena and mock infections; Isoforms_LV_vs_Lena: differentially expressed isoforms between the LV and Lena infections. TSS_Mock_vs_LV: differentially expressed TSSs between the LV and mock infections; TSS_Mock_vs_Lena: differentially expressed TSSs between the Lena and mock infections; TSS_LV_vs_Lena: differentially expressed TSSs between the LV and Lena infections. Promoters_Mock_vs_LV: differentially used promoters between the LV and mock infections; Promoter_Mock_vs_Lena: differentially used promoters between the Lena and mock infections; Promoter_LV_vs_Lena: differentially used promoters between the LV and Lena infections. Each page is composed of 7 columns: test_id, gene_id, gene, gene-map, locus, sample_1, sample_2, status, value_1, value_2, log2.fold_change, test_stat, p_value, q_value and significant.

(XLS)

Click here for additional data file.