Introduction

Human cytomegalovirus (HCMV) is nearly ubiquitous in the human population. While primary infection of immunocompetent individuals is typically mild or asymptomatic, the virus can cause significant disease in developing neonates and immunocompromised individuals ^1,2. HCMV establishes a lifelong latent infection following primary infection that can periodically reactivate throughout the patient’s lifetime, leading to shedding of infectious virus or promotion of virus-associated morbidity in immunocompromised individuals ³. HCMV can infect a broad range of cells including fibroblasts, epithelial cells, and hematopoietic progenitor stem cells ⁴. Fibroblasts and epithelial cells are typically supportive of a lytic infection, while hematopoietic cells are the reservoir of latent or persistent infection. Fibroblasts have been widely used to study the cellular changes that occur in response to viral infection and for viral propagation in cell culture ^4–6. Like all other viruses, HCMV is an obligate intracellular pathogen that relies on the host cell machinery for its replication, assembly, and egress. These complex host-pathogen interactions vary among different cell types, leading to viral infection manifesting itself in primarily a lytic or latent fashion ^7,8. Infection of permissive cells is accompanied by widespread changes in host cell gene expression. Indeed, several previous studies have examined pathways altered by HCMV infection of in-vitro cultured cells at the level of gene expression ^8–12.

Regulation of gene expression involves intricate interplay between the binding of transcription factors and cofactors, chromatin accessibility, histone marks, and long-range chromatin looping interactions ¹³. A hallmark feature of viruses is their dependence on the host nuclear environment and host transcriptional processes. For example, some viruses, such as Epstein-Barr virus and Kaposi sarcoma-associated herpesvirus extensively alter human chromatin architecture in order to replicate ^14,15. HCMV uses host-derived histones to chromatinize its genome for temporal regulation of the Immediate Early protein IE1 ¹⁶. However, the mechanisms affecting HCMV driven regulation of host gene expression are not well defined, thus limiting our understanding of the complex interplay between virus and host.

To determine the human gene regulatory mechanisms impacted by HCMV, we examined the changes to chromatin accessibility, transcription factor occupancy, chromatin looping, and gene expression resulting from HCMV infection. We show that HCMV infection in human fibroblasts and retinal epithelial cells induces large-scale global changes in the human chromatin landscape. Accessible chromatin regions that closed upon HCMV infection were highly enriched for predicted TEAD sites and depleted of CTCF (CCCTC-Binding Factor) sites in both human fibroblasts and retinal epithelial cell lines. Chromatin-immunoprecipitation followed by sequencing (ChIP-seq) experiments confirm that TEAD1 binding sites are significantly depleted at sites that become inaccessible with HCMV infection, with concomitant reduction of chromatin looping interactions and H3K27ac levels.

TEAD1 is a downstream effector of the Hippo signaling pathway, which regulates cell proliferation and cell fate to control organ growth and regeneration ^17,18. TEAD1 activity is primarily controlled by its co-activator YAP1 ^18,19. We show that HCMV-induced loss of TEAD1 binding is mediated by four distinct mechanisms: 1. Extensive closing of human chromatin that is normally occupied by TEAD1; 2. Reduction of YAP1 and pYAP1 protein levels; 3. Reduction of TEAD1 transcript and protein levels; 4. Exclusion of TEAD1 exon 6. Consistent with these observations, pathway enrichment analysis of differentially expressed genes upon HCMV infection reveals extensive perturbation of the Hippo/TEAD signaling pathway. GWAS-based enrichment analysis reveals that these TEAD1 binding loss events are specifically enriched for genetic variants associated with ear and eye development. Collectively, our data provide novel insights into the mechanisms employed by HCMV upon infection of human cells, with important implications in HCMV-induced human growth defects.

Results

Human cytomegalovirus infection extensively re-organizes human chromatin upon infection

Several previous studies have reported extensive changes in human gene expression profiles in HCMV infected cells ^20–24. However, the regulatory mechanisms underlying these vast changes remain largely unknown. Gene expression programs defined by cis-acting DNA elements such as enhancers and promoters typically occur in open or accessible regions of chromatin. We thus sought to determine if HCMV infection leads to genome-wide alterations in chromatin structure. To this end, we used ATAC-seq ²⁵ to measure HCMV-induced changes to human chromatin accessibility in two widely used HCMV infection models: human fibroblasts (HS68 cells) and retinal epithelial cells (ARPE-19 cells) (Figure 1A). ATAC-seq experiments were performed in duplicate in uninfected conditions or 48 hours post infection (hpi) (Figure 1B). For each cell line, we obtained >200,000 peaks in both uninfected and infected cells, with Fraction of Reads in Peaks (FRiP) scores greater than 0.5 and Transcription Start Sites (TSS) enrichment scores >30 (average of 35.1) (Supplemental Dataset 1), all of which greatly exceed the recommendations set forth by the ENCODE Project consortium ^26–29. Likewise, our ATAC-seq peaks align strongly with relevant publicly available datasets (Supplemental Dataset 1), and we observed very strong agreement between replicates, with samples tightly clustering first by cell type and then by infection status (Supplemental Figure S1A). Collectively, these results highlight the high quality and internal consistency of our datasets.

Figure 1.

We next sought to systematically identify regions of the human genome with differential chromatin accessibility between infected and uninfected cells. To this end, we used DiffBind ³⁰ (see Methods), which identified thousands of chromatin regions (ATAC-seq peaks) with statistically significant changes in signal (Supplemental Dataset 1). As expected, most accessible chromatin is common to both uninfected and HCMV infected cells in both cell types – we designate such regions ‘unchanged’ (Figure 2A and 2B). In fibroblasts, 38,651 peaks were unique to uninfected cells, and hence were closed following infection with HCMV (Figure 2A and 2C). Likewise, 49,003 regions were newly accessible in HCMV infected fibroblasts. We also observed large HCMV-dependent changes to chromatin accessibility in retinal epithelial cells (Figure 2B and Figure 2D). Collectively, these data indicate that HCMV infection of human cells has widespread effects on the human chromatin accessibility landscape.

Figure 2.

Regions of open chromatin largely reflect regulatory regions such as enhancers, which can interact with promoters at large genomic distances spanning many megabases through chromatin looping interactions ³¹. To measure the impact of HCMV infection on functional chromatin looping across the human genome, we next performed HiChIP with an antibody against H3K27ac in uninfected and infected fibroblasts (see Methods). Analysis of the resulting HiChIP data revealed 143,882 and 97,815 chromatin looping interactions in uninfected and infected cells, respectively (Supplemental Dataset 2). QC analyses using HiC-Pro ³² indicate that the data are of high quality: the final set of unique valid interaction pairs was 67% of the total sequenced pairs for the uninfected cells and 52% for the infected cells. The number of trans interactions was 21% of the sequenced pairs for the uninfected cells and 11% for the infected cells (Supplemental Dataset 2), similar to or better than the results obtained in the original HiChIP study ³³. Peaks were called from the HiChIP data and compared to peaks called from H3K27ac ChIP-seq data (see Methods), revealing between 82% and 92% HiChIP peak overlap with H3K27ac ChIP-seq peaks (Supplemental Dataset 2). Collectively, these results indicate that our HiChIP data are of high quality.

We next used the HiChIP data to identify HCMV-dependent differential chromatin looping events (see Methods), identifying 74,797 loops that open and 120,864 loops that close following HCMV infection (Supplemental Dataset 2). Of the 23,018 loops present in both infected and uninfected cells, only 10 are differential at a 2-fold cutoff and a false discovery rate (FDR) < 0.01. Because the number of loops unique to the infection state were so striking, we repeated loop calling and differential analysis with FDR values less than 0.05, 0.01, and 0.001 (Supplemental Dataset 2). For all three cutoffs, the percentage of loops specific to infected cells remains at least 70%. We also randomly downsampled the number of input pairs used for calling loops to verify that our results were not due to a difference in read depth (Supplemental Dataset 2). For the three smaller subsets of data, the number of unique loops changed by only one percent. HCMV-dependent losses in chromatin accessibility significantly correspond with HCMV-dependent losses of chromatin looping (enrichment: 2.1-fold, adjusted p-value: 9.6×10⁻¹⁹²), but not HCMV-dependent gains of chromatin looping (enrichment: 0.6-fold, adjusted p-value: 1) (Supplemental Dataset 3). The full quantification of each chromatin looping event and comparisons of events between conditions are provided in Supplemental Dataset 2.

Next, we examined HCMV-altered regions of chromatin accessibility for enriched transcription factor binding site (TFBS) motifs (see Methods). We use unchanged accessibility regions as a baseline for motif enrichment analysis, in order to identify motifs specific to chromatin regions that are uniquely accessible in either uninfected or infected cells. In human fibroblasts, all three categories of regions (unchanged, closed with infection, open with infection) were highly enriched for AP-1 motifs (Figure 3A, purple dots), reflecting the important role played by this TF family in virtually every cell type regardless of infection status ^34,35. Intriguingly, we found substantial differences in motif enrichment for other TF families in HCMV-altered chromatin regions. In regions closed upon HCMV infection, CCAAT-enhancer-binding protein (CEBP), E-box (bound by bHLH TFs), and TEA domain (TEAD) binding motifs were uniquely enriched (Figure 3A, left: blue, orange, and maroon dots along the Y-axis, respectively). In contrast, newly closed chromatin in HCMV infected cells was depleted of CTCF binding sites (Figure 3A, left: dark green dots along the X-axis), and newly opened chromatin was enriched slightly more strongly for CTCF sites (Figure 3A, right: dark green dots approaching the Y-axis).

Figure 3.

We next examined transcription factor binding site motif enrichment in retinal epithelial cells. Similar to human fibroblasts, we found that TEAD motifs were specifically enriched in peaks that were closed with infection (Figure 3B, left). Unlike in fibroblasts, we did not identify strong enrichment of CEBP or E-box motifs, presumably due to cell-type differences. In addition to differences in enrichment levels, we also observed substantial and consistent changes in the overall predicted occupancy of TEAD TFs and CTCF in both cell types (Figures 3C and 3D). Taken together, these results indicate that HCMV infection extensively re-organizes the human genome in both fibroblasts and retinal epithelial cells. Further, this reorganization alters the accessibility of predicted binding sites for the TEAD TF family, and avoids closing the binding sites of CTCF, in both cell types.

TEAD1 transcription factor genomic binding is substantially depleted upon HCMV infection

Because chromatin re-organization was more pronounced in fibroblasts (Figure 2), we focused on this cell type for subsequent chromatin immunoprecipitation experiments and analyses. We prioritized the two TF families (TEAD and CTCF) with consistent enrichment patterns across cell types within HCMV-dependent chromatin accessibility. We selected TEAD1 to represent the TEAD family because it is the only member whose expression levels are significantly altered by HCMV infection (see Figure 5D). We also examined the genome-wide distribution of H3K27ac, a histone mark correlated with active enhancers.

ChIP-seq experiments for TEAD1, CTCF, and H3K27ac were performed in duplicate in uninfected and HCMV infected conditions (48 hours post infection). The resulting data were of high quality, meeting or exceeding ENCODE data quality standards (Supplemental Dataset 4), with strong agreement between replicates (Supplemental Figure S2). In particular, all H3K27ac datasets had > 65,000 peaks and all TEAD1 and CTCF ChIP-seq datasets had > 15,000 peaks, with FRiP scores ranging from 0.02 to 0.13 and very strong motif enrichment results (CTCF and TEAD motifs rank #1 in every respective experiment, with p-values < 10⁻³⁰⁰⁰) (Supplemental Dataset 4). As predicted, the number of TEAD1 ChIP-seq peaks was substantially diminished in infected cells (75,554 in uninfected compared to 17,567 peaks in infected cells), along with a substantial drop in H3K27ac peaks (110,308 vs. 66,644) (Supplemental Dataset 4), suggesting that loss of TEAD1 binding might have a strong impact on enhancer functionality. In contrast, the total number of CTCF peaks was largely consistent with and without infection (54,697 vs. 57,936 peaks).

To further quantify these changes, we performed differential ChIP-seq peak analysis using the DiffBind software package ³⁰ (see Methods). For TEAD1, we identified 30,740 peaks specific to the uninfected condition, with only 237 peaks specific to infected cells. Likewise, we identified 40,656 H3K27ac peaks specific to uninfected cells, with 8,542 peaks specific to infected cells (Supplemental Dataset 4). As predicted by the motif enrichment analysis of differentially accessible chromatin, regions where the chromatin closed with infection did not contain CTCF ChIP-seq peaks. Instead, we found that losses of TEAD1 occupancy strongly coincide with losses in chromatin accessibility (enrichment: 9.1-fold, adjusted p-value: 1.4×10⁻²¹⁵), losses in chromatin looping events (enrichment: 2.1-fold, adjusted p-value: 3.5×10⁻¹⁹²), and losses in H3K27ac levels (enrichment: 6.1-fold, adjusted p-value: 8.7×10⁻²¹⁴) (Figure 4A and Supplemental Dataset 3). For example, we observe significant HCMV-induced loss of TEAD1 binding, chromatin accessibility, and H3K27ac levels proximal to the promoters of the Hippo pathway genes FRMD6 and RASSF2 (Figure 4B), both of which have significantly diminished gene expression subsequent to infection (see next section). In the case of FRDM6, which involves a likely enhancer, we also observe loss of a chromatin looping interaction (blue loop) between the promoter (dashed box) and two TEAD1 binding sites lost upon infection (solid box). In the case of RASSF2, a TEAD1 binding site within the promoter (solid box) interacts (blue loop) with a downstream enhancer (dashed box). Both the TEAD1 binding site and the interaction are lost upon infection. These examples indicate that loss of TEAD1 binding rather than CTCF binding correlates with loss of chromatin interactions. Collectively, these observations confirm our motif-based predictions that loss of chromatin accessibility co-occurs with loss of TEAD1 binding, in addition to loss of enhancer function (H3K27ac) and enhancer looping interactions (HiChIP).

Figure 4.

HCMV infection induces large-scale changes to human gene expression, including perturbation of the Hippo signaling pathway

Our data indicate substantial changes to human gene regulatory features upon HCMV infection. Therefore, we next systematically assessed HCMV-dependent alterations to human gene expression levels. To this end, RNA was extracted from uninfected and infected human fibroblasts at 48 hours post infection and analyzed by RNA-seq (see Methods). The resulting data were of high quality (Supplemental Dataset 5) and displayed strong agreement between replicates (Supplemental Figure S1B). As expected, we detected HCMV-encoded genes only in infected cells (Supplemental Dataset 5). We identified 2,352 differentially regulated human genes (1,155 upregulated with infection and 1,197 downregulated) at a 2-fold cutoff with an FDR < 0.01 (Figure 5A and Supplemental Dataset 5).

As expected, pathway enrichment analysis of differentially expressed genes revealed strong enrichment for pathways generally involved in the infection response, including “viral protein interaction with cytokine and cytokine receptor” and “human papillomavirus infection” (Figure 5B). We also observe strong enrichment for many of the major developmental pathways, including WNT, TGF-beta, and Hippo (Figure 5B). Strikingly, many of the genomic regions with extensive (five or more) TEAD1 binding loss events encode major regulators of these pathways, including WNT5B and FZD7 (WNT pathway), SMAD3, TGFB2, and TGFBR2 (TGF-beta pathway), and TEAD1 itself (Hippo pathway), all of which have HCMV-altered expression levels (Supplemental Dataset 5). The extensive loss of TEAD1 binding (with concordant alterations to gene expression levels) suggests that TEAD1 might be a key HCMV-targeted regulator. TEAD transcription factors are direct effectors of the Hippo signaling pathway ¹⁷, and our data support TEAD being at the top of a hierarchy involving many other key developmental pathways.

Overall, we observe extensive HCMV-induced alterations to Hippo gene expression levels. Of the 163 Hippo pathway genes annotated in the KEGG database, 21 were significantly downregulated (including TEAD1) and 20 were upregulated (including PARD6B, a negative regulator of Hippo ³⁶) (Supplemental Dataset 5 and Figure 5C). Notably, only TEAD1 was differentially expressed among the four TEAD family members, with a 3.7-fold decrease upon infection (adjusted p-vale: 1.9×10⁻²⁹) (Figure 5D and Supplemental Dataset 5). We confirmed alterations in expression levels at the protein level for the classic TEAD1 targets Thrombospondin 1 (THBS1) and Cellular Communication Network Factor 1 (CCN1), both of which were substantially downregulated upon HCMV infection. Specifically, THBS1 protein was not detected in infected cells, while CCN1 was reduced by over 3-fold (Figure 5E).

Figure 5.

HCMV infection diminishes TEAD1 transcription factor activity through four distinct mechanisms

We next aimed to further explore the virus-induced mechanisms underlying the reduction in TEAD1 activity. The transcriptional regulatory activity of TEAD1 is primarily controlled by the phosphorylation-dephosphorylation cycle of its co-activator YAP1 ^18,19. We thus examined the protein expression levels and phosphorylation status of YAP1 (pYAP1) with and without infection, along with the expression of HCMV proteins IE1/2, TEAD1, and the H3K27ac histone mark in uninfected and infected cells (Figure 6A, Supplemental Figure S3). Consistent with diminished TEAD1 activity, we observed a reduction in levels of total YAP1 protein (1.8-fold, p<0.03) and pYAP1 levels (1.5-fold, p<0.0005) (Figure 6A, Supplemental Figure S3). It is likely that the reduction of pYAP1 is a direct function of reduced YAP1 expression. These results are consistent with a large-scale HCMV infection proteomics screen performed in human fetal foreskin fibroblasts ³⁷. Likewise, TEAD1 protein levels decreased (3.7-fold, p<0.0009) with HCMV infection (Figure 6A, Supplemental Figure S3). Consistent with the H3K27ac ChIP-seq results (Figure 4B), global H3K27ac was also substantially reduced (1.9-fold, p<0.19) in infected cells. Since YAP1 activity is primarily controlled by phosphorylation of its serine residues (S127 and S381), which leads to its cytoplasmic translocation ³⁸, we next assessed the cytoplasmic-nuclear shuttling of pYAP1. These experiments revealed increased YAP1 and pYAP1 in the nuclear fraction of uninfected cells relative to HCMV infected cells (Figure 6B). Taken together, these results indicate that HCMV-induced loss of TEAD1 activity can also partially be accounted for by loss of active YAP1.

Figure 6.

Recent reports have shown that alternative splicing of TEAD1 also regulates TEAD1 activity ³⁹. We thus used the AltAnalyze software package ⁴⁰ to systematically examine alternative splicing changes between HCMV infected and uninfected cells within our RNA-seq data (see Methods). This analysis identified two differential exon use events within the TEAD1 gene that are impacted by HCMV infection: 1) skipping of exon 6 with HCMV infection; and 2) diminished inclusion of an upstream intronic region just upstream (∼90 bases) of exon 5 with HCMV infection (Supplemental Figure S4 A and B). A recent study carefully examined the functional effect of the alternative exon 6 usage event, concluding that skipping of exon 6 diminishes TEAD1 activity by disrupting the intramolecular interaction between its DNA-binding domain and the YAP1-interacting domain ³⁹. We confirmed HCMV-induced differential splicing of TEAD1 exon 6 by RT-PCR of total RNA isolated from uninfected and HCMV infected fibroblasts (Figure 6C and Supplemental Figure S4C).

Taken together, our results reveal multiple interrelated mechanisms by which HCMV infection impairs TEAD1 activity (Figure 6D): 1) Exclusion of TEAD1 exon 6 (Figure 6C); 2) Reduction of TEAD1 gene and TEAD1 protein levels (Figures 5D and 6A); 3) Reduction of YAP1 and pYAP1 protein levels (Figure 6A); and 4) Extensive closing of human chromatin that is normally occupied by TEAD1 (Figure 4).

HCMV-induced TEAD1 binding loss coincides with genetic variants associated with ear and eye growth defect phenotypes

Next, we examined the enrichment of GWAS hits within regions of HCMV-induced loss of TEAD1 binding. To this end, we used our RELI tool ⁴¹ to estimate GWAS-associated risk locus enrichment for 1337 diseases and phenotypes. This analysis produced three phenotypes with specific enrichment for TEAD1 binding loss events compared to unchanged TEAD1 binding events (see Methods and Supplemental Dataset 3): lobe attachment ^42,43, optic cup area ^44–47, and percent mammographic density ^48–50. Among these, enrichment replicated in multiple GWAS for two phenotypes: lobe attachment and optic cup area (Figure 6E). Notably, HCMV infection of fibroblasts has previously been shown to contribute to hearing impairment ^51–54 and MCMV infection can lead to hearing loss in mice ⁵⁵. Likewise, Hippo signaling plays established roles in ear growth and hearing loss ^56–59. Of note, the thrombospondin (THBS1) gene, which plays important, well-established roles in hearing loss ⁶⁰, has extensive proximal HCMV-induced TEAD1 binding loss and HCMV-induced gene (Supplemental Dataset 5) and protein level reduction (Figure 5E). Meanwhile, the “optic cup area” phenotype refers to the overall size of the optic nerve, a phenotype influenced by HCMV infection of the eye ^61,62. Hippo signaling plays well-established roles in eye development ⁶³. Collectively, these results implicate HCMV-induced TEAD1 binding loss in both hearing and eyesight growth defects at multiple specific genomic loci, providing mechanistic insights into the well-established roles of both HCMV and Hippo signaling in these disorders.

Discussion

HCMV is a leading cause of birth defects in developing fetuses, including ear and eye developmental disorders. Upon infection, HCMV alters the expression of thousands of human genes. A better understanding of the mechanisms underlying these changes is critical for understanding the processes underlying HCMV-associated diseases. In this study, we used several complementary functional genomics approaches along with biochemical validation to show that HCMV employs multiple interrelated molecular mechanisms that diminish the activity of the human TEAD1 transcription factor: closing of chromatin in regions normally bound by TEAD1 (with concomitant reduction of H3K27ac levels and chromatin looping interactions), lowering of TEAD1 protein expression levels, preferential exclusion of TEAD1 exon 6, and lowering of protein expression of the TEAD1 co-activator YAP1. This widespread reduction in TEAD1 activity results in substantial alteration to the expression of many human genes, including targets of Hippo and other key developmental signaling pathways.

The human genome encodes four TEAD transcription factors (TEAD1-4), all of which are mediators of the Hippo signaling pathway. The Hippo pathway regulates key aspects of development, including control of organ size, stem cell identity, and lineage specification ^17,18. The activity of TEAD transcription factors (TFs) is regulated by the availability of their co-factors YAP1, TAZ (WWTR1), and members of the vestigial-like protein family (VGLL1-4). While YAP1 and TAZ function as transcriptional co-activators of TEAD TFs, VGLL proteins function as co-repressors ^17,64,65. The activity of TEAD1 has previously been shown to be affected by differential use of exon 6, with exclusion of exon 6 reducing activity due to an impaired ability to interact with YAP1 ³⁹. The YAP1-TEAD complex exerts its activity by binding to enhancers of genes involved in extracellular matrix organization, actin cytoskeleton organization, and cell adhesion ^64,66, all of which are pathways for which we see significantly altered gene expression levels subsequent to HCMV infection (Figure 5B).

A previous report has shown that HCMV infection inhibits the proliferation and invasion of extravillous cytotrophoblasts (EVT) through mRNA downregulation of core Hippo pathway components, including downregulation of YAP1/TAZ and all four TEAD TFs ⁶⁷. Here, we observe a 3-fold depletion in TEAD1 protein levels and a substantial reduction in YAP1 protein levels, in addition to other mechanisms impacting TEAD1 activity. Additionally, the expression of key TEAD1 targets, including Connective Tissue Growth Factor (CTGF), Fibronectin (FN), α-smooth muscle actin (ACTA2) and Collagen 1a (COL1A) have all previously been found to be downregulated with HCMV infection of trabecular meshwork cells, leading to an increase in intra-ocular pressure, a leading cause of glaucoma ⁶⁸. Although the exact mechanisms seem to differ depending on the cell type, HCMV infection of permissible cells clearly perturbs the Hippo pathway in multiple cellular contexts by diminishing TEAD TF activity.

Recent studies have shown that other viruses perturb the Hippo pathway as part of their replicative cycle ⁶⁹. For example, oncogenic viruses, such as Kaposi’s sarcoma-associated herpesvirus (KSHV) and Human papillomavirus (HPV) promote dephosphorylation and nuclear translocation of YAP1 ⁷⁰. In the case of KSHV, the viral GPCR (vGPCR; ORF74) has been shown to inhibit the Hippo kinases LATS1/2, thereby activating YAP1/TAZ ⁷⁰. Likewise, HPV16 oncoprotein E6 induces cell proliferation by dephosphorylating a serine residue of YAP1 (S379) and preventing its degradation by the proteasome complex ⁷¹. More recently, it has been discovered that Epstein-Barr virus subverts the YAP/TAZ pathway in lytic reactivation in epithelial cells ^72,73.

The results of this study represent one of likely many examples of a single virus targeting the same host molecule through multiple mechanisms. For example, Human Papillomavirus targets TP53 through independent mechanisms involving both the E6 and E7 proteins ⁷⁴. Likewise, poxviruses target caspase activity through multiple mechanisms ⁷⁵. Given the vast amount of evolutionary time that viruses have at their disposal due to their highly elevated mutation rates, and the reliance of viruses on host-encoded pathways, it is likely that many other viruses also target a single molecule through multiple mechanisms.

Several previous studies have examined changes to human gene expression induced by HCMV infection in a variety of contexts (reviewed in ⁷⁶). To date, only two studies have employed more than one type of genome-scale measurement to compare HCMV infected to uninfected cells. In the first study, multiple single cell approaches were employed (including scATAC-seq and CITE-seq) to compare NK cells with or without infection, revealing a possible role for AP-1 family TFs ⁷⁷. In the second study, integration of Pol2, H3K27ac, and H3K27me3 ChIP-seq data with ATAC-seq data revealed epigenetic reprogramming of the virus and host genomes in myeloid progenitor cells ⁷⁸. Interestingly, this study found only limited changes to chromatin accessibility between infected and uninfected cells, likely due to differences in cell type (Kasumi-3 vs. HS-68 and ARPE-19 in our study), time point (24 vs. 48 hours), infection strength (MOI of 1 with <50% infection vs MOI of 5 and 10 with >90% infection), and replication efficiency (∼30 vs. ∼250 viral genomes per cell). To our knowledge, our study is the first to employ a systematic approach to identify important host transcriptional regulators through unbiased analysis of virus-altered chromatin accessibility regions, with subsequent ChIP-seq-based validation. Future studies will be needed to systematically identify important host-encoded regulators in the context of other cell types and other viruses.

HCMV infects up to 86% of the world’s population ⁷⁹. For most, this infection is largely benign. However, infection can cause complications in newborn infants and the immunocompromised. Our comprehensive analyses reveal that HCMV targets the Hippo signaling pathway through TEAD1, setting off a cascade of effects on other key developmental pathways. Importantly, these TEAD1 binding loss events are highly enriched for genetic variants associated with eye and ear phenotypes, providing possible mechanistic insights into the well-established role of HCMV infection in eye and ear disorders. Collectively, the results of our study offers new avenues to investigate HCMV-dependent mechanisms in healthy and disease states.

Methods

Experimental design, cell culture, and viral infections

Experiments were performed in human foreskin fibroblasts (HS68 cells). ATAC-seq experiments were also performed in human retinal epithelial cells (ARPE-19 cells). HCMV infections were performed using the TB40/E clinical isolate. Fibroblasts were infected at a multiplicity of infection (MOI) of 5 and cells were harvested 48 hours post infection for RNA-seq, ATAC-seq, ChIP-seq, and HiChIP (Figure 1). Retinal epithelial cells were infected at an MOI of 10. The rate of infection was monitored by both FACS and microscopy.

Gene expression (RNA-seq)

Total RNA was isolated from 5×10⁶ cells, either uninfected or infected (48 hrs), using the RNeasy RNA Isolation kit (Qiagen# 74104) as per instructions provided by the manufacturer. A total of 500 ng of RNA per sample was used as input material for ribosomal-RNA depleted RNA-sequencing. First, ribosomal RNA (rRNA) was removed, and rRNA free residue was cleaned up by ethanol precipitation. Subsequently, sequencing libraries were generated using the rRNA-depleted RNA by the Directional RNA Library Prep Kit. Briefly, fragmentation was carried out using divalent cations under elevated temperature in first strand synthesis reaction buffer (5X). First strand cDNA was synthesized using random hexamer primer and M-MuLV Reverse Transcriptase (RNaseH-). Second strand cDNA synthesis was subsequently performed using DNA Polymerase I and RNase H. In the reaction buffer, dNTPs with dTTP were replaced by dUTP. Remaining overhangs were converted into blunt ends via exonuclease/polymerase activities. After adenylation of 3’ ends of DNA fragments, adaptors with hairpin loop structure were ligated to prepare for hybridization. In order to select cDNA fragments of preferentially 150∼200 bp in length, the library fragments were purified with the AMPure XP system (Beckman Coulter, Beverly, USA). Then 3 µL USER Enzyme (NEB, Ipswich, MA, USA) was used with size-selected, adaptor-ligated cDNA at 37 °C for 15 minutes followed by 5 minutes at 95 °C before PCR. PCR was performed with Phusion High-Fidelity DNA polymerase, Universal PCR primers, and Index Primer. Finally, products were purified (AMPure XP system) and library quality was assessed on an Agilent Bioanalyzer 2100 system (Agilent Technologies, Inc; Santa Clara, CA, USA). The libraries were sequenced on an Illumina NovaSeq 6000 at Novogene (paired-end, 150 bp read length).

RNA sequencing data were processed using the nf-core/rnaseq pipeline (version 3.8.1) ^80–82. Initial quality control of the raw sequencing data was performed using FastQC (version 0.11.9) ⁸³. Low-quality bases and adapter sequences were then trimmed and filtered from the reads using Cutadapt (version 3.4) ⁸⁴ and Trim Galore (version 0.6.7) ⁸⁵, respectively. Ribosomal RNA sequences were subsequently removed from the aligned data using SortMeRNA (version 4.3.4) ⁸⁶ to eliminate potential contamination from non-target RNA species. STAR aligner (version 2.7.10a) ⁸⁷ was used to align the trimmed reads to a custom reference genome of hg19 and Human herpesvirus 5 strain TB40/E clone TB40-BAC4 sequence (GenBank: EF999921.1). The resulting alignments were sorted and indexed using SAMtools (version 1.15.1). To further evaluate the quality of the RNA sequencing data, several tools were employed, including RseqQC (version 3.0.1) ⁸⁸, Qualimap (version 2.2.2-dev) ⁸⁹, dupRadar (version 1.18.0) ⁹⁰, and Preseq (version 3.1.1)⁹¹. These tools provide comprehensive assessments of various quality metrics, such as read distribution, GC content, duplication rates, and library complexity, ensuring reliable data for downstream analysis. For transcript quantification, Salmon (version 1.5.2) ⁹² was utilized to estimate the expression levels of transcripts. Finally, differential gene expression analysis was performed using DESeq2 (version 1.30.1) ⁹³. Genes were considered differentially expressed if they had a 2-fold change and an adjusted p-value threshold of less than 0.01. Pathway enrichment analyses for the differentially expressed genes were performed using Enrichr ^94,95. An adjusted p-value threshold of 0.05 was used to generate the pathway enrichment figure (Figure 5B).

Alternative splicing analysis

Exon-exon and exon-intron spanning reads were identified with the software AltAnalyze (version 2.1.4) using the Ensembl-72 human database, along with splicing event calculation and annotation with the MultiPath-PSI algorithm (see: http://altanalyze.readthedocs.io/en/latest/Algorithms for algorithm details and benchmarking). To detect high-confidence splicing changes between HCMV infected cells and uninfected cells, splicing events that were not detected in all replicates in each biological group were excluded from the differential splicing analysis. Significant splicing changes were defined as splicing events with a change in percent spliced in (PSI) value between the two groups > 10% (ΔPSIL> |L0.1|), with p-value <0.01.

Chromatin accessibility (ATAC-seq)

Omni-ATAC-Seq was performed in replicates as previously published ²⁵. Briefly, HCMV infected and uninfected control cells were harvested by trypsinization for 5 minutes. Approximately 70,000 cells were transferred to a microfuge tube and washed once with PBS. The cell pellet was resuspended in 50 µL of ice-cold lysis buffer (10 mM Tris-HCl [pH 7.5], 10 mM NaCl, 3 mM MgCl2, 0.1% NP-40, 0.1% Tween-20, and 0.01% Digitonin) and incubated on ice for 3 minutes, followed by centrifugation at 500g for 5 minutes at 4 °C. For the transposition reaction, nuclei were resuspended in 50 µL of Nextera transposition reaction mix consisting of 25 µL 2× TD Buffer, 2.5 µL Nextera Tn5 transposase (Illumina # FC-121-1030), and 22.5 µL of nuclease free water. The reaction mixture was mixed 5-10 times by gentle pipetting and incubated at 37 °C for 45 minutes. The transposed DNA was then purified by Qiagen MinElute kit (Qiagen # 28004) and eluted in 11 µL of elution buffer. 1 µL of the eluted DNA was used to assess quality on an Agilent Tapestation 4510. For library preparation, 10 µL of tagmented DNA was PCR amplified for 14 cycles in a 50 µL reaction volume using NEBNext High-Fidelity 2X PCR Master Mix (NEB # M0541S) and Nextera primers. The amplified library was purified and size-selected using AMPure XP beads (Beckman # 63880) in a two-step protocol. First, DNA fragments >1000 bp were removed from the PCR reaction mixture by adding 22.5 µL (0.5× volume) of AMPure XP beads followed by incubation at room temperature (RT) for 10 minutes. The beads were magnetized, and the supernatant was transferred to a new tube. In the second size selection step, 58.5 µL (1.2× volume) of AMPure XP beads were added to the supernatant, mixed by pipetting, and incubated for 10 minutes at room temperature. Beads were washed twice with 80% ethanol, air-dried and DNA eluted in 30 µL of 10 mM Tris-HCl [pH 8.0]. Library DNA concentrations and quality were assessed by a Qubit dsDNA HS assay kit (Q32854) and Agilent 4510 TapeStation system, respectively. Libraries were sequenced on an Illumina NovaSeq 6000 at the Cincinnati Children’s Hospital Medical Center Genomics Sequencing Facility (paired-end, 150 bp read length).

ATAC-seq data were processed and aligned to the hg19 genome using the ENCODE ATAC-seq pipeline (V2.0.0) ^26,27,29,96. Peaks were called within the pipeline using MACS2 ⁹⁷. Differential chromatin accessibility analysis was performed using DiffBind 3.6.5 ^30,98 in R 4.2.1 ⁹⁹. The “conservative overlap” peak sets generated from the individual replicates by the ENCODE ATAC-seq pipeline were used for differential analysis. Peaks were considered to be differentially accessible if the FDR was less than 0.01 and the fold change was 2-fold or greater. A modified version of HOMER ¹⁰⁰ using a log base 2 likelihood scoring system was used to calculate motif enrichment statistics for a large library of human position weight matrix (PWM) binding site models contained in build 2.0 of the CisBP database ¹⁰¹. DeepTools (v2.0.0) ¹⁰² was used to generate heatmaps of signal tracks across differentially accessible chromatin for each set of comparisons.

Chromatin immunoprecipitation sequencing (ChIP-seq)

2 x 10⁷ cells were seeded in 150 mm culture dishes and were either infected with HCMV at an MOI of 5 or kept uninfected. After 48 hours, the medium was removed, and cells were harvested by trypsinization followed by two washes with PBS. Cells were crosslinked with 10 mL of formaldehyde cross-linking solution (1x PBS, and 1% formaldehyde) for 10 minutes at RT. Crosslinking reaction was quenched by adding 2.5 M Glycine to a final concentration of 125 mM and incubated for 5 minutes at RT. Cells were washed twice with ice-cold 1x PBS. Crosslinked cells were then transferred to a microfuge tube and frozen at −80 °C until further use.

For chromatin preparation, crosslinked cells were resuspended in L1 buffer (50 mM Hepes-KOH [pH 8.0], 140 mM NaCl, 1 mM EDTA, 10% glycerol, 0.5% NP-40, 0.25% Triton X-100, and 1x protease inhibitors) and incubated at 4 °C on a rotator for 10 minutes. Nuclei were pelleted by centrifugation and resuspended in L2 buffer (10 mM Tris-HCl [pH 8.0], 200 mM NaCl, 0.5 mM EGTA, and 1x protease inhibitors) with rotation for 10 minutes at RT. Isolated nuclei were then resuspended in sonication buffer solution (10 mM Tris-HCl, 1 mM EDTA, and 0.1% SDS). To obtain chromatin fragments of 200-500 bp length, nuclei were sonicated using an S220 ultrasonicator (Covaris, LLC., Woburn, MA, USA) at 10% duty cycle, 175 peak power, 200 burst/cycle for 7 minutes at 4 °C. A portion of the sonicated chromatin was run on an agarose gel to verify fragment sizes. The chromatin solution was centrifuged to pellet the debris and the supernatant was collected in a fresh microfuge tube. Since the sonicated chromatin solution is devoid of any detergents or salts, the following were added at indicated final concentrations: Triton X-100 (1%), Sodium deoxycholate (0.1%), Glycerol (5%), NaCl (150 Mm), and Protease inhibitor cocktail (1x). Chromatin was pre-cleared with Protein G Dynabeads for 45 minutes at 4 °C on a rotator. Chromatin immunoprecipitations were performed in an SX-8X IP-Star Compact automation system (Diagenode. LLC.,Denville, NJ CA, USA) with 200 µL of pre-cleared chromatin (approximately 3 million cells), 21 µL of Protein A or G Dynabeads (Thermo Fisher Scientific, Waltham, MA USA), and 5 µg of antibody. The following antibodies were used for ChIP: CTCF (ActiveMotif # 61311, lot 17118005), TEAD1 (ActiveMotif # 61643, lot 34614001) and H3K27ac (Abcam # ab4729, lot GR3357415-1). Immunoprecipitation of chromatin was carried-out for 8 hours, after which the Dynabeads were sequentially washed with Wash Buffer 1 (10 mM Tris-HCl [pH 8.0], 150 mM NaCl, 1 mM EDTA, 0.1% SDS, 0.1% sodium deoxycholate, and 1% Triton X-100), Wash Buffer 2 (10 mM Tris-HCl [pH 8.0], 250 mM NaCl, 1 mM EDTA, 0.1% SDS, 0.1% sodium deoxycholate, and 1% Triton X-100), Wash Buffer 3 (50 mM Tris-HCl [pH 8.0], 2 mM EDTA, and 0.2% N-Lauroylsarcosine sodium salt), and Wash Buffer 4 (TE + 0.2% Triton X-100) for 5 minutes each. For library preparation, immunoprecipitated chromatin was eluted in elution buffer (1x TE, 250 mM NaCl and 0.3% SDS). Chromatin proteins and RNA were digested with Proteinase-K and RNase A respectively and DNA purified with the Qiagen MinElute kit. A ChIP-DNA library was constructed using the NEBNext Ultra II DNA library preparation kit (E7645S) as per manufacturer’s instructions and purified with 0.6x volume of AMPure XP beads. Purified library DNA quality and quantity were assessed by an Agilent TapeStation 4150 and sequenced using an Illumina NovaSeq 6000 at the Cincinnati Children’s Hospital Medical Center Genomics Sequencing Facility (Single end, 100 bp read length).

ChIP-seq data were processed and aligned to the hg19 genome using the ENCODE ChIP-seq pipeline (V2.0.0) ^26,27,29,103. Peaks were called within the pipeline using MACS2 ⁹⁷. Differential peak analysis was performed using DiffBind 3.6.5 ^30,98 in R 4.2.1 ⁹⁹. The “conservative overlap” peak sets generated from the individual replicates by the ENCODE ChIP-seq pipeline were used for differential analysis. Peaks were considered differentially enriched if the FDR was less than 0.01 and the fold change was 2-fold or greater.

Chromatin looping interactions (HiChIP)

HiChIP libraries were prepared following the protocol from Mumbach et al ³³. Experiments were performed using biological triplicates in both uninfected and infected cells. In brief, 10 million cells were cross-linked with 1% formaldehyde, followed by a cell lysis with Hi-C lysis buffer (10 mM Tris-HCl [pH 8.0], 10 mM NaCl, 0.2% NP-40, 1x Roche protease inhibitors) for 30 minutes at 4 °C. Following cell lysis, 375 U of MboI (NEB, R0147) was used to digest, in-situ, the cross-linked chromatin for 2 hours at 37 °C. After filling-in the DNA ends by biotin-dTAP (Thermo # 19524016), dCTP, dGTP, and dTTP with 5 U/µL DNA polymerase (NEB # M0210), the DNA was ligated by T4 DNA ligase (NEB # M0202) at 4 °C overnight. The ligated DNA was sonicated using a Covaris S220 platform (Covaris, LLC., Woburn, MA, USA) in nuclear lysis buffer (50 mM Tris-HCl [pH 7.5], 10 mM NaCl, 0.2% NP-40, 1x Roche protease inhibitors) at 4

°C. The fragmented DNA was 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 the samples were precleared with Protein A Dynabeads at 4 °C for 1 hour, followed by immunoprecipitation with 8 µL H3K27Ac antibody (Abcam # ab4729, lot GR3374555-1) at 4 °C overnight. DNA-protein complexes were captured by Protein A beads with rotation at 4 °C for 2 hours. The protein A beads were then washed sequentially with low salt, high salt, and LiCl wash buffer, followed by elution twice with 100 µL freshly prepared DNA elution buffer at 37 °C (50 mM NaHCO₃, 1% SDS). The eluted chromatin was reverse cross-linked and purified by a PCR purification Kit (Qiagen). 5 μL of Streptavidin C-1 beads (Thermo Fisher # 65001) was used to capture the Biotin Datp labeled DNA. The captured DNA was transposed with 2.5 μL Tn5 transposase (Illumina). The beads were then sequentially washed by Tween wash buffer (5 mM Tris-HCl [pH 7.5], 0.5 mM EDTA, 1 M NaCl, 0.05% Tween-20), 50 mM EDTA, and 10 mM Tris. After washing, PCR was performed by re-suspending the beads with 23 μL H₂O, 25 μL 2X Phusion HF (New England Biosciences), 1 μL Nextera forward primer (Ad1_noMX), and 1 μL Nextera reverse primer (Ad2.X) at 12.5 μM. The PCR was run at 8 cycles of (1) 72 °C for 5 minutes, (2) 98 °C for 1 minute, (3) 98 °C for 15 seconds, (4) 63 °C for 30 seconds, followed by an extension at 72 °C for 1 minute. The post-PCR size selection was performed by two-sided selection with AMPure XP beads to capture the fragments between 300 and 700 bp. Samples were then sequenced on an Illumina NovaSeq 6000 at Novogene (paired-end, 150 bp read length).

HiC-Pro (version 2.11.4) was used to align and filter read pairs and identify the contact map ³². MboI restriction sites and default parameters were used to align the reads to the hg19 genome. After alignment, the read pairs were filtered to remove those that mapped to multiple locations, were not in a valid orientation, or were duplicated. Quality control metrics were generated at each step. HiC-Pro also reports the number of trans pairs, short-range cis pairs and long-range cis pairs, as well as the number of pairs in each orientation (Supplemental Dataset 2). The overlap conservative peaks generated by the ENCODE pipeline for the H3K27ac ChIP-seq data were used as input for calling loops using FitHiChIP (version 11.0) ¹⁰⁴. The parameters for FitHiChIP were the following: interaction type of “peak to all”; bin size of 10 kb; lower loop distance of 20 kb; upper loop distance of 20 mb; background model of loose; bias correction of coverage bias regression; merge filtering enabled; and FDR of 0.01 (Supplemental Dataset 2). Peaks were called from the HiChIP data using HiChIP-Peaks (version 0.1.2) ¹⁰⁵ using the MboI restriction sites and an FDR of 0.01. The peaks were compared to the H3K27ac ChIP-seq peaks using bedtools (version 2.30.0) ¹⁰⁶ (Supplemental Dataset 2). We observed strong agreement between experimental replicates (Supplemental Figure S5), so we pooled the reads for the replicates and analyzed the resulting data as described above (Supplemental Dataset 2). We identified differential looping events using the DiffAnalysisHiChIP.r script included in the FitHiChIP software with a default bcv parameter value of 0.4 and cutoffs of >2-fold change and FDR < 0.01 (Supplemental Dataset 2). We used the merged loops from the combined replicates as input. Replicate comparisons (Supplemental Figure 5) were performed using the multiBamSummary and plotCorrelation programs from the deepTools software package ¹⁰². The parameters used for multiBamSummary were –-binSize 10000 and –-distanceBetweenBins 0. The parameters used for plotCorrelation were –-corMethod pearson, –-skipZeros, and –-removeOutliers.

Assessment of protein expression levels (Western blots)

Uninfected and HCMV infected cells were harvested, in triplicates, by trypsinization and washed twice with PBS. The cell pellet was then resuspended in 1x RIPA buffer supplemented with Halt Protease inhibitor cocktail (Sigma) and kept on ice for 30 minutes. Cell lysates were cleared by centrifugation at 12,000 RPM for 10 minutes and supernatant was collected in a new tube. Protein concentration was measured by the BCA method and 30 µg of protein loaded on 4-12% Bis-tris novex gel (Invitrogen). After electrophoresis, proteins were transferred on to a Nylon membrane using Invitrogen’s semi-dry transfer method. Proteins were blocked with Intercept blocking buffer (LI-COR) and membrane incubated with HCMV IE1/2 (Millipore # MAB810R), TEAD1 (Cell Signaling # 12292), YAP1 (Cell Signaling # 14074), pYAP1 (Phospho-YAP1; Ser127; Cell Signaling # 4911), THBS1 (Novus Biologicals # NB100-2059), CCN1 (Novus Biologicals # NB100-356), H3K27ac (Abcam # ab4729), β-Actin (Cell Signaling # 8457) or GAPDH (Invitrogen # ma5-15738) antibodies at 4 °C for overnight. Following antibody binding, membrane was washed 3 times with 1x TBST and incubated either with IRDye® 800CW or IRDye® 680RD conjugated secondary antibodies for 45 minutes at RT. For target protein quantification against a loading control (GAPDH and β-Actin), secondary antibodies against both target protein and loading control were added to the blot. The membrane was washed 3 times with 1x TBST and imaged on an Odyssey Dlx system (LI-COR, Biosciences LLC., Lincoln, NE, U.S.A). Scanned images were then used for target protein quantification using Emperia Studio software (LI-COR, Biosciences LLC., Lincoln, NE, U.S.A).

Assessment of YAP1 localization

Uninfected and HCMV infected cells were harvested by trypsinization and washed twice with PBS. Cytosolic and nuclear protein extracts were prepared using Pierce NE-PER Nuclear and Cytoplasmic Extraction Reagent Kit (Thermo Scientific # 78833), as per manufacturer’s instructions. Protein concentration was measured by the BCA method and 15 µg of protein loaded on 4-12% Bis-tris novex gel (Thermo Fisher Scientific, Waltham, MA USA). Western blotting and imaging were performed as given above with the following antibodies: YAP1 (Cell Signaling # 14074), pYAP1 (Phospho-YAP1; Ser127; Cell Signaling # 4911), Histone H3 (Abcam # ab1791) and GAPDH (Invitrogen # ma5-15738).

Validation of TEAD1 splicing (RT-PCR)

Total RNA was extracted from uninfected and HCMV infected fibroblasts 48 hours post infection (hpi) using the RNA mini kit from Qiagen as per manufacturer’s instructions. 2.5 µg of total RNA was reverse transcribed using SuperScript™ IV VILO™ Master Mix (Thermo Fisher Scientific, Waltham, MA USA) and 3 µL of cDNA was used in PCR to amplify the regions immediately flanking TEAD1 exon-6 using forward and reverse primers. Forward primer: 5’-ATCTCGTGATTTTCATTCCAAGC. Reverse Primer: 5’-TGAGGACATGGCCGCCATGTGC.

Data Visualization

Visualization tracks for each dataset were created for the UCSC Genome Browser ¹⁰⁷. Signal tracks (in bigWig format) were created using the bamCoverage program from the deepTools software package ¹⁰² with the parameters –-normalizeUsing BPM and –-binSize 10. For the RNA-seq signal tracks, the command was run twice in order to maintain strandedness. For the forward strand track, the parameter –-filterRNAstrand forward was included. For the reverse strand track, the parameters –-filterRNAstrand reverse and –-scaleFactor –1 were included. The .hic files for the HiChIP data were created using the hicpro2juicebox.sh script from HiC-Pro (version 2.11.4) ³². Juicer (version 1.22.01) ¹⁰⁸ was used for KR normalization.

TEAD1 binding event phenotype enrichment analysis

To determine the significance of the overlap between TEAD1 binding events and GWAS-derived disease– and phenotype-associated genetic variants, we generated a custom GWAS catalogue for each phenotype in an ancestry specific manner. To this end, we downloaded the Genome Wide Association Studies Catalogue (https://www.ebi.ac.uk/gwas/) v1.0.2, as queried on June 20^th 2019. Independent risk loci for each disease/phenotype were identified based on linkage disequilibrium (LD) pruning (r2 < 0.2). Risk loci across these independent genetic risk variants were identified by linkage disequilibrium expansion (r2>0.8) based on 1000 Genomes Data using PLINK (v.1.90b). This created a list of disease risk loci, along with the corresponding genetic variants within the LD block. This list of variants was then used for RELI analyses as previously published ⁴¹. Phenotypes were considered “TEAD1 loss specific” if: (1) they had 3 or more overlaps with TEAD1 binding event losses; (2) this overlap was significant according to RELI (corrected p<0.01); (3) this overlap was not significant for unchanged TEAD1 peaks (p>0.01), and (4) TEAD1 binding loss events were at least 2-fold enriched according to RELI. The three phenotypes meetings these criteria are presented in the Results, with further data included in Supplemental Dataset 3.

Data access

All raw and processed sequencing data generated in this study have been submitted to the NCBI Gene Expression Omnibus (GEO; https://www.ncbi.nlm.nih.gov/geo/) and will be made available upon peer review.

Acknowledgements

This research was funded by National Institutes of Health (NIH) R01 HG010730, R01 GM055479, U01 AI130830 to M.T.W.; R01 AR073228, R01 NS099068, and R01 AI024717 to M.T.W. and L.C.K.; R01 AI148276, U01 HG011172, and U19 AI070235, and P30 AR070549 to L.C.K.; T32 ES007250 to A.V.H.; R21 DE026267, R01 AI121028 to W.E.M; F32 AI172329 and T32 AI007245 to L.A.M-N; R01 AI164709 to B.E.G; US Department of Veterans Affairs, I01 BX003850 to K.M.K; US Department of Veterans Affairs, I01 BX001834 and I01 BX006254 to J.B.H. We would like to thank Jaimie Miser for her assistance with culture work, Phillip Dexheimer for his help with the virus expression database (https://vexd.cchmc.org/), the Cincinnati Children’s Hospital Medical Center Genomics Sequencing Facility (RRID SCR_022630) for high-throughput sequencing, the support of the Informatics Shared Facility in Information Services for Research (IS4R) at Cincinnati Children’s Hospital Medical Center (RRID:SCR_022622), and Kevin Ernst for help with database management and computational support.

Additional information

Author contributions

Project conception: M.T.W, L.C.K, J.B.H, and K.S; W.E.M provided TB40/E stocks, M.J.B, W.E.M, J.W., and K.S performed cell cultures and initiated viral infections; K.S conducted most experiments; M.H. and C.F. provided experimental support; S.P performed all initial informatics and QC analysis; K.S, S.P, L.A.M.-N., L.C.K, and M.T.W analyzed and interpreted data; S.P, A.V, L.E.E, O.D, and M.T.W designed and wrote informatics pipelines used for data analysis; M.M.L.L and B.Z performed the H3K27ac HiChIP experiment; L.E.E performed HiChIP analysis; O.A.D quantified and pooled the sequencing libraries for ATAC-seq, ChIP-seq and HiChIP-seq; A.C performed alternative splicing analysis; S.R assisted with qRT-PCRs; K.K provided critical insights at the early stages of the project and helped with data analysis; K.S, L.C.K, and M.T.W wrote the first draft of the manuscript; K.S, S.P, L.E.E, L.C.K, and M.T.W finalized the manuscript with comments from all other authors.

Declaration of interests

All authors declare no conflicts of interests related to this manuscript.

Supplemental Figures

Supplemental Figure S1.

Supplemental Figure S2.

Supplemental Figure S3.

Supplemental Figure S4.

Supplemental Figure S5.