An evolutionarily conserved constellation of functional cis-elements programs the virus-responsive fate of the human (epi)genome

Abstract

Human health depends on perplexing defensive cellular responses against microbial pathogens like Viruses. Despite the major effort undertaken, the (epi)genomic mechanisms that human cells utilize to tailor defensive gene expression programs against microbial attacks have remained inadequately understood, mainly due to a significant lack of recording of the in vivo functional cis-regulatory modules (CRMs) of the human genome. Here, we introduce the virus-responsive fate of the human (epi)genome as characterized in naïve and infected cells by functional genomics, computational biology, DNA evolution, and DNA Grammar and Syntax investigations. We discovered that multitudes of novel functional virus-responsive CRMs (vrCRMs) compose typical enhancers (tEs), super-enhancers (SEs), repetitive-DNA enhancers (rDEs), and stand-alone functional genomic stretches that grant human cells regulatory underpinnings for layering basal immunity and eliminating illogical/harmful defensive responses under homeostasis, yet stimulating virus-responsive genes and transposable elements (TEs) upon infection. Moreover, extensive epigenomic reprogramming of previously unknown SE landscapes marks the transition from naïve to antiviral human cell states and involves the functions of the antimicrobial transcription factors (TFs), including interferon response factor 3 (IRF3) and nuclear factor-κB (NF-κB), as well as coactivators and transcriptional apparatus, along with intensive modifications/alterations in histone marks and chromatin accessibility. Considering the polyphyletic evolutionary fingerprints of the composite DNA sequences of the vrCRMs assessed by TFs-STARR-seq, ranging from the animal to microbial kingdoms, the conserved features of antimicrobial TFs and chromatin complexes, and their pluripotent stimulus-induced activation, these findings shed light on how mammalian (epi)genomes evolved their functions to interpret the exogenous stress inflicted and program defensive transcriptional responses against microbial agents. Crucially, many known human short variants, e.g. single-nucleotide polymorphisms (SNPs), insertions, deletions etc., and quantitative trait loci (QTLs) linked to autoimmune diseases, such as multiple sclerosis (MS), systemic lupus erythematosus (SLE), Crohn's disease (CD) etc., were mapped within or vastly proximal (±2.5 kb) to the novel in vivo functional SEs and vrCRMs discovered, thus underscoring the impact of their (mal)functions on human physiology and disease development. Hence, we delved into the virus-responsive fate of the human (epi)genome and illuminated its architecture, function, evolutionary origins, and its significance for cellular homeostasis. These results allow us to chart the "Human hyper-Atlas of virus-infection", an integrated "molecular in silico" encyclopedia situated in the UCSC Genome Browser that benefits our mechanistic understanding of human infectious/(auto)immune diseases development and can facilitate the generation of in vivo preclinical animal models, drug design, and evolution of therapeutic applications.

Graphical Abstract

Figure 1.

The human epigenome grants cells with regulatory underpinnings fundamental for virus-stimulated gene expression programs establishment. (A) Delineation of common and cell-type-specific virus-stimulated gene expression programs in human epithelial cells (HL) and B-lymphocytes (NM); Left upper and lower panels: RNA-seq heatmaps depict the transcriptional changes acquired in NM and HL [SVI; Naïve (mock-infected) (0 h) and virus-infected cells (3 and 6 h)] and revealed 678 vrDEGs in NM (490 upregulated and 188 downregulated) and 581 in HL (466 upregulated and 115 downregulated). Central upper panels: Intersection of the NM and HL vruDEGs records common (n = 167) and cell-type-specific spectrums (n = 299 in HL; n = 323 in NM). Central lower panel: GOs (clusterProfiler), sharply highlight the specificity of the 167 common vruDEGs for antiviral/defensive cellular processes. Right panels: Heatmaps of gene expression depicting the kinetics of transcriptional upregulation of striking examples of common C/I-III vruDEGs (e.g. secreted effectors, PRRs, and TFs). Additional data are described in Supplementary Figs S1–S3 and Supplementary Table S1. (B) “ON” and “OFF” transcriptional states profiling of common C/I-III vruDEGs; Box plots depict the distribution of the average log₂ normalized expression values of each vruDEG per time point and the average of log₂ (FC) as computed in biological replicates. The boxplots illustrate the interquartile range, their upper and lower limits represent the 75th and 25th percentile, respectively, and the whiskers extend to the minimum (lowest) and maximum (highest) non-outlier values. The central horizontal line depicts the median of the values, and the outliers are shown as dark black dots. The red dot depicts the value of the mean. Two-sided Wilcoxon rank-sum test was applied to the total common vruDEGs of each cluster (C/I: 43 in NM and 47 in HL; C/II: 63 in NM and 64 in HL; C/III: 61 in NM and 56 in HL) (Supplementary Table S1). The results depict statistically significant differences in basal levels of expression between C/I and C/II in NM (P-value = 2.81e⁻⁷), and in HL (P-value = 9.24e⁻⁷); C/I and C/III in NM (P-value = 7.22e⁻¹¹), and in HL (P-value = 2.94e⁻¹⁰); C/II and C/III in NM (P-value = 0.028) and in HL (P-value = 0.0259). These results unveil significant heterogeneity in basal levels of vruDEGs’ expression and underscore that higher levels predominantly correlate with medium or minimal transcriptional induction upon virus-infection. Additional data for common- and cell-type-specific clusters of vruDEGs are described in Supplementary Figs S1–S3 and Supplementary Table S1. (C and D) Massive whole-epigenome- and regional-scale epigenomics studies in naïve and virus-infected human cells; Shown are chromatin states accessibility profiling by DNaseI-seq, and mapping of histone modifications, antimicrobial TFs, coactivators, and transcriptional apparatus by ChIP-seq. (C) The in vivo reconstitution of CHs in human cells: (i) concentric cycles demonstrate the epigenome-wide, sequential, overlaid intersections of the NGS-peaks obtained from DNaseI-seq, H3K27ac-, MED1-, and RNA pol II-ChIP-seq assays in HL and NM (SVI, 0 and 6 h). In virus-infected cells, IRF3, p65, and CBP incorporation within CHs is also highlighted. Genomic localization assessments mapped common and cell-type-specific vruDEGs in cis proximity (±2 kb) to established CHs. Bar graphs depict the number of CHs-associated vruDEGs and highlight that in naïve cells, C/II and C/III predominantly reside proximal to CHs, while upon virus-infection all clusters become enriched for CHs-associated vruDEGs, in HL and NM, in compatibility with the corresponding transcriptional states, as exchanged in vivo (additional numerical data for CHs are illustrated Supplementary Fig. S9 and Supplementary Table S2). (ii and iii) DNA Grammar assessments highlight an enriched lexicon of TFBMs legible by antimicrobial TFs in total and in IRF3-, p65-, and CBP-ChIPed, “pre-printed” and “newly established” CHs. Genomic localization assessments mapped the “pre-printed” CHs more excessively within promoter regions compared to “newly established” CHs. (D) QRTMs [Aggregation plots (signals’ average) and heatmaps (signals’ coverage)] applied within the TSS-proximal (±2 kb) chromatin microenvironments of residence of C/I-III common vruDEGs in HL (SVI, 0, 3, and 6 h). In naïve cells, C/II and C/III reside within “epigenetically triggered” whereas C/I in epigenetically fingerprinted chromatin microenvironments that lack hallmarks of transcriptional activation in line with their enhanced and constrained basal levels of expression, respectively. Upon virus-infection, C/II and C/III vruDEGs maintain or slightly enhance their “archetype” epigenetic features, which become divergently (co)-occupied by antimicrobial TFs and CBP, whereas C/I gradually acquire “open” chromatin, MED1, H3K27ac, H3K4me3, and RNA pol II accompanied by robust binding of antimicrobial TFs and CBP. QRTMs in NM verified the above epigenetic states (Supplementary Fig. S9). QRTMs for cell-type-specific genes are illustrated in Supplementary Fig. S12. (E) Virus-stimulated epigenomic fingerprints hallmark the transition from naïve-to-antiviral states in human cells and correlate with defensive cellular responses; Comparative GOs highlight tremendous specificity of the virus-stimulated epigenomics signals obtained from the full spectrum of DNaseI- and ChIP-seq assays for antiviral/defensive cellular processes in HL. The same analysis in NM is illustrated in Supplementary Fig. S13.

Figure 2.

Common-Cores of human vruDEGs hallmark the establishment of defensive virus-stimulated gene expression programs in diverse human cell-types infected with distinct Viruses: Viral-stimulated gene expression programs mounting in human cells of epithelium, blood, lung, respiratory tract, and renal/kidney origins, and in human organoids, are comparatively investigated. Computational analyses were applied on the datasets derived from HL, NM, and MRC-5 experiments and publicly available RNA-seq datasets retrieved from Human Calu-3, male lung epithelial adenocarcinoma cells, naïve or infected with Severe Acute Respiratory Syndrome Coronavirus 2 (SARS-CoV-2; 0 and 24 h); Human A549, male lung epithelial cells, naïve or infected with: (i) Influenza A Virus (IAV; 0 and 6 h), (ii) Respiratory Syncytial Virus (RSV; 0 and 24 h), (iii) Dengue Virus 2 (DENV2; 0 and 24 h), and (iv) parainfluenza virus type 3 (HPIV3; 0 and 24 h); Human cell cultures established from upper respiratory tract specimens, naïve, or infected with enterovirus (EV-D68; 0 and 48 h); Human HEK293 female embryonic kidney cells, naïve or infected with Rift Valley Fever Virus (RVFV; 0 and 48 h); Human CD4⁺ T cells naïve or infected with Ebola Virus (EBOV; 0 and 24 h); Human Pluripotent Stem Cells (HPSCs)-derived Airway organoids naïve or infected with SARS-CoV-2 (0 h, 2-day post infection); HPSCs-derived brain organoids naïve or infected with Zika Virus (ZIKV 0 h, 3-day post infection); and Primary Human macrophages isolated from healthy donors naïve or infected with West Nile Virus (WNV; 0 and 24 h). Heatmaps depict the relative expression of the 167 common vruDEGs of NM and HL in the various human cell systems investigated prior to and upon infection with distinct Viruses. The results depict substantial correlation between the gene expression profiles assembled and delineate the mounting of broadly applied profiles and fates of virus-stimulated defensive/antiviral transcriptional responses within human tissues. Additional data are described in Supplementary Figs S3 and S4, and Supplementary Table S1.

Figure 3.

SEs’ “decommission”, “maintenance”, and “genesis” hallmark the “naïve-to-antiviral” exchange of human cell states. (A) Upper panels: Ranking plots from ROSE-analyses applied on MED1-ChIP-seq experiments depict the in vivo reconstitution of SEs and tEs, in naïve and virus-infected epithelial cells (HL) and B-lymphocytes (NM) (SVI, 0 and 6 h). Striking examples of vruDEGs residing in cis proximity to SEs are highlighted. ROSE-analyses applied to H3K27ac-ChIP-seq experiments are illustrated in Supplementary Fig. S15. Bottom panels: QRTMs [Aggregation plots (signals’ average) and heatmaps (signals’ coverage)] illustrate prolonged epigenomic entities corresponding to in vivo reconstituted SEs epigenetically marked by H3K27ac and MED1 in naïve and virus-infected (SVI, 0 and 6 h) NM (left panels) and HL (right panels). (B) High-resolution IGVs snapshots of genomics tracks illustrate MED1- and H3K27ac-ChIP-seq signals intensities and distribution across striking examples of pdSEs’, pmSEs’, and viSEs’ genomic loci of residence, in NM (upper panels) and HL (bottom panels). Horizontal lines and rectangles demarcate SEs genomic coordinates, as mapped by the ROSE algorithm. (C) Upper and lower panels: SEs median lengths [total and individual classes (pdSEs, pmSEs, and viSEs)]; Naïve cells: ∼40 kb in NM and ∼20.5 kb in HL. Infected cells: ∼36 kb in NM and ∼18.6 kb in HL. Median lengths are also plotted for each cell type and state. Boxplots are constructed as detailed in “Materials and methods” section and Fig. 1 legend. (D) Upset plots demonstrate the (co)-occupation of IRF3, p65, CBP, and RNA pol II across SEs genomic coordinates. Naïve cells (SVI, 0 h); CBP: 811 (∼62%) in NM, 505 (∼48%) in HL. RNA pol II: 1193 (∼92%) in NM, 997 SEs (∼96%) in HL; virus-infected cells (SVI, 6 h); CBP: 733 (∼46%) in NM, 537 (∼40%) in HL. RNA pol II: 1442 (∼92%) in NM and 1250 (∼94%) in HL. IRF3: 490 (∼31%) in NM and, 798 (∼60%) in HL. p65: 379 (∼24%) in NM, and 299 (∼22%) in HL. In NM and HL, 26.3% and ∼25.6% of SEs, respectively, are co-occupied by IRF3 and p65. Rectangles indicate the spectrum of “fully-equipped” SEs with the above markers examined (>200 SEs, in both cell types). (E) 3D investigations for SEs in vivo reconstitution proximal to vruDEGs in naïve HeLa cells captured several pairs of vruDEGs and SEs to coinhabit (neighbor, reside proximal, or entirely coincide) within chromosomal domains in naïve HL cells (e.g. chr20: B4GALT5/SE; chr20: CEBPB/SE; chr5: PTGER4/SE). Rectangles indicate SEs. Additional data are described in Supplementary Fig. S22 and Supplementary Table S3.

Figure 4.

SEs genomic localization highlights extensive epigenome reprogramming phenomena that match the gene expression characteristics as steady-stated prior to and exchanged upon virus-infection in HL. (A) Chromosomal Ideograms (e-karyotypes) illustrate SEs in vivo reconstitution in naïve (upper panel) and virus-infected (bottom panel) HL (SVI, 0 and 6 h). Horizontal lines within chromosomes depict gene density. Squares label pdSEs, circles label pmSEs, and arrows/triangles label viSEs (on the right). Boxes (on the left) label vruDEGs expressed at the cell state of examination that reside within ±0.5 Mb from SEs. Striking examples of SE-associated-expressed vruDEGs are highlighted above each human chromosome. Upper and bottom right bar graphs: The percentage of allocation (%) of the 110 and 236 SEs-associated expressed vruDEGs within C/I-III in naïve and virus-infected HL, respectively. The anatomical patterns of SEs assembly match the gene expression programs as established, in vivo (SVI, 0 and 6 h). The same analysis for NM is illustrated in Supplementary Fig. S16. (B) viSEs double-marked by H3K27ac and MED1 in HL; Left panel: GOs depict substantial specificity for negative regulation of viral processes, signal transduction, immune system processes etc. Right panel: QRTMs [Aggregation plots (signals’ average) and heatmaps (signals’ coverage)] illustrate the epigenomics signals in naïve (SVI, 0 h, blue line) and virus-infected cells (SVI, 6 h; red line) corresponding to H3K27ac, MED1, DNaseI, RNA pol II, IRF3, p65, and CBP as distributed across the genomic coordinates. Predominantly, enhanced signals across expanded epigenomics entities that follow a virus-inducible mode of enrichment are captured for the markers/factors examined. The incorporation of antimicrobial TFs, coactivators, and transcriptional apparatus further indicates the functional conductance of viSEs in defensive cellular mechanisms. (C) SEs genomic localization and transcriptional regulation of linearly associated vruDEGs in naïve and infected HL. Left panel: A comprehensive description of the workflow applied. Flowchart showing the steps followed for the assignment of SEs to genes leading to the construction of chromosomal Ideograms (e-karyotypes) in naïve and infected HL. Numerous neighboring, overlapping, or entirely embedded expressed vruDEGs genes of C/I-III clusters were mapped in cis association with at least one SE in naïve and infected HL. Right panel: Boxplots depict the average log₂ normalized expression values of SEs-associated and SEs-far-distant vruDEGs allocated in C/I-III in naïve (SVI, 0 h) and virus-infected cells (SVI, 6 h) as computed in biological replicates. Boxplots are constructed as detailed in “Materials and methods” section and Fig. 1 legend. Two-sided Wilcoxon rank-sum test computed the statistical significance in the divergence of expression levels between the SEs-associated and SEs-far-distant vruDEGs of each cluster [naïve HL: C/I P-value = 0.001, C/II P-value = 4.82e⁻⁹, C/III P-value = 0.0056; virus-infected HL: C/I P-value = 0.29 (nonsignificant), C/II P-value = 0.004, C/III P-value = 0.0058]. Numbers next to boxes highlight the number of SEs-associated or SEs-far-distant expressed vruDEGs, respectively. Significantly, higher diversification in basal expression levels between SEs-associated and SEs-far-distant vruDEGs for C/I-III in naïve cells, is recorded. Importantly, C/II and C/III are enriched in SEs-associated-expressed vruDEGs compared to C/I. In virus-infected cells these diversifications are more balanced. The same statistical analysis validates similar findings in B-lymphocytes (NM), as illustrated in Supplementary Fig. S16.

Figure 5.

Massive-in-parallel in vivo functional authentication by TFs-ChIP-STARR-seq unravels “bona fide” IRF3- and NF-κB-targeted vrCRMs. (A) (i) The principles of TFs-ChIP-STARR-seq methodology: STARR-seq is based on the capability of functional cis-elements to activate the transcription of an intervening reporter-gene (e.g. GFP) when they are located several kbs downstream (3′ prime) of a promoter sequence, thus becoming self-transcribed as 5′-GFP:cis-element:poly (A)-3′ mRNA hybrids. TFs-ChIP-STARR-seq examination platforms are established in naïve and virus-infected HL cells. In brief, upon recombination of NGSed ChIPed DNA fragments in STARR-vectors, the generated multiplexed plasmid libraries are back-transfected in HL cells that are then subjected to virus-infection (SVI, 6 h) or remain uninfected (mock-infected; SVI, 0 h). Total mRNA is isolated and processed to generate STARR-transcripts-specific-cDNA-libraries for NGS (“Materials and methods” section). (ii) Genome-wide correlation analyses of the mapped NGS reads of each ChIP-seq experiment that was utilized to generate the respective input multiplexed-STARR-plasmid library for the IRF3-ChIP-STARR-seq (left) and the p65-ChIP-STARR-seq (right) assays. In these analyses, the human genome was fragmented in silico into 5000 bp bins, and the log1p (natural log + 1) reads per bin are displayed. A strong statistically significant positive correlation is demarcated, in both assays (Pearson correlation coefficient: 5Aii left: R = 0.76, P-value < 2.2e⁻¹⁶; 5Aii right: R = 0.8, P-value < 2.2e⁻¹⁶). (B) Left upper and lower panels: Aggregation plots (signals’ average) and heatmaps (signals’ coverage) illustrate the ChIP-STARR-seq RNA-seq signals identified in naïve (SVI, 0 h, green line) and virus-infected HL (SVI, 6 h, red line). Blue line corresponds to the signals derived from direct NGS of the multiplexed-STARR-plasmid library. In both IRF3- and p65-ChIP-STARR-seq assays, aggressive transcriptional upregulation is detected upon virus-infection, that indicates the function of vrCRMs. Computational processing resulted in the in vivo functional authentication of 1949 IRF3- and 1601 p65-SHAe. Central upper and lower panels: QRTMs meta-profiling of epigenomics signals across endogenous SHAe loci of residence illustrate that during the transition from naïve (SVI, 0 h) to antiviral (SVI, 3 and 6 h) cell states, IRF3- and p65-SHAe exhibit gradual enrichment of H3K27ac labeling, antimicrobial TFs binding, recruitment of coactivators and RNA pol II, yet stability of H3K4me3 labeling; whereas “open” chromatin is gradually enriched across IRF3-SHAe yet pre-established and enriched across p65-SHAe. Right upper and lower panels: Genomic localization assessments (% distribution) map IRF3- (upper bar graph) and p65-SHAe (lower bar graph) within human SEs, tEs, and “safGs”, thus validating their functional fitness/conductance for antiviral cellular processes. (C) (i) and (ii) left panels: IRF3- and p65-SHAe neighbor in cis vruDEGs: Binomial tests followed by Wilcoxon rank sum test (within the first 10 kb) verified statistically significant in cis proximity of IRF3- and p65-SHAe to the vruDEGs (continuous lines), compared to STARR-active-non-SHAe (dashed lines) or randomly sampled loci of equal length tracked from the rest of the human genome (square dotted lines) (P-value < 2.2e⁻¹⁶, for both Ci and Cii). These results show that the virus-responsive fate of the human epigenome harbors in vivo functional vrCRMs embedded within chromatin microenvironments linearly associated with the vruDEGs. Striking examples are illustrated in Fig. 6 and the entire catalog of SHAe-associated vruDEGs is charted in Supplementary Fig. S23B. (i) and (ii) right panels: GOs validate that IRF3- and p65-SHAe are circuited in the cellular mechanism of defensive transcriptional responses; Defense response to virus, negative regulation of viral genome replication, type I IFN signaling, etc., are highlighted among other cellular functions. (D) Integrative analyses of the SHAe genomic coordinates in alternative human systems challenged with immunogenic or inflammatory stimuli. The results highlight that several SHAe loci work as IFN-α responsive cis-elements in K562 cells. In addition, extended fractions of SHAe genomic coordinates host p65-binding in Detroit 562 cells upon TNF-a- or LPS-stimulation, and exhibit outstanding specificity for defensive and immune cellular responses, such as those against other organisms, LPS, cytokines, etc. For these analyses, the genomic coordinates of SHAe were extended bidirectionally (±2.5 kb) to avoid exclusion of biologically meaningful coincidences.

Figure 6.

Antimicrobial TFs cooperate at the SHAe level. (A) (i) Central upper panel: Genomic localization assessments capture similar patterns of distribution between IRF3- and p65-SHAe within intra- and inter-genic regions. Upper left and right panels: Analyses of the lexicon of TFBSs of SHAe elucidate DNA grammar characteristics legible, among others, by members of the microbial-activated families of TFs e.g. IRFs, RELA, STATs, FOS/JUN, etc. Central middle panel: Diagram depicts a spectrum of 219 IRF3/p65-SHAe (common) within the cohorts of IRF3- and p65-SHAe, that emerged upon intersection of the respective genomic coordinates. Below: Aggregation plots (signals’ average) and heatmaps (signals’ coverage) illustrate the ChIP-STARR-seq signals captured in naïve (SVI, 0 h; green line) and virus-infected HL (SVI, 6 h; red line) upon assessments of IRF3/p65-SHAe in individual IRF3-ChIP-STARR-seq (left panel) and p65-ChIP-STARR-seq (right panel) assays. Blue line corresponds to the signals derived from direct NGS of the multiplexed-STARR-plasmid library. In both IRF3- and p65-ChIP-STARR-seq assays, aggressive upregulation is depicted upon virus-infection for IRF3/p65-SHAe. Below: Scatter plot depicts that the vast majority of IRF3/p65-SHAe (∼90%) executes strong (FC ≥ 2) virus-stimulated upregulation of STARR-transcripts (average value of FCs derived from the individual assays). Tables summarize striking examples of SHAe-associated vruDEGs (GREAT analysis; The entire catalog is charted in Supplementary Fig. S23B). (ii) Topographic maps of human SHAe endogenous loci. High-resolution IGVs snapshots of genomics tracks illustrate the epigenetic and transcriptional states of in vivo reconstituted IRF3-SHAe (TNFRSF8 locus), p65-SHAe (OPTN locus), and IRF3/p65-SHAe (IL7R and OASL loci) in HL (SVI, 0, 3, and 6 h). IL7R locus is also discussed in Fig. 10. CHs assembly, histone modifications, and distinct patterns of antimicrobial TFs binding and CBP recruitment verify the specificity of the results. (iii) Bar graph depicts the allocation (absolute number) of SHAe-associated vruDEGs. (B) Bar graph depicts the distribution (%) of IRF3/p65-SHAe within SEs, tEs, and “safGs.” (C) (i) GOs validate statistically significant correlations of IRF3/p65-SHAe with defensive cellular functions, e.g. response to virus, type I IFN signaling, response to organism. (ii) QRTMs meta-profiling of epigenomics signals across the endogenous IRF3/p65-SHAe loci: During the transition from naïve (SVI, 0 h) to antiviral (SVI, 3 and 6 h) cell states IRF3/p65-SHAe exhibit gradual aggressive enrichment of “open” chromatin, and MED1 recruitment, H3K27ac labeling, strong occupation by antimicrobial TFs, CBP, and RNA pol II, while a slight change in H3K4me3 labeling is also detected. Hence, IRF3-, p65-, and IRF3/p65-SHAe work as virus-inducible transcriptional enhancers, in vivo.

Figure 7.

Tracing the evolutionary origins of SHAe. (A) The phylogeny of the human SHAe Repeatome; Upper left and right panels: Sankey plots depict the distribution of the 1415 Repetitive-IRF3-SHAe (left) and the 1052 Repetitive-p65-SHAe (right) in between the diverse classes of repetitive DNA. DRs, TRs, and hybrids of those including c(T:D)Rs composite the initial nodes (entities) of the diagrams that are then subclassified in additional linked nodes including class I retrotransposons (LINEs, SINEs, and LTRs, class II DNA transposons, STRs etc.). In both Repeatomes, DRs are the predominant class followed by hybrids and STRs. Central bottom panels: Monitoring of the evolutionary history of repetitive-SHAe within the cohort of the 100 vertebrates examined traced the vast majority (∼90%) of those upon the emergence of mammalian species (e-mammals), whereas a limited spectrum displays a more prolonged evolutionary history and captured in the genomes of nonmammalian vertebrates (b-mammals). Left and right bottom panels: Density plots depict statistically significant differences between the evolutionary conservation of repetitive- versus nonrepetitive- IRF3- and p65-SHAe, as computed by the average PhyloP score of each element [two-sided Wilcoxon rank sum test (P-value = 3.33e⁻¹⁶ for IRF3-SHAe and P-value = 5.28e⁻⁸ for p65-SHAe)]. In principle, repetitive-SHAe are “newly acquired” in the human genome compared to nonrepetitive. (B) Upper panels: Bar graphs illustrate the evolutionary distribution of human IRF3- and p65-repetitive-SHAe between the classes of e-primates, e-mammals, and b-mammals accompanied by representative examples of DNA sequence conservation within the cohort of 100 vertebrates. These patterns validate the evolutionary hierarchies of IRF3- and p65-Repeatomes. (C) The evolutionary history of human higher primates-enriched SHAe (HHPe-primates); Side-by-side monitoring of the conservation of vrCRMs within the 100 vertebrates’ (left panels) and the 241 placental mammals’ (right panels). The entire DNA sequence of this SHAe is not traced in the genomes of other organisms in both cohorts and is also classified in G3 cCREs. Part of Fig. 7 was created in BioRender: Agelopoulos, M. (2025) https://BioRender.com/x93p313

Figure 8.

DNA Grammar and Syntax assessments and phylogenomics investigations mapped the evolution of hundreds of homotypic clusters of TF-binding sites (HCTFBSs) within composite sequences of the IRF3- and p65-SHAe. (A) The pyramid of the >180 human SHAe that encompass extended composite DNA sequences (>150 bp) highly conserved in one or multiple viral genomes. Central panels: Analyses of the distribution among the classes of SHAe. Right panels: A striking example of pairwise alignment of the DNA sequence of the SHAe and BeAn 58058 virus genome. TFBMs analyses applied exclusively on the segment of SHAe that exhibits conservation with viral genome(s) highlights that CTCF, STAT2, NFKB2, and CREB1 TFBMs are among the top-ranked. (B) Striking examples of IR- (left) and κBR-HCTFBSs (right) composite DNA sequences which illustrate the density of perfect, degenerate, and half TFBSs of IRF3 and NF-κB, respectively. Magenta (IRF3) and blue (NF-κB/Rel) letters label half sites, green letters perfect IRF3 binding sites, and red letters NF-κB/Rel-binding sites. (C) (i) and (iii): Sunburst charts illustrate the distribution of HCTFBSs-SHAe within distinct classes of repetitive DNA. In principle, IR-HCTFBSs (n = 406) and to a lesser extent κBR-HCTFBSs (n = 343) reside predominantly within the repetitive-SHAe spectrums, consistent with their concatemerized structures. IR-HCTFBSs predominantly bear hybrid sequences, while within κBR-HCTFBSs, the DRs are vastly represented. (ii) and (iv): Monitoring of the evolutionary history underscores that repetitive-IR-HCTFBSs (n = 406) and repetitive-κBR-HCTFBSs (n = 343) are “newly acquired” compared to nonrepetitive-IR-HCTFBSs (n = 76) and nonrepetitive-κBR-HCTFBSs (n = 245), respectively. Statistically significant differences of PhyloP scores between the spectrums of repetitive and nonrepetitive HCTFBSs were identified [two-sided Wilcoxon rank sum test (P-value = 9.4e⁻¹⁵ for IR-HCTFBSs and P-value = 0.014 for κBR-HCTFBSs)]. This analysis is complemented in Supplementary Fig. S28. (D) (i) and (iii): DNA Grammar and Syntax analysis on the composite sequences shows that the core elements of both IR-HCTFBSs and κBR-HCTFBSs exhibit “Head-to-Tail” and to a lesser extent “Head-to-Head” orientation of their embedded perfect TFBSs, a feature that reminisces TRs expansion within the genomes. (ii) and (iv): GOs on IR-HCTFBSs and κBR-HCTFBSs show substantial specificity for defensive and immune-response cellular processes. (E) (i) and (ii): Left and middle panels: IR-HCTFBSs and κBR-HCTFBSs host robust/stronger virus-inducible binding of IRF3 and p65, respectively, and execute stronger/wealthier virus-inducible STARR-transcripts upregulation in vivo compared to non-IR- and non-κBR-HCTFBSs, respectively [two-sided Wilcoxon rank sum tests (Di left: P-value < 2.2e⁻¹⁶; Di middle: P-value < 2.2e⁻¹⁶; Dii left: P-value = 3.08e⁻⁶; Dii, middle: P-value = 9.54e⁻¹¹)]. The Boxplots are constructed as detailed in materials and methods and Fig. 1 legend. Right panels: Genomic localization assessments map the on-cis-elements distribution (%) of HCTFBSs: IR-HCTFBSs: 58 are embedded within 54 unique SEs (10 elements in pmSEs and 48 elements in viSEs), 260 within tEs, and 164 “safGs”. κBR-HCTFBSs: 69 are embedded within 66 unique SEs (15 elements in pmSEs, and 54 elements in viSEs), 131 within tEs, as verified by the ROSE algorithm, and 388 within “safGs”. (F) Chromosomal Ideograms (e-karyotypes) illustrate the broad in vivo reconstitution of functional virus-responsive IR-HCTFBSs, κBR-HCTFBSs, and IR//κBR-HCTFBSs in virus-infected HL (SVI, 6 h), across chromosomal telomeres, intra-chromosomal domains, and centromeres. Significant linear on-genome associations are marked between HCTFBSs and proximal vruDEGs (GREAT analysis) and SEs. Green circles: IR-HCTFBSs; Red circles: κBR-HCTFBSs; Orange circles: IR//κBR- HCTFBSs; Magenta circles: pmSEs; Red arrows/triangles: viSEs. The above findings validate the functional conductance of clustered cis-elements of diverse phylogeny that are recognized by IRF3 and/or NF-κB in virus-stimulated gene expression in human cells. Part of Fig. 8 was created in BioRender: Agelopoulos, M. (2025) https://BioRender.com/d37w614

Figure 9.

Tracing IRF3 HCTFBSs fingerprints in evolutionary recent, old, and ancient genomes, and characterization of vruDTTEs. (A) HOMER2 software tool (v4.11) and the de novo function computed the most enriched motif of an adjusted 25 bp length within the 482 IR-HCTFBSs which is termed “cIR-motif,” encompasses two perfect IRF3-binding sites (“RRAARGGAAAGGAAAGGAAAGGAAA”), and efficiently probes in silico “IRF3-clustered binding.” Lower panel: The 97 genomes investigated (retrieved either as default assemblies (PWMScan) or manually curated from the NCBI RefSeq database) belong to organisms with divergent evolutionary origins: (i) Metazoans: vertebrates (e.g. primates and mammals) and invertebrates (e.g. worms, insects, and mollusk), (ii) Plants: Arabidopsis thaliana, (iii) unicellular eukaryotes (e.g. fungi and protozoa), (iv) Prokaryotes (bacteria and archaea), and (v) and Human infectious Viruses. Percentage (%) of genome coverage (bar graph) and number of “cIR-motif” instances (right lane) are displayed. (B) Highlighted in dark are the 102 IRF3-target genes orthologous in mice and human (HL vruDEGs; SVI, 6 h). Upper bar graph: The spectrum of 102 IRF3-target genes is associated with 280 virus-inducible IRF3-ChIP-seq signals in NIH/3T3 (SVI, 6 h), 187 (66%) of which are orthologous in human according to lift-over analysis. Bottom bar plot: The 102 IRF3-target genes within the full spectrum of mouse genes that reside in cis proximity (1 Mb; GREAT analysis) to virus-inducible IRF3-ChIP-seq signals. (C) Genome innovation episodes mark the evolution of human vrCRMs: (i) The CCL5-Ccl5 paradigm: Three distinguished IRF3-SHAe compose a repertoire of vrCRMs that regulate the human CCL5 gene and reside within the modular architecture of its viSE (rectangle and horizontal line). All three of those are IR-HCTFBSs, one TSS-proximal (+237 bp), one TSS-distal (-9687 bp), and one intervening (-3866 bp) that was captured by the “cIR-motif.” The orthologous mouse locus bears two mIR-HCTFBSs (-68 bp, -4341 bp). The segment of the human intervening sequence that is missing from the mouse genome hosts the -3 866 IR-HCTFBSs, is classified in higher Primates, and harbors multiple insertions of TEs, including SINEs (Alu, MIR) and ERVs (LTR35A, LTR8). This genome innovation episode might happen via horizontal transmission of viral DNA and grants human cells an additional vrCRM that might facilitate the in vivo reconstitution of the viSE. (ii) Human CHEK2 harbors an IR-HCTFBS that is located ∼11.3 kb downstream of the TSS (GREAT analysis), inhabits the viSE (rectangle and horizontal line) of the gene, and coincides with its third exon (all classified in e-primates). It is structured by TR and Alu elements, and has viral origins. This genome innovation episode leads to the generation of another vrCRM, presumably the “birth” of a viSE, and the encoding of additional amino acids that lead to a new isoform of the human protein (exonization). (D) Dissecting the human virus-stimulating TEs-encoded transcriptome; (i) TEtranscripts tool analyses identified 60 distinct families of vruDTTEs in virus-infected HL (SVI, 6 h) e.g. MER107, LTR26, and MER57B1. (ii) Binomial test followed by Wilcoxon rank sum test within the first 20 kb highlight statistically significant enrichment of in cis proximity between vruDTTEs and vruDEGs, in HL, when total (continuous line) or distal intergenic copies (dotted line) were assessed, compared to an equal number of nonvirus-upregulated TEs (dashed line) (P-value < 2.2e⁻¹⁶). (iii) vruDTTEs inhabit genomic loci more proximal to SHAe (yellow line) compared to the total STARR-Active elements (IRF3; green line, p65; red line) [Dunn’s test in 0.1 cutoff in the respective relative distances (P-value < 2.2e⁻¹⁶)]. (iv) Bidirectional demarcation of vruDTTEs across ∼6 Mb flanking distance was conducted by step-wise increase (doubling) of the flanking sequences and utilizing the center of each SHAe as the reference point of the genome. This analysis highlights linearly increased inhabitance across ∼300 kb, which is validated in (v) by increasing the resolution of analysis within ± 40 kb flanking distances (Spearman correlation coefficient, R = 1, P< 2.2e⁻¹⁶), and mapped 222 SHAe topologically associated with core antiviral genes vruDTTEs (Supplementary Table S5). Part of Fig. 9 was created in BioRender: Agelopoulos, M. (2025) https://BioRender.com/k72m711.

Figure 10.

Integrative analyses of SHAe loci that are enriched in SNPs and other short variants linked to human (auto)immune diseases. (A) Upper left panel: A description of the workflow. Upper center panel: The transcriptional and epigenomics states of the IRF3/p65-SHAe and its flanking ±2.5 kb sequence. These genomic coordinates are encompassed within the IL7R genomic locus and are constituents of the viSE of the gene (horizontal line). A CH is assembled prior to (SVI, 0 h) and becomes enhanced upon virus-infection (SVI 3 and 6 h) and is accompanied by the binding of the antimicrobial TFs (IRF3 and p65) and the recruitment of CBP, thus leading to the in vivo reconstitution of functional vrCRMs. RNA pol II travels across the gene body and transcribes its coding region, which apart from protein-coding sequences harbors multiple vruDTTEs (orange lines) that generate TEs-encoded RNAs (lower orange rectangle). Several autoimmune-disease-associated SNPs are mapped within the above genomic coordinates, and the majority of those are linked to MS. Interestingly, rs10063294 inhabits the IRF3/p65-SHAe sequence within an IRF3-binding site, as shown in (B) (sequence panel). Upper right panel: A catalog of the SNPs and their associated Human autoimmune diseases. Middle left panel: A chart of striking examples of characterized deletions and duplications that reside within the IL7R locus, accompanied by their associated phenotypes’ descriptions. The entire analysis is charted in GitHub Repository. (B) Upper panel: The evolutionary origins of the SHAe and its in cis associated SNPs within the cohort of 100 representative vertebrate species (phylop100way); In principle, the chr5:35874430–35880176 locus of the human genome is not traced beyond mammalian species, and its sequence conservation is significantly increased from lower mammals such as Aardvark, to primates. Middle panel: The sequence architecture of the locus [TFBSs for IRF1/3 (green letters) and NF-κB (red letters); IRF3/p65-SHAe (blue shadow); SNPs (brown letters)]. Lower panel: Linkage disequilibrium (LD) analyses performed by the calculation of the squared correlation (r²) between SNPs. A pair-wised approach was conducted for the total of SNPs and the results highlighted significant correlations with the rs6881706–rs6897932 pair that is classified at the top of the ranking, reaching an absolute correlation. Both SNPs are associated with MS. Part of Fig. 10 was created in BioRender: Agelopoulos, M. (2025) https://BioRender.com/b54k071.