Time Series Transcriptomic Analysis of Bronchoalveolar Lavage Cells from Piglets Infected with Virulent or Low-Virulent Porcine Reproductive and Respiratory Syndrome Virus 1

ABSTRACT Porcine reproductive and respiratory syndrome virus (PRRSV) has evolved to escape the immune surveillance for a survival advantage leading to a strong modulation of host’s immune responses and favoring secondary bacterial infections. However, limited data are available on how the immunological and transcriptional responses elicited by virulent and low-virulent PRRSV-1 strains are comparable and how they are conserved during the infection. To explore the kinetic transcriptional signature associated with the modulation of host immune response at lung level, a time-series transcriptomic analysis was performed in bronchoalveolar lavage cells upon experimental in vivo infection with two PRRSV-1 strains of different virulence, virulent subtype 3 Lena strain or the low-virulent subtype 1 3249 strain. The time-series analysis revealed overlapping patterns of dysregulated genes enriched in T-cell signaling pathways among both virulent and low-virulent strains, highlighting an upregulation of co-stimulatory and co-inhibitory immune checkpoints that were disclosed as Hub genes. On the other hand, virulent Lena infection induced an early and more marked “negative regulation of immune system process” with an overexpression of co-inhibitory receptors genes related to T-cell and NK cell functions, in association with more severe lung lesion, lung viral load, and BAL cell kinetics. These results underline a complex network of molecular mechanisms governing PRRSV-1 immunopathogenesis at lung level, revealing a pivotal role of co-inhibitory and co-stimulatory immune checkpoints in the pulmonary disease, which may have an impact on T-cell activation and related pathways. These immune checkpoints, together with the regulation of cytokine-signaling pathways, modulated in a virulence-dependent fashion, orchestrate an interplay among pro- and anti-inflammatory responses. IMPORTANCE Porcine reproductive and respiratory syndrome virus (PRRSV) is one of the major threats to swine health and global production, causing substantial economic losses. We explore the mechanisms involved in the modulation of host immune response at lung level performing a time-series transcriptomic analysis upon experimental infection with two PRRSV-1 strains of different virulence. A complex network of molecular mechanisms was revealed to control the immunopathogenesis of PRRSV-1 infection, highlighting an interplay among pro- and anti-inflammatory responses as a potential mechanism to restrict inflammation-induced lung injury. Moreover, a pivotal role of co-inhibitory and co-stimulatory immune checkpoints was evidenced, which may lead to progressive dysfunction of T cells, impairing viral clearance and leading to persistent infection, favoring as well secondary bacterial infections or viral rebound. However, further studies should be conducted to evaluate the functional role of immune checkpoints in advanced stages of PRRSV infection and explore a possible T-cell exhaustion state.

ABSTRACT Porcine reproductive and respiratory syndrome virus (PRRSV) has evolved to escape the immune surveillance for a survival advantage leading to a strong modulation of host's immune responses and favoring secondary bacterial infections. However, limited data are available on how the immunological and transcriptional responses elicited by virulent and low-virulent PRRSV-1 strains are comparable and how they are conserved during the infection. To explore the kinetic transcriptional signature associated with the modulation of host immune response at lung level, a time-series transcriptomic analysis was performed in bronchoalveolar lavage cells upon experimental in vivo infection with two PRRSV-1 strains of different virulence, virulent subtype 3 Lena strain or the low-virulent subtype 1 3249 strain. The time-series analysis revealed overlapping patterns of dysregulated genes enriched in T-cell signaling pathways among both virulent and low-virulent strains, highlighting an upregulation of co-stimulatory and coinhibitory immune checkpoints that were disclosed as Hub genes. On the other hand, virulent Lena infection induced an early and more marked "negative regulation of immune system process" with an overexpression of co-inhibitory receptors genes related to T-cell and NK cell functions, in association with more severe lung lesion, lung viral load, and BAL cell kinetics. These results underline a complex network of molecular mechanisms governing PRRSV-1 immunopathogenesis at lung level, revealing a pivotal role of co-inhibitory and co-stimulatory immune checkpoints in the pulmonary disease, which may have an impact on T-cell activation and related pathways. These immune checkpoints, together with the regulation of cytokine-signaling pathways, modulated in a virulence-dependent fashion, orchestrate an interplay among pro-and anti-inflammatory responses. IMPORTANCE Porcine reproductive and respiratory syndrome virus (PRRSV) is one of the major threats to swine health and global production, causing substantial economic losses. We explore the mechanisms involved in the modulation of host immune response at lung level performing a time-series transcriptomic analysis upon experimental infection with two PRRSV-1 strains of different virulence. A complex network of molecular mechanisms was revealed to control the immunopathogenesis of PRRSV-1 infection, highlighting an interplay among pro-and anti-inflammatory responses as a potential mechanism to restrict inflammation-induced lung injury. Moreover, a pivotal role of co-inhibitory and co-stimulatory immune checkpoints was evidenced, which may lead to progressive dysfunction of T cells, impairing viral clearance and leading to persistent infection, favoring as well secondary bacterial infections or viral rebound. However, further studies should be conducted to evaluate the functional role of immune checkpoints in advanced stages of PRRSV infection and explore a possible T-cell exhaustion state.
severe interstitial pneumonia and getting to the highest score in the Lena group for the additional presence of extensive areas of pulmonary consolidation in cranial and middle lobes from 6 dpi onwards, which was confirmed by histopathological evaluation. Virulent Lena strain induced an earlier (2 out of 5 piglets PRRSV-1 positive at 1 dpi) and higher replication in the lung compared to low-virulent 3249 strain (P , 0.01 at 3 dpi; P , 0.05 at 6 and 8 dpi), displaying a lung viral load peak at 6 (Cq 18.9 6 0.9) and 8 dpi (Cq 22.9 6 2.5), respectively (Fig. 1A). Control animals did not show clinical signs and minimal lung lesion, remaining as RT-qPCR-negative for PRRSV-1 throughout the study.
CD163 + cells decreased within live PAMs together with a mixed influx of neutrophils, monocytes, and lymphocytes in BAL cells of Lena-infected pigs. FCM analysis in BAL cells were thoroughly described by Rodríguez-Gómez et al. (2019). Briefly, CD163 expression was studied only in live BAL cells from either control or PRRSV-infected pigs. Control animals showed a homogeneous and stable subset of cells within live PAMs, compatible with PAMs because of size and granularity properties as well as CD163 labeling, throughout the study (Fig. 1B, red circle, and Fig. 1C). By contrast, both subsets, PAMs and CD163 1 cells, decreased alongside the study in PRRSV-1-infected pigs (Fig. 1B, red arrows, and 1C). This drop was more marked and occurred earlier (3 dpi) in Lena-infected pigs than in 3249 group (8 dpi) and was accompanied by a mixture of neutrophils, monocytes, and, to a lesser extent, lymphocytes (Fig. 1B, green circle) from 6 dpi onwards in Lena-infected pigs and at 13 dpi in 3249-infected animals.
Time-series differential gene expression pattern related to PRRSV-1 infection. A total of 1,935 and 8,352 significant variable genes were detected by MasigPro during the time course for 3249-and Lena-infected pigs, respectively. After clustering, these significant variable genes were subdivided into 5 different clusters according to their expression profiles, thus, the genes included within a cluster presented a similar expression pattern. Fig. 2 illustrates the median expression level of clusters 3, 4, and 5 at each time point. Clusters 1 and 2 showed an inconsistent trend among experimental groups. Cluster 3 consisted of 98 genes whose median expression level increased in 3249-infected piglets from 8 dpi throughout the study, but this increase was higher and earlier, from 6 dpi, in Lena-infected piglets. GO enrichment-functional analyses were performed to gain biological understanding of variable genes included in cluster 3, which were mainly enriched in GO involved in "regulation of immune system process," "leukocyte differentiation," "regulation of cytokine production," "negative regulation of immune system process," "lymphocyte activation," "response to cytokine," and "positive regulation of apoptotic process" (Fig. 2A). Genes in cluster 4, a total of 59, followed a similar kinetics to the one described above for cluster 3 in 3249-infected piglets, but for Lena-infected pigs, the expression of these genes had a marked and rapid increase from 1 to 3 dpi remaining roughly constant throughout the study (Fig. 2). The significantly enriched GO terms related to cluster 4 were "cytokine-mediated signaling pathway," "response to bacterium," "regulation of cytokine production," "IFN-I signaling pathway," and "cellular response to tumor necrosis factor" (genes and pathways are listed in Fig. 2B). Cluster 5 comprised 45 genes whose expression rose from 8 to 13 dpi in both infected groups (Fig. 2). Genes in cluster 5 were involved in the GO terms "lipid transport" and "carboxylic acid transport" (Fig. 2C).
DEGs profile in BAL cells of PRRSV-1-infected piglets. DEGs were identified in response to 3249 or Lena infection by comparing the gene expression levels of PRRSV-1-infected animals to control animals at each time point using an FDR , 0.05 and a log2FC $ 1.0 (upregulated genes) or log2FC # -1.0 (downregulated genes) as the cutoff criteria. Thus, a total of 14,368 DEGs in 3249 data set and 11,656 DEGs in Lena data set were found to be differentially expressed along the study ( Fig. 3A and B).
Similar dynamic changes in the gene expression level were observed in both 3249and Lena-infected piglets, with an increasing number of DEGs along the study. As expected, in both infected groups, the lowest number of DEGs was observed at 1 dpi (data not showed). At 3 dpi, 38 DEGs (downregulated: 5; upregulated: 33) and 863   3A and B).
Identification and functional classification of conserved DEGs in response to 3249 strain infection. A total of 14,368 DEGs were identified in response to 3249 infection during the study. Most of the DEGs (9,218) were expressed at only one of the evaluated time points, hence, we focused on those DEGs that were conserved in two or more time points, to achieve a better understanding of the infection. Thus, whereas no DEGs were concomitantly identified during the early time points of infection (1 or 3 dpi), 96 overlapped DEGs were found at 6-8-13 dpi, of which 13 genes were downregulated and 83 were upregulated (Fig. 4A). To gain further insight into the potential functions of these 96 DEGs, a functional-enrichment analysis based on the GO database was performed to explore BPs and ISPs. The significantly enriched GO terms with a P value , 0.05 are illustrated in Fig. 4B to D. The results of the GO analysis reported that most DEGs were mainly involved in the "innate immune response" (17 DEGs, 28.80%), "MAPK cascade" (11 DEGs, 18.64%), "defense response to virus" (6 DEGs, 7.55%), "regulation of JNK cascade" (5 DEGs, 8.47%), "cellular response to interferon g (IFN-g)" (5 DEGs, 8.47%), and "regulation of a/b T-cell activation" (4 DEGs, 6.78%), among others. The patterns of expression of representative genes for low-virulent 3249 strain are illustrated in Fig. 4C. Fig. 4B-4D list the symbol and number of genes enriched in each GO term.
Moreover, the number of conserved DEGs increased at 8-13 dpi, with 854 DEGs identified at these time points (Fig. 4A). Among these genes, only 16 genes were downregulated, whereas 838 were upregulated. In this case, the GO functional-enrichment analysis (BPs and ISPs categories) revealed that T-cell-related pathways were significantly altered with 7 different GO terms ("T-cell co-stimulation," "regulation of a/b T-cell differentiation," "T-cell selection," "regulation of T-cell receptor signaling pathway," "regulatory T-cell differentiation," "CD8 1 a/b T-cell activation," and "CD4 1 CD25 1 a/b T-cell differentiation") and an occurrence of enriched genes over 50%. In addition, other highly enriched GO terms were related to "regulation of interleukin-2 (IL-2) production" (17 DEGs, 14.41%), "regulation of lymphocyte chemotaxis" (9 DEGs, 7.63%), or "regulation of tolerance induction" (6 DEGs, 5.08%). Fig. 5A and B list the symbol and number of genes enriched in each GO term.
PPI network construction and Hub genes identification for 3249 infection. The PPI network was constructed with the 854 conserved DEGs found at 8 and 13 dpi by the STRING database, with 402 nodes and 2,400 edges. The Hub genes were selected from the whole PPI network using the MCC and DMNC algorithm of the CytoHubba plugin. According to the MCC and DMNC scores, the top 10 highest-scored genes found for each approach include PDCD1, LAG3, TNFRSF18, ICOS, CD8A, CD6, CTLA4, CD28, GZMB, TBX21, SELL, IL7R, IL2RA, IL2RB, PRF1, FASLG, KLRK1, CD226 and CD40LG ( Fig. 6A and B). Also, the kinetic of expression of these genes is illustrated in Fig. 6C-6D. Because of the complexity, the whole PPI network was analyzed using the MCODE plugin to identify essential PPI network modules, making the understanding of molecular mechanisms more approachable. The top 3 significant clusters with a k-core . 6 were named as A, B,  Because of the low number of DEGs not data was showed at 1 dpi. and C, respectively, and selected as sub-networks (Fig. 7A). Interestingly, all the Hub genes identified for the 3249 strain were included in cluster A. "Positive regulation of leukocyte activation" (11 out of 20 Hub genes), "regulation of T-cell activation (11 out of 20 Hub genes) and differentiation (9 out of 20 Hub genes)," "lymphocyte migration" (5 out of 20 Hub genes), "positive regulation of a/b T-cell activation" (2 out of 20 Hub genes), "IL-2 receptor activity" (3 out of 20 Hub genes), and "tolerance induction" (3 out of 20 Hub genes) were the main terms involved after screening of GO enrichment analyses (BPs and ISPs categories) of genes included within cluster A (Fig. 7B, Table 1).
Identification and functional classification of conserved DEGs in response to virulent Lena strain infection. A sum of 11,656 DEGs were found in the Lena data set (Fig. 8A), with most of them (6,760 genes) being expressed only at one single time point during the infection. The same approach described above was followed to analyze the 3249 data set; the analysis was addressed in DEGs identified at two or more time points.
At 3-6-8-13 dpi, 82 overlapped DEGs were identified among which only 3 were downregulated, whereas 79 were upregulated (Fig. 8A). In the next step, a GO functional-enrichment analysis was conducted (BPs and ISPs categories), with most of the DEGs enriched in terms associated with "IFN-I signaling pathway" (8 DEGs, 19.05%), "neutrophil migration" (7 DEGs, 16.78%), "negative regulation of T-cell mediated cytotoxicity" (3 DEGs, 7.14%) and "negative regulation of CD8 1 a/b T cell activation" (3 DEGs, 33.33%). Fig. 9 lists the symbol and number of genes enriched in each GO term, as well as the patterns of expression of representative genes for the virulent Lena strain (Fig. 9B).
PPI network construction and Hub genes identification for virulent Lena infection. In this case, the 1,227 overlapped DEGs detected at 8-13 dpi were utilized to construct the PPI network by the STRING database, finding 560 nodes and 4,461 edges. Afterwards, according to the MCC and DMNC algorithm of the CytoHubba plugin, the top 10 highest-scored genes found for each approach from the whole PPI network were identified as Lena strain Hub genes, highlighting CTLA4, CCR7, SELL, CD28, CD274, IL-10, GZMB, TBX21, CD40LG, CD27, LAG3, PDCD1, HAVCR2, ICOS, TNFRSF18, TNFRSF4, TNFRSF9, GPR18, CD1D and CD226 ( Fig. 12A and B). The kinetics of expression Hub genes disclosed for the virulent Lena strain was displayed in Fig. 12C and D. Of note, several of these Lena strain Hub genes (11 out of 20) were previously identified for the 3249 strain. After analyzing the whole PPI network by MCODE plugin, the top 3 more significant sub-networks were selected and named as cluster A, cluster B, and cluster C (Fig. 13A). Cluster A, which included most of Lena strain Hub genes, was subsequently subjected to a GO functional-enrichment analysis (BPs and ISPs categories). "Regulation of lymphocyte apoptotic process" (4 out of 20 Hub genes), "T-cell selection" (7 out of 20 Hub genes), "positive regulation of IL-4 production" (5 out of 20 Hub genes), "T-cell co-stimulation" (5 out of 20 Hub genes), "regulation of tolerance induction" (4 out of 20 Hub genes), "regulation of IL-2 production" (4 out of 20 Hub genes), "regulation of regulatory T-cell differentiation" (4 out of 20 Hub genes), "regulation of lymphocyte chemotaxis," and "regulation of IFN-a production" (2 out of 20 Hub genes) were reported as highly enriched GO terms (Fig. 13B, Table 3).

DISCUSSION
PRRSV have evolved different mechanisms to evade the host's immune response (18,20), favoring also secondary infections that result in severe disease such as PRDC (23). Understanding the mechanisms involved in the dysregulation of the immune response at lung level is a cornerstone for unraveling the immunopathogenesis of PRRS, one of the most economically devastating diseases for swine production.
Time-series analysis pointed out an early and sustained upregulation of coinhibitory receptors related to T-cell and NK cell functions over time points. In the present study, RNA-seq was used to perform a time course analysis in BAL cells upon PRRSV in vivo infection with strains of different virulence to gain further insight into the potential pathways involved in modulating the host's immune response. Clusters 3, 4, and 5 consisted of genes whose expression levels increased during the time course, although with different patterns and intensity for each strain. Genes in cluster 3, which exhibited a similar kinetics to lung lesion, were enriched in the GO terms "regulation of immune system process," "leukocyte differentiation," "regulation of cytokine production," "lymphocyte activation," and "response to cytokine and positive regulation of apoptotic process," among others. These terms related to viral infection and immune response signaling pathways were upregulated earlier and higher in virulent Lena-infected than in low-virulent 3249-infected piglets, which may explain the severe lung injury observed. In addition, a group of genes was enriched for "negative regulation of immune system process" (Fig. 2A), including CD96, CTLA4, RFN128, TIGIT, and SAMHD1, which play an inhibitory role in T-cell and NK cell functions (33)(34)(35)(36)(37)(38), potentially hampering viral clearance and the onset of an effective adaptive immune response.
The interplay among IFN-I/IFN-II orchestrated an antiviral response resulting in IFN-stimulated genes unable to control efficiently viral replication, probably due to PRRSV countermeasures. IFN-I leads to a signal transduction cascade of hundreds of IFN-stimulated genes (ISGs), a powerful instrument that interferes with viral replication (44,45). In our case, a screen of ISGs, enriched in IFN-I signaling pathway, was upregulated in BAL cells from both virulent Lena-infected (BST2, IFI6, IFITM1, ISG15, ISG20, MX1, OASL, RSAD2, STAT1) and low-virulent 3249-infected piglets (EIF2AK4, GBP1, MX1) from 3 and 6 dpi onwards, respectively. Interferon regulatory factors (IRF) have been also found to be overexpressed in PAMs infected in vitro with Lena and Lelystad virus strains (29), as well as during the early stage of PRRSV-2 infection (24,46,47). Many of these upregulated genes could point to a continuous antiviral state (45,48), since ISG20, MX1, and RSAD2 have the capacity to inhibit viral replication (49-51), nevertheless, PRRSV has been reported to inhibit IRF3 as well as several of these ISGs such as STAT1, ISG15, and IFITM1 (52,53), by viral non-structural proteins (nsp1a, nsp1b, nsp2, nsp4, and nsp11) (53)(54)(55)(56). Thus, unlike the 3249 strain, the virulent Lena strain induced an upregulation of DEGs enriched in negative regulation of IFN-I and type II IFN (IFN-II) at 8-13 dpi in our study. Altogether, these results draw an interplay among IFN-I/IFN-II induced antiviral response and PRRSV countermeasures resulting in ISGs unable to control efficiently viral replication, leading to an incomplete viral clearance and potential persistent infection (57).
Early upregulation of chemokines and cytokines in response to virulent Lena infection led to an acute lung injury. Virulent PRRSV-1 strains trigger an exacerbated and early pro-inflammatory response compared with low-virulent strains, leading to severe lung injury, secondary bacterial infection, and hence, acute respiratory symptoms (12,39,58,59). A considerable number of DEGs were enriched in GO terms involved in macrophages activation, pro-and anti-inflammatory cytokines, as well as chemokines for the recruitment of neutrophils, monocytes, and lymphocytes to the site of infection in both PRRSV-1-infected groups during the study. Interestingly, Lenainfected pigs showed 5 DEGs, which were upregulated from 1 to 13 dpi, with three of them (CSF1, CCL4, and EDN1) showing striking inflammatory properties, such as monocyte and macrophage proliferation and differentiation (CSF1) (60), vasoconstriction and pro-inflammatory response associated with acute lung injury (ALI) (EDN1) (61), or neutrophil migration (CCL4) (62). Although, PAMs are the key players during PRRSV infection, the GO analysis suggests that chemotaxis of neutrophils and their likely degranulation are continued over the time points (25,46), leading to alveolar-capillary barrier damage, increased vascular permeability, and hence, a higher extent of lung damage, as has been observed in this study.
According to the dynamic changes observed in BAL cells from PRRSV-1-infected piglets, the "positive regulation of lymphocyte migration" term was observed in virulent Lena-infected piglets from 6 dpi onwards, including upregulated expression of ADAM8, CCL11, CCL20, CCL3, CCL4, and CXCL10. By contrast, in the 3249 group not only were these genes overexpressed later (8-13 dpi), but also, genes enriched in "negative regulation of lymphocyte migration" were found (CCL2, GCSAM, KLRK1, PADl2, RIPOR2). CCL20/MIP-3a is a strong chemotactic factor for lymphocytes, while a weak factor for neutrophils, that is upregulated by IFN-g, tumor necrosis factor (TNF) and lipopolysaccharide, inducing a marked lung damage upon viral infection (63). CCL2, CCL3, CCL4, and CCXL10 are produced in the lung in the early phase of lung inflammation, working as a chemotactic factor for macrophages, NK cells, and T cells during viral infection (64). In addition, there was an upregulation of genes associated with "positive regulation of apoptotic cell clearance" (C3, C4, CD300LF, CCL2) in Lena-infected piglets from 6 to 13 dpi, probably, due to the severe activation of regulated cell death phenomena at lung level induced by the virulent Lena strain (42), suggesting an attempt to remove dead cells and cellular debris to restore tissue damage.
Interestingly, many cytokine-signaling pathways were found to be regulated in Lenacompared with 3249-infected piglets, highlighting an interplay among pro-and antiinflammatory responses that are likely to be activated in a virulence dependent fashion and hence influencing the profile of innate cytokines as well as the development of adaptive immunity. IFN-g release pathway was observed in the 3249 group from 6 dpi onwards, whereas SOCS1 and SOSC3, potent IFN-g inhibitors (65), were early overexpressed in the Lena group (6 dpi) followed by an interface of upregulated genes that may promote IFN-g production (8-13 dpi). IL1B (IL-1b) and TNF genes were upregulated in the Lena group, unlike the 3249 group. While classical PRRSV strains are reported to show a delayed and weak systemic production of these cytokines (66)(67)(68)(69), this proinflammatory profile seems to be a hallmark of virulent PRRSV strains (24,25,39,58,59), resulting in being useless in controlling viral replication but exacerbating clinical signs and lung damage.
Moreover, the virulent Lena strain induced a strong regulation of adaptive-cytokines-related pathways such as "regulation of IL-2, IL-4, and IL-12 production," as well as an upregulation of IL-10 (IL-10) and IL18BP (IL-18BP, IL-18 binding protein) from 8 to 13 dpi, in contrast to the 3249 group, in which only "regulation of IL-2 production" was found. Most of the genes enriched in "regulation of IL-12 production" (7 out of 12), such as IL-10 and IDO1, show a negative regulation of this pathway. Likewise, upregulation of IL-10 expression has been reported to impair Th1 response during PRRSV infection, inducing a niche for viral-specific regulatory T cells (Tregs) in a strain-dependent fashion, though Tregs induction after PRRSV infection is controversial (70)(71)(72). Although no changes in the serum concentration of IL-10 were observed in Lena-and 3249-infected pigs (41), IL-10 could be locally secreted in the lung, exerting a paracrine in situ effect. IL18BP gene encodes a soluble receptor that hijacks IL-18, blocking the engagement of IL-18 to its receptor, which in turn inhibits IL-18-induced IFN-g production and alters Th1 response (73,74). The anti-inflammatory mediators herein reported may play a double-edged sword role, limiting the inflammatory response during the early stage of the infection but also, when sustained in time, weakening the host's immune response against PRRSV and hence impairing viral clearance.
Activation of T-cell signaling pathways was overlapped in virulent and lowvirulent PRRSV-1 strains in a time-dependent manner. Activation of T-cells is crucial for anti-PRRSV adaptive immunity since T cells are key players in tackling viral infection, particularly, during the first stages of infection (17,75). The virulent Lena strain caused a prompt upregulation of genes enriched in "negative regulation of T cell mediated cytotoxicity," "negative regulation of CD8 1 a/b T-cell activation," "negative regulation of T-cell activation via T-cell receptor contact with antigen bound to MHC molecules on antigen presenting cells," and "regulation of T cell apoptotic process," which would impair the induction of CD8 1 cytotoxic T lymphocytes (CTLs) needed to eliminate virus-infected cells, release of pro-inflammatory cytokines, and establish memory cells (18,75). Besides, genes associated with "regulation of chronic inflammatory response" were found to be overexpressed, suggesting that virulent Lena strain would be able to overcome the host's antiviral strategies, inducing a persistent infection. Thus, many of these DEGs (SOCS1, LILRB1 (CD85J), CD274 (PD-L1), HAVCR2 (TIM3), and IDO1) enriched in the above-mentioned pathways are immune checkpoints associated with other persistent viral infections inducing Tcell exhaustion (33,76,77). On the other hand, a substantial number of genes were regulated from 8 dpi onwards in either Lena-or 3249-infected groups, and these DEGs were enriched in quite similar T-cell-related pathways, highlighting "T-cell co-stimulation," "T-cell selection," "regulation of a/b T-cell differentiation," and "CD8 1 a/b T-cell activation." Tcell stimulation is an essential mechanism that relies on co-stimulatory and co-inhibitory receptors to achieve an effective T-cell activation (78). Several viruses involved in persistent infections exploit these adaptive mechanisms to evade immune-mediated viral clearance, leading to viral persistence and, hence, a functionally inferior T-cell response (33,79). Thereby, either Lena or 3249 elicited an upregulation of co-stimulatory receptors (CD28, ICOS, CD40LG (CD154), CD27, TNFRSF4 (OX40), TNFRSF9 (CD137) or TNFRSF18 (GITR)) but also inhibitory receptors (CD5, CTLA4, CD274 (PD-L1)/PDCD1, HAVCR2 (TIM3), LAG3, TIGIT, TOX), some of which were also overexpressed in the thymus of these animals (80), suggesting an attempt to hamper T-cell activation. The transcription factors GATA3 and TBX21 (Tbet), overexpressed in both low-virulent 3249 and virulent Lena infected groups, have been reported to be involved in T-cell maturation and Th1 polarization of porcine alphabeta T-cells (81); however, EOMES, upregulated as well, and TBX21 could also cooperate to sustain exhausted CD8 1 T-cell subsets (78,82,83). Besides, we found a profile of upregulated DEGs enriched in "regulation of chronic inflammatory response" (ADORA2B, CCL5, 1001, IL-10, TNF, TNFAIP3) and "regulation of tolerance induction" (CD274, CD3E, HAVCR2, 1001, IL2RA, MARCHF7, PDCD1). Considering all the above-mentioned, two different scenarios could be drawn. First, these pathways regulate self-tolerance, minimizing bystander tis-   sue damage and hindering the inflammatory response observed mainly in Lena-infected piglets. A second possibility is that PRRSV might induce T-cell anergy if upregulation of inhibitory receptors is sustained over time and associated with a high level of viral antigen, as it has been described for other arterivirus, such as Equine arteritis virus (84). Although multiple elements are responsible for T-cell exhaustion (79), these findings, together with the well-known persistence of PRRSV in lymphoid organs (57), point out that PRRSV-1 could induce adaptive immune tolerance, making viral transmission in the farm easier. Immune checkpoints were revealed as Hub genes conserved in both virulent and low-virulent PRRSV-1 infection. A Hub gene analysis was conducted by CytoHubba, since Hub genes set up more complex interactions compared with other genes in a PPI network and hence play a central role in the mechanisms of disease (85). Consequently, exploring Hub genes and their associated key pathways would be essential for a better understanding of PRRSV-1 infection at lung level. Eleven of 20 Hub genes were shared in both 3249 (Fig. 6) and Lena (Fig. 12) strains. Most of the Hub genes followed a similar kinetic for both strains, Lena and 3249, increasing their expression from 8 to 13 dpi. Nevertheless, several of them, SELL, GZMB, and CCR7 (at 3 dpi), and HAVCR2 (TIM-3), TNFRSF9 (CD137), and GPR18 (at 6 dpi) were promptly upregulated for the virulent Lena strain, probably associated with the marked clinical signs, the earlier and stronger peak of replication, and severe pneumonia. Remarkably, 13 of 20 Hub genes are considered as co-stimulatory immune checkpoints that are gaining a special interest nowadays, because of their role in cancer and chronic viral infections (86,87). Mainly the Lena, but also the 3249 strain, induced an upregulation of co-stimulatory receptors such as CD28, ICOS, and several members of the TNF receptor family (CD40LG (CD154), CD27, TNFRSF4 (OX40), TNFRSF9 (CD137) and TNFRSF18 (GITR)), which are considered as key players in CD4 1 and CD8 1 T-cell differentiation and expansion in acute and chronic viral infections (88)(89)(90). By contrast, co-inhibitory receptors, such as CTLA4, CD274 (PD-L1), PDCD1 (PD-1), LAG3, and HAVCR2 (TIM-3), working together with IL-10, have been reported in acute viral infection, such as SARS-CoV-2 or influenza virus, playing a double role, one limiting inflammation induced tissue damage and ALI, as can particularly be the case for the virulent Lena strain, and another one, inducing a progressive dysfunction of T cells, impairing viral clearance and leading to chronic infection (33,77,79). A sustained upregulation of these Hub genes has been closely associated with CD8 1 and CD4 1 T-cells dysfunction during persistent viral infections together with high viral loads and antigen levels (33, 79, 87). Subsequently, a MCODE analysis was conducted to evaluate the interaction of Hub genes with other DEGs, clustering most of them in cluster A in both Lena-and 3249-infected groups. Therefore, "regulation of lymphocyte apoptotic process," "T-cell selection," "positive regulation of IL-4 production," "T-cell co-stimulation," and "regulation of tolerance induction" seem to be significant pathways related to virulent Lena infection, whereas "positive regulation of leukocyte activation," "regulation of T-cell activation and differentiation," "lymphocyte migration," and "tolerance induction" would be the essential ones for the low-virulent 3249 strain. These results highlight that the mechanism of PRRSV infection with strains of different virulence could follow common pathways with a greater activation of pro-and anti-inflammatory response, and indeed, severe lung damage in Lena-infected piglets, probably associated with a superior replication rate (41).

CONCLUSION
The present study dissects the panoply of molecular mechanisms that govern the immunopathogenesis of PRRSV-1 infection with strains of different virulence at lung level, involving many conserved molecular pathways for virulent Lena and low-virulent 3249 strain. These included the activation of co-inhibitory and co-stimulatory immune checkpoints, revealed as Hub genes, in the pulmonary disease, which may have an impact on activation of CD8 1 and CD4 1 T cells and other T-cell related-pathways. Nevertheless, virulent Lena infection resulted in more pronounced and earlier

MATERIALS AND METHODS
Porcine reproductive and respiratory syndrome virus strains. The low-virulent 3249 strain (subtype 1 PRRSV-1) was isolated from the serum of a piglet with pneumonia from a PRRSV-positive herd located in Spain in 2005 (91). The virulent Lena strain (subtype 3 PRRSV-1), considered as the prototype of PRRSV-1 virulent strains, was isolated from the lung of weak born piglets from a PRRSV outbreak characterized by high mortality rate, reproductive failure, and respiratory disorders in a Belarusian farm in 2007 (15). Both strains were propagated in porcine alveolar macrophages (PAMs), then the viral stocks were produced from the fourth passage of each strain, titrated by means of immunoperoxidase monolayer assay and expressed as tissue culture infectious doses 50 (TCID 50 )/mL (Lena strain: 10 5.66 TCID 50 /mL; 3249 strain: 10 5.79 TCID 50 /mL).
Animals and experimental design. Animals used in this study belong to a large project carried out in order to investigate the pathogenesis of PRRSV-1 strains of different virulence (40). Briefly, 65 4-weekold piglets (Landrace Â Large White crossbreed), obtained from a historically PRRSV-negative farm, were randomly assigned to three different groups and housed in separate pens in biosafety level III containment facilities (IRTA-CReSA, Cerdanyola del Vallès, Barcelona, Spain): 3249 group (n = 25), Lena group (n = 25), and control group (n = 15). At the beginning of the study, all piglets were ELISA and PCR-negative for porcine circovirus type 2 (PCV2), PRRSV, and Mycoplasma hyopneumoniae (92,93).
At necropsy, gross lung lesions were evaluated and recorded by the same pathologist as previously described by Halbur et al. (94). Parallel samples from cranial, middle, and caudal lobes of the right lung were collected and immediately frozen at 280°C for the quantification of lung viral load or fixed in 10% neutral buffered formalin (Fisher Scientific, Ltd., Loughborough, UK) for histopathological evaluation.
The histopathological findings were blindly evaluated and graded by two different pathologists. The severity of lung lesions for the interstitial pneumonia was scored as previously described (94), whereas the score for suppurative bronchopneumonia was previously described by Rodríguez-Gómez et al. (2019). Overall, the final score for each piglet was calculated as the sum of both the interstitial pneumonia and the bronchopneumonia scores.
A lung tissue homogenate was made previously to carry out RNA isolation and purification by using TRIzol LS Reagent (Thermo Fisher Scientific) followed by NucleoSpin RNA virus columns kit according to manufacturer's protocols (Macherey-Nagel, Düren, Germany). According to subgenomic copies of the virus, results for viral load were expressed by changes in quantification cycle (Cq), to not overestimate the number of PRRSV viral particles in the lung, as previously described (95). The ORF7 RT-PCR product from both 3249 and Lena strains was firstly precipitated in ethanol and purified using ExoSAP-ITTM (Thermo Fisher Scientific). The purified products were quantified using Nanodrop 2000 (Thermo Fisher Scientific). Serial 10-fold dilutions of 3249 or Lena ORF7 RT-PCR products (ranging from 10 8 to 10 2 genomic copies/ mL) were used as standards to determine the limit of detection (1 copy/mL) and PCR efficiency (99%). Then, viral load for either Lena-or 3249-infected piglets was determined by RT-qPCR using VetMAX PRRSV EU/NA 2.0 kit (Thermo Fisher Scientific).
Flow cytometry (FCM) staining and analysis. BAL samples for FCM analysis and sequencing were collected at the time of necropsy from the left lung of each piglet (40). Freshly isolated cells from BAL were adjusted to 1.5 Â 10 6 cells per sample in a final volume of 200 mL. Then, BAL cells were stained for CD163 (clone 2A10/11, IgG1, 10 mg/mL; Bio-Rad Laboratories, S.A., Alcobendas, Madrid, Spain). After washing, a second incubation step with a fluorochrome-labeled isotype-specific secondary antibody (Alexa Fluor 647 goat anti-mouse IgG1, 7 mg/mL; Invitrogen, Carlsbad, CA, USA) in combination with Live/Dead Fixable Aqua Dead Cell Stain (Invitrogen) was performed. Following surface labeling, cells were fixed and permeabilized with methanol overnight at 220°C (VWR International, Llinars del Vallès, Barcelona, Spain). FCM analysis was performed on a FACSCanto II (BD Biosciences, NJ, USA) recording BAL cells, RNA isolation, and sequencing. RNA seq analyses were performed on representative piglets (3 animals/group/time point) over time points (n = 45). These animals were selected according to their clinical signs, gross and microscopic lung lesion scores, as well as viral load. Total RNA was isolated from 1.5 Â 10 6 cells per sample re-suspended with 1 ml of TRIzol reagent (Thermo Fisher Scientific), and then treated with Turbo DNA-free Kit (Thermo Fisher Scientific) to remove traces of contaminating genomic DNA (gDNA) following the manufacturer's guidelines. RNA integrity number (RIN) was assessed in the Bioanalyzer 2100 system (Agilent Technologies, CA, USA). All RNA samples used in the present study showed absence of gDNA and RIN values over 8.5.
Library preparation and sequencing were performed at the Functional Genomics Core of the Institute for Research in Biomedicine (IRB Barcelona, Spain). The cDNA libraries were prepared with 500 ng of total RNA for each individual sample using the Illumina TruSeq Stranded mRNA Sample Prep Kit (Illumina, Inc., San Diego, CA, USA). Libraries were quantified with Qubit dsDNA HS assay (Thermo Fisher Scientific) and quality was assessed by Agilent 2100 Bioanalyzer (Agilent Technologies). Subsequently, the indexed libraries were sequenced on a HiSeq2000 device (Illumina Inc.), and paired-end reads (2 Â 75 base pair) were generated. The raw RNA-seq data set was deposited at NCBI Sequence Read Archive database with the accession number PRJNA704925.
RNA sequencing data processing, time-series, and DEG analysis. Bioinformatics analysis was carried out by the Andalusian Bioinformatics Platform of the University of Málaga. Raw data were firstly processed using the in-house customizable pre-processing pipeline SeqTrimNext (96). In this step, clean data (clean reads) were obtained by removing contaminants, sequencing adapters, PolyA/PolyT tails, and short (, 17 nucleotide) and bad quality reads (Phred score , 20); thus, all the downstream analyses were based on this high-quality clean data set. Then, clean reads were mapped to the reference porcine transcriptome (Sus_scrofa. Sscrofa11.1.cdna, with 46.076 unigenes) using Bowtie (v2.2.3) and Samtools (v1.9). The reads count of each transcript was extracted using the python script sam2counts.py (v0.91) (Buffalo, 2010).
Since our data were collected at different time points, a time-series analysis was conducted on clean data using MaSigPro R package (97) to identify clusters of genes with significant temporal expression changes related to 3249 and Lena strain infection over time. Afterwards, DEgenes Hunter R pipeline and DeSeq2 R package (98) were used to identify differentially expressed genes (DEGs) between both infected groups and control group at each time point (1, 3, 6, 8, and 13 dpi). The P values were adjusted using Benjamini and Hochberg's correction for controlling the false discovery rate (FDR). Genes with an FDR , 0.05 and an absolute log 2 fold change (log2FC) $ 1 were assigned as DEGs. The DEGs of Lena-and 3249-infected groups' data set for each time point were visualized as a volcano plot by using GraphPad Prism (GraphPad Prism software v8.0), underlining the top 5 DEGs with a higher log2FC. Overlapped DEGs at 1, 3, 6, 8, and 13 dpi from Lena-or 3249-infected animals were identified and presented as a Venn diagram using an online tool (http://bioinformatics.psb.ugent .be/webtools/Venn/).
Gene ontology (GO) and pathway enrichment analysis. In order to explore GO of DEGs (http:// geneontology.org), functional enrichment analyses were conducted using ClueGO (v2.3.3) and CluePedia (v1.3.3), plugins for Cytoscape (v3.8, http://cytoscape.org/) describing biological processes (BPs) and immune system processes (ISPs) GO categories. ClueGO determines the distribution of the DEGs for various GO terms and pathways, generating a functionally grouped GO annotation network (99,100). The P value was calculated using right-sided hypergeometric tests and the Benjamini and Hochberg's correction for multiple testing (FDR , 0.05). This FDR threshold, together with a high kappa value (0.4), enabled us to precisely select significantly enriched and highly connected GO terms. The most significant term defines the name of the group. In order to avoid over-interpretation of data, a minimum of 3 genes were considered to evaluate the relevance of selected pathways.
Construction of protein-protein interaction (PPI) network and screening of Hub genes. The Search Tool for the Retrieval of Interacting Genes (STRING) database was applied to predict PPI and construct a PPI network of selected DEGs (101). Using the STRING database, DEGs with a score $ 0.4 were chosen to build a network model visualized by Cytoscape (102). Molecular complex detection (MCODE) is a Cytoscape plugin used to identify the finest PPI sub-network modules (102). The Hub genes are defined as genes with the highest degree of connectivity in the key module. Since the biological networks usually are heterogeneous, it seems to be reasonable to use more than one method for identifying Hub nodes (103); thus, Maximal Clique Centrality (MCC) and Density of Maximum Neighborhood Component (DMNC) algorithms were calculated for each node by CytoHubba plugin in Cytoscape (103). The genes with the top 10 MCC and DMNC values were considered as Hub genes.
Real-time quantitative reverse transcriptase PCR (RT-qPCR). To verify the major results drawn from RNA-seq and Hub genes analysis, the expression levels of a panel of 12 identified Hub genes were performed by RT-qPCR. cDNA was synthesized by reverse transcription using iScript cDNA Synthesis Kit (Bio-Rad) from 1 mg of the original RNA sample used for the RNA-seq according to manufacturer's guidelines. Subsequently, amplifications were run in triplicate using iTaq Universal SYBR Green Supermix (Bio-Rad) on the MyiQ2 Two-Color Real-Time PCR Detection System (Bio-Rad). For each reaction well, 50 ng of cDNA from each animal and 0.5 mM each primer were used. GAPDH was chosen as reference gene, and the relative expression level of each Hub gene was calculated by the 2 -DDCt method (104). The primers set used for RT-qPCR are listed in Table 5.