Transcriptional response to West Nile virus infection in the zebra finch (Taeniopygia guttata).

West Nile virus (WNV) is a widespread arbovirus that imposes a significant cost to both human and wildlife health. WNV exists in a bird-mosquito transmission cycle in which passerine birds act as the primary reservoir host. As a public health concern, the mammalian immune response to WNV has been studied in detail. Little, however, is known about the avian immune response to WNV. Avian taxa show variable susceptibility to WNV and what drives this variation is unknown. Thus, to study the immune response to WNV in birds, we experimentally infected captive zebra finches (Taeniopygia guttata). Zebra finches provide a useful model, as like many natural avian hosts they are moderately susceptible to WNV and thus provide sufficient viremia to infect mosquitoes. We performed RNAseq in spleen tissue during peak viremia to provide an overview of the transcriptional response. In general, we find strong parallels with the mammalian immune response to WNV, including upregulation of five genes in the Rig-I-like receptor signalling pathway, and offer insights into avian-specific responses. Together with complementary immunological assays, we provide a model of the avian immune response to WNV and set the stage for future comparative studies among variably susceptible populations and species.

Introduction

West Nile virus (WNV) is a single-stranded RNA flavivirus that exists in an avian-mosquito transmission cycle, where birds (typically Passeriformes) act as the primary amplification hosts. In addition to birds, nearly 30 other non-avian vertebrate species have been documented as hosts [1]. Although many WNV-infected hosts are asymptomatic, WNV infection can cause severe meningitis or encephalitis in those that are highly susceptible. Avian species for the most part exhibit low to moderate susceptibility. That is, individuals become infected and develop sufficient viremia for transmission via mosquito blood meal, but the hosts recover and avoid significant mortality (reviewed in [2]). First described in 1937, WNV has not resulted in widespread avian decline throughout its historical range [3], perhaps due to host–parasite coevolution. However, the emergence of WNV in North America in 1999 has negatively impacted a wide range of populations [4,5]. Surveys of North American wild birds have shown a variety of competent WNV hosts, with varying degrees of susceptibility, morbidity and pathogenicity [2]. American robins (Turdus migratorius) appear to be the main host in spreading WNV infection in North America [6], but infection appears most detrimental to members of Family Corvidae [7]. Despite great variation in susceptibility, the mechanisms underlying this variation are primarily unknown [2].

Largely due to interest in human health implications, most work describing the host immune response to WNV infection has been performed in mammalian systems [8]. From these studies, we know that in mammals, both the innate and adaptive arms are critical for virus detection and clearance [9,10]. Within the innate immune response, the retinoic acid-inducible gene 1 (Rig-I)-like receptor (RLR) pathway appears to play a key role in viral clearance. This pathway recognizes viral products and initiates type I interferon expression [11]. Mice lacking the viral recognition RLR genes in this pathway, DDx58 (Rig-I) and IFIH1 (MDA5), become highly susceptible to WNV infection [12]. In the adaptive immune system, a broad range of components appear to play important roles in mounting a response, including antibody and CD4+ and CD8+ T cells [9,13,14]. Interestingly, major histocompatibility complex (MHC) class I genes are upregulated post-infection [15,16]. Viruses typically evade MHC class I detection [17,18], as MHC class I molecules bind and present viral peptides to CD8+ T cells. However, the purpose of WNV-induced MHC expression is unclear.

While the mammalian immune response to WNV infection has been extensively studied, the avian immune response remains mostly unknown. Of the studies in birds, many involve experimentally infecting wild caught birds (reviewed in [2]), or domestic chickens (Gallus gallus) [19]. These studies primarily focus on viral detection, tissue tropism, antibody production or lymphocyte counts [2,19,20]. Little is known about the molecular mechanisms driving the immune response to WNV infection (but see [21]). Furthermore, current avian WNV studies suffer many challenges. Wild caught birds may be co-infected with other parasites (e.g. avian malaria) and are difficult to maintain in captivity for experimental infection studies. Chickens, although an avian model species, are uncommon hosts and highly resistant to WNV infection [22]. Therefore, chickens are not ideal to describe the avian immune response to WNV infection. Passeriformes and Galliformes are also highly divergent bird lineages, with distinctive immune gene repertoires and architecture [23].

As passerine birds are the main hosts for WNV, we have sought to develop a passerine model to study the impacts of WNV infection on a taxonomically appropriate host [24]. We have recently shown that zebra finches, Taeniopygia guttata, are moderately susceptible hosts for WNV [25]. That is, WNV rapidly disseminates to a variety of tissues and is detectable in most samples by 4 days post-inoculation (dpi). Despite rapid development of sufficient viremia for arthropod transmission, zebra finches develop anti-WNV antibodies, clear WNV by 14 dpi, and avoid significant mortality [25]. This moderate disease susceptibility is similar to what is observed in many natural WNV hosts. Zebra finches are also an established biomedical model system with a suite of genetic and genomic tools available [26].

In this study, we experimentally infected zebra finches and performed RNAseq to describe their transcriptional response up to the point of peak viremia. In doing so, we characterize the zebra finch immune response to WNV infection, explore expression of the avian RLR pathway in response to WNV, gain insights into the avian immune response to this widespread infectious disease, and uncover conserved evolutionary responses in avian and mammalian systems.

Results

Experimental infection

We challenged six individuals with 105 plaque-forming units (PFU) WNV and sequenced RNA (Illumina RNAseq) isolated from spleens, an organ critical to the avian immune response. Three birds served as procedural controls and on day 0 were injected subcutaneously with 100 µl of BA1 media, as previously described [27]. Peak viremia occurs at 4.6 ± 1.7 dpi as quantified via RT-PCR [25] and thus, we characterized the transcriptional response leading to (2 dpi, n = 3) and at peak viral load (4 dpi, n = 3) in the present study. WNV RNA was detected by culture in lung and kidney RNA pools of two out of three birds sampled at day 2, and all three birds sampled at 4 dpi. These findings were verified by semi-quantitative RT-PCR. WNV is rarely detected in spleen by 2 dpi, but all birds previously inoculated at 105 PFU developed WNV antibodies [25]. We therefore treated all six birds inoculated with WNV as being infected.

Sequencing results and read mapping

We obtained 18–30 million paired-end, 100 bp reads for each sample and removed 0.57–1.24% of the total bases after adapter trimming (electronic supplementary material, table S1). On average, 79.0–80.8% total trimmed reads mapped to the zebra finch reference genome (electronic supplementary material, table S2), corresponding to 18 618 Ensembl-annotated genes [28]. Of these, 14 114 genes averaged at least five mapped reads across all samples and were used for differential expression (DE) analyses.

Sample clustering and differential expression

We tested for DE in two ways: as pairwise comparisons between treatments with DEseq2 [29] and as a time-course grouping genes into expression paths with EBSeqHMM [30]. To visualize patterns of expression variation among samples, we conducted principal component analysis (PCA) and distance-based clustering (electronic supplementary material, figures S1 and S2). The first three principal components explained 93.04% of the variance in gene expression, but none of the PCs were significantly correlated with treatment (ANOVA, PC1: p = 0.288, PC2: p = 0.956, PC3: p = 0.202).

Although clustering analyses suggest that across the genome, much of the variation in expression was independent of the experimental treatment, pairwise comparisons revealed many genes that were regulated in response to infection (electronic supplementary material, table S3). When comparing Control versus 2 dpi, we found 161 differentially expressed genes (false discovery rate (FDR) < 0.10, average log2 fold-change (FC) = 1.74). This gene list includes several immune-related genes associated with the innate (e.g. IL18) and adaptive (e.g. MHC IIB) immune system (table 1, figure 1). Sixty-five genes were differentially expressed between Control and 4 dpi (average log2FC = 1.61), also with several immune-relevant genes including five genes in the RLR pathway (table 1, figures 2 and 3). Lastly, we observed 44 DE genes between 2 dpi versus 4 dpi individuals (average log2FC = 1.56). Three of these have described functions in immunity. We also combined 2 dpi and 4 dpi cohorts and compared with Control, but due to high variation in gene expression between days 2 and 4 dpi, we only found 16 DE genes (average log2FC = 1.64) between Control and Infected cohorts, one of which was associated with immunity.

Table 1.

Candidate immune genes differentially expressed in the present study and comparisons with mammals.

Ensembl ID	gene name	log2 fold change	FDR	regulation pattern observed	regulation pattern in mammals	references
Control versus Infected
ENSTGUG00000013615	NFKBIZ	0.73	0.064	up	up	[31]
Control versus 2 dpi
ENSTGUG00000000297	IL18	1.01	0.010	up	no change	[32]
ENSTGUG00000000678	TIM1	1.49	7.99 × 10−05	up	up	[33,34]
ENSTGUG00000001485	IRF6	−2.09	0.037	down	up	[35]
ENSTGUG00000003354	NKRF	−2.35	4.85 × 10−05	down	unknown
ENSTGUG00000005295	C-C motif chemokine	2.08	0.007	up	up	[31,35,36]
ENSTGUG00000008638	UCHL1	−1.91	0.029	down	unknown
ENSTGUG00000008991	APOD	1.76	0.053	up	up	[37]
ENSTGUG00000009454	IFITM10	1.24	2.71 × 10−04	up	unknown
ENSTGUG00000009769	TNFRSF13C	0.89	0.010	up	unknown
ENSTGUG00000015634	Novel gene (MHC IIB)	2.23	0.001	up	up	[38]
ENSTGUG00000016383	SIGLEC1	1.39	0.046	up	up	[31]
ENSTGUG00000017149	Novel gene (MHC IIB)	1.57	0.099	up	up	[38]
Control versus 4 dpi
ENSTGUG00000001516	DDx58	1.50	1.45 × 10−08	up	up	[31]
ENSTGUG00000002144	IRF4	1.39	0.022	up	up	[35]
ENSTGUG00000002305	LY86	−1.07	1.70 × 10−05	down	unknown
ENSTGUG00000002516	DHx58	1.67	4.05 × 10−06	up	up	[31]
ENSTGUG00000004105	ADAR	1.13	1.70 × 10−05	up	up	[39]
ENSTGUG00000006914	IFIH1	0.95	0.093	up	up	[31]
ENSTGUG00000007454	TNFRSF13B	1.61	0.010	up	unknown
ENSTGUG00000008354	IFIT5	2.97	1.07 × 10−09	up	unknown
ENSTGUG00000008788	EIF2AK2	2.15	9.86 × 10−07	up	up	[35,40]
ENSTGUG00000009162	PXK	1.19	0.011	up	unknown
ENSTGUG00000009536	TRIM25	1.24	0.001	up	up	[31]
ENSTGUG00000009838	IRF7	0.80	0.047	up	up	[31]
ENSTGUG00000011784	ZC3HAV1	1.31	7.09 × 10−08	up	up	[31]
ENSTGUG00000017534	MOV10	1.29	0.017	up	up	[39]
2 dpi versus 4 dpi
ENSTGUG00000003354	NKRF	1.70	0.033	up	unknown
ENSTGUG00000005206	ADA	−1.48	0.001	down	unknown
ENSTGUG00000011784	ZC3HAV1	0.82	0.058	up	up	[31]

Figure 1.

Immune genes differentially expressed between day 2 post-inoculation and Control. (a) Heatmap of expression levels (log transformed read counts) across all treatments of immune genes differentially expressed at 2 dpi relative to Control. (b–d) Expression values (normalized read counts) for three key immune genes and their regulation pattern classification by EBSeqHMM. Asterisks represent statistical significance in DEseq2 analysis after FDR correction (*p < 0.10, **p < 0.05, ***p < 0.01).

Figure 2.

Immune genes differentially expressed between day 4 post-inoculation and Control. (a) Heatmap of expression levels (log transformed read counts) across all treatments of immune genes differentially expressed at 4 dpi relative to Control. (b–d) Expression values (normalized read counts) for three key immune genes and their regulation pattern classification by EBSeqHMM. Asterisks represent statistical significance in DEseq2 analysis after FDR correction (*p < 0.10, **p < 0.05, ***p < 0.01).

Figure 3.

Regulation of the zebra finch RLR pathway. Colour represents log2 fold change between Control and 4 dpi. Asterisks represent statistical significance in DEseq2 analysis after FDR correction (*p < 0.10, **p < 0.05, ***p < 0.01).

When analysed for DE as a time course in EBSeqHMM, 686 genes showed evidence of differential expression (posterior probability >0.99, FDR < 0.01). Most DE genes (n = 561) were suppressed relative to Controls on days 2 and 4 post-infection (‘Down-Down’). Seventy-five genes were ‘Up-Down’, 49 were ‘Down-Up’ and one was ‘Up-Up’. As expected, we found overlap of several immune genes between the two analyses. For example, IL18, APOD and IFITM10 are ‘Up-Down’ and this trend is reflected in the DEseq2 Control versus 2 dpi analysis (figure 1).

Functional annotation of differentially expressed genes

To place differentially expressed genes into groups based on their biological function, we performed a gene ontology (GO) analysis using the GOrilla tool [41,42]. GOrilla uses the rank ordered of genes from DEseq2 based on FDR-adjusted p-values. An enrichment score is calculated based on the number of genes in the top of the list that belong to a particular GO category relative to the expected number based on the frequency of functionally related genes in the total list. As above, we conducted GO analyses based on multiple pairwise analyses of gene expression. We found five significantly enriched GO categories between Control versus Infected (2 and 4 dpi) cohorts, of which ‘response to virus’ is the most significant (FDR = 0.008, Enrichment = 5.34) (table 2). We observed the strongest evidence of functional enrichment in the Control versus 2 dpi (n = 120) and Control versus 4 dpi (n = 36) contrasts (FDR < 0.05) (electronic supplementary material, table S4). Many enriched GO terms in the Control versus 2 dpi contrast are involved in membrane components, metabolism and cellular processes. Four GO categories were immune relevant, including ‘inflammatory response’ and ‘positive regulation of cytokine biosynthetic process’ (electronic supplementary material, table S4). The immune response manifests itself most strongly in the Control versus 4 dpi contrast with many enriched GO terms being immune-related categories (table 2; electronic supplementary material, table S4) and a broad range of differentially expressed immune genes (n = 14, table 1). Only two GO categories are enriched between 2 and 4 dpi: ‘inner mitochondrial membrane protein complex’ and ‘mitochondrial protein complex’ (table 2; electronic supplementary material, table S4).

Table 2.

Top five statistically enriched gene ontology (GO) categories, FDR-adjusted p-value, and GOrilla enrichment score, among DEseq2 pairwise comparisons. Enrichment is calculated as (b/n)/(B/N), where N is the total number of genes, B is the total number of genes associated with a specific GO term, n is the number of genes in the top of the input list, and b is the number of GO-term-associated genes in the top of the list [41,42].

GO ID	description	FDR	enrichment
Control versus WNV
GO:0009615	response to virus	0.008	5.34
GO:0051276	chromosome organization	0.019	1.97
GO:1903047	mitotic cell cycle process	0.035	1.69
GO:0034723	DNA replication-dependent nucleosome assembly	0.039	14.27
GO:0006335	DNA replication-dependent nucleosome organization	0.048	14.27
Control versus 2 dpi
GO:0044425	membrane part	1.77 × 10−07	1.29
GO:0098800	inner mitochondrial membrane protein complex	1.15 × 10−06	3.73
GO:0044459	plasma membrane part	2.66 × 10−06	1.98
GO:0031224	intrinsic component of membrane	4.04 × 10−06	1.35
GO:0098798	mitochondrial protein complex	4.06 × 10−06	3.32
Control versus 4 dpi
GO:0009615	response to virus	0.001	19.23
GO:0051607	defence response to virus	0.002	51.27
GO:0060337	type I interferon signalling pathway	0.002	24.25
GO:0051707	response to other organism	0.002	11.73
GO:0098586	cellular response to virus	0.003	70.99
2 dpi versus 4 dpi
GO:0098800	inner mitochondrial membrane protein complex	0.01	3.09
GO:0098798	mitochondrial protein complex	0.04	2.69

We also conducted a similar analysis of genes identified as DE by EBseqHMM, which revealed one (Up-Up), 199 (Up-Down), 69 (Down-Up) and 527 (Down-Down) significantly enriched GO categories (FDR < 0.05) (electronic supplementary material, table S5). Interestingly, Up-Down GO categories had the strongest representation of immune-related GO terms, including ‘immune response’ (FDR = 4.85 × 10−4) and ‘negative regulation of immune system process’ (FDR = 6.01 × 10−4). Among Down-Down genes, we observed enrichment of many metabolic and membrane processes and only one immune-related category (‘positive regulation of innate immune response’, FDR = 0.01). We found enrichment of mitochondrial components and processes among Down-Up genes, similar to the 2 dpi versus 4 dpi contrast in the DEseq2 analysis. Additionally, 10 categories involved in immunoglobulin processes were significantly enriched among ‘Down-Up’ genes, driven by the presence of the joining chain of multimeric IgA and IgM (JCHAIN) gene. Lastly, as in the DEseq2-based analysis, we also detected a strong enrichment signature of membrane proteins. Genes annotated as ‘plasma membrane part’ were highly enriched among those showing an Up-Down pattern (FDR = 1.61 × 10−12, electronic supplementary material, table S5). Combined, we found broad overlap in GO representation between the EBseqHMM and DEseq2 approaches.

In addition to placing genes into broad systematic functions in the GO analysis, we were also interested in placing our gene expression results in the context of immune pathways of interest. The RLR antiviral pathway is critical to WNV clearance in mammals [12] and appears important in mounting an immune response to avian influenza in ducks [43-45]. Using Pathview v1.8.0 [46], we find that WNV infection induces the RLR pathway. Five genes, including the two RLR genes, DDx58 and IFIH1, which encode the Rig-I and MDA5 viral detection molecules, are significantly upregulated (table 1, figure 2; electronic supplementary material, figure S3). We detect expression of 36/37 genes in the pathway, many of which are also upregulated, though not always significantly (figure 3).

Discussion

We have characterized the zebra finch transcriptional response to WNV infection. Overall, we find that as in mammalian systems, components of both the adaptive and innate immune pathways are activated following infection. While WNV is primarily an avian-specific infectious disease, most work describing the host immune response to infection has been performed in mammals. Despite genomic, physiological and evolutionary differences between birds and mammals, the host immune response shows broad similarity between taxa (table 1).

We were particularly interested in the role of the innate RLR pathway. This pathway mounts an antiviral innate immune response and is critical for WNV detection and clearance in mammals [12]. We have shown here that the RLR pathway in zebra finches is induced by WNV infection. Five genes in this pathway are significantly upregulated at 4 dpi (figure 3; electronic supplementary material, figure S3), including DDx58 and IFIH1 (figure 2b,d), which encode molecules that recognize WNV particles in mammals [47]. This results in a corresponding over-representation of genes in the interferon signalling and regulation GO categories (table 2; electronic supplementary material, table S4). While no studies have investigated the role of the RLR following WNV infection in birds, this pathway appears important for avian influenza clearance in ducks [43-45], Buggy Creek virus clearance in house sparrows [48], and probably for the broad avian antiviral immune response, including WNV. Interestingly, chickens (Gallus gallus), which are often used as sentinels for WNV, have lost the gene encoding the DDx58 RLR during their evolution [43] yet do not develop disease post WNV infection [22]. This suggests that chickens respond to WNV using a Rig-I independent mechanism and highlights the importance of future work targeting the evolution of avian innate immunity.

We observed other parallels with mammals as well (table 1). For example, T-cell immunoglobulin mucin receptor 1 (TIM1) is upregulated at 2 dpi in zebra finches (figure 1a,c). In human cell lines, expression of TIM1 promotes infection of WNV virus-like particles (VLPs) [33,34], suggesting that the upregulation of TIM1 seen in zebra finches may promote viral entry as well. Similarly, C-C motif chemokine (ENSTGUG00000005295) is upregulated in our study at 2 dpi and in previous human cell line and mouse experiments, suggesting a conserved role in chemokine production following WNV infection [31,35,36]. Apolipoprotein D (APOD), a gene typically involved in brain injury and potentially responding to the neurodegenerative nature of WNV, is upregulated in WNV-infected mice [37], as well as in our study. Two interferon-stimulated genes (ISGs), ADAR and MOV10 are both significantly upregulated at 4 dpi relative to Control. Schoggins et al. [39] showed ADAR expression to enhance WNV replication and MOV10 expression to have antiviral activity. While further testing of these genes is needed to validate their roles in avian WNV infection, they nonetheless offer insights into a broad range of conserved responses between mammals and birds.

Within the adaptive immune response, the role of the MHC in the host response to WNV is also particularly interesting. The MHC plays a key role in antigen processing and presentation. The MHC comprises two main gene families (Class I & II) and both are upregulated in mammals following WNV infection [15,16,38]. Similarly, two genes encoding MHC class IIB proteins are significantly upregulated in zebra finches at 2 dpi (figure 1). Unlike mammals, however, we found that MHC class I is not significantly DE in any comparison (e.g. C versus 4 dpi, log2FC = 0.001, FDR = 0.99). In mammals, upregulation of MHC class I may not be adaptive for the host, as upregulation may be a mechanism by which the virus evades Natural Killer (NK) cell detection by the innate immune system [15]. It has also been suggested that MHC upregulation is a by-product of flavivirus assembly [15]. Interestingly, at 2 dpi, interleukin-18 (IL18) is significantly upregulated (table 1, figure 1a,b). IL18 can enhance NK cell activity [49] and is potentially a mechanism by which the immune system can counteract WNV evasion strategies via NK cell activation, although further testing is needed to quantify NK cell activity in zebra finches to support this hypothesis.

Despite many similarities, several immune genes differentially expressed in our analyses have not been previously reported in the mammalian WNV literature or are expressed differently in zebra finches (table 1). For example, at 2 dpi, the proinflammatory cytokine IL18 was significantly upregulated in zebra finches (figure 1), contrasting a previous study in human cell lines, which showed no difference in IL18 expression following WNV infection [32]. Furthermore, interferon regulatory factor 6 (IRF6) was downregulated at 2 dpi, but upregulated in human macrophages following infection [35]. Another significantly downregulated gene at 2 dpi, ubiquitin carboxyl-terminal hydrolase L1 (UCHL1), has been previously associated with pattern recognition receptor (PRR) pathway (e.g. RLR) function in human cell lines infected with high-risk human papilloma virus [50]. When upregulated, UCHL1 supresses PRR expression leading to viral evasion of the host immune response. However, downregulation of UCHL1 restores functional PRR pathways [50]. Thus, the downregulation of UCHL1 2 dpi in zebra finches may be associated with the upregulation of the PRR RLR pathway in this study (table 1, figure 3). Interferon-induced protein with tetratricopeptide repeats (IFIT) and interferon-inducible transmembrane protein (IFITM) gene families are known innate antiviral proteins and have been shown to restrict WNV entry in human cells lines [39,51]. Both IFIT5 and IFITM10 are upregulated (figure 2a,c) in our study and yet, to our knowledge, neither have previously been implicated in the WNV immune response. This potentially reveals an avian-specific function of IFIT5 and IFITM10. Lastly, several genes involved in metabolic and mitochondrial processes were DE in our analyses. Viral alteration of host metabolism typically benefits viral replication [52,53] and highlights the need for future work investigating the role of WNV on host physiology.

Functional enrichment of immune-related GO terms primarily appears in Up-Down path defined by EBseqHMM (electronic supplementary material, table S5), as many genes in the immune system are upregulated post-infection (table 1, figures 1 and 2). In both the EBseqHMM and DEseq2 analyses, most of the significant immune GO categories are innate immune responses, although adaptive immune categories involved in immunoglobulin complexes and B & T cell proliferation appear in the EBseqHMM analysis (table 2; electronic supplementary material, tables S4 and S5). Similar to the mammalian model, broad organismal processes, encompassing both innate and adaptive immunity, are represented in the zebra finch response to WNV.

Like many passerine birds infected in nature, zebra finches are moderately susceptible to WNV, developing sufficient viremia to serve as competent hosts, but generally resist mortality due to infection [25]. While there are clear differences among treatments in terms of differentially expressed genes (table 1), the modest effect of treatment on overall expression profile (electronic supplementary material, figures S1 and S2) may be a reflection of this moderate susceptibility. Most zebra finches are able to clear WNV infection by 14 dpi [25]. WNV infection intensity varies among tissues [20], but due to the spleen's important role in the avian immune system [54,55] we expect the results presented here to be representative of the overall immune response. Although we expect to have missed some genes that are regulated in response to infection, DEseq2 has been shown to perform very well (low false positive rate) in experiments with a sample size of three [56]. Further studies will also be required to document more subtle, and tissue-specific patterns of gene regulation in response to infection. We note that we only sampled our control group at 4 dpi and thus, do not have a direct procedural control at 2 dpi. Changes in gene expression at 2 dpi therefore could be in part due to the injection itself. Pronounced DE of immune-related genes at 2 dpi, however, suggests that changes in gene expression were driven by WNV infection rather than by the injection, which might be predicted to trigger a more general stress response.

We have begun to develop the zebra finch as an avian model for the host response to WNV infection. We show here that in terms of gene expression, the zebra finch immune response is largely conserved with that seen in mammalian-based studies (table 1). Additionally, we identify many components of the immune system that have not been previously implicated in the host immune response to WNV. This potentially reveals an avian-specific immune response and highlights avenues for future research. Combined with our recent immunological characterization [25], we have broadly described the immune response of a moderately susceptible avian host for WNV. This sets the stage for future comparative work to uncover the genetic basis of variable avian susceptibility to WNV infection.

Methods

Experimental set-up

All animal use was approved by the USGS National Wildlife Health Center Institutional Animal Care and Use Committee (IACUC Protocol: EP120521) and this study was performed in accordance with USGS IACUC guidelines. The experimental infection set-up is described in detail in [25]. Briefly, nine female zebra finches were randomly divided into three cohorts, one unchallenged and two challenged (n = 3 each). Birds were challenged subcutaneously with 100 µl BA1 media containing 105 plaque-forming units (PFU) of the 1999 American crow isolate of WNV (NWHC 16399-3) and sacrificed at 2 and 4 dpi, corresponding to peak viremia. Uninfected individuals were injected with 100 µl BA1 media and sacrificed at 4 dpi. WNV infection was confirmed by RT-PCR, as previously described [26], in lung and kidney pooled tissue [25]. Due to the critical role of the spleen in the initiation of the immune response, and its common use in experimental infection gene expression studies [57-60], we focused our study on gene expression in the spleen. Spleens from each individual were removed, placed into RNAlater (Qiagen, Valencia, CA, USA), and frozen at −80°C until RNA extraction.

RNA extraction and sequencing

Whole spleen tissue was homogenized in Tri-Reagent (Molecular Research Company) and total RNA was purified with a Qiagen RNeasy (Valencia, CA, USA) mini kit following the manufacturer's protocol. RNA was DNAse treated and purified. Purified RNA was quality assessed on a Bioanalyzer (Agilent, Wilmington, DE, USA) to ensure RNA quality before sequencing (RIN = 6.6–8.1). All library prep and sequencing was performed at the University of Illinois Roy J. Carver Biotechnology Center. A library for each sample was prepared with an Illumina TruSeq Stranded RNA sample prep kit. All libraries were pooled, quantitated by qPCR, and sequenced on one lane of an Illumina HiSeq 2000 with a TruSeq SBS Sequencing Kit producing paired-end 100 nt reads. Reads were analysed with Casava 1.8.2 following manufacturer's instructions (Illumina, San Diego, CA). Sequencing data from this study have been deposited in the NCBI Sequence Read Archive (BioProject: PRJNA352507).

Adapter trimming and read mapping

We removed Illumina adapters from reads with Trim Galore! v0.3.7 (http://www.bioinformatics.babraham.ac.uk/projects/trim_galore/) which makes use of Cutadapt v1.7.1 [61]. Reads were then mapped to the zebra finch genome (v3.2.74,26) using TopHat v2.0.13 [62], which uses the aligner Bowtie v2.2.4 [63]. We specified the library type as fr-firststrand in TopHat2. Successfully mapped reads were converted from SAM to BAM format with SAMtools View v1.2 [64,65] and counted in htseq-count v0.6.0 specifying ‘-s rev’ [66]. This assigned zebra finch Ensembl gene IDs and we only retained genes with >5X mapping across each sample.

Differential expression

Gene counts were then normalized for read-depth and analysed for DE in DEseq2 v1.8.1 [29]. We analysed DE across four comparisons: Control versus Infected, Control versus 2 dpi, Control versus 4 dpi, and 2 dpi versus 4 dpi. We visualized expression profiles in R v.3.3.0 [67] by PCA with the R package pcaExplorer [68], and hierarchical clustering heat maps with the ggplot2 library [69] following the DEseq2 manual. DEseq2 tests for DE with a Wald test and genes were considered differentially expressed if the Benjamini & Hochberg [70] FDR correction for multiple testing p-value less than 0.10. We chose this significance threshold as DEseq2 is generally conservative in classifying DE [71]. Furthermore, this cut-off is used by the DEseq2 authors [29] and has been used in other RNAseq experimental infection studies [72]. We plotted genes of interest individually with the plotCounts function in DEseq2 and clustered expression profiles of these genes with the pheatmap R library to view expression levels across samples and treatments.

We tested DE genes for enriched gene ontology (GO) categories with GOrilla [41,42]. GOrilla does not perform analyses with zebra finch Ensembl IDs, so we converted zebra finch Ensembl IDs to human Ensembl IDs using BioMart [73]. We used this set of 10 152 genes for analysis. For each pairwise comparison, we used the FDR ranked order DE genes from DEseq2. Statistical significance was determined with p-values corrected for multiple hypothesis testing (p < 0.05) using the Benjamini & Hochberg method [70]. To visualize DE results in the context of the RLR pathway, we used Pathview v1.8.0 [46] to plot the log fold change of each gene detected in our dataset into the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway (KEGG ID = 04622) [74,75].

Time-course gene expression

In addition to the pair-wise comparisons performed in DEseq2, we were interested in understanding how clusters of genes are differentially expressed over the time course of infection. Thus, we performed DE analyses in EBSeqHMM [30]. EBSeqHMM uses a Bayesian approach with a hidden Markov model to identify DE between ordered conditions. Genes are then grouped into expression paths (i.e. ‘Up-Down’, ‘Down-Down’), in which DE occurs when expression paths change between at least one adjacent condition. For example, a gene upregulated at both 2 dpi relative to Control and 4 dpi relative to 2 dpi would be classified as ‘Up-Up’. We included three time points, with Control individuals classified as t1, 2 dpi as t2 and 4 dpi as t3. Genes were considered DE at posterior probability >0.99 and FDR <0.01. We chose a more stringent cut-off in this analysis as EBseq can be liberal in classifying differential expression [71] and based on visual inspection of expression profiles. We ordered genes based on posterior-probability for each expression path and performed the GO analysis described above.