Specimens were obtained from the Fred Hutch Research Cell Bank, which prospectively collects and cryopreserves peripheral blood mononuclear cells (PBMCs) from donors and recipients. The University of Washington Institutional Review Board approved use of the iciHHV-6 specimens from the Fred Hutchinson Cancer Research Center and use of anonymized excess HHV-6-positive samples submitted for testing at the University of Washington Virology lab.

SMC cells harboring integrated HHV-6A virus (iciHHV-6A cells) were immortalized following transduction with a lentiviral vector (pLenti CMV/TO SV40 small + large T (W612-1), a gift form Eric Campeau and obtained from Addgene (plasmid #22298)) and expressing the SV40 T antigens. SMCs were cultured in Dulbecco’s modified Eagle’s medium (DMEM) supplemented with 20% fetal bovine serum (FBS), 1% penicillin streptomycin (Pen/Strep), and 1% L-Glutamine. Human epithelial kidney 293 T (293 T, ATCC CRL-11268) cells were cultured in the same medium but supplemented with 10% FBS. All cells were maintained in 10 cm² flasks as a monolayer culture in a humidified 5% CO2 air incubator at 37°C. Bacterial artificial chromosome (BAC)-derived HHV6A (strain U1102) expressing green fluorescent protein (GFP) under the control of the HCMV major immediate early (IE) promoter (293-HHV-6A virus) was propagated in JJHan cells as described previously (Tang et al., 2010b). 293 T cells were infected with 293-HHV-6A virus and GFP positive cells were isolated using a FACS AriaIII cell sorter (BD Biosciences). Clones harboring the integrated 293-HHV-6A genome (ciHHV-6A cells) were identified by quantitative PCR (qPCR) and HHV-6A integration was confirmed by fluorescent in situ hybridization (FISH). iciHHV-6B+ human lymphocytes (NCBI GenBank: KY315552) were isolated from an infected patient as described by Aswad et al. (Aimola et al., 2020; Aswad et al., 2020).

FISH was performed with digoxigenin-labeled bacterial artificial chromosome (BAC) probes against both the HHV-6A genome and individual chromosomal reads essentially as described previously (Kaufer et al., 2011; Kaufer, 2013). Slides were mounted using DAPI Vectashield (Vector Laboratories) and images were taken with an Axio Imager M1 (Zeiss).

RNA was extracted from an HHV-6A infected 293T cell line using Trizol Reagent (Life Technologies) and purified using Direct-zol RNA MicroPrep Kit (Zymo Research #R2060) following the manufacturer’s instructions. One microgram of total RNA was depleted of ribosomal RNA (rRNA) using the KAPA RiboErase Kit (#KR1142) and libraries were prepared using the KAPA Stranded RNA-seq Library Preparation Kit (#KR0934). Libraries were quantified using Qubit (Life Technologies), and quality was assessed using the Agilent Bioanalyzer High-Sensitivity DNA kit (Agilent Technologies). Barcoded libraries were pooled and sequenced on an Illumina HiSeq 2000 to obtain 100-bp paired-end reads. RNA-seq data were processed using STAR as described previously (Dobin et al., 2013) using human (hg38) and a custom HHV-6A reference genomes (NC_001664.4) modified to contain GFP transgenes. BAM files were sorted and converted to SAM files using SAMtools (Li et al., 2009) and reads were counted with featureCounts (Liao et al., 2014) against the corresponding viral GTF file and normalized in DESeq2 (Love et al., 2014). Alignment of viral reads across the HHV-6A genome was visualized using the “seqsetvis” R library (Boyd, 2020) and in-house scripts.

In order to identify the higher order chromatin structure and virus-host DNA contacts, 4C-seq libraries were prepared using both ciHHV-6A and iciHHV-6A samples as described previously (Krijger et al., 2020), using HindIII/DpnII restriction enzymes. Cells were cultured as previously described (Saviola et al., 2019) and fixed with 1% formaldehyde prior to being snap frozen. For each 4C library, ten million cells were used unless otherwise noted. Inverse PCR was performed in 50 µl reactions using 200 ng 3C library template, 25 µl Q5 Hot Start High-Fidelity DNA Polymerase 2X Mastermix (New England Biolabs), 1.5 µl of 10uM forward and reverse inverse PCR primers and water under the following conditions: 30 s 98°C, 10 s 98°C, 30 s 55°C, 2 m 72°C, 5 m 72°C for a total of 25 cycles. PCR products were cleaned up using the PCR purification kit (Thermofisher Scientific) using the B3 reagent to exclude fragments less than 300 base pairs or with Ampure XP beads (Beckman Coulter) and eluted in water. Library indexing was performed using the DNA HT Dual Index Kit (Illumina) or with custom dual indexing primers with an additional eight cycles of PCR. Indexing of nanopore libraries was performed with custom 12-nucleotide dual index primers with eight additional PCR cycles. All primers sequences are listed in Supplementary Table 1 . Libraries were pooled and sequenced with either HiSeq2000 or MinION platforms. The experiments were performed with samples from separate cell culture dates. For the Illumina viewpoint comparison libraries, two replicates were performed for two separate viewpoints. For the ONP titration and viewpoint comparison studies, one and two replicates were performed, respectively.

We developed a 4C-seq primer design tool ¹ to facilitate the generation of the inverse PCR primers necessary for 4C-Seq assays in a genome-agnostic manner for use with experiments involving custom genomes (i.e., chromosomally integrated HHV-6A). A fasta file containing the HHV-6A genome (NC_001664) was obtained from NCBI and uploaded into the design tool, selecting HindIII and DpnII restriction endonucleases, and minimum fragment sizes 500 bp for the first fragment and 300 bp for the second fragment were used (van de Werken et al., 2012).

4C-seq trans interactions were determined using a “window” method as described by Kim et al., (2020). Briefly, the virus viewpoint regions at the beginning of the reads were trimmed off using Cutadapt (Martin, 2011) software and the remaining human portions were aligned to the human UCSC hg38 reference genome using bowtie2 (version 2.3.5) (Langmead et al., 2009; Langmead and Salzberg, 2012). P-values for each mapped base were then calculated for each 10kb window across all chromosomes using a Poisson formula. Final, significant peak regions are then called using the MACS2 (version 2.2.7) bdgpeakcall function (Zhang et al., 2008). Significant Interacting trans peak regions were identified using the above window method for each replicate, for each viewpoint, in both the ciHHV-6A and iciHHV-6A samples. These regions were plotted as circos plots using the Circlize R library (Gu et al., 2014). Bedtools (Quinlan and Hall, 2010) was used to find the intersection of the trans regions and the 127 epigenomes 25-state imputation based chromatin state model created with ChromImpute (Ernst and Kellis, 2012; Ernst and Kellis, 2015) and made available from the Wang Lab ² (Jin Wook Lee). For ciHHV-6A cells, chromatin state annotation was also performed using fetal kidney and fetal adrenal gland tissues – also made available from the Wang lab². 4C trans peaks and signals were also compared to signals for the repressive histone marks H3K9me3, NCBI GEO: GSE66530 (Hattori et al., 2016), and H3K27me3, NCBI GEO: (Lamb et al., 2019), in HEK293 cells and compared to the same marks, NCBI GEO: GSE121984, in iciHHV-6A cells (Saviola et al., 2019).

4C-seq cis interactions were determined using the program peakC (Geeven et al., 2018). Telomere content comparison was assessed by TelomereHunter (Feuerbach et al., 2019). Illumina data was simulated using ART (Huang et al., 2012) set to HiSeq2500 data. Reads were produced to cover the entire UCSC hg38 reference genome at 1X coverage. Then, a number of simulated reads equivalent to the number of reads sequenced between 4C data were randomly sampled from the simulated data. TelomereHunter was run on all samples and a non-parametric Wilcoxon test was performed to compare 4C-seq reads to simulated reads.

To identify the chromosomes harboring HHV-6A integration, a scoring system was developed and coded using the R programming language. Accordingly, trimmed reads were aligned to the hg38 genome using bowtie2, or BWA-MEM (Li, 2013) for MinION reads and each read that maps within 500 kb of any autosome terminus was scored. The total scores for reads within these 500 kb regions are tallied and then compared to each other in one of two ways. If replicates are available, a two-way ANOVA is first performed with chromosome end sum scores as the outcome variable and individual chromosome ends as the categorical predictor variable. TukeyHSD pairwise-post hoc tests are then performed and the chromosome end with the highest mean -log₁₀ p-value is determined to be the most likely candidate. A significance of 0.05 was used for this approach. If replicates were not available for a specific viewpoint/cell type combination, the chromosome end scores were all compared to each other using Wilcoxon nonparametric tests. The chromosomal end with the largest -log₁₀ p-value adjusted for multiple test comparison was chosen as the most likely candidate and an alpha of 0.05 was used as a cutoff of significance. Furthermore, statistical ranking via the aforementioned pairwise post-hoc tests from the algorithm score sums (for each chromosome end) to identify the ends of each chromosome as most significant candidate for HHV-6A integration.

We recently reported that the chromosomally integrated HHV-6A genome exists in a compacted transcriptionally silent state (Saviola et al., 2019). To study the kinetics of HHV-6A gene silencing upon HHV-6A infection, we performed a 7-day time course in human 293T cells. 293T cells are susceptible to HHV-6A infection and were previously established as a model for studying virus integration (Arbuckle et al., 2010; Gravel et al., 2017; Saviola et al., 2019). Cells were infected with recombinant HHV-6A virus expressing GFP under the control of the major immediate-early (IE) HCMV promoter (Tang et al., 2010a; Wallaschek et al., 2016). Infected GFP-positive cells were isolated by FACS sorting 16 h post infection and subsequently cultured. The presence of the HHV-6A genome was confirmed by qPCR and integration was validated by FISH at day seven post-infection. Infected cells were collected daily until 7 days post-infection (dpi) and processed for gene expression profiling via RNA-sequencing ( Figure 1A ). The expression of all HHV-6A genes progressively decreases over the 7-day period as visualized via a clustered heatmap ( Figure 1B ). The observed gene silencing over this period course can be grouped into three clusters based on expression levels of all viral genes (early, mid and late). The collective expression of genes is reduced over a 3-day period, followed by near complete silencing from 5–7 dpi ( Figure 1C ). Thus, HHV-6A gene expression upon infection is silenced within 7 days in 293T cells harboring the integrated virus genome.

The epigenetic mechanisms that regulate gene silencing in chromosomally integrated HHV-6A/B remain poorly characterized. We previously demonstrated that the integrated viral genome forms repressed heterochromatin domains (Saviola et al., 2019). We further hypothesized that higher-order chromatin interactions play a role in virus gene silencing. To investigate higher-order chromatin interactions between chromosomally integrated HHV-6A, we designed 4C-seq assays using distinct viewpoint regions designed against distinct HHV-6A regions ( Figure 2A ). To obtain the inverse PCR primers required for 4C analysis, we developed a genome agnostic 4C primer design tool (see methods). Two distinct viewpoints (vp1 vs. vp2) were designed based on optimal fragment length and restriction enzyme sequence along the length of the ~160 kb HHV-6A genome. The HHV-6A genomic regions targeted by viewpoint vp1 and vp2 primers are located adjacent to the U39 gene-encoding glycoprotein B (gB) (position: 64 kb) and the U95 gene (position: 148 kb), respectively (Nicholas and Martin, 1994). With this approach, we assessed the HHV-6A chromatin conformation in a previously described 293T cells line harboring the chromosomally integrated HHV-6A (ciHHV-6A) (Saviola et al., 2019) and in iciHHV-6A+ patient-derived cells (iciHHV-6A). We obtained approximately 2 hundred thousand to 2 million reads per viewpoint replicate across 2 independent 4C-seq assays, of which 85.93% of reads mapped to the combined reduced human+HHV-6A genome (hg38 + NC_001664.4 or recombinant GFP genome, depending on the sample) ( Supplementary Table 2 ). Significant ‘trans’ interacting regions between the HHV-6A genome and the human genome were identified for both in vitro-derived (ciHHV-6A) and patient-derived iciHHV-6A cells ( Figure 2B ). For ciHHV-6A, we found a total of 6 high-confidence interactions that intersected between viewpoints, and a total of 3 interactions were between viewpoints for iciHHV-6A cells. In addition, a number of cis-interacting regions were identified for both viewpoints ( Figure 2C ; Supplementary Table 2 ). In line with other 4C-seq results (Krijger et al., 2020) there were relatively high densities of interactions that occur near and between viewpoints. Overall, there were a number of cis interactions identified in each cell types suggesting that virus genome may reside in a highly folded compartment similar to a topologically associated domain.

To assess the chromatin features of the identified interaction regions, we annotated the identified interacting human regions using available chromatin state segmentations compiled from multiple human tissue types (Ernst and Kellis, 2015). This annotation dataset is derived from 12 epigenetic features and from 127 original reference epigenomes. We found that the human chromatin that interacted with HHV-6A is largely composed of quiescent and heterochromatin states, whereas enhancer and transcriptional activation annotations occur to a lesser degree ( Figure 2D ). Similarly, chromatin state annotations derived from single representative cell types exhibit similar chromatin state patterns at the HHV-6A interacting genomic regions ( Supplementary Figure 3 ). We further inspected the enrichment of the repressive histone modifications H3K9me3 and H3K27me3 at both ciHHV-6A and iciHHV-6A 4C-seq peaks using ChIP-seq data from each respective cell-line ( Supplementary Figure 2 ). In both ciHHV6-A and iciHHV6-A, the majority of the significant 4C-seq peaks overlap directly with these repressive histone marks. These results indicate that the integrated HHV-6A genome forms higher order chromatin structures within the viral genome as well as with regions of the human genome that are enriched with repressive chromatin. Overall, these results provide insight into HHV-6 chromatin structures within the host cell nucleus.

HHV-6A/B integration sites were first detected by FISH (Nacheva et al., 2008), and the virus-host junction of three different integrations was sequenced using a PCR-based approach. Although PCR amplification and Sanger sequencing provided sequence information, this approach requires previous knowledge of the chromosomal location of the virus genome (Huang et al., 2014; Tweedy et al., 2016) and is prone to amplification biases and sequencing errors due to the repetitive nature of the region (Arbuckle et al., 2010; Arbuckle et al., 2013; Huang et al., 2014; Ohye et al., 2014; Gulve et al., 2017). Genome-wide mapping of short read sequence data to human subtelomeres is challenging due the repetitive nature of subtelomeric regions. It is well-documented that higher order chromatin interactions largely occur in cis in large megabase pair (Mbp)-sized topologically associated domains (TADs) (Szabo et al., 2019). Because HHV-6A integrates into a host chromosome, the extra-chromosomal trans 4C-seq interactions that are identified should behave like endogenous intrachromosomal cis interactions, i.e., 4C-seq should show distinct interacting regions between the host and HHV-6A genomes. With this in mind, we hypothesized that 4C-seq data could identify HHV-6A integration sites, and visualization of 4C alignments to the human genome reveals particular telomeric ends that are enriched with 4C signal. Figure 3A shows enrichment at the distal end of chromosome 15q relative to the rest of the genome in ciHHV-6A cells. In addition, clustering of reads at the distal ends of chromosome 15q and chromosome 19q in ciHHV-6A and iciHHV-6A, respectively, can be seen in Figure 3B . Further, we identified significantly more telomeric sequences in our 4C-seq data relative to simulated reads of equal size and GC content ( Supplementary Figure 1 ). To enable systematic analysis of HHV-6A integration sites using 4C-seq, we established a scoring procedure to assess the chance of integration across all chromosomes based on clustered read mapping (see methods). The ends of chromosome 15q and chromosome 19q rank as the highest scoring chromosome ends for ciHHV-6A and iciHHV-6A, respectively, in terms of read mapping density and integration probability ( Figures 3C, D ). This indicates that these two chromosome ends were the most likely candidates for harboring integrated HHV6A in the ciHHV-6A and iciHHV-6A cell models, respectively. We confirmed that these chromosomal regions harbor the integrated virus by FISH using probes specific for the virus genome and the respective human chromosomes ( Figure 3E ). The integrated HHV-6A genome indeed was present in the respective chromosomal loci, highlighting that the 4C-seq analysis is an unbiased method to identify HHV-6 chromosomal integration sites.

Next generation sequencing technologies have made it possible to generate large quantities of sequence data required for a variety of genomic assays including 4C-seq. However, 4C-seq library sizes derived from inverse PCR are relatively large and often difficult to cluster on Illumina flowcells. Long read sequencing platforms, including the Oxford Nanopore Technologies (ONT) minION, have become common an option for long read sequencing and have been applied for structural genome mapping approaches (De Coster et al., 2019). We therefore evaluated minION sequencing on the 4C-seq libraries as a method to identify integration sites from HHV-6 samples. We sequenced the same 4C-seq libraries with minION flowcells and obtained 150,646,076 bp of demuxable reads with a mean read length of 604.0 bp for the ciHHV-6A and iciHHV-6A sequencing run. For the cell counts titration sequencing run we obtained 419,121,238 reads with a mean read length of 730.4 bp. Alignment of reads that pass quality filters to the human genome reveals a similar grouping of telomere proximal reads at 15q and 19q for ciHHV-6A and iciHHV-6A samples, respectively ( Figure 4A ). Current 4C-seq peak calling software are not well-suited for long-read with low quality scores. We thus adapted our own method based off clustering of adjacent 4C-seq reads (methods). Using this approach, we found an overall comparable number of trans peaks called between replicates ( Figure 4B ). Thus, minION sequencing data represents a sequencing method to generate 4C-seq read data.

We applied the minION 4C-seq reads to investigate the required input material required for identifying HHV-6A integration sites by titrating the number of cells for 4C library construction. We compared a range of input cell quantities for 4C library construction (10⁴ to 10⁷ cells) using the iciHHV-6A cells combined with minION sequencing ( Figure 4C ). These results confirm that 19q harbors the integrated HHV-6A genome for iciHHV-6A SMC cells, and this can be reliably detected in as few as 100,000 cells for the 4C-seq assay. However, using just 10,000 cells 19q also scores highest, however there appears to be some ambiguity as 14q also scores highly. Finally, we also applied minION 4C-seq to detect HHV-6A integration sites using frozen patient-derived B cells from iciHHV-6A+ individuals for this assay. Like the iciHHV-6A SMC samples, 19q also scores as the likely integration site for HHV-6 integration in additional patient-derived lymphocyte samples ( Figure 4D ). FISH validation using chromosome 19q probes of the patient samples confirms HHV-6A integration into chromosome 19q (Aswad et al., 2020; Wight et al., 2020). Thus, 4C-seq libraries can be generated as little as 10,000 cells and using cryopreserved lymphocyte cell samples. In summary, this workflow allowed us to efficiently identify the chromosomal ends that most likely harbor the integrated virus genome.