An Acute Stress Model in New Zealand White Rabbits Exhibits Altered Immune Response to Infection with West Nile Virus

The immune competence of an individual is a major determinant of morbidity in West Nile virus (WNV)-infection. Previously, we showed that immunocompetent New Zealand White rabbits (NZWRs; Oryctolagus cuniculus) are phenotypically resistant to WNV-induced disease, thus presenting a suitable model for study of virus-control mechanisms. The current study used corticosteroid-treated NZWRs to model acute “stress”-related immunosuppression. Maximal effects on immune parameters were observed on day 3 post dexamethasone-treatment (pdt). However, contrary to our hypothesis, intradermal WNV challenge at this time pdt produced significantly lower viremia 1 day post-infection (dpi) compared to untreated controls, suggestive of changes to antiviral control mechanisms. To examine this further, RNAseq was performed on RNA extracted from draining lymph node—the first site of virus replication and immune detection. Unaffected by dexamethasone-treatment, an early antiviral response, primarily via interferon (IFN)-I, and induction of a range of known and novel IFN-stimulated genes, was observed. However, treatment was associated with expression of a different repertoire of IFN-α-21-like and IFN-ω-1-like subtypes on 1 dpi, which may have driven the different chemokine response on 3 dpi. Ongoing expression of Toll-like receptor-3 and transmembrane protein-173/STING likely contributed to signaling of the treatment-independent IFN-I response. Two novel genes (putative HERC6 and IFIT1B genes), and the SLC16A5 gene were also highlighted as important component of the transcriptomic response. Therefore, the current study shows that rabbits are capable of restricting WNV replication and dissemination by known and novel robust antiviral mechanisms despite environmental challenges such as stress.


Introduction
West Nile virus (WNV) is a mosquito-borne flavivirus that has the potential to cause fatal encephalitis in humans, horses, birds and several other animals. WNV is a pathogen of global

Virus Kinetics in Untreated & Dexamethasone-Treated WNV-Infected Rabbits
Contrary to expectations, the viremia in the dexamethasone-treated group was significantly lower on day 1 p.i. compared to the mock-treated group. This was consistent across the two methods of virus titration, plaque assay and tissue culture infective doses at 50% cut-off (TCID50) (Figure 2). There was no statistically significant difference in viral load in the left popliteal lymph node draining the site of inoculation (hereafter referred to as the draining lymph node) between the two groups ( Figure 3). This was also consistent across the three methods of virus titration performed (plaque assay, TCID50, and qRT-PCR). By immunohistochemistry (IHC), we assessed the presence of WNV-NS1 antigen in the draining lymph node and other lymphoid and non-lymphoid organs (Table 1) as an indication of the extent of virus dissemination to distant sites. While in most rabbits viral antigen was largely restricted to the draining lymph node, one dexamethasone-treated rabbit (2105) had detectable antigen in several of the distant organs (Table 1). Another rabbit of the same group also

Virus Kinetics in Untreated & Dexamethasone-Treated WNV-Infected Rabbits
Contrary to expectations, the viremia in the dexamethasone-treated group was significantly lower on day 1 p.i. compared to the mock-treated group. This was consistent across the two methods of virus titration, plaque assay and tissue culture infective doses at 50% cut-off (TCID 50 ) ( Figure 2). There was no statistically significant difference in viral load in the left popliteal lymph node draining the site of inoculation (hereafter referred to as the draining lymph node) between the two groups ( Figure 3). This was also consistent across the three methods of virus titration performed (plaque assay, TCID 50 , and qRT-PCR). By immunohistochemistry (IHC), we assessed the presence of WNV-NS1 antigen in the draining lymph node and other lymphoid and non-lymphoid organs (Table 1) as an indication of the extent of virus dissemination to distant sites. While in most rabbits viral antigen was largely restricted to the draining lymph node, one dexamethasone-treated rabbit (2105) had detectable antigen in several of the distant organs (Table 1). Another rabbit of the same group also had detectable antigen Pathogens 2019, 8, 195 5 of 34 in the spleen, in addition to the draining lymph node (Table 1). The majority of NS1-positive cells were large pleomorphic leukocytes, with morphology consistent with dendritic cell and macrophage-like cells ( Figure 4). By day 7 p.i., no infected cells were detected by IHC in the draining lymph node, despite a low viral titer on virus isolation and qRT-PCR ( Figure 3). This is likely explained by the difference in the sensitivity of the assays, where the amount of tissue used in the homogenates for virus isolation/qRT-PCR provided higher sensitivity than that represented by the 4 µm thick tissue sections used for IHC. had detectable antigen in the spleen, in addition to the draining lymph node (Table 1). The majority of NS1-positive cells were large pleomorphic leukocytes, with morphology consistent with dendritic cell and macrophage-like cells ( Figure 4). By day 7 p.i., no infected cells were detected by IHC in the draining lymph node, despite a low viral titer on virus isolation and qRT-PCR ( Figure 3). This is likely explained by the difference in the sensitivity of the assays, where the amount of tissue used in the homogenates for virus isolation/qRT-PCR provided higher sensitivity than that represented by the 4 µ m thick tissue sections used for IHC.   had detectable antigen in the spleen, in addition to the draining lymph node (Table 1). The majority of NS1-positive cells were large pleomorphic leukocytes, with morphology consistent with dendritic cell and macrophage-like cells ( Figure 4). By day 7 p.i., no infected cells were detected by IHC in the draining lymph node, despite a low viral titer on virus isolation and qRT-PCR ( Figure 3). This is likely explained by the difference in the sensitivity of the assays, where the amount of tissue used in the homogenates for virus isolation/qRT-PCR provided higher sensitivity than that represented by the 4 µ m thick tissue sections used for IHC.   . Viral load in the left popliteal lymph node (draining the footpad injection site) in untreated and dexamethasone pre-treated rabbits, measured by three different approaches: (A) plaque assay, (B) TCID 50 assay, (C). qRT-PCR. LOD = limit of detection. No statistical significance was detected when groups were compared by two-way ANOVA. Each point represents the viral load of each of the rabbits in each group. The error bars represent the standard error of the mean and the middle bar represents the mean. Note that titers below the LOD are represented as zero.  the rabbits in each group. The error bars represent the standard error of the mean and the middle bar represents the mean. Note that titers below the LOD are represented as zero.
Dex, dexamethasone-treated group; PBS, mock-treated group; L. Ln Popl, left popliteal lymph node; R. Ln Popl, right popliteal lymph node; Ln. Mes, mesenteric lymph node; 0, no antigen detected; +, low antigen load; ++, moderate antigen load; +++, heavy antigen load.   Samples from lymphoid organs (thymus, spleen and lymph nodes), injection site in the footpad, brain (multiple sites) and liver were examined by microscopy. The most notable changes were present in the left popliteal lymph node-the lymph node draining the injection site-where very prominent hypertrophy of the post-capillary venules was evident from day 1 p.i. in the mock-treated animals and from day 3 p.i. in the dexamethasone-treated animals. Development of secondary follicles was also evident by day 7 p.i. for both treated and untreated animals. In contrast, a decrease in the thymic cortex, similar to that seen in the preliminary study (Supplementary Figure S3) was evident at day 1 and 3 p.i. in the dexamethasone treated animals, but had resolved by day 7 p.i. This change was not evident in any of the mock-treated animals at any time point p.i. Similarly, a slight decrease in the white pulp of the spleen was evident in the day 1 and 3 p.i. dexamethasone-treated animals, but this was followed by reconstitution and hypertrophy by day 7 p.i. Mild inflammation in the injection site was only detected in one mock-treated animal on day 7 p.i., but it cannot be excluded that similar reactions were present in other animals, but missed as only a few 4 µm sections were examined for each animal.
Immunolabeling for T lymphocytes (CD3), B lymphocytes (CD79a) and myeloid-lineage cells (macrophages and heterophils; Mac387) did not reveal any major or consistent differences in relative numbers and distribution between the mock-and dexamethasone-treated, WNV-infected animals for any of those cell populations at any time point (data not shown).

Antibody Responses
Sera from animals culled on day 1 and 7 p.i. were tested by blocking-ELISA and virus neutralization test (VNT) for antibodies specific for WNV. All day 1 p.i. sera were negative in both assays, while all five (three dexamethasone-treated and two PBS-treated) of the day 7 p.i. sera were positive for flavivirus E-protein specific antibodies in the blocking-ELISA. Of these, three had WNV-neutralizing activity: one dexamethasone-treated animal had a titer of 1:40 and the two PBS-pretreated animals had titers of 1:20 and 1:40, respectively. All day 7 p.i. sera were negative for WNV-NS1 specific antibodies in the blocking-ELISA, suggesting that virus replication was quickly contained before release of substantial amounts of NS1-protein.

Differential Gene Expression Analyses of RNA Extracted from Fraining Lymph Nodes
In our previous studies, we demonstrated the draining popliteal lymph node was the first tissue site with detectable tissue viral load [20,21], and thus, a suitable site for study of early antiviral responses. Therefore, RNA was extracted from the draining lymph node and RNAseq was performed in order to examine potential differences in the early immune response between dexamethasone and mock-treated rabbits. Table 2 summarizes the denotation assigned to each of the groups for downstream differential gene expression and gene ontology (GO) analyses. Two mock-treated, non-infected animals from the pilot study were also included (denoted C0 in Table 2). Of the 16 WNV-infected rabbits (nine dexamethasone-treated and seven mock), one RNA sample (from animal 2103) was excluded for RNAseq, due to an unworkably low RNA integrity number (RIN score = 3.3). The RIN score for the remaining samples ranged from 4.9-6.8, and an adjustment to the fragmentation time was used during the sequencing library preparation in order to account for the generally low RIN score (see Methods section for more details).
RNAseq of the 17 RNA samples yielded an average of 29067807 reads (range: 22019057-37133333), with all fastq files passing sequence quality check using the Galaxy application, FastQC (Galaxy version 0.72). Pairwise comparisons of transcriptome were performed between dexamethasone-treated WNV-infected groups (T) and untreated-uninfected controls (C); between untreated WNV-infected groups (N) and C; and between T and N. We employed a highly stringent protocol for differential gene expression analysis that used three R-based packages (edgeR, DEseq2, and Voom-limma) under the Galaxy tool "Differential _Count models using BioConductor packages" (Galaxy Version 0.28). Only differentially expressed genes (DEGs) that intersected across all three packages were used for downstream analyses.
The number of intersecting DEGs was highest on day 1 p.i., irrespective of the comparisons ( Figure 5A,D,G), suggesting early transcriptional activity is associated with WNV NSW2011 infection. The direct comparisons between treated and untreated WNV infected groups (T-N) showed only low numbers of DEGs, the highest number being observed on day 1 p.i. (132 DEGs, Figure 5G). The transcriptome of the draining lymph nodes of rabbits treated with dexamethasone and infected with WNV NSW2011 was generally comparable with that of mock-treated rabbits, at least for day 3 and 7 p.i. In order to investigate the transcriptional dynamics of genes associated with WNV NSW2011 restriction, pairwise differential gene expression analyses were performed for each of the groups, comparing against the mock-treated uninfected controls. Consistent with the T-N comparisons, the number of DEGs between T-C and N-C comparisons was generally comparable for day 1 and 3 p.i. (1410 versus 1937 on day 1 p.i.; 698 versus 462 on day 3 p.i.; Figure 5D,E,G,H). A notable difference was observed on day 7 p.i., however, where only the N-C comparison yielded DEGs (530 DEGs; Figure 5I). No DEG was identified in the T7C0 comparison ( Figure 5F).

Gene Ontology Analyses
The intersecting DEGs from each analysis were assessed for enrichment of functional annotation, using the gene ontology (GO) tool, DAVID 6.8. An adjusted p-value threshold of 0.05 using the Benjamini correction was set as the cut-off for GO term enrichment. Enriched GO terms from each of the pairwise comparisons were then analyzed for commonalities and differences using Venn diagrams (Supplementary Figure S4). Across all the comparisons against the control group that yielded DEGs (i.e. T1C0, T3C0, N1C0, N3C0, and N7C0, Figure 5), two GO terms "defense response to virus" (GO: 0051607) and "extracellular exosome" (GO: 0070062) were consistently enriched, . The days in the column headings refer to the day p.i. The numbers within circles represent the number of DEGs that were common (intersects) between each of the R package analyses or unique to each of the analyses (non-overlapping regions). The number external to the three circles represents the number of genes that were not differentially expressed. Each of the comparisons is denoted by a code (in brown) where the alphabets represent the treatment/infected status, and the number represents the day post-infection. T, dexamethasone treated WNV-infect groups; N, mock-treated WNV-infected groups; C, untreated uninfected group.

Gene Ontology Analyses
The intersecting DEGs from each analysis were assessed for enrichment of functional annotation, using the gene ontology (GO) tool, DAVID 6.8. An adjusted p-value threshold of 0.05 using the Benjamini correction was set as the cut-off for GO term enrichment. Enriched GO terms from each of the pairwise comparisons were then analyzed for commonalities and differences using Venn diagrams (Supplementary Figure S4). Across all the comparisons against the control group that yielded DEGs (i.e. T1C0, T3C0, N1C0, N3C0, and N7C0, Figure 5), two GO terms "defense response to virus" (GO: 0051607) and "extracellular exosome" (GO: 0070062) were consistently enriched, highlighting the importance of these two annotations in the virus challenged groups regardless of dexamethasone treatment and time post-infection.
On day 1 p.i. and across the treated and mock comparisons against controls, other commonly enriched GO terms relevant to the context of virus control included "negative regulation of viral genome replication" (GO:0045071), "positive regulation of JNK cascade" (GO:0046330) and terms associated with various cellular compartments (Table 3, Supplementary Figure S4). Unique to the DEGs from the T1C0 comparison, there were a further 20 enriched GO terms, many of which were associated with immune response, in particular ones relating to response to lipopolysaccharide, ATP binding, innate immunity, tumor necrosis factor production, macrophage/monocyte chemotaxis, nuclear factor-κB (NF-κB) activity, interleukin-1β (IL-1β) secretion, integrin signaling and IL-8 production ( Table 3).
On day 3 p.i., intersecting GO terms across treated and mock-treated WNV-infected rabbits that were within the theme of antiviral activity and immune response included "double-stranded RNA binding" (GO:0003725), "immune response" (GO:0006955), and "ATP binding" (GO:0005524), as well as "defense response against virus" and "extracellular exosome" as mentioned above. On day 7 p.i., only the comparison N7C0 yielded DEGs. GO terms enriched within the N7C0 DEG list that were relevant for further analyses included "defense response to virus" and "extracellular exosome", which were enriched in all of the other comparisons as mentioned above.
It should be noted that for the following downstream analyses, we excluded several enriched GO terms relating to translation and various intracellular compartments, as their links to antiviral mechanisms seems to be indirect, and it is not expected that they be relevant in resistance to infection. Table 3. Summary of selected enriched gene ontology (GO) terms and the time post-infection at which they are enriched. Terms are selected on the criteria that they are relevant to antiviral processes and are enriched on two or more occasions.

Antiviral DEGs in Rabbits Infected with WNV NSW2011
In order to identify and compare which DEG contributed most to enriched GO terms associated with antiviral immune responses on day 1 and 3 p.i., DEG lists from each of the GO terms outlined in Table 3 (excluding ones within the "extracellular exosome" term which were analyzed separately) were compared for similarities and differences ( Figure 6A). A total of 209 DEGs were associated with the antiviral GO terms; 14 of which were common to the day 1 and 3 p.i. analyses regardless of dexamethasone-treatment ( Figure 6A,B); 22 were uniquely common to only the comparisons on day 1 p.i. (Figure 6A,C); and 31 were uniquely common to only the comparisons on day 3 p.i. (Figure 6A,D). Notably, there were 86 antiviral DEGs that were only expressed in the T1C0 comparison, indicating a higher number of unique DEGs in dexamethasone-treated rabbits on day 1 p.i., as compared to the number of unique DEGs in the remaining groups ( Figure 6E). DEGs unique to the comparisons, N1C0, T3C0 and N3C0, are illustrated in Figure 6F,G,H, respectively.
To assess the magnitude of up/downregulation of the antiviral DEGs, heatmaps of the DEG fold change (log 2 FC) were generated ( Figure 6B-H, and Supplementary Figure S5) based on subsets of the various intersects of interest in the Venn diagram ( Figure 6A). Hierarchical clustering of the logFC was performed for each heatmap in order to identify top upregulated and downregulated clusters of DEGs, which are summarized in Tables 4-6, and Supplementary Table S1. To assess the magnitude of up/downregulation of the antiviral DEGs, heatmaps of the DEG fold change (log2FC) were generated ( Figure 6B-H, and Supplementary Figure S5) based on subsets of the various intersects of interest in the Venn diagram ( Figure 6A). Hierarchical clustering of the logFC was performed for each heatmap in order to identify top upregulated and downregulated clusters of DEGs, which are summarized in Tables 4-6, and Supplementary Table S1.      Figure 6B).   Figure 6C,E,F).

Day 1 p.i. Antiviral Genes Affected by Dexamethasone Treatment Antiviral DEGs Irrespective of Dexamethasone Treatment Expressed in Treated Group Expressed in Mock-Treated Group
Up   Figure 6D,G,H).

Antiviral Genes Across Day 1 and 3 p.i. Unaffected by Dexamethasone Treatment
Amongst the 14 antiviral DEGs that were similarly differentially expressed across all comparisons on day 1 and 3 p.i., the chemokine genes CXCL9 and CXCL10, and the IFN-stimulated gene (ISG), RSAD2 (also known as viperin), were clustered as DEGs with the highest logFC on day 1 and 3 p.i. ( Figure 6B). The second cluster to follow this included the ISGs: OASL, OAS2 and IFIT5 ( Figure 6B). Of lower magnitude but nevertheless upregulated in their expression across the day 1 and 3 p.i. groups were the pattern recognition receptors, Toll-like receptor (TLR)-3 and transmembrane protein (TMEM)-173 (also known as stimulator of IFN genes [STING]); the IFN signaling molecule, IFN regulatory factor (IRF)-1 and NOD-like receptor (NLR) family CARD domain containing 5 (NLRC5); and further ISGs such as OAS3, and RNaseL.
The persistent expression of these genes across the groups during these early time-points highlights their significance in a conserved, early antiviral response to WNV NSW2011 infection, irrespective of dexamethasone treatment.
Unique to dexamethasone treated rabbits on day 1 p.i., there were upregulated expression of the T-helper type I response cytokine IL-23A; chemokine receptor CCR1 and ligand CCL2; and a transmembrane ATP-binding protein, multidrug resistance-associated protein 1 (LOC100346553; also known as ABCC1). Notably, also unique to treated rabbits on day 1 p.i., there were several downregulated genes, in particular ones with nucleotide/ATP binding function (ABCA8, ABCC5, LOC100358984, NPR1, CNNM2, MYO6, KATNAL1, FGFR1, and EPHB3); some of which are transmembrane receptors (e.g. TGFBR3, FGFR1 and EPHB3) ( Table 5 and Supplementary Figures  S3 and S4). Unique to mock-treated rabbits on day 1 p.i., the expression of growth arrest and DNA damage inducible gamma (GADD45G) gene was upregulated (Table 5).

Antiviral Genes Uniquely Expressed on Day 3 p.i.
On day 3 p.i., common across the treated/mock-treated analyses, the interferon-inducible T cell chemokine, CXCL11, was the top upregulated gene ( Figure 6D), followed by few genes that encode proteins with double-stranded RNA binding function, such as DExH-box helicase 58 (DHX58), DExD/H-box helicase 58 (DDX58; also known as retinoic acid-inducible gene I [RIG-I]), and IFN induced with helicase C domain 1 (IFIH1; also known as melanoma differentiation-associated protein 5 [MDA5]). The latter two are well-known pattern recognition receptors (PRRs) that detect cytosolic dsRNA. The tumor necrosis factor superfamily member 10 (TNFSF10, also known as TRAIL) and the ISG, transporter 1 ATP binding cassette subfamily B member (TAP1) were also in the second top cluster of DEGs (Table 6). Concurrently, all day 3 p.i. infected groups had downregulated expression of the chemokine CCL17 and its receptor CCR8, and ATP-binding proteins Alpha Kinase 2 (ALPK2) and Kinesin Family Member 5C (KIF5C; Table 6).
Unique to the dexamethasone treated group on day 3 p.i., expression of IFN-γ was upregulated to the highest extent ( Figure 6G and Table 6). Chemokine genes, lymphotactin and CCR2, were also upregulated along with HLA class I histocompatibility antigen B-7 alpha chain. Most prominent downregulated expression was observed for HLA class II histocompatibility antigen DRB1-4 beta chain, and the chemokine ligand CCL22 and its receptor CCR4 (Table 6).
On the other hand, the top upregulated DEGs in the mock-treated rabbits included the GTP-binding tubulin alpha-1B chain-like molecule, the pore-forming cytotoxic perforin 1 (PRF1), and heat shock protein family A (Hsp70) member 8 (HSPA8; Table 6). To a lesser extent, the second cluster of top upregulated genes included the ATP-binding, actin-related protein 3-like protein (LOC100358336), the immunomodulatory FKBP Prolyl Isomerase 4 (FKBP4), and the ATP-binding cassette transporter ABCA1 (Table 6). Specific to the mock-treated group, the expression of DNA topoisomerase II beta (TOP2B) and diacylglycerol kinase eta (DGKH) genes were downregulated to the highest extent (Table 6).
2.2.11. Antiviral Genes Expressed on Day 7 p.i. On day 7 p.i., the draining lymph node transcriptome of the dexamethasone-treated group was comparable to that of the baseline control, yielding no DEG ( Figure 5F). On the other hand, 530 DEGs were identified in the untreated group to control comparison (N7C0; Figure 5I). This suggests that the treated group returned to baseline gene expression faster than mock-treated rabbits, despite similar levels of viral load in the lymph node by this time (Figure 3). Similar to the comparisons in the earlier time-points, the two GO terms "defense response to virus" and "extracellular exosome" were enriched (Supplementary Figure S4C). Within the top upregulated cluster of DEGs in the former term, DEGs that have been featured in earlier comparisons included ISG15 and CXCL9 (Supplementary Figure S5). The chemokine CCL5 and cytidine/uridine monophosphate kinase 2 (CMPK2) were also among the top clustered DEGs (Supplementary Figure S5). Notably, ATP-binding cassette subfamily A member 3 (LOC100353012) was the top downregulated DEG (Supplementary Figure S5).

Extracellular Exosome
Belonging in the GO category of cellular compartment, the commonly enriched GO term "extracellular exosome" was analyzed separately to the antiviral GO terms as outlined above (Figure 7). A total of 386 DEGs were associated with this cellular compartment GO term across the day 1 and 3 p.i. comparisons, 36 of which were common across these comparisons ( Figure 7A). The top upregulated DEG across all groups was complement 2 (C2; Figure 7B,C). The second cluster of top upregulated DEGs included tissue inhibitor of metalloproteinases 1 (TIMP1), galectin 3 binding protein (LGALS3BP), tumor necrosis factor ligand superfamily member 10 (TNFSF10) and carbonic anhydrase 2 (CA2; Figure 7B,C). Of note, both TIMP1 and TNFSF10 are genes involved in cell proliferation and apoptotic processes, while LGALS3BP has been implicated in NK and lymphokine-activated killer (LAK) cell mediated immune responses. TNFSF10 was also highlighted as one of the top upregulated DEGs on day 3 p.i., irrespective of dexamethasone treatment ( Table 6). The top cluster of downregulated DEGs included myosin 1B (MYO1B), centlein (CNTLN), sorting nexin 29 (SNX29) and Fc fragment of immunoglobulin E (IgE) receptor II (FCER2; Figure 7B,C).   Figure 7B).

Unsupervised Analysis of Gene Expression
While gene ontology analyses allowed us to identify groups of genes of interest by using gene functional annotations in a curated database, for example the DAVID 6.8 database in our case, a major limitation to this approach is that genes that are not yet included in the database may not be highlighted as an integral part of the transcriptome. While an extensive unsupervised analysis of gene expression is beyond the scope of this study, we have conducted two simple and common approaches to identifying further genes of interest without a priori assumptions. Firstly, we conducted hierarchical clustering of DEGs from each group analysis based purely on their logFC, followed by principle component analysis (PCA; Figure 8A-C).
For hierarchical clustering, we first subsetted the matrix of DEG logFC to only genes that were differentially expressed across all the pairwise comparisons (i.e. N1C0, N3C0, N7C0, T1C0 and T3C0). This allowed us to identify DEGs that were highly up-or downregulated across the comparisons. A total of 78 genes were differentially expressed across all the comparisons. Of these, a notable top upregulated cluster is identified, containing only one gene, solute carrier family 16 member 5 (SLC16A5; ENSOCUG00000016231), which was highly upregulated across all five comparisons, irrespective of dexamethasone treatment and time post-infection ( Figure 8A).
The second top cluster included some genes already identified as key genes across all the comparisons from the GO analyses in Table 4. These included CXCL10 and CXCL9 ( Figure 8A). In addition to this, within the second cluster, there were also several genes that have not been previously highlighted in the GO analyses. These included ubiquitin specific peptidase 18 (USP18), complement factor B (CFB), IFIT-1-like protein (LOC100342826), 2'-5'-oligoadenylate synthase-like protein (LOC100356242), and a novel unannotated gene (ENSOCUG00000004197 with closest annotated orthologue being IFIT1B in Gibbon [ENSNLEG00000018636; 68.03% query % identity; 66.88% target % identity], and human orthologue IFIT1B [ENSG00000204010]; Figure 8A). On the other hand, the top downregulated gene cluster contained one gene, the chemokine receptor, CX3CR1, which was most downregulated on day 1 p.i., and to lesser extent, on the subsequent days ( Figure 8A).  The subsetted matrix containing 78 DEGs used for hierarchical clustering was also used for PCA, which aimed at identifying DEGs that have the largest effect in explaining the variation of logFC across the different comparisons, N1C0, N3C0, N7CO, T1C0 and T3C0. This allowed identification of DEGs that were the "most different" across the comparisons. The PCA plot in Figure 8B shows that the different comparisons cluster most closely within their corresponding time p.i., with the greatest degree of separation between treatment versus mock-treatment on day 1 p.i. Again, similar to the unsupervised hierarchical clustering, some DEGs already highlighted as important from GO analyses were also present in the top 10 DEGs of the first principal component (PC1), which explained 97.3% of the variation of logFC in the data matrix ( Figure 8C and Table 7). These genes included CXCL10, RSAD2 and DDX58 (Table 7). Common to the top clusters from the unsupervised hierarchical clustering, IFIT2, CFB, IFIT1-like protein (LOC100342826) and the novel gene ENSOCUG00000004197 (human orthologue of IFIT1B) were also identified as some of the top ten contributors to the first principal component (Table 7). Unreported in the analyses prior to this, guanylate-binding protein 1 (GBP1; LOC100358539), IFIT3 (ENSOCUG00000029154) and an additional novel gene (ENSOCUG00000015849; Human orthologue: HERC6 [ENSG00000138642]) were identified as important contributing genes to principal component 1 ( Table 7).

Discussion
We have presented data for our studies examining the effect of acute "stress", in the form of a single systemic injection of the synthetic corticosteroid drug dexamethasone, on viral replication and the antiviral innate immune response following intradermal WNV infection in young rabbits. While virus challenge was done at the time of maximal effect of dexamethasone treatment on immune parameters, the infection, contrary to our hypothesis, resulted in significantly lower viremia in treated rabbits relative to that in the untreated controls. This was likely explained by a different initial innate immune response at the first tissue site of virus replication, the draining lymph node, as demonstrated by the difference in the composition of DEGs within commonly enriched antiviral GO terms in a treatment and time dependent manner.

• Gene ontology
Through the GO analysis, we identified several GO terms with relevance to antiviral immune response that were enriched in each of the pairwise comparisons against control. Enriched biological process and molecular function GO terms "defense response to virus", "ATP binding", "negative regulation of viral genome replication", "positive regulation of JNK cascade", "double-stranded RNA binding" and "immune response" highlighted the early activation of immune processes in response to the virus challenge. Some of these terms were enriched in a time-dependent manner, in particular, the term "double-stranded RNA binding" was commonly enriched only on day 3 p.i. This coincided with the peak of viral load in the draining lymph node, and therefore, suggests that the expression of molecules for "double-stranded RNA binding" was likely important for restricting virus replication, consistent with the reported literature on WNV-induced immune response [28].
• IFN-I At the gene level, we pooled DEGs that were associated with enriched GO terms relevant to antiviral immune responses and examined their level of fold change from control. This analysis showed that the IFN response featured prominently, in particular on day 1 p.i., where IFN-β1 was commonly upregulated regardless of dexamethasone treatment. Furthermore, depending on whether the rabbits were treated with dexamethasone or not, different repertoires of IFN-α-21-like and IFN-ω-1-like genes at different loci were expressed on day 1 p.i. While it could be expected that IFN-α and -β would be an important part of the antiviral response, consistent with previous studies on the innate immune response of WNV, flavivirus and viruses in general [28], the finding of concurrent expression of IFN-ω subtypes alongside IFN-α/β is notable. IFN-ω has been reported to possess a comparable level of antiviral activity against flaviviruses such as hepatitis C virus in vitro [29], with increased in vitro antiviral activity when glycosylated, in particular against WNV and other flaviviruses [29]. Our study supports an important role for IFN-ω, perhaps in combination with IFN-α/β, in restricting early WNV replication in rabbits. The role of the various rabbit IFN-α and -ω subtypes in the antiviral action also requires further characterization. Understanding this may help elucidate the complex redundancy of the IFN system available for ancillary support against viral pathogens, when perhaps a particular arm of the immune system may be disturbed by environmental factors such as acute stress.
• IFN-stimulated genes The IFN system, most notably the type I IFNs, IFN-α and -β, initiates a first-line antiviral response upon viral infection, mediated by the induction of >300 ISGs, leading to the establishment of an antiviral state [30,31]. In the current study, a number of ISGs were highlighted as top upregulated DEGs either in a time/treatment-independent or dependent manner. Regardless of dexamethasone treatment or time post-infection, expression of ISGs of the OAS family, specifically OAS2, OAS3 and OASL, and its downstream RNaseL, were highly upregulated. Generally, ISGs of the OAS family has been well studied and are frequently highlighted in their ability to activate the latent RNase, RNaseL, which in turn degrades viral RNA [32]. Recently, OASL, which lacks the 2'-5' oligoadenylate catalytic activity, has been described as possessing the ability to boost the antiviral IFN-I response via enhancement of the sensitivity to activate RIG-I [32]. In other words, OASL has been shown to enhance the ability of RIG-I to detect low levels of viral RNA, levels which may be reflective of very early phase of viral replication. In the infected rabbits of the current study, expression of RIG-I (referred to as DDX58 in this manuscript) was upregulated on day 3 p.i., regardless of dexamethasone treatment. Given the overall limited viremia and virus spread in infected rabbits, along with high degree of upregulated expression of OASL and RIG-I, we propose that one of the major mechanisms of early restriction on viral replication in the draining lymph node was mediated by OASL and its enhancement of RIG-I's sensitivity to low levels of viral RNA. This in turn induced a robust IFN-I response, which was also evident.
Another family of ISGs highly featured in the antiviral response of the rabbits is the IFIT family. IFIT5 was consistently upregulated across all comparisons (Table 4), and from the unsupervised analyses, IFIT1-like, IFIT2, IFIT3 and an unannotated gene with putative IFIT1B function were identified as key genes in the antiviral defense against WNV NSW2011 . Individually, IFIT1 has been reported to disrupt the interaction of MITA, MAVS and TBK1, which then negatively regulates the cellular antiviral response. IFIT2 interacts with MITA, and induces apoptosis via the mitochondrial pathway that is induced by the innate immune response. IFIT3 bridges TBK1 to MAVS in mitochondria, which synergizes the activation of IRF3 and NF-κB to activate the immune response [33]. We propose that these characterized functions apply also to rabbits, although further investigation into the antiviral mechanisms of the novel IFIT1-like and the putative IFIT1B gene should be conducted. RSAD2, also known as 'virus-inhibitory protein, endoplasmic reticulum-associated, interferon-inducible' (viperin), is another important ISG that was consistently upregulated in infected rabbits in this study. RSAD2/viperin has many diverse functions, antiviral activity being one of these [34], and mediated by multiple pathways, depending on virus species [35,36]. In case of flaviviruses, including WNV, viperin catalyzes conversion of cytidine triphosphate, thereby generating 3'deoxy-3',4'-didehydro-CTP, which acts as a chain terminator for the RNA-dependent RNA polymerase [35,36]. Viperin also targets the flavivirus NS3 for proteasomal degradation [37] and may interfere with the assembly of the flavivirus capsid [38]. Specifically for WNV, viperin knockout mice were shown to be more susceptible to WNV challenges [39]. The same study also showed tissue and cell specificity for its antiviral function [39].
In a time-specific manner, ISG15 was commonly upregulated on day 1 p.i. regardless of dexamethasone treatment. ISG15, a 15kD protein containing two ubiquitin-like domains [40], forms conjugates with a diverse array of proteins, the specificity of which is in part derived from the temporal expression of the target proteins, in a process referred to as ISGylation. However, there is no evidence that ISGylation leads to proteasome mediated degradation of the target proteins, and thus, the exact molecular mechanism(s) of ISG15 antiviral effect remain inadequately understood and may vary for various viruses and hosts [40]. While ISG15 expression is upregulated in the context of in vivo infection with several flaviviruses [17,41], the antiviral role of ISG15 in WNV infection in vivo has been called into doubt: ISG15 -/mice did not appear to have increased susceptibility to WNV infection despite apparent in vitro antiviral activity of ISG15 towards WNV [40,42]. These conflicting findings need further dissection. Given its time-specific expression in our WNV NSW2011 -infected rabbits, perhaps ISG15's role in the overall antiviral defense may also be temporally specific.
Similarly, the expression of a recently characterized ISG, TAP1, was also upregulated in a time-specific manner, except this was observed on day 3 p.i. Traditionally, TAP1 was known for its role as a transporter of proteasome-degradation products to major histocompatibility complex-I (MHC-I) receptors inside the endoplasmic reticulum lumen [43,44]. Hence, it was a regulator of the adaptive immune response to viral infection. However, a study into the antiviral mechanisms against vesicular stomatitis virus infection suggested that TAP1 was an ISG with post-translational processing and degradation of viral protein [44]. Its role in WNV control is poorly understood. Our finding of its common expression on day 3 p.i. highlights a potentially important role and further characterization of this ISG is required.
• IFN-signaling It is well characterized that the IFN-I response is initiated by the pattern recognition receptors (PRRs) that detect pathogen associated molecular patterns (PAMPs), followed by intracellular signaling. In the current study, consistent with the findings of a robust IFN and ISG response, there was high expression of a variety of PRRs: TLR-3, TMEM173 (also known as STING), NLRC5, DDX58 (RIG-I) and IFIH1 (MDA5). These were upregulated on day 1 and/or 3 p.i. irrespective of dexamethasone treatment. TLR-3 also remained upregulated on day 7 p.i. in the untreated group.
TLR-3 binds double-stranded RNA within endosomes and has been heavily emphasized in its role in initiating an anti-WNV IFN-I response via either the IRF3 or IRF7 signaling pathway [45]. In our previous studies examining the expression of TLR genes in rabbit and equine peripheral blood mononuclear cells (PBMCs) infected with WNV, we observed similar upregulated expression of TLR-3 along with type I IFN and ISGs [46,47]. Therefore, observing upregulated TLR-3 expression in the rabbits of the current experiment further affirmed the importance of this PRR; however, more importantly, that its upregulated expression was maintained even when treated with dexamethasone. This suggests that despite acute stress, the rabbit innate immune system retains its ability to rapidly and persistently deal with WNV infection via TLR-3.
TMEM173 (also known as STING), like TLR-3, is a known initiator of the type I IFN response but unlike TLR-3, it is a PRR that detects cytosolic DNA [48]. Despite this, it has been implicated in the antiviral control of both RNA and DNA viruses [48,49]. The finding of TMEM173/STING transcript upregulation in our current study is consistent with published literature for WNV-induced innate immune response [48]. Using STING knockout mice, it was shown that STING was required for bone-marrow derived macrophages to control WNV replication and to initiate both an effective innate immune response and an appropriate T cell response in vivo in order to restrict WNV-associated neuropathology [48]. Moreover, degradation of STING by flavivirus NS2B3 protease may be one of the immunoevasive mechanisms displayed by these viruses, with the caveat that STING is not subject to cleavage in all species [49]. However, the question why and how a PRR, which senses cytosolic DNA, can be an integral part of the antiviral response to an RNA virus like WNV remains to be elucidated.
NLRC5 was another commonly upregulated gene in the infected rabbits, regardless of dexamethasone treatment. This NOD-like CARD domain-containing intracellular protein has recently been identified as a major histocompatibility complex class I (MHC-I) transactivator [50], but its association with flavivirus infection has apparently not been described. Guo et al. [51] previously showed that respiratory syncytial virus (RSV) infection upregulates MHC-I expression through the induction of NLRC5. In our study, the HLA-I antigen, B7-alpha chain, was concurrently upregulated on day 3 p.i. The link between NLRC5 and MHC-I expression may, therefore, be present in rabbits infected with WNV NSW2011 , and this interaction may play a key role in virus control, perhaps via the induction of an appropriate cell-mediated adaptive response by cytotoxic T cells. Notably, upregulated expression of PRF1 was observed on day 3 p.i. for untreated rabbits, which may support an important role of the cytotoxic T cell response. NLRC5 may also induce antiviral response against WNV NSW2011 via RIG-I signaling pathway, as seen in influenza virus study [52].
In the current study, both DDX58 (RIG-I) and IFIH1 (MDA5) expression were upregulated on day 3 p.i. The importance of DDX58 (RIG-I) was further highlighted by the unsupervised PCA. These two PRRs are well-characterized sensors of dsRNA and mediate antiviral IFN-I response via IRF3 and 7 [50]. However, it is notable that the expression of the regulator of both the RIG-I and MDA5 pathway, DHX58 (also known as LGP2), was also upregulated at the same time. This suggests that a fine-tuning regulatory mechanism within the RIG-I and MDA5 pathways is present in the rabbit's antiviral response against WNV NSW2011 .
Throughout the analyses of day 1 and 3 p.i., there were indications that adaptive immune response involving specific T cell subsets may be integral to the restricting WNV infection. As mentioned already, upregulated expression of PRF1 on day 3 p.i. of untreated rabbits hinted at a cytotoxic T cell or NK cell response. The early upregulation of IL23A, as well as IFNγ on day 3 p.i. in dexamethasone treated rabbits suggest a CD4 + T lymphocytic response (helper 17 and helper 1, respectively) was preferred in this group. The kinetics of leukocytes through the lymph node was also likely driven by the temporal pattern of chemokine ligand and receptor expression [53].

• Chemokines ligand and receptor expression
Mobilization of effector cells to sites of active infection is driven by, amongst others, chemokines. CXC type chemokines that lack the ELR-sequence, including CXCL9 and CXCL10, are potent attractants for activated T and B lymphocytes, monocytes, dendritic cells and NK cells [53,54]. Both of these chemokines bind to the chemokine receptor CXCR3, the signaling through which has been shown to be very important in many viral infections [54][55][56][57]. We previously demonstrated upregulation of CXCL10 in brain and draining lymph node following WNV infection in rabbits [20,21] and this was validated for the draining lymph node in this study in both dexamethasone treated and untreated rabbits. We also demonstrated equally potent upregulation of CXCL9 in a similar time-and treatment-independent manner ( Figure 6B, Table 4). CXCL11 was also highlighted as a top upregulated gene specifically on day 3 p.i. (Figure 6D), emphasizing the importance of CXC type chemokine response to restricting WNV replication in vivo.
Despite the upregulation of the CXC type chemokines, it is noteworthy that genes encoding other mononuclear leukocyte chemokines like CCL17 (irrespective of dexamethasone treatment) and CCL22 (specific to dexamethasone treated rabbits) were prominently downregulated, alongside their receptors CCR8 and CCR4, respectively. This potentially hints at a complex regulation of chemoattraction of specific mononuclear leukocytic subsets with specific activation state to the lymph node in order to achieve virus control.
• Extracellular exosome Alongside the use of biological process and molecular function GO terms, we included terms in the category of cellular compartment for the GO analyses. Notably, the "extracellular exosome" compartment was repeatedly enriched across all comparisons. Analysis into the composition of the top DEGs within this term revealed that the genes encoding C2, as well as TIMP1, LGAL3BP, TNFSF10 and CA2, were highly upregulated on day 1 and 3 p.i. regardless of dexamethasone treatment. Upregulation of C2 was the most consistently pronounced across these comparisons. This suggests an important role of complement in the antiviral response against WNV, potentially via extracellular exosomes. TIMP1 and TNFSF10 have been implicated in the control of cell proliferation and apoptotic processes [58,59]. These molecules in exosomes may potentially drive the formation of germinal centers in the lymph nodes for achieving a competent antibody response, which was apparent in these rabbits (Section 2.2.4). This exciting field of exosome-mediated antiviral mechanisms should be further investigated.
• Unsupervised analyses identified genes with little-known antiviral activity Apart from validating the importance of antiviral genes, such as the ISGs, DDX58 and the chemokines CXCL9 and 10, the unsupervised analyses also highlighted genes with little-known antiviral activity. Notably, SLC16A5 was consistently upregulated to the highest extent across all groups, regardless of dexamethasone treatment or time post-infection. This protein, also known as monocarboxylate transporter 6, aids in the transport of monocarboxylates (e.g., lactate, pyruvate), ketone bodies and branched-chain oxo acids across cell plasma membrane [60]. Nothing is known about its association with antiviral activity.
Another gene identified as important by our unsupervised approach was the unannotated gene with putative HERC6 function. HERC6 is a ubiquitin ligase and has been implicated in MHC class I mediated antigen processing [61]. Further characterization of this novel rabbit gene with putative HERC6 function is warranted.

• Corticosteroids and immunosuppression
The current study showed that the altered immune function associated with a single bolus injection of dexamethasone was insufficient for induction of severe disease upon WNV infection. This challenges the perception that corticosteroid administration or acute stress will compromise immune function to the point of increasing susceptibility to infectious diseases. It is evident from our study that the action of corticosteroids on the immune function is complex and the outcome of its use in the context of WNV infected patients may be dose-, frequency-and timing-dependent. Further dissection into these aspects by varying the corticosteroid dose, frequency and time of administration in relation to WNV infection is required.
In conclusion, acute stress as modeled by a single dose dexamethasone injection in rabbits does not induce a more severe disease phenotype when infected with WNV NSW2011 , contrary to our initial hypothesis. We have shown that despite acute stress, WNV-infected rabbits retained their ability to activate an early antiviral response, primarily via type I IFNs and in particular IFN-β1, IFN-α-21-like and IFN-ω-1-like subtypes, as well as the induction of a range of ISGs (OASs, IFITs, RSAD2/viperin, and ISG15). The involvement of IFN-ω in the anti-WNV immune response is a unique highlight of the rabbit model. We propose that this robust IFN-I response was likely initiated and sustained by ongoing expression of TLR-3 and TMEM173/STING, with DHX58-regulated IFIH1 (MDA5)-and DDX58 (RIG-I)-mediated signaling modulating the IFN-I response. Despite these commonalities, a lower viremia and faster return of the draining lymph node transcriptome to baseline levels were observed in the dexamethasone treated rabbits. This highlights a slightly different immune response, as supported by the expression of the different repertoire of IFN-α-21-like and IFN-ω-1-like subtypes on day 1 p.i., which in turn may have driven the different chemokine response between treated and untreated rabbits on day 3 p.i. Future studies in this model should be extended to transcriptomic analysis of lymphoid organs distal to the draining lymph node, e.g., spleen and other lymph nodes, with a comparison of the transcripts in the draining versus distal lymphoid organs. We have also identified few genes upregulated to high magnitude, but with little known association with WNV restriction. These included SLC16A5 and an unannotated gene with putative HERC6 function. We have also identified another unannotated gene with putative IFIT1B function, paving potential new research into the profile of novel ISGs in rabbits. The acute stress model as presented in the current study therefore shows that rabbits are capable of restricting WNV replication and dissemination by a robust innate immune system that can retain its antiviral function despite environmental challenges such as stress.

Virus
The WNV NSW2011 stock used in this study is the same as that previously used in the initial characterization study of the rabbit WNV-challenge model [20]. The stock has the same passage history and same titer as previously described [20]. Briefly, the WNV NSW2011 stock originally derived from a brain isolate of a horse during the 2011 outbreak, propagated in C6/36 cells, then passaged once in Vero cells, and twice further in C6/36 cells, before being used in the current study [3,4,20].

Animals and Experimental Design
All experiments were performed at the University of Queensland, Gatton Campus, in accordance with the ethical standards of The Australian Code for the Care and Use of Animals for Scientific Purposes and were approved by the University of Queensland Production Animal Ethics Committee before commencement of the study (Permit No. SVS/555/17). Rabbits used in this study were 4-5 weeks old NZWRs, sourced from Nanowie (Victoria, Australia) and the University of Queensland Biological Resources (UQBR). The study had two phases: In an initial pilot experiment, four weanling NZWRs were injected with 6.5 mg/kg dexamethasone sodium phosphate (Dexapent, 5 mg/mL; Troy Laboratories, Glendenning, NSW, Australia) intramuscularly in the semimembranosus muscle, split into two doses of 0.3-0.4 mL between the left and right hind legs. This dose of dexamethasone was based on that used by Jeklova et al. [26], but modified to combine the multiple injections to a single bolus injection. The effects of this modified dexamethasone dosing regime on immune parameters were assessed, and this subsequently informed the optimal day post-dexamethasone treatment for WNV inoculation. Two age-matched NZWRs were used as mock-treatment control animals for comparison. Post-dexamethasone treatment, all rabbits were monitored twice daily and bled daily to assess hematological changes and blood cytokine transcript profile. Two treated and one control rabbit were culled on day 3 and 7 post-treatment for tissue collection for histopathology and downstream molecular assays. For the euthanasia, all animals were sedated with xylazine (Xylazil, 20 mg/mL, Troy Laboratories, Glendenning, NSW, Australia), followed by an overdose of pentobarbital (Lethabarb, Troy Laboratories), as previously described [20].
This pilot experiment suggested the optimal day post-dexamethasone injection for WNV inoculation was day 3 pdt (see Section 2.1). For the challenge experiment, nine NZWRs were treated with dexamethasone, as described above, 3 days prior to challenge with 10 5 TCID 50 of WNV NSW2011 in 50 µL media intradermally in the left footpad, as per our previous studies [20]. Seven rabbits were used as mock-treatment controls, receiving a volume of PBS similar to that of the dexamethasone intramuscularly. These rabbits were challenged with the same dose of WNV NSW2011 by the same route as the dexamethasone-treated group. The animals were assessed clinically twice daily and bled once daily. The termination schedule was as indicated in Table 8.
At termination, a complete post-mortem examination was performed on each rabbit with collection of a range of organs and tissues for histopathology and downstream virological and molecular assays.
In order to reduce the amount of environmental variability that may influence the stress level of the rabbits, the ambient temperature and day/night cycle were maintained constant throughout the trial, and husbandry practices and blood sampling were performed at approximately the same time daily. The two animal trials were also performed in the same room in the same animal holding facility under the same conditions, thus justifying the use of the tissues from untreated uninfected rabbits 2001 and 2002 from the pilot trial as baseline controls for the differential gene expression analysis between dexamethasone-treated and untreated WNV-infected rabbits. The use of these controls was also required, as the number of animals used in the second trial reached maximum capacity for the animal room (at n = 16), and no additional clean animals could be included in the trial as baseline controls.

Hematology
Immediately after blood sample collection into EDTA tubes (Multivette ® 600, Sarstedt Ag & Co.), 3 µL of whole blood was used to prepare a routine blood smear on single frosted glass slides. Three blood smears were prepared from each animal at each bleed. After air drying by waving slides in the air for few seconds, smears were fixed immediately in 100% methanol for 5 min. Smears were then air dried, before staining with Wright-Giemsa stain in a Siemens Hematek ® 1000 autostainer (Siemens Corp., Tarrytown, NY, USA). All blood smears were examined using a Nikon Eclipse Ci microscope at 10× or 20× objectives for manual cell counting and 20× to 40× objectives for differential counts.

Histopathology
Tissue samples were harvested immediately after euthanasia and fixed in 10% neutral buffered formalin solution for 48 h before being transferred into 70% ethanol for storage until trimming and routine processing for paraffin embedding. 4 µm thick sections were stained with hematoxylin and eosin and examined on a Nikon Eclipse 51 E microscope. Digital microphotographs were taken using a Nikon DS-Fi1 camera with a DS-U2 unit and NIS elements F 4.60 software. Images are reproduced without manipulations other than cropping and adjustment of light intensity.

Immunohistochemistry
Serial sections (4 µm) were cut from formaldehyde-fixed, paraffin-embedded samples from the popliteal lymph nodes, spleen, thymus, mesenteric lymph node, the virus injection site (left footpad) and brain and subjected to IHC-labeling as previously described in detail [17,21]. Briefly, following deparaffinization, antigen retrieval and several blocking steps, the sections were incubated with primary antibody (details listed in Table 9) followed by visualization of the binding using the Envision kit from DAKO, specific for either mouse or rabbit immunoglobulin as appropriate. The sections were counterstained with Meyer's hematoxylin and examined on a Nikon Eclipse 51 E microscope. A semi-quantitative score was applied that took into account both number of positive cells and labeling intensity for each marker [17]. Digital micrographs were generated as described above.

Cell culture and Virus Isolation
Virus titration was performed by both tissue culture infective dose 50% end-point (TCID 50 ) and plaque assay on serum (viremia) and tissues (10% tissue homogenate) using monolayers of Vero cells (African green monkey kidney cells) as previously described in detail [20].

RNA Extraction and Viral Transcript Quantitation
Half of each draining lymph node harvested during post-mortem was submerged into liquid nitrogen immediately after harvest in order to preserve RNA integrity. RNA extraction was performed by first homogenizing the tissue in 1 mL of Qiazol (Qiagen, Chadstone Centre, Vic, Australia), followed by DNase I digestion in solution and RNA cleanup with RNeasy columns, as previously described [20,21]. Prior to homogenization in Qiazol, the tissue was weighed, in order to allow expression of virus quantitation per gram of input tissue. WNV RNA in samples was quantified by first reverse transcription into cDNA using the qScript first strand synthesis kit (Quanta Bioscience, Gaithersburg, MD, USA) and then qPCR using the Rotorgene SYBR green kit (Qiagen) [20,21]. Using 10-fold dilution series of cDNA derived from RNA extracted from WNV stock with known TCID 50 titer, a standard curve was generated along each qPCR run for absolute quantitation (R 2 of 0.9386 and a limit of detection of 1.3 TCID 50 ). RNA extraction from whole blood was achieved by lysing up to 100 uL of whole blood in 1mL of Qiazol by vortexing. The downstream processing and cytokine transcript qRT-PCR was as described above, using primers as previously described [21].

RNAseq
RNA quality was initially screened by Agilent 2100 Bioanalyzer. Samples with a RIN score of greater than or equal to 4.9 were used for RNA-Seq library preparation, using the Illumina TruSeq Stranded Total RNA LT (Ribo-Zero Gold) Sample Prep Kit (Illumina, RS-122-2301/RS-122-2302), according to the standard manufacturer's protocol (Illumina, 15031048 Rev. E October 2013). Briefly, to enrich for mRNA, 1 µg of total RNA was depleted of rRNA using Ribo-Zero Gold. The enriched mRNA was then subjected to a heat fragmentation step aimed at producing fragments between 130-290 bp (average 185 bp). This step was adjusted to 7 minutes in order to account for the generally low RIN score. cDNA was synthesized from the fragmented RNA using SuperScript II Reverse Transcriptase (Invitrogen, 18064014) and random primers. The resulting cDNA was converted into dsDNA in the presence of dUTP to prevent subsequent amplification of the second strand and thus maintaining the 'strandedness' of the library. Following 3' adenylation and adaptor ligation, libraries were subjected to 15 cycles of PCR to produce libraries ready for sequencing. The libraries were quantified on the Perkin Elmer LabChip GX with the DNA High Sensitivity Reagent kit (Perkin Elmer, CLS760672). Libraries were pooled in equimolar ratios, and the pool was quantified by qPCR using the KAPA Library Quantification Kit-Illumina/Universal (KAPA Biosystems, KK4824) in combination with the Life Technologies Viia 7 real time PCR instrument. Sequencing was performed using the Illumina NextSeq500 (NextSeq control software v2.1.0/Real Time Analysis v2.4.11). The library pool was diluted and denatured according to the standard NextSeq protocol (Document # 15048776 v02), and sequenced to generate paired-end 76 bp reads using a 150 cycle NextSeq500/550 High Output reagent Kit v2 (illumina, FC-404-2002). After sequencing, fastq files were generated using bcl2fastq2 (v2.18.0). Library preparation and sequencing was performed at the Institute for Molecular Bioscience Sequencing Facility (The University of Queensland). Note that RNA samples extracted from the two uninfected control rabbits were spiked with RNA extracted from the WNV NSW2011 inoculum used for the challenge study. This was performed for the possibility to investigate viral genomic dynamics in the draining lymph node in the future.

Differential Gene Expression Analysis
All RNAseq data preparation and differential gene expression analyses were performed using tools within the web-based Galaxy Australia platform [62]. Paired-end reads for each sample, including controls, were mapped to the rabbit genome (OryCun2.0, GenBank assembly ID: GCA_000003625.1) using the gapped alignment tool TopHat (Galaxy version 2.1.1). A matrix of gene read counts for each sample was generated using the tool, "SAM/BAM to count matrix using HTseq code" (Version 0.5), with read quality filter set at MAPQ score of 5. Three R-based differential gene expression analysis packages under the Galaxy tool "Differential _Count models using BioConductor packages" (Galaxy Version 0.28) were used: 1) EdgeR, 2) DEseq2 and 3) Voom-limma. DEGs from each package were determined by an adjusted p-value ≤ 0.05. Only genes that were commonly differentially expressed across all three tools were then used for downstream GO analyses using the web tool, DAVID 6.8. Enriched GO terms and their associated genes were analyzed and graphically visualized using the R packages, GOplot (version 1.0.2) and Vennerable (Version 3.1.0.9000). Heatmaps were generated using heatmap.2 from the R package, gplots (Version 3.0.1.1), with rows ordered based on dendrogram and hierarchical clustering using the hclust function (R) under "complete" mode. Principle component analysis was conducted using the R package, mixOmics (version 6.7.0).

Serology
Quantification of anti-WNV specific antibodies was performed by two methods: blocking enzyme-linked immunosorbent assay (blocking-ELISA) [63] and microneutralization assay [64]. The blocking ELISA protocol has been described extensively in previous publications [2,63,64]. The current assay used lysate from C6/36 cells infected with KUNV strain MRM61C as the coating antigen. Monoclonal antibodies used for these competitive assays were either an anti-flavivirus envelope monoclonal antibody, 6B6C-1, or an anti-WNV NS1 specific monoclonal antibody, 3.1112G [63,64]. A cut-off of 30% inhibition was used to determine positive samples [63].
Supplementary Materials: The following are available online at http://www.mdpi.com/2076-0817/8/4/195/s1, Figure S1: Effect of bolus administration of dexamethasone on hematologic profile in uninfected rabbits. Figure  S2: Effect of dexamethasone treatment on mRNA expression for a select set of cytokines in whole blood. Figure S3: Effect of bolus administration of dexamethasone on thymic morphology in young rabbits. Figure S4: Enriched GO terms from the pairwise comparisons. Figure S5: DEGs associated with GO terms relevant to antiviral immune responses on day 7 p.i. Table S1: Top clusters of antiviral DEGs on day 7 p.i.