Interaction sites of the Epstein-Barr virus Zta transcription factor with the host genome in epithelial cells.

Epstein-Barr virus (EBV) is present in a state of latency in infected memory B-cells and EBV-associated lymphoid and epithelial cancers. Cell stimulation or differentiation of infected B-cells and epithelial cells induces reactivation to the lytic replication cycle. In each cell type, the EBV transcription and replication factor Zta (BZLF1, EB1) plays a role in mediating the lytic cycle of EBV. Zta is a transcription factor that interacts directly with Zta response elements (ZREs) within viral and cellular genomes. Here we undertake chromatin-precipitation coupled to DNA-sequencing (ChIP-Seq) of Zta-associated DNA from cancer-derived epithelial cells. The analysis identified over 14 000 Zta-binding sites in the cellular genome. We assessed the impact of lytic cycle reactivation on changes in gene expression for a panel of Zta-associated cellular genes. Finally, we compared the Zta-binding sites identified in this study with those previously identified in B-cells and reveal substantial conservation in genes associated with Zta-binding sites.

Introduction

Infection with Epstein–Barr virus (EBV) results in a life-long association with the virus. EBV infects B-lymphocytes, NK/T lymphocytes and some epithelial cells. EBV infection of B-lymphocytes results in a pre-latent phase of viral gene expression and the activation of many host genes [1-4]. Within 4 days infected cells start to proliferate and, potentially through epigenetic reprogramming, a sub-set of latency genes are expressed, and long-term latency is established [2, 5, 6]. The site for EBV latency in healthy people is the non-proliferating memory B-lymphocyte pool [7, 8]. Here, the viral genome resides in the cell’s nucleus, but there is little viral gene expression and no viral replication. EBV is associated with both lymphomas and carcinomas: Burkitt’s lymphoma [9, 10]; Hodgkin’s disease [11]; NK/T lymphomas [12]; post-transplant lymphomas [13]; diffuse large B-cell lymphomas [14]; nasopharyngeal carcinoma (NPC) [15, 16] and gastric carcinomas [17]. The viral genome is latent in these cells, with between 1–10 % of the viral genes expressed [18] and no viral replication.

The viral genome consists of 170 kb of double-stranded DNA. Epigenetic modifications to the EBV genome occur during viral latency, resulting in establishing a repressive chromatin environment [19] that restricts viral gene expression. The epigenetic changes are associated with a combination of DNA methylation of CpG motifs in non-transcribed regions of the genome [19, 20] and tri-methylation of histone H3 at lysine 27 and lysine 9 [21, 22] at viral gene promoters. The virus can replicate in B-lymphocytes and epithelial cells as they differentiate. Reactivation of EBV from latency results in the activation of an ordered cascade of the expression of around 90 viral genes, coordinated with lytic replication of the viral genome, assembly, packaging of virions and viral egress [23-25]. Histone deacetylase inhibitors, DNA methylation inhibitors and activators of signal transduction can induce viral reactivation. An important step in viral reactivation is the expression of the viral gene BZLF1, encoding the Zta transcription factor [26]. Zta then plays a direct role in driving the expression of other EBV genes, as shown by the ability of Zta to activate viral genes following ectopic expression in cells harbouring latent virus [25]. The viral gene BRLF1 encodes the transcription factor Rta, which can initiate the lytic replication cycle in epithelial cells and co-operates with Zta to activate many viral genes [25, 27]. Zta binds directly to a DNA sequence motif that resembles the cellular AP1 motif [28-31]. Zta possesses the unusual ability of preferentially binding to related DNA-binding sites that contain an internal CpG motif when they are methylated [32-37]. In addition, Zta binds to and activates the expression of some cellular genes [35, 38–42]. Binding sites for Zta on the cellular cell genome in several B-cell lines have been identified [43, 44]. It was recently revealed that Zta binding to the human genome results in global changes to the architecture of cellular chromatin, opening some areas and making other areas inaccessible [44]. The high number of Zta-binding sites has been postulated to contribute a buffer zone for Zta function; to prevent low levels of Zta initiating virus lytic cycle erroneously [44]. During the viral lytic cycle in B-cells, the expression of some cellular genes located near Zta-binding sites is activated while others are downregulated [43, 44].

Here we focussed on another cell type infected by EBV, epithelial cells. We asked which binding sites Zta recognizes on the cellular genome during the reactivation of EBV latency in epithelial cells and whether these relate to the sites that we previously identified in Burkitt’s lymphoma cells undergoing lytic reactivation [43]. Integration of the data reveals the conservation of a common sub-set of cellular genes associated with Zta binding in both cell types.

Methods

Cell culture and induction of EBV lytic replication

HONE1-EBV cells were received from G. Tsao (University of Hong Kong). They were derived from poorly differentiated squamous cell carcinoma of the nasopharynx and re-infected with recombinant Akata EBV genomes [45]. HONE1-EBV was since shown to contain part of the genome of the human epithelial cell line HeLa [46]. HONE1-EBV cells were maintained in RPMI medium supplemented with 10 % (vol/vol) foetal bovine serum, 100U of penicillin ml−1, 100 g of streptomycin ml−1 2 mM l-glutamine (Invitrogen) at 37 °C with 5 % CO2. Transfection of HONE1 cells with EBV is not stable, but the EBV plasmid contains a G418 resistance gene, and EBV-positive HONE1 cells are continually selected with 600 µg ml−1 G418. For EBV lytic induction, HONE1-EBV cells were induced with 10 µM suberoylanilide hydroxamic acid (SAHA) for 48 h at 70 % confluence. Only the adherent cell population was used in experiments.

Antibodies and Western blot

Cells were harvested and boiled in a protein sample buffer (Sigma). The protein was resolved on a NuPAGE 12 % Bis-Tris polyacrylamide gel (Life Technologies) in 1× morpholinepropanesulfonic acid buffer (Invitrogen) at 200 V for 50 min. The proteins were transferred onto nitrocellulose membranes (Santa Cruz Biotechnology) and incubated overnight with the indicated antibodies at 4 °C. The following day, the membranes were incubated with IRDye secondary antibodies (LiCor), and the proteins were detected by fluorescence at 680 and 800 nm. In addition, the BZ1 mouse monoclonal antibody to Zta [47] and a rabbit polyclonal recognizing beta-Actin (Sigma) were used to detect proteins by Western blotting. Intracellular analysis of Zta expression was detected using BZ1 mouse monoclonal antibody to Zta [47] and a Fix and Perm Cell fixation and permeabilization system (Life Technologies) and subsequently detected using a FACS Canto analyzer (Beckton Dickinson).

Chromatin precipitation

Chromatin was fixed by applying 1 % (v/v) formaldehyde to cells for 10 min at 20 °C and was extracted as described in Bark-Jones et al. [48]. A Vibra-Cell (SONICS) was used to shear the chromatin. HONE1-EBV cell lysates were subject to sonication in a volume of 100 µl on ice for ten rounds, with 10 s pulses and a 21 % amplitude output. The efficacy of sonication was assessed by analysing the sample on a 1 % (w/v) agarose gel. Chromatin from 5×106 cells was used to immunoprecipitate Zta associated-DNA using 10 ug EBV ZEBRA antibody (sc-53904) (Santa Cruz). Immune complexes were collected with 50 % pre-blocked protein A/G-sepharose beads. Following wash steps and protein digestion of immune complexes, eluted DNA was purified using a PCR purification kit (Qiagen) and eluted in 55 µl. For the chip-sequencing, the chromatin-precipitation was scaled up and undertaken twice (ChIP1 and ChIP2), with chromatin from 3.5×107 cells using EBV ZEBRA antibody (sc-53904) (Santa Cruz) and the DNA was eluted into 45 µl.

PCR analysis

RNA was extracted from cells exposed or not to 10 uM SAHA for 48 h with RNeasy Plus (Qiagen). In addition, first-strand cDNA synthesis (Roche) was undertaken with random primers, and relative changes in gene expression were determined using q-PCR Taqman assays as described in [43] using assays for SCL6A7 (Hs00204454_m1), ANO1 (Hs00216121_m1), FOSB (Hs00171851_m1), SCIMP (Hs010294_m1;), RASA3 (Hs00183698_m1) and FSCN1 (Hs00602051_mH). For chromatin, loci were analysed with GoTaq qPCR master mix (Promega) with primers for SCIMP described in [37, 43].

The EBV genome load was determined using QPCR analyses undertaken using Sensimix SYBR (Bioline) on an ABI 7500 real-time QPCR system (Applied Biosystems). The absolute quantitation method was used with dilutions of input DNA to generate the standard curve. The signal resulting from amplification with primers specific for the EBV genome (DNA polymerase gene) AGTCCTTCTTGGCTAGTCTGTTGAC and CTTTGGCGCGGATCCTC were compared with signal resulting from amplification with the human genome GGCAACCCTAAGGTGAAGGC and GGTGAGCCAGGCCATCACTA as described previously [21, 49].

ChIP DNA was analysed using primer sets that include the binding peak and primer sets that flank the peak. The DNA sequences for the primers that amplify the EBV sequences from Orilyt and its flank are Orilyt CAGCTGACCGATGCTCGCCA and ATGGTGAGGCAGGCAAGGCG and Flank Orilyt GCGCAACAGTGCCACCAACC and CAGGACCTGGCGGTAGTGCAG [21].

ChIP library preparation

The DNA libraries from ChIP1 and ChIP2 were prepared according to the NEBNext ChIP-Seq Library Prep Reagent Set for Illumina kit (New England Bioscience), switching the PCR and the size selection steps. Index four was used for the input, index six for ChIP1 and index 12 for ChIP2. Size selection was performed on a 2 % (w/v) agarose gel with SYBR Gold dye to exclude the adaptor from the sequencing process and purify DNA fragments in the 175–225 bp size range. Quality control was performed using the Agilent 2100 Bioanalyzer.

ChIP-sequencing and data analysis

DNA sequencing was undertaken on an Illumina HiSeq sequencer (Eurofins). Base calls were performed using RTA 1.18.61. The sequencing library size ranged between 3760 and 4393 Mbp. Mapping of reads to the reference sequences was performed using BWA-BACKTRACK version 0.6.2-r126 [50]. The sequencing data includes the human genome and the Epstein–Barr virus genomes. The data are deposited in GEO: GSE83354.

Peaks were called using deepTools Galaxy MACS-1.4.2 with the following setting: tag size 30 bp, bandwidth 200 bp, P-value 10−9, arbitrary shift size 100 bp, format. wig, MFOLD range 10, 30, Dynamic lambda 1000–10000.

Peaks common to ChIP1 and ChIP2 were determined using deepTools Galaxy’s ‘Operate on Genomic Intervals’ tool (overlapping intervals for at least 1 bp). Overlapping pieces of intervals were selected as return format. The binding of Zta to the human genome was visualized in the UCSC genome browser (http://genome.ucsc.edu/) [51]. For this purpose, the .bam files were normalized to fragments per kilobase per million reads and converted to .bigwig format with the help of the bamCoverage tool in deepTools Galaxy (http://deeptools.ie-freiburg.mpg.de/) [52] using the following parameters: 200 bp fragment size and 30 bp bin size .bigwig files were then converted to .bed files by UCSC for downstream analysis. MACS peaks were determined using the MACS-1.4.2 tool in deepTools Galaxy with a cut off P-value of 1×10−9. Overlaps in binding were identified by intersecting MACS peaks within deepTools Galaxy [52].

The nearest transcription start sites, gene lists and gene ontology associated with neighbouring genes associated with Zta-associated DNA peaks were determined using Great 4.0.4 [53]. Finally, motif analysis was performed with the help of the MEME-ChIP suite using default parameters (http://meme.nbcr.net/meme/tools/meme-chip) [54].

Results and discussion

Zta association with cellular DNA in epithelial cells

We undertook chromatin precipitation coupled with DNA sequencing (ChIP-seq) to identify loci that the viral transcription factor Zta interacts within the cellular genome of EBV-infected epithelial cells (HONE1-EBV [55]). Cells were exposed to the histone deacetylase inhibitor suberoylanilide hydroxylamine (SAHA, Vorinostat) to promote Zta expression and reactivation of the EBV lytic replication cycle [45] (Fig. 1a). SAHA addition resulted in intracellular detection of Zta in 20 % of cells (data not shown) in good accordance with the initial characterization of the cells and induction method [45]. We also identified a substantial increase in EBV genome abundance (Fig. 1b), confirming the onset of the EBV lytic replication cycle. We undertook two independent Zta chromatin precipitation experiments, and ChIP-seq libraries were prepared and evaluated for Zta binding. For each chromatin precipitation, we found Zta bound to a previously characterized ZRE in the Orilyt region of the EBV genome, but not to a flanking site 3 Kb upstream (P≤0.001) (Fig. 1c); this validated the Zta-binding specificity of the library. Sequence analysis and alignment to the human genome resulted in identifying 21 694 and 16 766 Zta-binding sites for each of the chromatin precipitations (P≤1×10−9). Of these, 14 666 locations showed specific interaction with Zta in both datasets (Tables S1–S3, available in the online version of this article) (Fig. 1d); we analysed this set of Zta-binding peaks in further experiments.

Fig. 1.

Zta interaction with the human genome in epithelial cells. HONE1-EBV cells were stimulated or not with SAHA for 48 h. (a) Total protein extracts were prepared and subject to Western-blot analysis with antibodies for Zta and actin. (b) DNA was prepared from the cells, and quantitative PCR was used to analyse the cellular and EBV genome amounts in duplicate. The relative EBV genome load was calculated, and average values are shown. * represents P≥0.05. (c) Two ChIP-seq experiments were undertaken with the Zta antibody; quantitative PCR was used to analyse the relative binding of Zta to a previously characterized site on the EBV genome at oriLyt and a region adjacent (in triplicate). *** represents P≤0.001; **** represents P≤0.0001. (d) Following alignment of the ChIP-Seq reads to the human genome, binding peaks were identified. The number of peaks identified with a P-value ≤10–9 is shown for each experiment. In addition, the DNA sequences associated with each peak were intersected, and the number of peaks present in both experiments is shown.

Characteristics of Zta-associated binding sites in cellular DNA

Interrogation of the mapping data revealed that the Zta-binding peaks lie both proximal and distal to the transcriptional start sites (TSS) of the closest genes, with the majority located more than 5 Kb from a TSS (Fig. 2a). This binding pattern resembles the overall pattern of genomic association that we identified in Burkitt’s lymphoma (BL) cells previously [43]. To further validate the ChIP-seq experiments, we investigated the interaction of Zta with a binding peak identified in one locus (SCIMP) using independent ChIP experiments (Fig. 2b). A greater association of Zta was observed at the ChIP-seq peak within the SCIMP locus (P≤0.05), compared to a region 2 kb from the observed binding peak (Fig. 2b). We then interrogated the Zta-binding peaks to identify any short DNA sequence motifs that were enriched by Zta binding. This approach determined that the most significant enriched DNA element (P=9.7×10−401) in the Zta-binding sites is an AP1 (Fos-Jun) motif. This discovery is likely to represent the direct interaction of Zta with the cellular DNA as many ZREs are known to be closely related to AP1 motifs [20, 36, 43] (Fig. 2c). We then searched for enrichment of common biological processes and discovered that the genes associated with Zta-binding sites are related to regulatory processes (Table 1).

Fig. 2.

Characterization of Zta-associated cellular genes. (a) The distance of Zta-binding sites identified from the ChIP-Seq experiment with the nearest two genes was calculated. The bar chart shows the distribution of genes at the indicated distances from the transcription start site (TSS). (b) HONE1-EBV cells were stimulated or not with SAHA for 48 h. A ChIP experiment was undertaken with the Zta antibody; quantitative PCR was used to analyse the relative binding of Zta to a peak identified from the ChIP-Seq analysis in the SCIMP locus and to an adjacent region (in triplicate). * represents P≤0.05. (c) Motif enrichment (MEME-CHIP) identifies conserved motifs within the Zta ChiP binding sites in epithelial cells (P=1.38×10−68).

Table 1.

Gene ontology (GO) enrichment analysis of the biological processes (BP) associated with the two closest genes within 100 kb of each Zta-binding site showing those terms with a false discovery rate (FDR) ≤1×10−15

GO BP term	BinomFDR Q-Val
Negative regulation of transferase activity	1.16E-38
Cellular component disassembly	3.13E-36
Negative regulation of kinase activity	2.93E-29
Negative regulation of protein kinase activity	5.02E-28
Negative regulation of sequence-specific DNA binding transcription factor activity	2.53E-24
Regulation of cell size	6.14E-19
Negative regulation of Notch signalling pathway	1.01E-18
Regeneration	3.26E-17
Vascular endothelial growth factor receptor signalling pathway	7.79E-17
Regulation of Notch signalling pathway	2.02E-15
Cellular response to reactive oxygen species	6.47E-15

Similarities and differences between Zta binding with cellular DNA in different cell lineages

We compared the Zta-binding peaks in cellular DNA identified in this study with those associated with a previously published set of Zta-binding sites to cellular DNA identified in BL cells [43]. The sequencing technologies have substantially improved since these data were generated. We selected a numerically equivalent sub-set of Zta-binding peaks from the epithelial data to directly compare. We then compared the Zta-binding peaks from the HONE1-EBV cells with the 5020 Zta-binding peaks identified in BL (Fig. 3). We found a remarkable coincidence of binding at over 40 % of sites. Although this represented fewer than half of the Zta-binding sites, we consider that the conservation of Zta-binding sites in these two very different cell lineages shows a robust programme of Zta-interaction with the cellular genome.

Fig. 3.

Common and unique Zta-binding peaks in epithelial and B-cells. An equivalent number of Zta-binding peaks from the HONE1-EBV cells were intersected with those previously identified in B-cells, and the proportions that are common or unique to one cell type are indicated.

We identified the genes within 20 kb from the common binding sites. The two cell types shared around 50 % of the genes (Table S4). We previously showed that about 11 % of genes nearby Zta-binding sites in the human genome are also transcriptionally regulated in BL cells [43]. We analysed the expression of a small panel of genes that (i) have Zta-binding peaks in both cell types and (ii) are regulated in Akata BL cells. We induced EBV lytic cycle by adding SAHA; this resulted in a significant up-regulation of RNA ranging from 1.5-fold to 6.8-fold for most of the genes evaluated (Fig. 4). Only around 20 % of the epithelial cells enter into the lytic cycle, which may relate to between 7.5–34-fold changes in abundance in the responsive cells. This subset of genes is regulated in common between two entirely different cell types using various means to induce the EBV lytic cycle (SAHA and anti-IgG). The use of different cell types and activators of the EBV lytic cycle, together with the proximity of Zta-binding sites to the regulated genes, supports the possibility that Zta binding may affect their regulation. Still, potential co-operation with other EBV lytic proteins or indirect routes of regulation may also mediate the effect.

Fig. 4.

Expression of a panel of genes associated with Zta-binding sites. HONE1-EBV cells were induced to initiate EBV lytic replication by exposure to SAHA for 48 h. RNA was prepared from cells exposed to SAHA or not and subject to analysis with Taq-man primer/probe sets for the indicated genes. The fold change in expression following SAHA addition is shown together with the standard deviation for the assay (*P≤0.05; **P≤0.01, ***P≤0.001).

Although we found more Zta-binding sites in HONE1-EBV than in Akata BL, we are cautious about suggesting these be epithelial-specific Zta-binding sites. The HONE1-EBV data set was generated at a higher depth of sequencing than the BL data set and may simply be more representative. Likewise, the recent publication by Buschle et al. did not draw any inferences from the uniqueness of the Zta-binding sites or gene regulation that they identified in the different B-cell data sets that they recently published [44].

Buschle et al. hypothesize that most Zta-binding sites on the cellular genome have a threefold role, to act as a sink to sequester low levels of Zta gene expression to prevent untimely entry into the lytic cycle in response to the following: (i) sub-optimal stimuli; (ii) during pre-latency; and (iii) to drive changes to chromatin structure to optimize the cellular environment for the conversion of the nucleus into a viral replication factory during the lytic cycle [44]. They present strong evidence documenting global changes in chromatin structure around cellular Zta-binding sites during EBV lytic cycle [44].

We propose that there is scope for the interaction of Zta with the cellular genome to have a different role in the cell; it could impact the transcription of specific cellular genes while also acting as a sink for Zta and a nucleus for chromatin change. Indeed, evidence of changes in the expression of some cellular genes located near Zta-binding sites is well established [35, 38–42]. This study identified a common set of 760 cellular genes (Table S4) that Zta binds near during the EBV-lytic cycle in both B-lymphocytes and epithelial cells. We propose that this dataset provides a core panel of cellular genes to investigate further to decipher the potential contribution of cellular genes to the two stages of the EBV life cycle where Zta is known to play a crucial role – the lytic cycle [27] and pre-latency [36].

Click here for additional data file.