Host Cell Transcriptome Profile during Wild-Type and Attenuated Dengue Virus Infection

Dengue viruses 1–4 (DENV1-4) rely heavily on the host cell machinery to complete their life cycle, while at the same time evade the host response that could restrict their replication efficiency. These requirements may account for much of the broad gene-level changes to the host transcriptome upon DENV infection. However, host gene function is also regulated through transcriptional start site (TSS) selection and post-transcriptional modification to the RNA that give rise to multiple gene isoforms. The roles these processes play in the host response to dengue infection have not been explored. In the present study, we utilized RNA sequencing (RNAseq) to identify novel transcript variations in response to infection with both a pathogenic strain of DENV1 and its attenuated derivative. RNAseq provides the information necessary to distinguish the various isoforms produced from a single gene and their splice variants. Our data indicate that there is an extensive amount of previously uncharacterized TSS and post-transcriptional modifications to host RNA over a wide range of pathways and host functions in response to DENV infection. Many of the differentially expressed genes identified in this study have previously been shown to be required for flavivirus propagation and/or interact with DENV gene products. We also show here that the human transcriptome response to an infection by wild-type DENV or its attenuated derivative differs significantly. This differential response to wild-type and attenuated DENV infection suggests that alternative processing events may be part of a previously uncharacterized innate immune response to viral infection that is in large part evaded by wild-type DENV.


Introduction
Dengue viruses 1-4 (DENV1-4) are the world's most prevalent arthropod-borne viruses [1]. DENVs are responsible for an estimated 50-100 million cases of debilitating or life-threatening infection every year and an estimated 2.5 billion people in over 100 endemic countries are at risk of infection [1,2]. The economic impact of DENVs has been estimated to be as high, if not higher than other major global health menaces such as malaria, tuberculosis, hepatitis, bacterial meningitis and others [3][4][5][6][7][8]. Despite the considerable health and economic impact, there are as yet no licensed vaccines or antiviral drugs to combat DENVs and an incomplete understanding of the biology of DENV infection has hampered progress on both of these fronts.
Given the limited coding capacity of their ,11 kb RNA genome, DENVs must parasitize the host cell machinery to complete their life cycle. At the same time, these viruses must effectively evade or suppress the host responses that act to restrict their replication [9][10][11]. This interplay between host and virus and the effect it has on host gene expression has been described previously [12][13][14][15][16][17][18][19][20][21][22][23][24][25][26][27][28][29]. Largely uncharacterized, however, is whether the transcriptional start site (TSS) and post-transcriptional variations of host RNA, leading to the production of different gene isoforms, may play a role in DENV infection. Differential RNA processing is known to be a major factor underlying cellular and functional complexity [30,31]. In order to interrogate TSS and post-transcriptional RNA variations across the entire genome in response to DENV infection, we harnessed the power of RNA sequencing (RNAseq). RNAseq is a recently developed approach to transcriptome profiling that permits a precise quantification of RNA levels and their alternatively processed variants by means of high throughput, massively parallel sequencing and subsequent mapping of the resultant short sequence fragments onto a reference genome [32,33].
We utilized two strains of DENV1 in our RNAseq study to identify strain-specific TSS and post-transcriptional variations in response to infection. The first strain, DENV1-16007, was isolated from the serum of a patient in Thailand in 1964. The second strain of DENV1 used in this study is an attenuated derivative of DENV1-16007. This attenuated virus, DENV1-PDK13 was passaged 13 times in primary dog kidney cells and was shown to be immunogenic but minimally reactogenic in human volunteers [34,35]. Our RNAseq data indicate that significant differences exist between these two strains of DENV1, not only at the transcript level but also at the level of alternative splicing. Similar trends were observed in RNAseq of an additional two lowpassaged DENV1 clinical isolates. These findings suggest that subversion of the host response includes TSS and posttranscriptional modification and is part of the mechanism of virulence. These findings also suggest that variations in the viral genome can have a profound effect in modifying host response to infection.

Cells and virus stock
HuH7, C6/36 and BHK-21 cells were purchased from the American Type Culture Collection (ATCC) and cultured according to ATCC recommendation. DENV1 strains 16007 and PDK13 were obtained from the Division of Vector-borne Diseases, Centers for Disease Control and Prevention. Sequence analysis in our laboratory indicates that these strains match the published sequences for these viruses (GenBank accession numbers AF180817.1 and AF180818.1, respectively). These viruses were amplified three times in C6/36 cells prior to use in the current study. Additionally, the DENV1 clinical isolates EDEN3300 and SL107, which were obtained from previously reported studies [36,37] and passaged in C6/36 cells for ,5 times were also included. The supernatant of all virus cultures were harvested six days post infection, clarified by centrifugation at 4506 g for 10 min at 4uC, filtered and concentrated by centrifugation at 30,0006 g for 3 hrs at 4uC. Virus pellets were re-suspended in DMEM medium with 2% FBS (Invitrogen) and stored at 280uC until use. Infectious titer was determined by plaque assay as described previously [38].

Virus infection in HuH7 cells
HuH7 cells were seeded at 3610 6 per flask in 25 mm flasks and incubated at 37uC for 24 hrs before infecting at a MOI of 20 with each of the DENV strains for 1 hr at 37uC/5% CO 2 , with gentle rocking every 15 min. The cells were then washed thoroughly and replaced with DMEM medium supplemented with 2% FBS and incubated for 20 hrs.

Immunofluorescence assay
Immunofluorescence assay was conducted according to a previously described method [39]. Briefly, the cells from the virus culture were washed once and re-suspended with PBS, and spotted onto a Teflon coated glass slide, air dried and then immersed in 80% acetone for 10 min. The slide was rinsed with PBS and airdried. 2 ml antibody against prM protein (2H2 monoclonal antibody) was added onto each well, incubated at 37uC for 45 min in a humidified chamber, and washed twice with PBS before drying. FITC-conjugated goat anti-mouse IgG were diluted 1:30 with 0.1% Evan's Blue and 2 ml was added onto each well. Slides were then incubated at 37uC for 45 min in the humidified chamber and then washed twice with PBS. Slides were dried and mounted with buffered glycerol before imaging under a fluorescent microscope.

Generation of whole-transcriptome cDNA library
Polyadenylated mRNA was isolated from HuH7 cells by three rounds of selection with the Dynabeads mRNA Direct Kit (Invitrogen) and assessed by electrophoresis on the Bioanalyzer 2100 (Agilent) for quality evaluation. For the RNAseq sample preparation, the NEBNext mRNA Sample Prep Master Mix Set 1 was used according to the manufacturer's protocol (NEB). Briefly, 0.5 ug mRNA was used for fragmentation and then subjected to cDNA synthesis using SuperScript III Reverse Transcriptase (Invitrogen) and random primers. The cDNA was further converted into double stranded cDNA and after an end repair process (Klenow fragment, T4 polynucleotide kinase and T4 polymerase), was ligated to Illumina paired end (PE) adaptors. Size selection was performed using a 2% agarose gel, generating cDNA libraries ranging in size from 275-325 bp. Finally, the libraries were enriched using 15 cycles of PCR and purified by the QIAquick PCR purification kit (Qiagen).

Author Summary
Dengue is the most common insect-borne viral disease globally. The continued absence of an effective therapy stems from an incomplete understanding of disease pathogenesis, of which the host response to infection is thought to play a central role. While previous studies have described the changes in total gene expression with dengue virus infection, they have not been able to provide any information on the subtle variations of the host RNA. These variations lead to the production of gene isoforms that can have a profound effect on gene function. In the current study, we have used the newly developed technique of RNA sequencing to more accurately interrogate the variations in the host RNA after infection with a wild-type dengue virus or its attenuated derivative. Findings from this study show that there is an extensive amount of previously uncharacterized variation in host RNA response to dengue infection. The response to infection with the wild-type dengue also differs significantly from infection with the vaccine strain. This suggests that variations in the host RNA comprise a part of the host response to viral infection that is in large part evaded by wild-type dengue viruses.
Validation of exon skipping by real-time PCR Total RNA derived from mock-infected and DENV-infected cells was used to synthesize cDNA using SuperScript III First Strand Synthesis System (Invitrogen) with random hexamers according to manufacturer's instructions. Quantitative real-time PCR was performed using LightCycler 480 Real-Time PCR System (Roche Diagnostics GmbH, Germany) and LightCycler 480 SYBR Green I Master (Roche Diagnostics GmbH, Germany). The reaction was carried out to simultaneously amplify exonskipped and exon-included isoforms using specific primers complementary to the exons flanking each target exon (Table  S13). Percent exon exclusion levels were calculated as the percentage of the isoform excluding an alternative exon divided by the total abundance of the isoforms including and excluding the alternative exon.

Statistical analysis
The statistical analysis employed to analyze biological triplicates for each condition are as previously described using the software Cufflinks and Mixture-of-Isoforms (MISO) [40,41]. Criteria used to define significance are according to the standard options in the Cufflinks and MISO programs. More detail is provided below.

Experimental design and RNAseq results
To investigate the effect of DENV1-16007 (wild-type strain) and DENV1-PDK13 (attenuated strain) infection on the transcriptome of the human host, we infected human hepatoma cells (HuH7) with each strain for 20 hours at a multiplicity of infection (MOI) of 20 ( Figure 1A). The MOI of 20 was chosen so that the subsequent RNAseq profiling would best reflect the infection-induced alterations to the host transcriptome and not be either masked by or derived from a large number of uninfected cells. We also measured viral RNA over the first 30 hours of infection to address the possibility that changes observed might be the result of a delayed replication cycle by one of the viruses. Although the absolute kinetics of the two viruses differ, the 20 hour time point was chosen as it represents the stage in the primary round of replication at which the genome copy numbers are most similar between the two viruses ( Figure 1B). Indirect immunofluorescence staining for the pre-membrane protein (prM) production with 2H2 monoclonal antibody also showed similar infection levels for both viruses at this time point ( Figure 1C). Infections with wild-type and attenuated strains were performed in three independent biological replicate experiments. Mock-infected HuH7 cells treated in the same way as the infected samples were also done in biological triplicate and served as the control for our experiment. At twenty hours post-infection, mRNA was extracted from all samples independently, and poly-A enriched cDNA libraries were constructed for 75-base, pair-end sequencing on an Illumina GAIIx (one sample for each condition) or Illumina HiSeq2000 machine (two samples for each condition). RNA sequencing was performed independently for each of these replicate experiments.
Mapping of reads (Bowtie, Tophat) and analysis of differential transcriptome response (Cufflinks) to infection was performed using the Tuxedo Suite of software [40,42]. Cufflinks utilizes sequence fragments mapped to the reference genome to estimate the abundance of each isoform arising from the gene. It then tests for differential expression between experimental conditions. Differential expression at the gene level is calculated and is defined as the sum of differential expression of all isoforms at a particular locus. Cufflinks also assesses differential splicing by comparing the relative abundance of isoforms using the same transcriptional start site [40]. As this definition of splicing encompasses all different types of splicing events (see below) [43], systematic analysis of these events for downstream validation work is exceedingly difficult. In order to attain more detail about the types of differential splicing in our samples, we used the MISO software following mapping [41]. MISO utilizes a fixed library of previously characterized splice events to predict alternative splicing and reports the number of occurrences for each type of splicing event: skipped exon (SE), mutually exclusive exons (MXE), alternative 39 splice site (A3SS), alternative 59 splice site (A5SS), alternative first exon (AFE), retained intron (RI), tandem untranslated regions (Tandem UTR) [43]. MISO also utilizes a different algorithm than Cufflinks to predict differential splicing between samples which, when compared with the results from Cufflinks, provides an additional level of stringency in selection of candidates for downstream analysis.
To assess the quality of our libraries and sequencing performance and to ensure that any differences observed between our samples was due to the biology and not bias in sequencing, we used the open source program RNAseqC to examine our data [44]. Results indicate that despite using two different Illumina machines to generate the sequences, the individual samples are highly comparable to each other across all the metrics interrogated (Dataset S1).
Wild-type and attenuated DENV-1 strains induce differential regulation of the host transcriptome Over 18,000 changes to the host transcriptome were observed in response to infection by the wild-type strain and .41,000 were observed in response to infection with the attenuated strain (Table 1). Differential isoform regulation is the largest category of response due to infection by both strains. Interestingly, there are over two-fold more differentially regulated isoforms following infection with the attenuated than with the wild-type strain. Similarly, infection with the attenuated strain also resulted in three-fold more differentially regulated genes than infection with the wild-type ( Figure 2). In order to determine whether this large differential response to infection was specific to these two strains or whether the 'quieter' response to the parental DENV1-16007 strain was typical of wild-type DENV1s, we repeated our experiment in HuH7 cells with two low passaged clinical isolates (EDEN3300 and SL107) and compared them to our uninfected control. RNAseq analysis for these isolates indicates even fewer transcriptomic changes ( Table 1), suggesting that attenuated virus triggers more host cell response than wild-type viruses. This observation is consistent with what has been reported for yellow fever virus and its attenuated derivative, YF17D [45].
Qualitative analysis of the RNAseq data also provides insights into the pathways that are essential to both strains or unique to one strain. Ingenuity Pathway Analysis (IPA) of isoforms indicates that the commonly regulated isoforms are enriched in pathways associated with viral infection and modulation of protein translation. Examples include EIF2 signaling, prolactin signaling, acute phase response signaling, regulation of EIF4 and p7056K and glucocorticoid receptor signaling. Differential expression of pathways associated with cellular growth and proliferation such as mTOR signaling and aryl hydrocarbon receptor signaling are also enriched following infection with both strains ( Figure 3A). Conversely, the pathways that are differentially regulated between the wild-type and attenuated strains include the innate immune response and cell cycle control. The top differentially regulated host pathways associated only with wild-type strain infection are involved with immunomodulation and cell cycle arrest, such as PPAR/RXR activation pathway, G2/M DNA damage checkpoint, ATM signaling and the PDGF signaling pathway ( Figure 3C).
Conversely, infection with the attenuated strain triggered pathways associated with inflammation, induction of apoptosis and stress such as the TNFR1 pathway, TWEAK signaling, NRF2-mediated oxidative stress response, IGF -1 signaling and ERK5 signaling ( Figure 3E). IPA also identified molecular and cellular functions associated with both or each of the strains ( Figures 3B, 3D and 3F). Taken collectively, the molecular and cellular functions in response to wild-type virus infection are associated with cell signaling and metabolism while those to attenuated virus are associated with transcriptional activation, cell cycle modification and post-translational modification.
Differential splicing in the host transcriptome is largely conserved between wild-type and attenuated DENV1 strains Next, we interrogated the data for alternative splicing within the isoforms sharing the same transcriptional start site. Interestingly, 80% and 74% of all differential splicing following infection with the wild-type strain and attenuated strain, respectively, are common to both viruses ( Figure 2). This degree of commonality in the splicing response to infection by each strain is significantly different than what was observed for differential isoform and gene regulation in response to infection. This novel finding suggests that mechanisms responsible for the specific regulation of host splicing may be critical for DENV1 propagation and thus remained relatively unchanged during the attenuation process.
The candidate list from our Cufflinks analysis of alternatively spliced transcripts for both strains of virus was then crossreferenced against the list of skipped exon (SE) events generated by MISO analysis. This cross comparison resulted in 79 total SE events. To assess the accuracy of our alternative splicing predictions, we performed qPCR on each of the predicted SE events. Of the 79 events tested by qPCR, 32 (40%) were differentially expressed in the wild-type strain, the attenuated strain or both in comparison to an uninfected control (Figure 4). The splicing patterns were similar for both viruses across all time points measured, indicating that our observations are not artifacts of temporal sampling bias ( Figure 4). Furthermore, the genes involved in these splicing events belong to many of the same pathways shown in Figure 3 suggesting that the virus is, at least in part, exerting its influence on these pathways through alternative splicing (Table S1). If this were true, the SE events should share mechanisms of regulation by utilizing common RNA motifs within and surrounding the identified SE's. Indeed, using the software RegRNA (http://regrna.mbc.nctu.edu.tw) [46] and setting an arbitrary boundary of 250 nucleotides upstream and downstream of the 39 and 59 splice-site of the SEs, respectively, we identified 74 predicted RNA regulatory motifs and elements which could be bound by 17 RNA binding proteins. By using GeneCards (http:// www.genecards.org) to convert these putative motif binding proteins into HUGO nomenclature, we observed that nearly two thirds (11 of 17) of genes encoding these motif-binding proteins are themselves differentially expressed following infection  Tandem UTR 0 0 Total 609 384 by one or both of the DENV1 strains (Table S2). These results suggest that DENV1 may be regulating alternative splicing of host mRNA by a hitherto unknown mechanism.
Briefly, we examined 11 canonical pathways that were enriched in our study and/or had been previously implicated in the DENV life cycle: apoptosis, autophagy, clathrin-mediated endocytosis, interferon signaling, lipid metabolism, oxidative phosphorylation, regulation of stress granules and P-bodies, splicing-related RNA post-transcriptional modification, ubiquitination, endoplasmic reticulum stress, virus recognition and interferon induction. The proportion of host factors for DENV or other flaviviruses identified through functional genomic studies that were differentially regulated following infection with the wild-type strain or the attenuated strain is indicated in Table 2. Differential isoform regulation of flaviviral host factors ranged from 15.5% to 60.6% and is uniformly higher for the attenuated strain than the pathogenic wild-type strain. Alternative splicing ranged from 7.3% to 24.2% and the rates of these events in the wild-type strain and its attenuated derivative were comparable ( Table 2). The intersection of specific host factors in these eleven canonical pathways that are either required for flaviviral propagation or involved in direct interaction with DENV proteins, which are differentially regulated are graphically depicted in Figure 5.

Discussion
Investigation of host gene expression to date has relied primarily on microarray technology [12][13][14][15][16][17][18][19][20][21][22][23][24][25][26][27][28][29]. This technology is insensitive to the regulation of genes through alternative RNA processing. Thus, detection of differential expression down to the level of alternative isoforms has not been examined and the subtleties of TSS and post-transcriptional modifications that can dramatically alter the function of the derived proteins have been largely ignored. These processes could be a mechanism by which DENV attains specific isoforms of required host factors while suppressing those that act to restrict its replication [27,52,53]. Indeed, the need to understand the host response beyond simple gene expression is underscored by the observation that DENV can modify the splicing pattern of an endogenous gene, XBP1, to its advantage [53]. Our findings indicate that there is an extensive amount of previously uncharacterized gene isoforms and alternative processing of host transcripts over a wide range of pathways and host functions in response to DENV infection. Interestingly, the DENV1-16007 and DENV1-PDK13 viruses only differ from each other by 14 nucleotide and 8 amino acids [54], yet the host transcriptional response to these viruses is pronounced. This suggests that infection with different strains of DENV can result in significantly different disease phenotypes despite few nucleotide differences.
We have attempted, in this study, to provide a comprehensive guide to the transcriptomic changes with DENV infection. By analyzing our RNAseq data using two different programs, Cufflinks and MISO [40,41], maximal information on RNA transcript regulation, could be gleaned. Specific splice events could hence be identified for subsequent mechanistic studies that clarify their role in the host response. In particular, the integrative analysis of this study with existing functional genomics data reveals previously undocumented expression and post-transcriptional regulation of required host factors that should serve as a road map for future mechanistic investigations. A caveat, however, is that our work only profiled the transcriptome at the end of a single round of DENV replication. As hinted by Figure 4, both quantitative and qualitative differences may exist in the transcriptome at different stages of the virus life cycle. Furthermore, the possibility exists that bystander uninfected cells may exert some influence on the observed transcriptional changes although we have attempted to minimize this by using a high MOI in our experiments. Future studies may need to take these possibilities into account.
The large difference in the number of alternative splicing events identified by Cufflinks and MISO also underscores the fledgling nature of RNAseq. While the former identifies all possible splicing events from the data de novo, the latter relies on a pre-defined library to map alternatively spliced transcripts. Additional studies are needed to clarify whether Cufflinks over-estimated the number of splice variants, or there are authentic variants absent from the library used by MISO. Regardless, for transcriptome analysis of host response to infection, RNAseq is superior to microarray in terms of the breadth of information derived.
Our results also indicate that the human transcriptome response to an infection by wild-type DENV or its attenuated derivative differs significantly (Table 1, Figure 2). These differences suggest that alternative processing events may be part of a previously uncharacterized innate immune response to DENV1 infection that is in large part evaded by wild-type strains. This second hypothesis is supported by the greater than two-fold increase in the number of differentially regulated transcripts when infected with the attenuated strain of DENV1 as compared to the parental wild type strain of DENV1, many of which belong to pathways associated with inflammation, induction of apoptosis and stress response (Table 1, Figure 3). This inability to escape the innate immune response achieved by the wild-type virus may explain the lack of  The reactogenicity. It may also explain why antibody titer engendered by vaccination is not as high as those observed following natural infection unless supplemented with an adjuvant that stimulates this innate immune response. This observation also suggests a mechanism of pathogenicity where DENV regulates host transcriptome changes by interacting with a group of RNA binding proteins to control multiple splicing events. The development of a live attenuated tetravalent vaccine for DENV1-4 has bedeviled researchers for the past 60 years. The less than optimal efficacy of the leading dengue vaccine candidate makes an improved understanding of the molecular basis of a good vaccine all the more critical [55]. The differential host transcriptome response to infection with DENV1-16007 and DENV1-PDK13 provides an insight into the characteristics of an attenuated virus which, is likely a complex phenotype [45]. A molecular understanding of the basis of attenuation could lead to a quantitative approach to balancing reactogenicity and immunogenicity, which presently remains a hit-or-miss finding made only after lengthy clinical trials.
In conclusion, we provide here a detailed view of the host cell transcriptome response to infection with wild-type DENV-1 and its attenuated derivative that could be useful for future studies on the genetic determinants of viral virulence and attenuation.

Supporting Information
Dataset S1 RNASeQC analysis. Table S1 Ingenuity pathway analysis of the 33 skipped exon events found to be differentially regulated by qPCR. 2log(p-value) reflects the statistical significance of the indicated observation while the ratio refers to the number of genes identified in the indicated pathway divided by the total number of genes belonging to that pathway.