Simian varicella virus causes robust transcriptional changes in T cells that support viral replication.

Varicella zoster virus (VZV) causes varicella (chickenpox) during acute infection. Several studies have shown that T cells are early and preferential targets of VZV infection that play a critical role in disseminating VZV in to the skin and ganglia. However, the transcriptional changes that occur in VZV-infected T cells remain unclear due to limited access to clinical samples and robust translational animal models. In this study, we used a nonhuman primate model of VZV infection where rhesus macaques are infected with the closely related Simian Varicella Virus (SVV) to provide novel insights into VZV-T cell interactions. RNA sequencing of bronchial alveolar lavage-resident T cells isolated from infected rhesus macaques show that SVV infection alters expression of genes important for regulation of gene expression, cell cycle progression, metabolism, and antiviral immunity. These data provide insight into cellular processes that may support viral replication, facilitate SVV dissemination, and evade host defense.

Introduction

Varicella Zoster Virus (VZV) is a highly contagious neurotropic alpha herpes virus that causes varicella (chickenpox) during acute infection. VZV is spread through the inhalation of virus particles or through direct contact with infectious fluid from skin lesions (Arvin, 2001, Grose, 1981). VZV establishes latency in sensory ganglia, from which it can reactivate causing herpes zoster (shingles), a painful disease most commonly affecting the elderly and immunocompromised individuals. Previous studies have shown that T cells are highly susceptible to VZV infection and may play a critical role in its dissemination to the skin and ganglia. Specifically, direct inoculation of fetal thymus/liver implanted under the kidney capsule of SCID-hu mice with VZV-infected fibroblasts results in the infection of CD4+ and CD8+ thymocytes (Moffat et al., 1995). In addition, in vitro experiments have shown that VZV has a high propensity to infect tonsillar T cells (Ku et al., 2002). Moreover, co-culture experiments showed that activated tonsillar CD4 T cells with skin homing markers were more likely to be infected with VZV (Ku et al., 2002). Importantly, human tonsillar CD4 T cells infected in vitro with VZV, but not fibroblasts, intravenously injected into SCID-hu mice were able to transport VZV to fetal human skin explant resulting in development of varicella rash (Ku et al., 2004) and fetal dorsal root ganglia xenografts (Zerboni et al., 2005).

Although these in vitro studies suggest that T cells play a critical role in VZV pathogenesis, none of these findings have been confirmed using T cells isolated from varicella patients. Moreover, the strict host tropism of VZV has precluded the development of animal models. An alternative is to use rhesus macaques intra-bronchially infected with the homologous simian varicella virus (SVV). This model mimics the key characteristics of VZV infection including the development of varicella, cellular and humoral immune responses, the establishment of latency in sensory ganglia, and reactivation (Messaoudi et al., 2009, Mahalingam et al., 2010, Mahalingam et al., 2007, Kolappaswamy et al., 2007). As described for VZV, we and others have demonstrated that SVV primarily infects T cells that traffic to the ganglia as early as 3 days post-infection (Ouwendijk et al., 2013, Arnold et al., 2016a). Moreover, we showed that T cells isolated from the broncho-alveolar lavage (BAL) during acute infection supported SVV replication (Arnold et al., 2016a). These data firmly establish the importance of T cells in SVV pathogenesis making this model ideal for investigating how VZV infection alters T cell behavior and function. In this study, we used this animal model to investigate the transcriptional changes that SVV infection induces within CD4 and CD8 T cells isolated from the BAL during acute infection. Our results show that SVV induces robust transcriptional changes involved with chromatin assembly, translation, cell cycle and cellular metabolism. In addition, several gene expression changes reveal possible mechanisms by which SVV may evade the host response.

Methods and materials

Animals and samples

CD4 and CD8 T cells were isolated from bronchial alveolar lavage (BAL) samples collected from 8 colony-bred Indian origin Rhesus macaques (Macaca mulatta, RM) 3–5 years of age previously described in (Haberthur et al., 2013) using Magnetic Cell Sorting (MACs). Four of these animals were inoculated intrabronchially with 4 × 105 PFU wild-type SVV as previously described and were euthanized 3, 7 and 10 days post infection (DPI) (Table 1 ) (Haberthur et al., 2013). Briefly, BAL cells were first incubated with CD4 microbeads (Miltenyi Biotec, San Diego) for 20 mins at 4 °C. Cells were then washed and applied to the magnetic column to capture CD4 T cells. The negative fraction was then stained with anti-CD8 PE (Beckman Coulter) for 20 mins in the dark at 4 °C followed by a wash and incubation with anti-PE beads (Miltenyi Biotec) for 15 mins in the dark at 4 °C. Cells were then washed and the cell suspension was applied to a second magnetic column to isolate CD8 T cells. To determine purity, cells were stained with antibodies directed against surface markers CD4 (Tonbo Biosciences, San Diego, CA) and CD8β (Beckman Coulter, Brea, CA). Samples were acquired using an LSRII instrument (Becton, Dickinson and Company, San Jose, CA) and analyzed using Flowjo software (TreeStar, Ashland, OR). All samples used had a purity of >90%.

Table 1

Samples used for RNA-Seq Analysis.

Animal ID	DPI	Cell Type	Viral Load
28081	3	CD8	411,559
28081	3	CD4	320,688
17036	7	CD8	10,674
27180	10	CD8	404
27916	10	CD8	517
Control #1	UNIF	CD8	0
Control #2	UNIF	CD8	0
Control #3	UNIF	CD8	0
Control #4	UNIF	CD8	0
Control #1	UNIF	CD4	0
Control #2	UNIF	CD4	0
Control #3	UNIF	CD4	0
Control #4	UNIF	CD4	0

RNASeq and bioinformatics

RNA was extracted from purified CD4 and CD8 T cells using the Ambion Purelink RNA Mini Kit extraction kit (Life Technologies, Carsbad, CA). RNA library preparation was done using the New England Biolab (NEB) Next Ultra Direction RNA Prep kit for Illumina (Ipswich, MA). We were unable to generate a library of the CD4 T cells isolated from the 7 DPI and 10 DPI animals. Therefore, we only have paired CD4/CD8 libraries from the animal euthanized 3 DPI and our control animals; and only a CD8 T cell library from the animals euthanized 7 DPI and 10 DPI (Table 1). Libraries were next multiplexed and sequenced on the Illumina NextSeq (Illumina, San Diego, CA) platform at single-ends 75bps. All data analysis steps were performed with the RNA-Seq workflow module of the systemPiperR package available on Bioconductor (Huber et al., 2015, Girke, 2015). Next generation sequencing quality reports were generated with the seeFastq function defined by the same package. Reads were mapped with the splice junction aware short read alignment suite Bowtie2/Tophat2 (Langmead and Salzberg, 2012, Kim et al., 2013) against the Macaca mulatta genome sequence downloaded from Ensembl (Cunningham et al., 2015). The default parameters of Tophat2 optimized for mammalian genomes were used for the alignments. Raw expression values in the form of gene-level read counts were generated with the summarizeOverlaps function(Lawrence et al., 2013). Only reads that overlapped with exonic regions of genes were counted, while reads that mapped to ambiguous regions of exons from overlapping genes were discarded.

Analysis of differentially expressed genes (DEGs) was performed with the GLM method from the edgeR package (Robinson et al., 2010, Anders et al., 2013). Differentially expressed genes (DEGs) were defined as those with a fold change of ≥2, a false discovery rate (FDR) of ≤0.05 and an average reads per kilobase per million (RPKM) value ≥5. These FC and FDR cut off values are well-established statistical criteria for transcriptomic experiments (Conesa et al., 2016). Functional enrichment analysis was performed using MetaCore™ software (GeneGo, Philadelphia, PA) to identify significant gene ontology (GO) and Network Processes. When conducting the longitudinal analyses, an FDR cut-off could not be applied for the 7 DPI data set due to the fact that only 1 sample was analyzed. Heatmaps, volcano plot, and PCA were generated using gplots, ggplot2 and DESeq2 package in R.

Validation of gene expression changes

Remaining RNA from samples was reverse transcribed using random hexamers and SuperScript ® IV RT in the SuperScript ® IV First-Strand Synthesis System (Invitrogen, Lithuania) to generate cDNA. Taqman gene expression assays (Thermo Fisher, Waltham, MA) of candidate genes and housekeeping gene (RPL32) were used with 50 ng of cDNA and carried out in duplicate on the ABI StepOne instrument (Applied Biosystems). mRNA expression levels were calculated relative to the housekeeping gene (RPL32) using 2−ΔCt calculations.

BAL cells collected 0 and 7 DPI during the course of another previous study (Arnold et al., 2016b) were stained with antibodies against CD8β (Beckman Coulter), CD4 (eBioscience, San Diego, CA), CD3 (Biolegend, San Diego, CA), CD2 (Biolegend), TLR4 (ThermoFisher, Waltham, MA), and IGFR2 (ThermoFisher). The samples were analyzed using the Attune NxT (ThermoFisher). Mean Fluorescence Intensity (MFI) was used to calculate FC of protein expression levels 7 DPI compared to 0 DPI.

Results

BAL-resident t cells isolated from SVV infected animals show robust transcriptional changes

We recently reported that SVV primarily infects both CD4 and CD8 T cells in the BAL, which can support viral replication (Arnold et al., 2016a). To determine SVV-induced transcriptional changes within T cells, we performed RNA-Seq on BAL-resident T cells shown to express viral transcripts (measured through AmpliSeq technology (Arnold et al., 2016a)) and uninfected BAL T cells (4 paired CD4/CD8 T cell samples isolated from 4 control animals) (Table 1). A principle component analysis (PCA) shows that infected and uninfected T cells have distinct transcriptional profiles (A).

We first identified global differences between the two groups of samples. Overall, a total of 1463 (739 upregulated and 724 downregulated) differentially expressed genes (DEGs) were detected, of which 1130 (586 upregulated and 544 downregulated) were characterized (Fig. 1 B, Supplemental Dataset 1). Changes in expression of 4 DEGs were confirmed using qRT-PCR and that of an additional 5 DEGs were confirmed at the protein level using flow cytometry (Supplemental Fig. 1). Upregulated DEGs enriched to network processes related to inflammation and cell cycle (Fig. 1C) and gene ontology (GO) processes related to immunity, metabolism and gene expression (Fig. 1D). Downregulated DEGs enriched to process networks involved with phagocytosis and signaling (Fig. 1C) and GO processes involved with intercellular transport and metabolism (Fig. 1D).

Fig. 1

SVV infection results in robust gene expression changes in T cells. (A) PCA plot of T cell samples subjected to RNASeq. Each dot represents a CD4 or CD8 T cell sample used for RNA-Seq analysis and each color represents a timepoint. (B) Volcano plot representing overall gene expression changes observed in infected T cells. Each gene is denoted by a dot with red dots representing genes with significant difference in expression in infected compared to naive T cells (FDR p-value <0.05 and fold change >2). Y-axis represents the FDR corrected p-value and X-axis represents fold change. (C-D) The most statistically significant upregulated (in blue) and downregulated (in red) (C) Process Networks and (D) GO terms to which differentially expressed genes (DEGs) enriched. Bars represent the number of DEGs enriching to each process network/GO term and lines represent the FDR −log10(p-value) associated with each enrichment.

SVV infection upregulates cell cycle processes and immune responses

Many of the upregulated genes, played a role in gene regulation, notably chromatin assembly (Fig. 2 A) and translation (Fig. 2B). DEGs important for chromatin assembly were primarily detected 10 DPI (Supplemental Fig. 2A) and include histone genes such as HIST1H2AJ (FC 13), HIST1H1B (FC 11), HIST2H2AA4 (FC 11), HIST1H3A (FC 10) and H2AFX (FC 9); and transcription factors C-FOS (FC 4), JUNB (FC 6) and JUND (FC 5) (Halazonetis et al., 1988) (Fig. 2A). DEGs that play a role in translation were primarily upregulated 3 and 7 DPI (Supplemental Fig. 2B) and include ribosomal protein L (RPL) and ribosomal protein S (RPS) genes that are important for 60S and 40S ribosomal assembly respectively (Fig. 2B, Supplemental Fig. 2B) (Provost et al., 2013); and eukaryotic translation initiation factor (EIF) subunits, which play critical roles in ribosome stabilization and the initiation of translation (Jackson et al., 2010) (Fig. 2B).

Fig. 2

Enrichment analysis of upregulated DEGs. (A–C) Heatmap of the upregulated DEGs that enriched to (A) chromatin assembly, (B) translation initiation, and (C) TCR Signaling. Each column of the heatmap represents the median reads per kilobase per million reads (RPKM) of the samples analyzed 3 DPI (D3), 7 DPI (D7), 10 DPI (D10) and uninfected (UN) T cells. (D) Network image of the DEGs that mapped to the GO process “viral process” and directly interact with each other.

In addition to regulation of gene expression/cell cycle progression, several upregulated DEGs play a role in immunity and inflammation (Fig. 1C,D), such as: cytolytic molecule granzyme B (GZMB, FC 20) (Lord et al., 2003); interferon stimulated genes ISG20 (FC 10), IFI17 (FC 4), MxA (FC 6) and OAS2 (FC 6) (Schoggins and Rice, 2011); and transcription factors IRF1 and IRF8 (FC 5) (Honda and Taniguchi, 2006). Moreover, several DEGs upregulated 3 and 7 DPI (Supplemental Fig. 2C) play a critical role in “TCR signaling” (Fig. 2C, Supplemental Fig. 2C), notably ZAP70 (FC 8), AKT (FC 4), ITGB2 (FC) CD3ε (FC 4), CD3δ (FC 3), and CD3γ (FC 3) (Ostermann et al., 2002, Paz et al., 2001, Xue et al., 2008, Yamazaki et al., 1999, Zhu et al., 1998) (Fig. 2C). Several of the 91 genes that mapped to the “viral process” (Fig. 2D) include several transcription factors that regulate antiviral immune processes such as NFκB (FC 3), STAT1 (FC 2) and STAT2 (FC2). Interestingly, some of the other DEGs in this network were previously reported to support viral replication such as HCF1 (FC 3, alpha herpesviruses (Kristie et al., 2010)), TGFβ (FC 3, Epstein-barr virus (Fahmi et al., 2000), PLCG1 (FC 4, simian virus 40 (Butin-Israeli et al., 2010)) and SMARCB1 (FC 2, HIV-1(Kalpana et al., 1994)).

SVV infection downregulates genes involved in antigen presentation and metabolic processes

Most of the downregulated genes were detected at 10 DPI and to a lesser extent 3 DPI, and mapped to the GO term “localization” (Fig. 1D). Fifty-four of these genes are part of a network regulated by transcription factors PPARγ (FC 8) and TCF/LEF (FC 10) (Fig. 3 A) that contained chemokines IL-8 (FC 29), CCL23 (FC 5) and CXCL16 (FC 3) (Fig. 3A). Enrichment analysis using network processes showed that several downregulated genes play a role in phagocytosis and antigen presentation (Fig. 1C and Supplemental Fig. 2D). These genes included pathogen recognition receptors (PRRs) CLEC7A (FC 4) (Brown, 2006), MSR1 (FC 5) (Dansako et al., 2013), TLR4 (FC 3) (Lu et al., 2008); genes involved in vesicle transport such as SPIRE1 (FC 8) (Kerkhoff et al., 2001), EXOC5 (FC 2) and EXOC6B (FC 3) (Heider and Munson, 2012); and receptors that mediate viral entry such as ITGAV (HSV-1, FC 4) (Gianni et al., 2013), SCARB1 (hepatitis C, FC 2) (Zeisel et al., 2007), and MRC1L1 (dengue virus, FC 6) (Miller et al., 2008) (Fig. 3B). A large number of downregulated DEGs were involved with metabolic processes (Fig. 3C and Supplemental Fig. 2E) such as glycerolipid biosynthesis (GPAM, FC 11), lipoprotein metabolism (APOC2, FC 12) (Bengtsson-Olivecrona and Sletten, 1990), insulin signaling (ENPP1, FC 14) (Kang et al., 2014) and calcium signaling (CACNB4, FC 6).

Fig. 3

Enrichment analysis of downregulated DEGs. (A) Network image of the genes that mapped to the GO process “localization” and directly interact with each other. (B-C) Heatmap of the DEGs that mapped to (B) phagocytosis/phagosome antigen presentation and (C) of the 30 most downregulated genes that mapped to “single organism metabolic process”. Each column of the heatmap represents the average reads per kilobase per million reads (RPKM) of the samples analyzed 3 DPI (D3), 7 DPI (D7), 10 DPI (D10) and uninfected (UN) T cells.

SVV-induced gene expression changes detected between 3 and 10 DPI in t cells show significant overlap

Since we used samples collected at different DPI, we next carried out a longitudinal analysis of gene expression in infected T cells. Overall, the smallest number of DEGs was detected at 3 DPI (887) while the largest number of DEGs was detected 10 DPI (1674) (Fig. 4 A). Enrichment analysis of the 417 upregulated DEGs detected 3 DPI revealed that these DEGs play a role in gene expression and immunity (Table 2 ). Specifically, DEGs enriched to process networks: “translation initiation” which included several ribosomal protein L (RPL) and S (RPS) genes; “cell cycle” (e.g CDK4 (FC 7) AKT1(FC 3) and MAP2K2 (FC 4)); “chromatin modification” (e.g. HDAC5 (FC 2), MBD3 (FC 3) and HIST1H1C (FC4)); TCR signaling (e.g. CD3 (FC 4), CD8 (FC 9), NF-AT (FC 4), and ZAP70 (FC 11)); and lymphocyte proliferation (e.g. CCL5 (FC 6), IL-2R (FC 9) and IL-4R (FC 2) (Table 2). Additional functional enrichment analysis revealed that several upregulated DEGS at 3 DPI mapped to GO processes related to the regulation of gene expression (e.g. SBNO2 (FC4), MAF1(FC 2), SMARCC2 (FC 2) and ERN1 (FC 4)) and antiviral immunity (such as IFITM1 (FC 4), IRF1 (FC4), GRMB (FC 17), and TNFα (FC3)). The 470 downregulated DEGs detected 3 DPI played a role in phagocytosis (e.g. CLEC7A (FC6), MSR1 (FC 6) and SYK (FC2)); antigen presentation (e.g. SIRPA (FC 3) ACTN1 (FC4) and FGR (FC3)); cell adhesion (e.g. CDH6 (FC 9), DAB2 (FC5) and PVR (FC 5); and cell localization (such as CD163 (FC 9), CD9 (FC 3), and IL-8 (FC 12)). A large number of downregulated DEGs mapped to the GO process “response to stress” and include genes that play a role in immunity as well as apoptosis (e.g. CARD8 (FC 4), CARD9 (FC 4) and HIPK2 (FC 6)) and cellular metabolism (e.g. HK3 (FC 5), GPAM (FC 11) and GK (FC 4)). The 86 DEGs detected only 3 DPI (Fig. 4A) include several that play a role in: 1) the immune response such as the chemokine CCL5 (FC6, upregulated); 2) chromatin modification such as HIST1H1AH (FC 6, upregulated); and the regulation of the cell cycle such as E4F1 (FC5, upregulated) (Supplemental Table 1).

Fig. 4

Longitudinal gene changes in SVV-infected T cells. (A) 3-way Venn diagram of the DEGs detected at 3, 7 and 10 DPI compared to uninfected samples. (B) Short Time-series Expression Miner (STEM) clustering analysis of the most statistically significant clusters of upregulated (B-C, Clusters 1 and 2) and downregulated (C-E, Clusters 3–5) DEGs detected 3, 7 and 10 DPI.

Table 2

Gene Enrichment analysis on the DEGs detected at 3, 7,and 10 DPI compared to uninfected T cells.

DEGs detected at 3 DPI
Upregulated
Process Networks	# Genes	FDR p value
Translation initiation	26	3.660E-08
Immune response-TCR signaling	26	3.660E-08
Cell Cycle-G1-S growth factor regulation	25	1.334E-06
Transcription-Chromatin Modification	17	6.508E-05
Lymphocyte Proliferation	20	7.464E-04
Go Processes
Negative Regulation of RNA Metabolic Process	118	5.663E-31
Negative Regulation of Gene Expression	128	3.391E-29
Viral Process	81	1.098E-18
Immune System Process	139	2.211E-18
Downregulated
Process Networks	# Genes	FDR
Immune Response-Phagocytosis	26	4.532E-07
Phagosome in Antigen Presentation	19	7.98E-03
Cell Adhesion-Cadherins	17	7.170E-03
Inflammation IFN-gamma Signaling	10	3.383E-02
GO Processes
Localization	231	3.567E-15
Vesicle Mediated Transport	92	3.649E-12
Endocytosis	52	4.411E-10
Response to Stress	169	1.740E-09

Functional enrichment analysis of the DEGs detected 7 DPI (1120 total DEGs) was comparable to that reported for 3 DPI (Table 2), in line with the large overlap in DEGs observed between these two timepoints (Fig. 4A). Most of the 452 upregulated DEGs play a role in the immune response and enriched to process networks/GO terms: “antigen presentation” (e.g. STAT1 (FC 3), and CD2 (FC 4)); “leukocyte chemotaxis” (e.g. CXCR4 (FC2)); “viral process” (e.g MX1 (FC 5) and BST2 (FC 4)); and “programmed cell death” (e.g. FAS (FC3), and TNF (FC 4)). The 668 downregulated genes at 7 DPI also enriched to similar process networks/GO terms as those reported 3 DPI: “phagocytosis”; “IFN-gamma signaling”; “localization”; and “single organism metabolic process”. Interestingly, some of the downregulated genes mapped to “leucocyte chemotaxis” such as CCL23 (FC5) and IL-8 (FC8); and “positive regulation of defense response” such as CD1D (FC 7) and IL-1 (FC 3). The 388 DEGs detected only 7 DPI (Fig. 4A) enriched to similar process networks and GO terms reported for the other timepoints including: antigen presentation, cell death, and cell cycle (Table 2).

As described for DEGs detected 3 and 7 DPI, the upregulated and downregulated genes detected 10 DPI enriched to process networks and GO terms: “immune response”, “cell cycle”, “translation initiation”, “chromatin organization”, “localization”, and “antigen presentation”. In addition, some DEGs mapped to the process networks associated with cytoskeleton rearrangement and apoptosis. In contrast to DPI 3 and 7, the 771 DEGs detected only 10 DPI did not enrich to GO or network processes related to immunity. However, as described for DPI 3 and 7, they enriched to process networks: cell cycle (e.g ANAPC10 (FC 4, downregulated), DNM2 (FC 2, upregulated) and NUMA1 (FC3, upregulated)); and chromatin remodeling (e.g. BRD9 (FC 4, upregulated) and BAZ2A (FC2 2, upregulated)). Additionally, DEGs downregulated 10 DPI enriched to GO terms “cellular metabolic process” (Table 2).

In addition, the 475 DEGs that were shared amongst all 3 timepoints enriched to similar process networks and GO terms as those described for the global and longitudinal analyses (Table 2). To identify clusters of DEGs with comparable changes across time, we performed a Short Time-series Expression Miner (STEM). We identified 5 statistically significant gene clusters. Two of these clusters contained up-regulated genes: 1) a cluster of 50 genes that play a role in chromatin organization and RNA metabolism (Fig. 4B, Cluster 1, Supplemental Dataset 2); and 2) a second cluster of 107 genes that regulate GTPase activity and transcription (Fig. 4C, Cluster 2, Supplemental Dataset 2). Three clusters of downregulated genes were also identified: 1) one cluster of 105 genes that play a role in apoptosis, cell adhesion and cellular homeostasis (Fig. 4D, Cluster 3, Supplemental Dataset 2); 2) a second cluster of 80 genes that play role in cell proliferation and lipid metabolism (Fig. 4E, Cluster 4, Supplemental Dataset 2); and 3) a third cluster of 52 downregulated genes that play a role in cell proliferation and intracellular transport (Fig. 4F, Cluster 5, Supplemental Dataset 2).

Discussion

Although many studies have demonstrated a role for T cells in varicella pathogenesis (Moffat et al., 1995, Ku et al., 2002, Ouwendijk et al., 2013, Schaap et al., 2005), the impact of VZV infection on T cell function remains poorly understood due to limited access to T cells infected in vivo. In this study, we leveraged a robust animal model where rhesus macaques are intrabronchially infected with the closely related SVV to address this question. We characterized SVV-induced changes in the trancriptome of T cells isolated from the BAL that were previously shown to harbor SVV transcripts (Arnold et al., 2016a). Our RNASeq analysis revealed large gene expression changes between T cells isolated from infected and naïve animals.

Several of the upregulated genes were involved with the regulation of gene expression. Several of these genes encode ribosomal components. Ribosomes play a critical role in protein synthesis and have indeed been found to support viral protein synthesis (Walsh et al., 2013). Specifically, RPL40 supports vesicular stomatitis virus translation (Lee et al., 2013) and EIF4 supports translation of herpes simplex 1 (HSV-1) (Walsh and Mohr, 2004) and human cytomegalovirus (HCMV) (Walsh et al., 2005). Many histone genes, which regulate chromatin assembly, were also upregulated after SVV infection. HIST1H3 and HIST1H1 in our dataset were also shown to support replication of other DNA viruses such as adenoviruses and HSV-1 (Komatsu and Nagata, 2012, Conn et al., 2008, Lieberman, 2008). Moreover, the inhibition of histone demethylase has been shown to block HSV and VZV replication (Liang et al., 2009) and VZV ORF 63 has been shown to regulate histone levels in infected cells (Ambagala et al., 2009). Several of the other upregulated transcription factors that were involved with chromatin assembly such as CFOS, JUNB and JUND were reported to be upregulated with other viral infections such as kaposi’s sarcoma-associated herpesvirus (Xie et al., 2008), Epstein-barr virus (Song et al., 2004) measles (Helin et al., 2002), JC virus (Sadowska et al., 2003) and echoviruses (Huttunen et al., 1997).

Expression of genes that regulate the cell cycle were also altered with SVV infection in T cells. Viruses have been shown to promote infected cells to enter the cell cycle in order to produce nucleotides that are used for viral replication (Bagga and Bouchard, 2014), while inhibiting cell cycle progression in order to maintain the cell in a state that supports viral replication and prevents apoptosis. For example, influenza A, severe acute respiratory syndrome coronavirus (SARS-CoV), and Epstein Barr virus induce GO/G1 arrest (Jiang et al., 2013, Yuan et al., 2006) while HIV-1 induces G2/M arrest (Goh et al., 1998). Furthermore, genes that mapped to “viral process” such as SMARCB1 and TGFβ have been shown to support replication of other viruses, such as HIV (Kalpana et al., 1994) and Epstein-barr virus (Fahmi et al., 2000) respectively. These data suggest that viruses that replicate in lymphocytes may manipulate similar cellular processes to support viral replication.

Functional enrichment revealed that a large portion of DEGs were involved with metabolic processes. For instance, genes important for calcium signaling such as CAMK4 and CACNP were upregulated in our dataset and were previously reported to enhance HSV1/2 and HIV infection by activating the inositol 1,4,5-triphosphate pathway (Cheshenko et al., 2003, Chami et al., 2006). More importantly, calcium treatment of keratinocytes has also been shown to increase VZV gene expression (Jones et al., 2014). In addition, we saw a downregulation of insulin signaling (e.g. ENPP1 FC14), IRS2 (FC3) and IGFR2 (FC2)). Similarly, Hepatitis C virus has been shown to reduce insulin signaling in infected cells (Gao et al., 2015). Moreover, as previously described for HIV-infected T cells (Rasheed et al., 2008), we observed a down-regulation of genes that regulate lipid metabolism (e.g. APOC2 (FC 12), PPARγ (FC 8)). These observations suggest that SVV, like other viruses, may dysregulate T cell metabolism in order to support successful viral replication. However, reduced expression of genes important for metabolic processes could be part of the host defense response to SVV. For instance, genes involved with nitric oxide regulation (e.g, EDF1, FC 3) known to inhibit viral replication, were upregulated in our data set (Reiss and Komatsu, 1998). In addition, fatty acid synthesis genes (e.g. ACP5 (FC3), FABP5 (FC2) and FAR1 (FC2)) were found to be downregulated in our dataset. Since fatty acid acylation has been shown to support the formation of VZV virions by facilitating production of VZV glycoproteins (Namazue et al., 1989), downregulation of these genes could be a host defense mechanism aimed at limiting viral egress.

SVV infection also resulted in the dysregulation of several immune genes within T cells. Most of the upregulated immune genes were involved in the antiviral innate immune response. For instance, several interferon stimulated genes (ISGs), which play a critical role in inhibiting various steps of the virus life cycle (Schoggins and Rice, 2011), were upregulated in our infected T cell population including ISG20 (FC 11), IFITM1 (FC 4) and IFIT2 (FC 3). Up-regulation of ISGs could be in response to the robust up-regulation of IFNα in BAL of SVV infected macaques, which we previously reported (Meyer et al., 2013). Other immune genes that were upregulated include cytolytic molecules and components of the TCR signaling complex. Our results are in line with those from previous in vitro studies that investigated how VZV infection affects the phenotype of human tonsil T cells using mass cytometry (CyTOF). Those studies reported increased expression of ZAP70, AKT, and CD3 (Sen et al., 2015). SVV-induced T cell activation may increase viral gene expression and enhance T cell trafficking to the ganglia/skin as previously reported for VZV (Ku et al., 2002, Ku et al., 2004, Zerboni et al., 2005).

SVV infection led to the down-regulation of immune genes that play a role in phagocytosis and antigen presentation. These included pathogen recognition receptors (PRR) such as CLEC7a and MSR1 which were also downregulated following influenza infection (Yu et al., 2011). The down-regulation of these genes may be indicative of viral immune evasion strategies employed by SVV. Indeed, down-regulation of antigen presentation has been reported in VZV infected T cells (Abendroth et al., 2000, Abendroth et al., 2001). Other downregulated genes in our dataset were reported to act as entry receptors for other viruses; specifically, SCARB1, MRC1L1 and ITGAV have been reported to facilitate entry of hepatitis C, dengue virus and HSV1 respectively. Their reduced expression could be a strategy employed by the host to block viral entry. Expression of several chemokines such as CCL5, IL-8, CCL23, and CXCL16 was also reduced in our infected T cell population, which could also indicate a strategy for SVV to inhibit the recruitment of other immune cells to the site of infection. Similarly, human herpesvirus-6A (HHV6A) has been shown to downregulate expression of CCL5 (Catusse et al., 2008).

We recently analyzed gene expression changes in lung biopsies collected longitudinally after intrabronchial SVV infection (Arnold et al., 2016b). Comparison of that dataset to the one included in this manuscript revealed very few common genes (1.9–5.5%) (Supplemental Fig. 3). This is most likely due to the fact that lung biopsies are primarily composed of type I and II pneumocytes (epithelial cells) in addition to infiltrating immune cells whereas the current study focused exclusively on purified BAL-resident T cells which represent a small fraction of the cells within a lung biopsy. In contrast, the current gene expression data set shares similarity with a previous gene expression study that investigated the impact of in vitro VZV infection on human T cells using microarray analysis and reported increased expression of genes involved with antiviral immunity, cell cycle, transcription and translation .e.g. MXA, JUNB and NRP2 (Jones and Arvin, 2003). The limited similarities of our study with this earlier study could be due to the use of the attenuated Oka vaccine strain and the in vitro infection of the earlier study, which excludes other environmental factors that our T cells were exposed to in vivo.

Enrichment analysis of gene changes at each individual timepoint showed significant overlap between the timepoints and the overall comparison between infected and naïve T cell populations strongly suggesting that SVV-infection exerts the largest influence over transcriptional changes rather than the DPI. At each timepoint, several upregulated genes enriched to immune response and play a critical role in anti-viral immunity. The presence of these DEGs is most likely due to the fact that we analyzed a heterogeneous T cell population that included both infected and antigen-specific responding T cells. Changes in expression of these genes are consistent with the development of an immune response against SVV beginning 3 DPI as we previously reported (Arnold et al., 2016b, Haberthur et al., 2014).

Our analysis also revealed reduced expression of several genes that play a role apoptosis at each timepoint. This may be a viral evasion strategy as both SVV (via ORF63 (Hood et al., 2006)) and VZV (via ORF12 (Liu et al., 2012)) have been shown to inhibit apoptosis in infected cells. In addition, genes encoding cell adhesion molecules were downregulated at each timepoint, which may lead to increased T cell motility normally associated with their activation (Andersson et al., 1994). In support of this possibility, decreased expression of cell adhesion molecules correlates with upregulated T cell activation genes in our data set. A down-regulation of cell adhesion molecules has been described for Epstein-Barr virus (Krishna et al., 2005) and human cytomegalovirus (HCMV) (Leis et al., 2004) infection. Interestingly, immune genes that mapped to antigen presentation were downregulated 3 and 10 DPI, but upregulated 7 DPI. This change in FC direction could be due to initial immune evasion strategies employed by SVV (Whitmer et al., 2015, Meysman et al., 2013) 3 DPI, followed by development of host defense 7 DPI, and again downregulation of these pathways as the immune response begins to resolve 10 DPI.

Conclusion

In summary, this study provides the first in vivo analysis of SVV-infection induced changes in T cells gene expression profile. Our data provides novel insight into how SVV may use the host cell machinery to support its replication while evading host defense responses. Overall, we found that SVV infection of T cells leads to increased expression of genes important for chromatin modification, TCR signaling, and decreases expression of genes that play a role in antigen presentation and cellular metabolism. Changes in expression of genes important for regulation of gene expression may be critical for supporting viral replication. Increased expression of T cell activation genes may enhance the ability of SVV infected T cells to traffic to the skin and ganglia. Since T cells play a critical role in VZV dissemination, our RNA-Seq analysis may guide future studies in the development of more efficacious antivirals. The major limitation of our study is that we were not able to isolate SVV-infected T cells from bystander (and potentially antigen-specific) T cells that precluded us from determining whether gene expression changes originated from the infected T cells or from T cells responding to SVV infection. This limitation will be addressed in future studies by utilizing single cell RNASeq analysis.

Funding

This work was supported by the National Institutes of Health (NIH)(RO1AG037042-06).