The Epstein-Barr Virus Enhancer Interaction Landscapes in Virus-Associated Cancer Cell Lines

3D genome organization of type III EBV latency gene expressing GM12878 LCL.

We first investigated the EBV epigenome and 3D genome organization in cells expressing type III latency genes. The encyclopedia of DNA elements (ENCODE) GM12878 H3K27ac ChIP-seq data were analyzed by first mapping the sequencing reads to the B95.8 strain of EBV genome using Bowtie. MACS2 was used to identify significant peaks. The most prominent H3K27ac peak was around the BILF2 gene, and was about 3 kb wide. Other peaks were found near LMP1, LMP2A/B, Cp, OriP, and BARF1 (Fig. 1A). RNA-seq reads were mapped using STAR. The analysis found that type III latency genes EBNAs and LMPs were expressed. Abundant transcripts were also found at BHRL1, OriLyt, BMRF1, BALFs, BARF1, and BILF2 (Fig. 1A). ENCODE CTCF and cohesin subunit RAD21 ChIP-seq data were mapped to the EBV genome. CTCF and RAD21 bound strongly at a site between LMP1 and BARF1. Another strong CTCF peak was near BXRF1. The region between BXRF1 and BARF1 may be within a functional domain shielded by the CTCF/cohesin sites.

We then reanalyzed published GM12878 LCL H3K27ac HiChIP data for interactions within EBV episomes. For HiChIP assays, cells were first cross-linked by formaldehyde and lysed. Mbo I was used to cut DNA. The DNA ends were filled in with biotinylated nucleotides and ligated in situ. H3K27ac chromatin immunoprecipitation (ChIP) was used to enrich the DNA fragments with H3K27ac marks. Avidin beads were used to enrich the DNA fragments with the ligation products. Paired-end sequencing was used to identify the interacting genomic regions. Bowtie was used to map the sequencing reads against EBV sequence to identify the ligation products derived from the EBV genome. HiChIP interaction reads were counted between pairwise 1-kb bins of EBV genome. The count table, consisting of data from the 7 cell lines in this analysis (GM12878, KAI3, NUNK1, SNT-13, SNT-16, BC1, and JSC1), was normalized using the Voom tool with the “quantile” option to account for the sparsity and dispersion variations in count-based data. Circos plots were generated based on a normalized interaction frequency of 4 or higher, which accounted for the top ~5% of interactions. Abundant interactions between different EBV genomic regions were identified. The strongest interaction was between the strong enhancer at the BILF2 and the BALF5 or BARF1. The distances between BILF2 and BARF1 were ~18 kb, excluding the B958 deletion. BARF1 and BILF2/BdRF1 enhancers also interacted with many regions of the rest of EBV genome (Fig. 1A and B). Less frequent interactions between BILF2 enhancer and the genes on the opposite side of the BILF2 enhancer, including Cp, BHLF1, and Zp, were also evident, (Fig. 1B, lower part).

CRISPRi and CRISPRa were used to perturb the BILF2 enhancer to determine its activity in regulating the genes linked by HiChIP, including BALF2 and BARF0. Lentiviruses expressing dCAS9-KRAB-MeCP2 or dCAS9-VP64 were transduced into GM12878 LCLs and selected with blasticidin. Single guide RNA (sgRNA) targeting BILF2 enhancer or non-targeting control was packaged into lentiviruses and transduced into LCLs expressing stable dCAS9 fusion proteins. After puromycin selection, cells were grown for 2 days. Total RNAs were extracted from these cells. Reverse transcription-quantitative PCR (qRT-PCR) was used to measure the expression levels of BILF2 and BALF2. CRISPRi significantly repressed BILF2 (P < 0.001) and BALF2 expression (P < 0.05 or P < 0.01) (Fig. 1C). The CRISPRi did not affect the expression of nearby BILF1, while some sgRNAs increased BALF4 and BALF5 expression (P < 0.05) (Fig. S1A in the supplemental material). CRISPRa significantly increased BILF2 and BALF2 expression compared with the control sgRNA. CRISPRa had no effect on LMP1 expression (Fig. 1D). These data indicated that the linkage between BILF2 enhancer and gene in the BALF2 locus is functional. Perturbation of BILF2 enhancer affected the expression of genes linked by H3K27ac HiChIP.

3D genome organization of type II EBV latency gene expressing T/NK lymphoma cells.

EBV-positive T/NK lymphomas are aggressive non-Hodgkin’s lymphoma. RNA-seq analyses determined that these lymphomas express type II latency genes, including EBNA1, LMP1, LMP2A/2B, and EBER. In addition, the EBV lytic genes BNRF1, BILF1, BALF5/4/3/2, and BNLF2b were also detected in these patients (39–42). Interestingly, in the same cohort of patients, whole-exome sequencing found fewer mutations than in other lymphomas that are EBV-negative (39). In another study, EBV immediate early gene BZLF1 can be detected by qRT-PCRs in some T/NK lymphomas (40). To identify model cell lines faithfully represent the primary tumor samples, we analyzed the expression profiles of EBV positive T/NK lymphoma cell lines, including two NK lineage lymphoma cell lines KAI3 (43) and NUNK1, and two T cell lineage lymphoma cell line SNT13 and SNT16 (44) by RNA-seq. For KAI3 cells, type II latency genes were expressed (Fig. 2A). In addition, BILF2, LF1/2, BALF5/4/3, BARF1, BNLF2a/b were also expressed (Fig. 2A and S).

H3K27ac ChIP-seq was used to identify enhancers in EBV genome in KAI3 cells. Multiple enhancers in the EBV genome were identified. The major H3K27ac peak was located around BdRF1/BVRF2/BILF1. Other major peaks were around LMP1 and EBERs. Minor peaks were found around oriLyt, BNLF2a/b, and oriP (Fig. 2A and B).

H3K27ac HiChIP was used to connect enhancers to their direct target genes. In KAI3 cells, BILF2 enhancer frequently interacted with LF1/2, BALF5/4/3/2, and BARF1 regions (Fig. 2A and B). Furthermore, BILF2 enhancers frequently interacted with the rest of the EBV genome, suggesting that BILF2 enhancer functioned at an interaction hub within the EBV genome in type II latency (Fig. 2A and B).

CRISPRi was used to evaluate the functional significance of the loop between BILF2 enhancer and BARF1 expression. sgRNAs targeting the BILF2 enhancer or non-targeting control sgRNAs were packaged into lentiviruses and transduced into KAI3 cells stably expressing dCAS9-KRAB-MeCP2 fusion protein. After puromycin selection, total RNAs were extracted from these cells. qRT-PCR was used to evaluate the effect of CRISPRi on BARF1 expression. CRISPRi at BILF2 enhancer significantly repressed BARF1 expression (P < 0.01). CRISPRi also repressed BILF2 expression (P < 0.01 or P < 0.05) (Fig. 2C). These data indicated that BILF2 enhancer can regulate the expression of its direct target gene BARF1 through long-range enhancer looping.

Similar expression and enhancer interaction patterns were observed for NUNK1, SNT13 and SNT16 cells (Fig. 3A and B).

3D genome organization of type I EBV latency gene-expressing primary effusion lymphoma cells.

The PEL cell lines BC1 and JSC express EBV type I latency genes. H3K27ac ChIP-seq of BC1 and JSC cells found peaks at the EBER loci, miRNA loci, and Qp, which drives the expression of EBNA1 in type I latency cells (45). H3K27ac HiChIP was used to evaluate the genome looping in these cells. Much less frequent interactions between BILF2 and the BALF gene cluster were found in BC1 and JSC cells. More frequent interactions were found between BILF2 and the genes on the opposite side (Fig. 4, bottom). We found no obvious enrichment of an anchor region that interacts with the rest of the genome (Fig. 4).

Cells.

EBV-positive NK/T cells KAI3, NUNK1, and SNT16 were cultured in Artemis medium 2 (Nihon Techno Service (NTS), Ushiku, Japan). NUNK1 cells were established from the peripheral blood of a 17-year-old boy suffering from NK-type chronic active EBV disease. NUNK1 cells were positive for CD56, HLA-DR, and CD2, but negative for CD20, CD21, CD3, CD4, and CD8. SNT13 was cultured in Artemis medium 2 supplemented with 700 U/mL IL-2. GM12878 was cultured in RPM1640 medium supplemented with 10% fetal bovine serum (FBS). 293T cells were maintained in Dulbecco’s modified Eagle’s medium with 10% FBS.

Ethics statement.

Written informed consent was obtained from the parents of patients. The Institutional Review Board of Nagoya University Graduate School of Medicine approved this study.

HiChIP.

HiChIP was performed according to described previously protocols (37). In brief, 10 million cells were harvested and cross-linked with 1% formaldehyde. Cells were lysed with Hi-C lysis buffer (10 mM Tris-HCl pH 7.5, 10 mM NaCl, 0.2% NP-40, 1× Roche protease inhibitors). DNA was digested by 375 U Mbo I (New England Biolabs [NEB] cat no. R0147M) at 37°C with shaking at 300 rpm for 2 h. DNA ends were filled in with biotin-dATP (Thermo Fisher Scientific cat no. 19524016), dCTP, dGTP, and dTTP mix, and DNA polymerase I, Klenow fragment (NEB, M0210) with shaking at 37°C for 1 h. Next, 4,000 U T4 DNA Ligase (NEB, M0202) was used for proximity DNA ligation with shaking at room temperature for 4 h followed by incubation at 16°C overnight. After ligation, DNA was sonicated using a Covaris M220. Fragmented DNAs were then diluted 10 times with ChIP dilution buffer (0.01% SDS, 1.1% Triton X-100, 1.2 mM EDTA, 16.7 mM Tris-HCl [pH 7.5], 167 mM NaCl) and samples were precleared with protein A Dynabeads at 4°C for 1 h. DNA-protein complexes were captured with 8 μL H3K27Ac antibody at 4°C overnight. DNA-protein complexes were captured with protein A beads with rotation at 4°C for 2 h. Beads were washed and DNAs were eluted twice with 150 μL of freshly prepared DNA elution buffer (50 mM sodium bicarbonate [pH 8.0], 1% SDS). The eluted DNA was reverse cross-linked. DNA was purified with a PCR purification kit (Qiagen, Hilden, Germany). Biotin dATP-labeled DNA was captured with 5 μL Streptavidin C-1 beads (Thermo Fisher, number 65001) and DNA was then transposed with 2.5 μL Tn5. Beads were then washed. After the final wash, beads were re-suspended with 23 μL double-distilled water, 25 μL 2× Phusion HF (New England Biosciences), and 1 μL each of Nextera forward primer (Ad1_noMX) and Nextera reverse primer (Ad2.X) at 12.5 μM. PCRs were run at: (i) 72°C for 5 min, (ii) 98°C for 1 min, (iii) 98°C for 15 s, (vi) 63°C for 30 s, and (v) repeat steps 1 to 4 for a total of 8 cycles, and final extension at 72°C for 1 min. After PCR amplification, a two-sided size selection with Ampure XP beads was performed for DNA size selection and purification. Samples were then sequenced on an Illumina NextSeq platform (2 × 75 bp).

ChIP-seq.

CHIP-seq was performed following the methods of Zhao et al. (13). Cells were cross-linked with 1% formaldehyde and DNA was sonicated into fragments with an average size of around 500 bp. Next, 4 μg H3K27ac antibody (Abcam cat no. ab4729) was used to immunoprecipitate histone H3 acetylated on lysine 27 and associated DNA. DNA was purified using a PCR purification column (Qiagen). ChIP-seq DNA libraries were prepared using a NEBNext Ultra II DNA Library Prep kit (E7645S). DNA libraries were sequenced on the Illumina platform.

CRISPRi and CRISPRa.

To silence the activity of BILF2 enhancer, CRISPRi was used. In brief, GM12878 and KAI3 cell lines were transduced with lentiviruses expressing dCAS9-KRAB-MeCP2 fusion protein, followed by selection with 5 μg/mL blasticidin for 5 days. Cells stably expressing dCAS9-KRAB-MeCP2 were then transduced by lentiviruses (PXPR-502) expressing sgRNAs (BILF2isg1-F: CACCGTAATTGCAGTGGACCCCGG, BILF2isg1-R: aaacCCGGGGTCCACTGCAATTAC; BILF2isg2-F: CACCGGAACACATATAACTAACGG, BILF2isg2-R: aaacCCGTTAGTTATATGTGTTCC; BILF2isg3-F: CACCGAGCAACATGGAGACCTCGAA, BILF2isg3-R: aaacTTCGAGGTCTCCATGTTGCTC; BILF2isg4-F: CACCGATTAGGCCTAAAGCCACCTG, BILF2isg4-R: aaacCAGGTGGCTTTAGGCCTAATC; BILF2isg5-F: CACCGCCCCTGCAGAGTAATTGCAG, BILF2isg5-F: aaacCTGCAATTACTCTGCAGGGGC). CRISPRa was used to activate BILF2 enhancer. GM12878 stably expressing dcas9-VP64 was made by transducing lentiviruses expressing dcas9-VP64 fusion protein. Lentiviruses expressing sgRNAs were used to transduce GM12878 cells expressing dcas9-VP64.

RT-qPCR.

After puromycin selection for 3 days, reverse transcriptase quantitative PCR (RT-qPCR) was used to measure gene expression. Total mRNAs were extracted using a PureLink RNA Mini Kit (Life Technologies, Carlsbad, CA). mRNAs were converted to cDNA using iScript Reverse Transcription SuperMix (Bio-Rad, Hercules, CA). cDNAs were diluted 20 times with RNase-free H₂O before amplified on an CFX96 Touch real-time PCR detection system (Bio-Rad). SYBR Green (Thermo Fisher Scientific) was used to detect cDNA amplification. GAPDH (glyceraldehyde-3-phosphate dehydrogenase) was used to normalize gene expression. RNA relative expression was calculated using the threshold cycle (2^–ΔΔCT) method. The value for the cells transduced with non-targeting sgRNA was set to 1. The primer pairs qRT-BARF1-F, GGCTGTCACCGCTTTCTTGG; qRT-BARF1-R, AGGTGTTGGCACTTCTGTGG; qRT-BILF2-F, TCTAAACGTCCCACGCATAAG; qRT-BILF2-R, CCACAGCACAAGTGATTAGGA; qRT-LMP1-F, GCCCTTTGTATACTCCTACTGATG; qRT-LMP1-R, AGACAAGTAAGCACCCGAAGAT; qRT-BALF2-F1, CCGGGCTTCAGCATCAAT; and qRT-BALF2-R1, GGCATTGTGGAACACGTAGA were used to measure BILF2, BALF2, BARF1, and LMP1 expression, respectively. The primer pairs qRT-BALF4-F, CGATGAGGGAGGTGTTTAGTG; qRT-BALF4-R, GATCCACGTCTACAACGACTAC; qRT-BALF5-F, GCAAAGAGGCGAGGAATCT; qRT-BALF5-R, GACTACAAGCTGGACACAGTAG; qRT-BILF1-F, TCAGGTACAGGCACACATTTC, and qRT-BILF1-R, CGTATGGACTCACTTGGGTATG were used to measure BILF2, BALF2, and BARF1.

Data analysis.

HiChIP paired-end reads were mapped to human (hg19) and EBV (Akata) genomes using HiC-Pro v2.11.133 (default settings with LIGATION_SITE set as GATCGATC for Mbo I). Due to the small genome size and densely located genes on the EBV genome, we examined EBV chromatin interactions based on genomic bins instead of anchors. In brief, the EBV genome was divided into 1-kb consecutive windows. HiChIP paired-end reads with both ends mapped to EBV genome were counted as bin interactions for two paired bins if the two read ends were located in the two bins separately. The interaction count table based on pairwise bins consisting of all cell lines was first filtered to keep only bin pairs with at least one interaction in one cell line, and was then normalized using the Voom tool with the “quantile” option to account for the sparsity and dispersion variations in count-based data. Circos plots were generated based on the normalized interaction frequency of 4 or above, which account for the top ~1% of pairwise bin interactions.