An atlas of the Epstein-Barr virus transcriptome and epigenome reveals host-virus regulatory interactions

Abstract

Epstein-Barr virus (EBV), which is associated with multiple human tumors, persists as a minichromosome in the nucleus of B lymphocytes and induces malignancies through incompletely understood mechanisms. Here, we present a large-scale functional genomic analysis of EBV. Our experimentally generated nucleosome positioning maps and viral protein binding data were integrated with over 700 publicly available high-throughput sequencing data sets for human lymphoblastoid cell lines mapped to the EBV genome. We found that viral lytic genes are coexpressed with cellular cancer-associated pathways, suggesting that the lytic cycle may play an unexpected role in virus-mediated oncogenesis. Host regulators of viral oncogene expression and chromosome structure were identified and validated, revealing a role for the B cell-specific protein Pax5 in viral gene regulation and the cohesin complex in regulating higher order chromatin structure. Our findings provide a deeper understanding of latent viral persistence in oncogenesis and establish a valuable viral genomics resource for future exploration.

Figure 1. Genomic profiling of Epstein-Barr virus in lymphoblastoid cell lines

(A) EBV infected cell lines were subjected to functional genomics assays coupled with high throughput sequencing, giving rise to reads that map to the human host and viral genomes. (B) Known functional elements of the 172kbp viral episome examined in this study. While we provide all coordinates with respect to the 172kbp reference, most LCLs were generated with the B95.8 strain, which contains a 12kbp deletion (139–151kbp). Unmappable loci have highly similar regions in the human genome, corresponding to reads that are of ambiguous origin. (C) Short reads from hundreds of ChIP, RNA, and other (DNase, FAIRE, MNase, HiC, etc) sequencing experiments were mapped to the viral and host genomes. For each experiment, percentage of aligned viral reads was determined as a fraction of aligned human reads. An average of more than 1% of reads map to the viral genome, though the percentage of viral reads can vary across 4 orders of magnitude. Acronyms: EBNA: EBV nuclear antigen; LMP: Latent membrane protein; oriP: Latent origin of replication; oriLyt: Lytic origin of replication; TR: Terminal repeats; BHLF1: Viral lytic non-coding RNA, BamHI H leftward fragment 1; BHRF1: Viral BCL2 homologue, BamHI H rightward fragment 1; BZLF1: Viral lytic activating gene, BamHI Z leftward fragment 1; RPMS1: Viral gene and promoter for BART miRNAs; RPM: Reads per million. See also Figure S1 and Table S1.

Figure 2. Transcribed and regulatory elements of the EBV genome

(A) The average RNA expression profile (transformed by square root) shows unexpectedly high RNA reads of the lytic cycle associated transcript BHLF1, in addition to the canonical genes associated with type III latency. Also shown are the number of transcription factor binding sites and average signal of activating histone modifications at known and candidate regulatory regions of the episome. Furthermore, CTCF segments the genome into many more domains than previously recognized. (B) The lytic origin of replication is highly enriched in transcription factor occupancy. The locus acts as a divergent promoter for the lytic-associated BHLF1 and BHRF1 transcripts separated by a central CTCF binding site. (C) Histone modifications at the highly transcribed and multiply processed RPMS1 promoter. A peak for H3K4me3 was more proximal, while H3K4me1 was distal to the transcription start site. The histone variant H2A.Z can be seen at the −1 and +1 nucleosomes, and downstream histones are acetylated up to the CTCF site. (D) Nucleosome occupancy correlates with promoter silencing in different latency types. Nucleosomes have higher occupancy at Cp and Wp in type I latent cells (MutuI), where these promoters are transcriptionally silenced, relative to type III latent cells (LCLs). (E) RNA-seq reveals transcript isoforms of BZLF1 (top) and EBNA-LP joined to BHRF1 (bottom). The BZLF1 transcript can be alternatively spliced to exclude the coding region for the DNA binding domain, and multiple exons of EBNA-LP can act as donor splice sites to multiple acceptor splice sites in BHLF1. Percentages are given as total of all BZLF1 transcripts and EBNA-LP transcripts, as estimated from reads that cross alternative and canonical splice junctions. BHRF1 is also transcribed canonically from OriLyt. (F) The BHLF1 RNA transcript differs from the DNA template, suggesting RNA-editing by ADAR-mediated deamination. Position 40080 is templated to be an A but RNA-seq reveals a significant portion of reads across 170 distinct RNA-seq experiments contain a G (p < 1e-71, t-test; error bars show standard deviation). See also Figure S2.

Figure 3. Lymphoblastoid cell lines cluster by viral reactivation propensity

(A) We analyzed gene expression from 201 independent profiles of across 143 lymphoblastoid cell lines and found clustering into type III latency (upper left) or lytic (lower right) gene expression programs. The lytic cluster has a higher percentage of reads mapping to EBV than the latent cluster (bar along top of heatmap). The clustering of genes reveals association with latent and lytic cycle, which is quantified by the first principal component (bar to the left of heatmap; Supplemental Material). While many of the early lytic genes are actively expressed, very few of the late lytic genes required for virion formation are present, suggesting that most reactivation may be abortive. (B) Propensity for lytic reactivation persists in subclones. Forty-six cell line subclones were assayed in independent labs and have significantly similar latent and lytic viral expression profiles (P < 0.0017, binomial test; Supplemental Material). (C) EBV copy number, as assayed by qPCR in a separate lab from the RNA-seq experiments, correlates with lytic reactivation clustering and confirms the heritability of expression program. See also Figure S3.

Figure 4. Identification of host factors and genes that correlate with spontaneous lytic reactivation

(A) Intercellular host cytokine concentrations, as measured by ELISA, are correlated with viral reactivation propensity (Holm corrected p-values, see Supplementary Material; error bars show standard deviation). (B) Specific EBV lytic and latent genes are expressed at similar levels relative to highly expressed human genes (error bars show standard deviation across the 201 samples). (C) The EBV transcriptome is expressed at similar levels to the human transcriptome. (D) The human WNT5B gene expression is highly correlated with the EBV lytic reactivation marker BHLF1. (E) Human genes that correlate with BHLF1 fall into distinct signaling pathways, including B cell receptor (BCR), interferon alpha (IFN-α), Wnt, and B cell chronic lymphocytic leukemia (BCLL) signaling pathways. These pathways represent an amalgamation of EBV-affected alterations and cellular anti-viral response mechanisms to increased viral load. (F) The viral transcription factor EBNA1 (EBV-encoded nuclear antigen 1) binds to genes that are activated during lytic reactivation. Correlation between BHLF1 and ChIP-seq derived high confidence targets of transcription factors is shown. (G) Host proteins that interact with EBV proteins have their RNA transcripts upregulated during lytic reactivation, as shown by increased correlation with BHLF1. See also Figure S4.

Figure 5. Pax5 binds the viral terminal repeats and regulates viral oncogene expression

(A) Pax5 and other transcriptional regulators bind the LMP control region. Pax5 binds the terminal repeats (TR) directly at two sites, as supported by ChIP-seq and sequence motif. In the promoter, EBF1, PU.1 and SP1 were previously documented to co-bind a GC-box element and regulate LMP1, while TCF12 has not been reported previously. (B) Pax5 binds with high occupancy at the TR locus compared with two known binding sites (EBV Wp promoter and human CD19 promoter) as determined by ChIP-qPCR. We also confirm that Pax5 binds the TR in multiple latency types and EBV strains. MUTU presents type I latency and does not express LMP1, whereas MUTU-LCL presents type III latency and expresses LMP1. (C) Pax5 protein expression is depleted by two independent lentiviral shRNA constructs. Western blot of LCL infected with shPax5-1 or shPax5-2 lentivirus and assayed 5 days post-infection followed by puromycin selection. (D) Western blot of BZLF1 and actin in shPax5-1 or shPax5-2 infected LCLs. (E) Real-time qPCR analysis of EBV genome copy number relative to cellular actin. Values are normalized to shControl infected LCLs. Error bars are standard deviation from the mean (n=3). (F) PFGE analysis of Raji cells infected with shControl, shPax5-1, or shPax5-2 for 5 days post-infection and then analyzed by Southern blot analysis. Circular (C), linear (L), and sub-linear (SL) forms of the viral genome are indicated. (G) Knockdown of Pax5 induces LMP1 and LMP2, while decreasing expression of EBNA1 and EBNA2. Pax5 depletion also leads to an increase in the BZLF1 lytic activator. All transcripts are quantified by qRT-PCR and normalized to cellular actin. Error bars are standard deviation from the mean (n=3). See also Figure S5.

Figure 6. Cohesin regulates EBV chromosome conformation and latent cycle gene expression

(A) ChIP-seq profiles for chromatin structure protein CTCF and Cohesin components (Smc3 and Rad21) demonstrate co-binding at several loci. Chromosome conformation capture (3C) anchor and acceptor primers (shown above) were used to determine DNA loop structures of the viral episome. (B) Co-occupancy of CTCF, Smc3, and Rad21 at LMP1 and OriP is accompanied by additional transcription factors JUND, BATF, TCF12, EBF1 and histone modifications. (C) Identification of a DNA loop between CTCF-cohesin at LMP1 (166.5kbp) and OriP (~6–9.5kbp). LCLs (top) subjected to 3C using MseI restriction digestion show specific interaction between oriP and 166.5kbp (anchor). Positive control (bottom) was performed on bacmid ligation products, where all permutations of acceptor sites and anchor are created. The assay is qualitatively specific to under 1kbp, as shown by the multiple negative regions proximal to OriP. (D) The long-range interaction is highly robust and specific. 3C-qPCR shows a highly specific quantitative interaction with the region surrounding OriP. Height of bars represents qPCR signal relative to positive control bacmid ligation products. Width of bars demonstrates the size of the MseI restriction enzyme fragment, which is also given at the top of the plot. Error bars show standard deviation from the mean (n=3). (E) Depletion of cohesin components through lentiviral shRNA results in loss of long-range interaction. 3C-qPCR as described in panel D, except LCLs were transduced with shControl (top panel), shSMC1 (middle panel), or shRad21 (lower panel). (F) Lentiviral shRNA constructs deplete protein expression for Smc1 and Rad21. (G) Depletion of cohesin components activates transcription of latency genes. Expression of EBNA1, EBNA2, LMP1, LMP2A, and BZLF1 mRNA was assayed by qRT-PCR in LCLs infected with shSmc1, shRad21, or shControl. mRNA expression is normalized to shControl. Error bars are standard deviation from the mean (n=3). (H) Model depicting a DNA-loop between the OriP and LMP loci, which is mediated by CTCF-cohesin binding sites and contains the Pax5 binding sites in the terminal repeats. Acronyms: OriP: latent origin of replication; FR: family repeats; DS: dyad symmetry repeats; EBER: EBV encoded RNA; Cp: BamHI C promoter. See also Figure S6.