Genomic and functional analysis of the host response to acute simian varicella infection in the lung

Varicella Zoster Virus (VZV) is the causative agent of varicella and herpes zoster. Although it is well established that VZV is transmitted via the respiratory route, the host-pathogen interactions during acute VZV infection in the lungs remain poorly understood due to limited access to clinical samples. To address these gaps in our knowledge, we leveraged a nonhuman primate model of VZV infection where rhesus macaques are intrabronchially challenged with the closely related Simian Varicella Virus (SVV). Acute infection is characterized by immune infiltration of the lung airways, a significant up-regulation of genes involved in antiviral-immunity, and a down-regulation of genes involved in lung development. This is followed by a decrease in viral loads and increased expression of genes associated with cell cycle and tissue repair. These data provide the first characterization of the host response required to control varicella virus replication in the lung and provide insight into mechanisms by which VZV infection can cause lung injury in an immune competent host.

Scientific RepoRts | 6:34164 | DOI: 10.1038/srep34164 In addition, levels of the potent antiviral cytokine IFNα were increased 3 and 7 DPI (Fig. 2b), while levels of IFNγ , IL-18 (induces IFNγ production by T cells), IL-1β , and its antagonist IL-1RA were increased only 7 DPI (Fig. 2b). Moreover, levels of growth factors VEGF-D and PDGF-BB were increased 3-14 DPI and those of the pro-angiogenic factor VEGF-A were increased 10 dpi (Fig. 2c). In contrast, the lung homogenates showed no significant change in protein levels of chemokines, cytokines or growth factors.

SVV infection leads to the recruitment of T cells and phagocytes.
Given the increased concentration of T cell, B cell and granulocyte chemokines in the BAL, we next characterized the immune cells that infiltrated the lungs during acute varicella (Fig. 3). As we previously reported 7 , CD20 cells were rare whereas CD4 and CD8 T cells were abundant pre-infection in BAL. After SVV infection, the frequency of CD8 T cells in the BAL increased 10 and 14 DPI (Fig. 3a). In contrast, frequencies of both B and T cells were low in naïve lungs, but increased as early as 3 DPI, with CD8 T cells showing the largest increase 10 and 14 DPI (Fig. 3b). T cell activation was indicated by the increased prevalence of the highly differentiated CD8 EM T cells and CD4 CM T cells 7 to 14 DPI and CD4 EM 7-10 DPI (Fig. 3c,d). We also observed an increased frequency of plasmacytoid dendritic cells (pDCs), myeloid dendritic cells (mDCs), and macrophages (MACs) 3-7 DPI in both BAL and lung (Fig. 3e,f). To further characterize the immune infiltrates, immunohistochemistry (IHC) was used to determine changes in distribution and frequency of CD3, CD20, CD68, Ki67 and granzyme B (GRZMB) positive cells (Fig. 4). As detected by flow cytometry (Fig. 3), intensity of CD3, CD20 and CD68 staining increased 3 DPI and peaked 7 DPI (Fig. 4b). Peak staining of GRZMB and Ki67 was also observed at 7 DPI and correlated with increased frequency of memory T cells detected by flow cytometry (Fig. 4b). By 14 DPI, CD20 and granzyme B staining was substantially reduced, while CD3, CD68, and Ki67 expression remained higher than baseline (Fig. 4c).
SVV infection induces robust changes in gene expression. RNA-Seq was used to characterize the host gene expression changes in the lungs during acute varicella. The number of differentially expressed genes (DEGs) correlated tightly with viral loads, with the highest number of DEGs detected 7 DPI (Fig. 5a and Supplemental Datasets S1-4). Changes in gene expression were validated for a subset of genes using qRT-PCR (Supplemental Fig. S1). The 54 DEGs that were differentially expressed 3, 7 and 10 DPI (Fig. 5b) segregated into 3 clusters. Cluster 1 contained 30 DEGs whose expression gradually increased until it peaked 10 DPI (Fig. 5c). Several of these genes are involved in cell cycle such as cyclin kinase subunit 2 (CKS2), Ribonucleoside-diphosphate reductase subunit M2 (RRM2), and topoisomerase II alpha (TOP2A). Cluster 2 contained 7 DEGs whose expression was substantially increased as early as 3 DPI and remained high until 10 DPI before decreasing back to   DEGs 3DPI play a role in innate immunity and lung development. Most of the 106 DEGs that were up-regulated 3 DPI mapped to Gene Ontology (GO) terms associated with antiviral immunity (Fig. 6a). Several genes that enriched to the GO process "immune system process" are part of a network regulated by Signal Transducers and Activators of Transcription 1 (STAT1; FC 4) (Fig. 6b). Most of these DEGs have antiviral function such as interferon-stimulated genes (ISGs) ISG15 (FC 8), MXA (FC 7), and GBP1 (FC 7) as well as dsRNA sensors DDX60 (FC 5) and RIG-I (FC 4) 9 (Fig. 6b)  (a) Bar graph shows the number of genes mapping to each of the 10 most statistically significant GO terms to which the 106 upregulated genes enriched. Line represents the − log(p-value) associated with each GO term. (b) Network of DEGs mapping to the GO process "immune system" that directly interact. (c) Bar graph shows the number of genes mapping to each of the 10 most statistically significant GO processes to which the 54 down-regulated genes enriched. Line represents the − log(p-value) associated with each GO term. (d) Heat map of the 30 downregulated DEGs that enriched to the GO terms described in (c) grouped by the GO term to which they mapped: ED = Embryo development, PSP = pattern specification process, BP = biosynthetic process and LD = Lung development.
Given the immune changes in immune infiltrates into the lung 3 DPI, we wanted to gain a better understanding of which immune cells are expressing these DEGs. To that end, we used the Immunological Genome Project (ImmGen), a database maintained by a collaborative group of immunologists and computational biologists that aim to determine the patterns of gene expression of specific immune cells in the mouse 17,18 . This analysis showed that a group of genes that play a role in the cell cycle progression (e.g. Top2a, BUB1, CEP55, SPC25 and CENPF) were highly expressed by B cells, dendritic cells and T cells (Supplemental Fig. S2, cluster 1). We also identified 2 additional clusters of genes that were mainly expressed by MACs (e.g. ADAMDEC1, IgJ, SIGLEC1, CXCL9 and CXCL10) (Supplemental Fig. S2, cluster 2a,b). Two small clusters are expressed by NK cells (e.g. GZMB, GZMA and SAMD3) (Supplemental Fig. S2 cluster 3a,b), and an additional cluster of antiviral genes was seen in the neutrophils (e.g. DDX60, OAS2, and RASD2) (Supplemental Fig. S2, cluster 4). These data correlate with our flow cytometry analysis (Fig. 3), which report significant increases in T cells, B cells, MACs and DC numbers at 3 DPI in the lung.
In line with our functional analysis, few of the DEGs down-regulated at 3 DPI mapped to immune cell populations. These include FSCN1 (pathogen recognition molecule 19 ) and ARC (apoptosis 20 ) highly expressed in dendritic cells, and SEPN1 (oxidation-reduction homeostasis 21 ), highly expressed in neutrophils (data not shown). DEGs detected 7 DPI primarily map to host defense. The 380 up-regulated DEGs detected 7 DPI primarily mapped to GO terms related to host defense (Fig. 7a). Over 100 DEGs involved in anti-viral immunity had a fold change > 10 including: 30), and IL8 (CXCL8, FC 28) (Fig. 7b). Additional immune related DEGs that didn't map to the 10 most significant GO processes include FAM26F (FC 92), involved in IFNγ production by NK cells 22 , Pentraxin 3 (PTX3, FC 71) involved in complement activation 23 , and tissue inhibitors of metalloproteinases 1 (TIMP1, FC 26). Other highly up-regulated DEGs play a role in respiratory diseases such as Hyaluronan synthase 2 (HAS2; FC 71) associated with asthma and tissue fibrosis 24 , and metallothionein 2A (MT2A, FC 24) involved in lung cancer 25 .
Of the 400 down-regulated DEGs 7 DPI, 300 enriched to GO terms "respiratory tract diseases" and "lung diseases" (Fig. 7c). The most down-regulated gene was keratin 5 (KRT5, FC 35), an intermediate filament protein involved in lung regeneration 26 (Fig. 7d). Down-regulated DEGs with a role in lung function include: insulin receptor substrate 2 (IRS2, FC 9), an anti-inflammatory agent in the lung 27 ; and SEC14-like 3 (SEC14L3, FC 33), a lipid packing sensor expressed by alveolar cells and involved in the prevention of lung collapse 28 .

Discussion
Although VZV is primarily transmitted via inhalation and replicates in the lungs 1 , our understanding of the impact of VZV infection on lung function and the host response to VZV in the respiratory tract remains incomplete. In this study, we leveraged a robust model of VZV infection wherein rhesus macaques are intrabronchially infected with the homologue SVV to carry out the first in depth analysis of the host-pathogen interactions during acute varicella infection in the lungs using a combination of immunohistochemistry, flow cytometry, multiplex ELISA, and RNA sequencing. Our analyses show for the first time that SVV spreads rapidly through the lungs replicating both within the lung parenchyma as well as BAL immune cells. We also report the down-regulation of genes important for lung function, indicative of lung injury during the acute phase of infection. Immune responses and changes in gene expression tightly correlate with viral loads with the up-regulation of pro-inflammatory cytokines, ISGs, chemokines and granzyme genes at peak viral loads. In contrast, genes associated with cell cycle and lung repair were up-regulated after cessation of viral replication. Immunohistochemistry staining and flow cytometry show that SVV infection results in significant infiltration of lymphocytes and MACs in the lungs as early as 3 DPI and before the appearance of varicella rash, typically detected 10-14 DPI. This is in agreement with previous reports of patients who developed varicella pneumonia before the development of the vesicular rash 3,4 . As previously reported for BAL 7 , we saw a higher accumulation of highly differentiated CD8 T cells compared to CD4 T and CD20 B cells within lung biopsies. CD8 T cells were also the most abundant lymphocyte subset in the lungs of mice intra-nasally infected with respiratory syncytial virus (RSV) 43 and rhesus macaques infected with H1N1 influenza 44 . Furthermore, CD8 T cells were shown to be essential for recovery in mice infected with Vaccinia virus 45 . In addition to lymphocytes, frequency of MACs within the lung biopsy also increased 3 and 7 DPI, returning to pre-infection levels 10 DPI. Bioinformatic analysis of DEGs detected 3-10 DPI using ImmGen also indicate that transcriptional changes are mediated by T cells, B cells, DCs, MACs, neutrophils and NK cells 3 and 7 DPI, with a bigger contribution of DCs, MACs, and neutrophils 7 DPI. In contrast, gene expression changes 10 DPI seem to be primarily originating from lymphocytes.
Concomitant with the influx of T cells, B cells, MACs and DCs into the BAL and lung tissue, we detected a significant increase in T cell attractants I-TAC and MIG, B cell chemokine BLC and granyulocyte attractant IL-8 and eotaxin in the BAL 3 and 7 DPI. Levels of pro-inflammatory cytokines IL-18, IL-1β , IL-6 and IFNγ also significantly increased in the BAL. These factors have been reported during other lung infections suggesting a common defense response to respiratory pathogens. For instance, IL-18 has been shown to activate NK cells and stimulate IFNγ production during influenza infection in the lungs 46 . Increases in IL-6 and IFNγ have also been observed in nasopharyngeal lavage samples taken from influenza patients, BAL samples taken from flu infected rhesus macaques, and BAL samples from rhesus macaques infected with pulmonary nontuberculous mycobacteria [47][48][49] . The sources of these factors could be alveolar macrophages and lung epithelial cells as previously reported during influenza infection and lung injury 50,51 . Indeed, ImmGen analysis also shows IL-18 and IL-1β to be highly expressed by MACs. Interestingly, no chemokines or cytokines were increased in the lung homogenate, suggesting that immune mediators are quickly secreted into the alveolar space and do not remain within the lung parenchyma. In support of this hypothesis, mRNA transcripts of CXCL9 (MIG), CXCL10, CXCL11 (ITAC), IL-1β, and IL-8 within lung homogenates were significantly increased. Protein levels of all these chemokines were increased in the BAL at the same time point.
The pDC infiltration 3 DPI correlated with a robust induction of IFNα in the BAL supernatant 3-7 DPI. Another source of IFNα could be epithelial cells, which, like pDCs, express the double stranded DNA sensor TLR9 52 , and could secrete type 1 interferons in a RIG-I dependent manner as described for RSV infection 53 . Alternatively, IFNα could also be secreted by alveolar macrophages, which were the primary producers of IFNα in mice infected with Newcastle disease virus 54 . Concomitant with the increased IFNα levels in the BAL, we detected a significant up-regulation of the transcription factor STAT1, induced by type 1 and 2 interferons 3 DPI 11 . STAT-1 in turn regulates the expression of interferon-stimulated genes (ISGs), which play a critical role in inhibiting viral replication. Indeed, ISGs were among the most highly up-regulated genes in our data set (average FC = 18). The increased expression of ISGs persists until 10 DPI and correlates with decreased viral loads 10 and 14 DPI. Our observations are similar to those described for influenza and RSV where increased expression of ISG15, RIGI, IRF1, CXCL10, CXCL11, IFIT1 and STAT1 have been shown to play a role in disease resolution [55][56][57] . Similarly, reduced ISG expression has been associated with severe cases of H5N1 infection in cynomologous macaques 58 .
STAT1 also regulates the expression of granzymes (GRZM) A, B and K 59 , which were amongst the most highly up-regulated DEGs in our dataset. NK cells and CD8 T cells produce these cytolytic molecules during viral infections to kill infected cells 60 . Increased expression of GRMZ genes correlated with increased frequency of EM CD8 T cells and was confirmed by IHC staining. Bioinformatic analysis using ImmGen also confirmed that GZMA and GZMB was highest in NK and T cells 3-10 DPI. We previously reported increased GRMZB levels within BAL-resident CD4 and CD8 T cells during SVV infection 7 . Elevated GRZMA and GRZMB expression have also been reported during RSV and influenza A infection 61,62 . Moreover, increased expression of GRZMK during acute lung inflammation can in turn induce IL-6, IL-8 and CCL2 production 63 . Both mRNA and protein levels of IL-6 Scientific RepoRts | 6:34164 | DOI: 10.1038/srep34164 and IL-8 were increased during acute varicella 3 DPI. Together with previous reports, data reported in this study suggest that the lung generates a similar response to respiratory viral pathogens lead by ISGs, granzymes, similar chemokines and cytokines followed by an infiltration of B and T cells.
Acute varicella is accompanied by lung injury as evidenced by the severe focal hemorrhaging and damaged alveolar walls as well as the down-regulation of genes involved in lung development and function. This damage is likely caused by both the anti-viral immune response as well as viral replication as previously described for influenza 64 . Indeed, the damage was most severe at the peak of viral loads and immune response (7 DPI). For instance, expression of IRX3 (3 DPI), important for lung development and airspace maintenance 12 and SEC14-like 3 (7 DPI), which plays a critical role in preventing lung collapse 28 , were significantly down-regulated by SVV infection. Similarly, down-regulation of TNNT1 (10 DPI) has been proposed as a cause of acute respiratory distress 65 . Moreover, genes that play a role in controlling lung inflammation were down-regulated such as SCGB3A1 (10 DPI), which regulates airway inflammation and tissue repair 66 . Interestingly, type 1 and 2 interferons, which were increased 3 and 7 DPI, have been shown to decrease SCGB3A1 expression 67 . Similarly, IRS2, which interacts with IL-4 and acts as an anti-inflammatory agent in the lungs 27 , was also down-regulated (7 DPI).
Several tumor suppressor genes were down-regulated, while genes that were involved in the cell cycle were up-regulated, which could have contributed to the increased T cell numbers seen 7-10 DPI. Indeed, analysis using ImmGen indicates that these genes are most highly expressed by T cells. Significant increase in expression of cell cycle genes was also observed in mice recovering from MRSA-induced pneumonia 68 . The up-regulated cell cycle genes could also indicate the beginning of the repair process in the lungs since young immune competent rhesus macaques resolve primary varicella infection without complication. In line with that hypothesis, we detected an increase in growth factors VEGF and PDGF-BB, and angiogenic factor VEGF-A 10 and 14 DPI, which have been shown to play a critical role in alveolarization in the lungs and repairing damaged alveolar walls 69,70 . Interestingly, several genes involved with lung function remained down-regulated 10 DPI and 14 DPI, suggesting that repair continues after cessation of viral replication. For instance, KRT5, essential for lung tissue repair 26 , remained one of the most significantly down-regulated genes at 10 DPI, and DNAH2, involved with cilium motility in the lungs 40 also remained down-regulated at 14 DPI.
In summary, data from this study provides the first kinetic analysis of host-pathogen interactions within the lungs during acute varicella infection. A robust immune response led by interferon-stimulated genes, granzymes, cytokines and chemokines is critical for the resolution of infection. These data also show that although SVV infection in young macaques is self-limiting, significant lung injury occurs during the acute phase as a result of the immune response and viral replication, a phenomenon that has not previously been described. Our results provide potential insights to both what happens in a self -resolving infection, and by extension, what may be happening when infection results in severe complications. Specifically, our RNA sequencing data show reduced expression of genes important for lung function during peak viral replication. In immune deficient individuals who lack the vigorous immune response, uncontrolled viral replication could lead to sustained reduction in the expression of these critical genes, which in turn results in severe lung damage and viral pneumonia. Furthermore, immune genes that are significantly up-regulated likely play a critical role in resolving infection; therefore, lack of increased expression of these genes may explain why the immunocompromised develop severe disease. On the other hand, adults may generate too vigorous of an immune response that results in severe inflammation and viral pneumonia. These data help design interventions to mitigate complications associated with VZV pneumonia and other infectious respiratory diseases.

Methods
Animals and sample collection. Fourteen colony-bred Rhesus macaques (Macaca mulatta, RM) 3-5 years of age and of Indian origin were used in these analyses. Animals were housed and handled in accordance with the Oregon National Primate Research Center (ONPRC) Institutional Animal Care and Use Committee (IACUC protocol #0779). The ONPRC Institutional Animal Care Use Committee approved all experimental protocols. The ONPRC has been continuously accredited by the American Association for Accreditation of Laboratory Animal Care since 1974 (PHS/OLAW Animal Welfare Assurance #A3304-01). Eleven RM were inoculated with 4 × 10 5 PFU SVV as previously described 6 and euthanized 3 (n = 3), 7 (n = 3), 10 (n = 2), and 14 (n = 3) days post-infection (DPI); and 3 were controls. Animals were either housed single or paired in caging that allowed for social interactions in a temperature and humidity controlled environment. Food and water were available ad libitum and enrichment was provided daily. All procedures were done under Ketamine anesthesia to minimize animal suffering. Necropsy was carried out in accordance with the recommendation of the American Veterinary Association guidelines for euthanasia. Blood, bronchial alveolar lavage (BAL) and lung tissues were collected at necropsy. Tissues were then flash frozen or stored in trizol at − 80°. mRNA Library Preparation. One microgram of RNA was used to generate RNA libraries using the using the New England Biolab (NEB) Next Ultra Direction RNA Prep kit for Illumina (Ipswich, MA) and AMPure XP beads (Beckman Coutler, Brea, CA). Each library was made with a unique index primer to allow for multiplexing and then sequenced on the Illumina HiSeq2500 (Illumina, San Diego, CA) platform at single-ends 100bps as previously described 71 .
RNA-Sequencing analysis. Sequence analysis was performed as previously described 71 . Due to the limited number of animals used in the this study, use used a stringent cutoff with differentially expressed genes (DEGs) being defined as those with a fold change ≥ 3 and a false discovery rate (FDR) ≤ 0.05 compared to naïve animals. Enrichment analysis was performed using MetaCore software (GeneGo, Philadelphia, PA). Gene clusters were predicted by STEM using k-means clustering of median normalized reads and maps were generated using ggplot package in R. GEO accession number SRP072466. Gene Validation. Gene validation was done using RNA from the same samples used for RNA-Seq. RNA was reverse-transcribed using 100 mM dNTPs (Applied Biosystems) and multiscribe reverse transcriptase (Applied Biosystems) to produce cDNA. Expression was determined using Taqman gene expression assays (Thermo Fisher, Waltham, MA) of BST2, IFITM1, IGJ, MX1 and a housekeeping gene (RPL32). Taqman gene expression assays with 50 ng of cDNA and a housekeeping gene (RPL32) was carried out in duplicates on the ABI StepOne instrument (Applied Biosystems,). Expression of mRNA was normalized to expression levels of RPL32 using Δ Ct values as previously described 71 . Statistical Analysis. Statistical analysis was conducted using GraphPad Prism software (GraphPad, Software, Inc., La Jolla, CA). Significant values were determined using one-way ANOVA with an alpha value of 0.05 or less.