Transcriptional response to West Nile virus infection in the zebra finch (Taeniopygia guttata), a songbird model for immune function

West Nile Virus (WNV) is the one of most widespread arboviruses worldwide. WNV exists in a bird-mosquito transmission cycle where 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 splenic RNAseq 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 up-regulation of five genes in the Rig-I-like receptor signaling 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.


West Nile Virus (WNV) is the one of most widespread arboviruses 17
worldwide. WNV exists in a bird-mosquito transmission cycle where passerine birds 18 act as the primary reservoir host. As a public health concern, the mammalian 19 immune response to WNV has been studied in detail. Little, however, is known 20 about the avian immune response to WNV. Avian taxa show variable susceptibility 21 to WNV and what drives this variation is unknown. Thus, to study the immune 22 response to WNV in birds, we experimentally infected captive zebra finches 23 (Taeniopygia guttata). Zebra finches provide a useful model, as like many natural 24 avian hosts they are moderately susceptible to WNV and thus provide sufficient 25 viremia to infect mosquitoes. We performed splenic RNAseq during peak viremia to 26 provide an overview of the transcriptional response. In general, we find strong 27 parallels with the mammalian immune response to WNV, including up-regulation of 28 five genes in the Rig-I-like receptor signaling pathway, and offer insights into avian 29 specific responses. Together with complementary immunological assays, we 30 provide a model of the avian immune response to WNV and set the stage for future 31 comparative studies among variably susceptible populations and species. 32 33

Introduction: 34
West Nile virus (WNV) is a single-stranded RNA flavivirus that exists in an 35 avian-mosquito transmission cycle, where birds (typically Passeriformes) act as the 36 primary amplification hosts. In addition to birds, nearly 30 other non-avian vertebrate 37 species have been documented as hosts (1). Although many WNV-infected hosts are 38 asymptomatic, WNV infection can cause severe meningitis or encephalitis in those that 39 are highly susceptible. Avian species for the most part exhibit low to moderate 40 susceptibility. That is, individuals become infected and develop sufficient viremia for 41 transmission via mosquito blood meal, but the hosts recover and avoid significant 42 mortality, reviewed in (2). First described in 1937, WNV has not resulted in widespread 43 avian decline throughout its historical range (3), perhaps due to host-parasite coevolution. 44 However, the emergence of WNV in North America in 1999 has negatively impacted a 45 wide range of populations (4,5). Surveys of North American wild birds have shown a 46 variety of competent WNV hosts, with varying degrees of susceptibility, morbidity, and 47 pathogenicity (2). American robins (Turdus migratorius) appear to be the main host in 48 spreading WNV infection in North America (6), but infection appears most detrimental to 49 members of Family Corvidae (7). Despite great variation in susceptibility, the 50 mechanisms underlying this variation are primarily unknown (2). 51 52 Largely due to interest in human health implications, most work describing the 53 host immune response to WNV infection has been performed in mammalian systems (8). 54 From these studies, we know that in mammals, both the innate and adaptive arms are 55 critical for virus detection and clearance (9,10). Within the innate immune response, the 56 retinoic acid-inducible gene 1(Rig-I)-like receptor (RLR) pathway appears to play a key 57 role in viral clearance. This pathway recognizes viral products and initiates type I 58 interferon expression (11). Mice lacking the viral recognition RLR genes in this pathway, 59 DDx58 (Rig-I) and IFIH1 (MDA5), become highly susceptible to WNV infection (12). In 60 the adaptive immune system, a broad range of components appear to play important roles 61 in mounting a response, including antibody and CD4+ and CD8+ T cells (9,13,14). 62 Interestingly, major histocompatibility complex (MHC) class I genes are up-regulated 63 post-infection (15,16 tropism, antibody production, or lymphocyte counts (2,19,20). Little is known about the 72 molecular mechanisms driving the immune response to WNV infection, but see (21)

Experimental infection 101
We challenged six individuals with with 10 5 plaque forming units (PFU) WNV 102 and sequenced RNA (Illumina RNAseq) isolated from spleens, an organ critical to the 103 avian immune response. Three birds served as procedural controls and on day 0 were 104 injected subcutaneously with 100 µL of BA1 media, as previously described (27). 105 Peak viremia occurrs at 4.6 ±1.7 dpi as quantified via RT-PCR (25) and thus, we 106 characterized the transcriptional response leading to (2dpi, n=3) and at peak viral load 107 (4dpi, n=3) in the present study. WNV RNA was detected by culture in lung and kidney 108 RNA pools of 2 out of 3 birds sampled at day 2, and all 3 birds sampled at 4dpi. These 109 findings were verified by semi-quantitative RT-PCR. Because WNV is rarely detected in 110 spleen by 2dpi, but all birds previously inoculated at 10 5 PFU developed WNV antibodies 111 When comparing Control vs. 2dpi, we found 161 differentially expressed genes 134 (adjusted p < 0.10, average log 2 fold-change (FC) = 1.74). This gene list includes several 135 immune related genes associated with the innate (e.g. IL18) and adaptive (e.g. MHC IIB) 136 immune system (Table 1, Figure 1). Sixty-five genes were differentially expressed 137 between Control and 4dpi (average log 2 FC = 1.61), also with several immune relevant 138 genes including five genes in the RLR pathway (Table 1, Figure 2, Figure 3). Lastly, we 139 observed 44 DE genes between 2dpi vs. 4dpi individuals (average log 2 FC = 1.56). Three 140 of these have described functions in immunity. The complete list of DEseq2 DE genes 141 (adjusted p < 0.10) across all comparisons can be found in Supplementary Table S3 and  142 immune relevant DE genes are listed in Table 1 along with their gene name, log 2 fold 143 change, padj value, and comparison to mammalian studies. We also combined 2dpi and 144 4dpi cohorts and compared with control, but due to high variation in gene expression 145 between days 2 and 4 dpi, we only found 16 DE genes (average log 2 FC = 1.64) between 146 Control and Infected cohorts, one of which was associated with immunity. 147 148 When analyzed for DE as a time course in EBSeqHMM, 686 genes showed 149 evidence of differential expression (posterior probability > 0.99, FDR < 0.01). Most DE 150 genes (n = 561) were suppressed following infection ("Down-Down") or varied between 151 conditions. Seventy-five genes were "Up-Down", 49 were "Down-Up" and 1 was "Up-152 Up". As expected, we found overlap of several immune genes between the two analyses. 153 For example, IL18, APOD and IFITM10 are "Up-Down" (Supplementary Table S5) Table S4) and a 165 broad range of immune functioning genes differentially expressed (n=14 , Tables 1 & 2). 166 When comparing genes differentially expressed between Control and 2dpi, only one GO 167 category was strongly enriched ("integral to membrane", expected = 15, observed = 35, p 168 = 0.00065). In this contrast we also observed several slightly enriched immune related 169 GO categories such as "inflammatory response" and "negative regulation of T cell 170 migration" but these fall outside of our significance threshold (Supplementary Table S4) 171 and reflect only few genes (e.g. APOD and IL18) that were DE in the Control vs. 2dpi 172 analysis. A set of 57 enriched GO categories describe genes DE between 2 and 4dpi. 173 Although some of these involved T and B cell regulation, these categories were also only 174 represented by a single gene each (ADA) ( Interestingly, Up-Down GO categories had the strongest representation of immune 182 related GO terms, including "cytokine receptor activity" (expected = 0, observed = 3 p = 183 0.015) and "positive regulation of interferon-gamma production" (expected = 0, observed 184 = 2, p = 0.028). However, we also observed a strong immune signature among Down-185 Down Genes (GO Term "Immune Response", expected = 3, observed = 11, p = 0.013), 186 which represents several innate immune genes including three complement genes and 187 interferon gamma. Twelve categories were significantly enriched among "Down-Up" 188 genes, including "double-stranded RNA binding" (expected = 0, observed = 2, adjusted p 189 = 0.049). One of the observed genes in this category is OASL, which has been 190 demonstrated to have antiviral activity towards WNV in chickens (21). Lastly, as in the 191 DEseq2-based analysis, we also detected a strong enrichment signature of membrane 192 proteins. Genes annotated as "integral to membrane" were highly enriched among those 193 showing an Up-Down pattern (expected = 6, observed =14, adjusted p < 0.028, 194 Supplementary Table S5). Combined, we find broad overlap in GO representation 195 between the EBseqHMM and DEseq2 approaches. To reveal any genes responding to infection missed by our genome based 209 analysis, we created a reference transcriptome with Trinity v2.1.1 (35) and annotated 210 against the NCBI non-redundant database with DIAMOND, a BLASTx-like aligner (37). 211 Our de novo transcriptome assembly resulted in 393,408 Trinity transcripts. We 212 visualized our annotated transcriptome with MEGAN (38), which places genes into 213 KEGG pathways. In doing so, we confirmed many expression results from the genome-214 based approach and importantly, this analysis revealed expression of two transcripts with 215 homology to interferon alpha, which was absent (but annotated) in the genome based 216 analysis described above. Combined, our genome-based and reference-free approaches 217 detect regulation of the complete RLR pathway in zebra finches. 218 219

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

228
We were particularly interested in the role of the innate RLR pathway. This 229 pathway mounts an antiviral innate immune response and is critical for WNV detection 230 and clearance in mammals (12)  We observed other parallels with mammals as well (Table 1) Figure 1). Unlike mammals, however, we found that MHC class I is not significantly DE 263 in any comparison. In mammals, upregulation of MHCI may not be adaptive for the host, 264 as upregulation may actually be a mechanism by which the virus evades Natural Killer 265 (NK) cell detection by the innate immune system (15). It has also been suggested that 266 interleukin-18 (IL18) is significantly up-regulated (Table 1, Figure1A,B). IL18 can 268 enhance NK cell activity (50) and is potentially a mechanism by which the immune 269 system can counteract WNV evasion strategies via NK cell activation, although further 270 testing is needed to quantify NK cell activity in zebra finches to support this hypothesis.  (Table 1). For example, at 2dpi, the 275 proinflammatory cytokine IL18 was significantly up-regulated in zebra finches ( Figure  276 1), contrasting a previous study in human cell lines, which show no difference in IL18 277 expression following WNV infection (50). Furthermore, interferon regulatory factor 6 278 (IRF6) was down-regulated at 2dpi, but up-regulated in human macrophages following 279 infection (43). Another significantly down-regulated gene at 2dpi, ubiquitin carboxyl-280 terminal hydrolase L1 (UCHL1), has been shown to suppress cellular innate immunity in 281 human cell lines infected with high-risk human papilloma virus (52). Down-regulation 282 restores functional pattern recognition receptor (PRR) pathways (e.g. RLR). The down-283 regulation of UCHL1 here thus is associated with the previously described regulation of 284 the PRR RLR pathway in this study (Table 1)  The zebra finch was the second bird to have its genome sequenced (26). However, 315 the zebra finch genome assembly and gene annotation remain incomplete. Immune genes 316 in particular are difficult to reconstruct in genomes due to their complex evolutionary 317 history (e.g. MHC, [23]). Our mapping results indicated that roughly 80% of our reads 318 mapped to the genome. While the unmapped 20% could contain poor quality reads or 319 non-avian sequences, it also represents zebra finch reads that could not be appropriately 320 mapped to the reference. Thus, we took a genome-independent approach and assembled a 321 de novo zebra finch transcriptome from our nine paired-end samples. Interestingly, RLR 322 induced interferon expression was not detected in the genome-based analysis. The 323 transcriptome analysis, however, revealed expression of two transcripts with homology to 324 interferon alpha, thus completing the pathway from virus detection to interferon 325 production. 326

327
We have begun to develop the zebra finch as an avian model for the host response 328 to WNV infection. We show here that the zebra finch immune response is largely 329 conserved with that seen in mammalian-based studies (Table 1). Additionally, we 330 identify many components of the immune system that have not been previously 331 We removed Illumina adapters from reads with Trim Galore! v0.3.7 369 (http://www.bioinformatics.babraham.ac.uk/projects/trim_galore/) which makes use of 370 Cutadapt v1.7.1 (56). Reads were then mapped to the zebra finch genome (v3.2.74,6) 371 using TopHat v2.0.13 (57), which utilizes the aligner Bowtie v2. 2.4 (58). We specified 372 the library type as fr-firststrand in TopHat2. Successfully mapped reads were converted 373 from SAM to BAM format with SAMtools View v1. 2 (59,60) and counted in htseq-count 374 v0.6.0 specifying '-s rev' (61). This assigned zebra finch Ensembl gene IDs and we only 375 retained genes that mapped an average of five times across each sample. 376 377

Differential expression 378
Gene counts were then normalized for read-depth and analyzed for DE in DEseq2  conditions. Genes are then grouped into expression paths (i.e. "Up-Down", "Down-405 Down"), in which DE occurs when expression paths change between at least one adjacent 406 condition. For example, a gene up-regulated at both 2dpi relative to control and 4dpi 407 relative to 2dpi would be classified as "Up-Up". We included three time points, with 408 control individuals classified as t1, 2dpi as t2 and 4dpi as t3. Genes were considered DE 409 at posterior probability > 0.99 and FDR < 0.01. We chose a more stringent cutoff in this 410 analysis as EBseq can be liberal in classifying differential expression (65)  The use of Ensembl gene annotations for read counting restricts analyses to 416 previously annotated genes. To test for potential regulation in any unannotated genes, we 417 created a de novo transcriptome from our trimmed, paired-end reads using Trinity v2.0.6 418 (36). We used default Trinity parameters with the exception of server specific settings 419 (e.g. --max_memory 100G) and strand specificity (--SS_lib_type RF). To annotate our 420 transcriptome, we used the BLASTX-like aligner DIAMOND v0.7.11 (37) Figure 1. Immune genes differentially expressed between day 2 post-inoculation and 678 control A) Heatmap of expression levels (log transformed read counts) across all 679 treatments of immune genes differentially expressed at 2dpi relative to control. B-D) 680 Expression values (normalized read counts) for three key immune genes and their 681 regulation pattern classification by EBSeqHMM. Asterisks represent statistical 682 significance in DEseq2 analysis after FDR correction (* p<0.10, ** p<0.05, *** p<0.01). 683 684 685 Figure 2. Immune genes differentially expressed between day 4 post-inoculation and 686 control A) Heatmap of expression levels (log transformed read counts) across all 687 treatments of immune genes differentially expressed at 4dpi relative to control. B-D) 688 Expression values (normalized read counts) for three key immune genes and their 689 regulation pattern classification by EBSeqHMM. Asterisks represent statistical 690 significance in DEseq2 analysis after FDR correction (* p<0.10, ** p<0.05, *** p<0.01). 691 692 693 694 Figure 3. Regulation of the zebra finch RLR pathway. Color represents log 2 fold 695 change between Control and 4dpi. Asterisks represent statistical significance in DEseq2 696 analysis after FDR correction (* p<0.10, ** p<0.05, *** p<0.01). 697 698 699