As previous studies demonstrated that non-B DNA motifs are enriched for somatic mutations,^11,14,15 we set out to analyze whether this enrichment also exists for germline variation. We took advantage of available whole-genome sequencing (WGS) datasets for thousands of individuals and analyzed them to determine whether non-B DNA sequences are enriched for variants. We measured the genome-wide distribution of 204,063,503 single-nucleotide polymorphisms (SNPs), including both rare and common variants as well as 25,925,202 small insertions and deletions (indels; <50 bp) derived from 15,496 genomes from the gnomAD project⁴⁰ relative to seven non-B DNA motifs: inverted repeats (IRs), direct repeats (DRs), mirror repeats (MRs), short tandem repeats (STRs), G4s, Z-DNA, and H-DNA motifs (Figures 1A–1G). To form a null distribution, we generated simulated SNPs, controlling for trinucleotide context and proximity to the original SNP (STAR Methods). We observed an excess of SNPs directly overlapping non-B DNA motifs (Figure S1A; Mann-Whitney U, p < 0.0001), but the magnitude of the effect was small, and the highly significant p value was due to the large sample size. Of note, H-DNA motifs and IRs showed the highest (1.56) and lowest (1.05) fold enrichments, respectively (Figures 1H, S1B, and S1C). Similarly, the proportion of indels overlapping non-B DNA motifs was substantially elevated relative to the simulated controls (2.26-fold, Mann-Whitney U, p < 0.0001; Figure S1D). The enrichment of genetic variants at individual non-B DNA motifs was higher for small indels than for SNPs, ranging from 2.44-fold for IRs to 13.68-fold for STRs (Figures 1H, S1E, and S1F). We further separated indels into insertions and deletions, finding differences depending on the non-B DNA motif category (Figure S1G). For example, STRs had a higher frequency of deletions, whereas G4s had a higher frequency of insertions.

Extending our analysis to 505,529 structural-variant breakpoints derived from the gnomAD project,⁴⁰ we found a strong association with non-B DNA motifs, with 14.61% of structural-variant breakpoints directly overlapping a non-B DNA motif versus 8.83% for simulated controls (Mann-Whitney U, p < 0.0001; Figure S1H), representing a 1.66-fold enrichment. For individual non-B DNA motifs, the enrichments ranged from 1.23-fold for G4s to 3.50-fold for H-DNA motifs (Figures 1H and S1I–S1K; Mann-Whitney p < 0.0001 for all non-B DNA motifs), consistent with previous reports finding an excess of non-B DNA motifs at structural-variant breakpoints.⁴¹ We separated structural variants into six categories: insertions, deletions, duplications, inversions, translocations, and complex.⁴⁰ We found that deletions, insertions, and duplications were the most enriched across non-B DNA motifs (Figure 1I). Taken together, these results suggest that non-B DNA motifs are hotspots of genetic variation in the human population across different categories of population variants.

To gain further insights regarding the regulatory potential of these variants, we investigated the relative frequency of variants overlapping non-B DNA motifs across six regulatory-element-associated sequences/functions defined by the Ensembl Regulatory Build:⁴² promoters, CTCF-binding sites, open chromatin regions, transcription factor binding sites, promoter flanking regions, and enhancers. The analysis was performed across twelve different cell lines (STAR Methods), finding that most non-B DNA motifs were enriched for SNPs, indels, and structural variants across the regulatory elements, but more so for indels than for SNPs and structural variants (Figures S2A–S2C). We also investigated the increase in mutagenicity for non-B DNA motifs across the seven annotated genic sub-compartments: genic, intronic, coding, and 5′ and 3′ UTRs as well as 1 kilobase (kb) upstream of the transcription start site (TSS) and 1 kb downstream of the transcription end site (TES). Most regions had elevated mutation rates, although the magnitude varied by mutation type and genic sub-compartment (Figures S2D–S2F). As expected, coding regions showed the lowest mutagenicity relative to other regions, most likely due to selection constraints and increased DNA repair.⁴³

To analyze whether variants in non-B DNA motifs could have a substantial impact on gene expression, we analyzed expression quantitative trait loci (eQTL). We examined the frequency of eQTLs, characterized by the GTEx consortium,⁴⁴ at each of the seven non-B DNA motifs genome wide. We found an enrichment of eQTLs across all non-B DNA categories relative to their flanking regions, with the most pronounced effect for G4s (Figure 1J). Although the excess of eQTLs in the vicinity of non-B DNA motifs can be explained by the higher background frequency of substitution and indel SNPs (Figure 1H), our results indicate that a subset of mutations overlapping non-B DNA motifs impact gene expression.

As G4s had the most pronounced effect on gene expression, we next analyzed G4 sequencing (G4-seq) and G4 chromatin immunoprecipitation (ChIP)-seq datasets for their overlap with population variants and eQTLs. We investigated the association between population variants and G4s using previously published G4-seq datasets from the HEK-293T cell line with Pyridostatin (PDS) and K⁺ treatments that provide in vitro evidence of G4 formation potential⁴⁵ and G4 ChIP-seq-derived peaks from the HaCat cell line that provide in vivo evidence of sites that form G4 structures.¹⁶ In accordance with the G4 motif analysis, we found that SNPs, indels, and structural variants were enriched at G4-seq and G4 ChIP-seq peaks (Figures S3A–S3F; Mann-Whitney U, p < 0.001). We considered the G4 ChIP-seq sites that overlapped both G4-seq K⁺ and G4-seq PDS peaks as the highest confidence, experimentally derived G4s (Figure S3G) and found consistent enrichments of 1.14-fold, 1.41-fold, and 1.36-fold for substitutions, small indels, and structural variants (Figures S3H and S3I). Next, we found that eQTLs are found more frequently than expected by chance in the experimentally derived G4 sites. In total, 20,310 eQTLs overlapped with the 8,955 ChIP-seq peaks, with 34% of the peaks having one or more eQTL (Figures 1K and S3J). Interestingly, the enrichment for the experimentally derived G4s was more pronounced than our results derived from the G4 motif analysis. This is likely the result of G4 formation occurring more frequently in open chromatin and transcribed regions.¹⁶

We also investigated if G4 ChIP-seq peaks overlapping genes display a preference for the template (non-coding) or non-template (coding) strands, using the G4 motif orientation within the peaks as proxy. After correcting for the background bias in the orientation of G4 motifs (Figure S4A), we found that G4 motifs on the non-template strand overlap G4 ChIP-seq peaks 1.71-fold more frequently than motifs on the template strand (binomial test, p < 1 × 10⁻¹²) (Figures S4B and S4C), suggesting significant bias in the formation of G4s, dependent on their orientation.

We next investigated the distribution of non-B DNA motifs across the six regulatory elements defined by the Ensembl Regulatory Build (promoters, CTCF-binding sites, open chromatin regions, transcription factor binding sites, promoter flanking regions, and enhancers). For most non-B DNA motifs, we found an enrichment at promoters and CTCF-binding sites relative to other regulatory elements (Figures 2A and S5A), in accordance with previous findings.⁴⁶ Next, we separated the gene body into six compartments: a 1 kb window upstream from the TSS, the 5′ and 3′ UTRs, coding exons and introns, and a 1 kb window downstream of the TES. Consistently, promoter regions displayed a higher density of non-B DNA motifs than the gene body for most non-B DNA motifs, with the enrichment ranging from 0.97-fold for IRs to 3.02-fold for G4s (Figures 2B and S5B). We also found a significant enrichment of G4-seq-derived peaks for both PDS and K⁺ treatments and for G4 ChIP-seq-derived peaks at promoters relative to other regulatory elements (Figure 2C). Across the gene body, we found the highest enrichments at promoters, coding regions, and 5′ UTRs (Figure 2D).

At promoters, for most non-B DNA motifs, the enrichment was higher upstream of the TSS than in the broader promoter region (Figure 2E). A close investigation of the distribution of non-B DNA motifs relative to the TSS showed an enrichment of peaks ∼50 bp upstream of the TSS ranging between 1.28- and 1.89-fold for DRs and G4 motifs, respectively (Figure 2F). Importantly, we observed a 5-fold enrichment approximately 100 bp upstream of the TSS for G4 ChIP-seq peaks, consistent with the literature.¹⁶ Interestingly, the ChIP-seq-derived enrichment was substantially larger than that of the G4 motif and the G4-seq datasets (Figure 2G), reflecting a preference in structure formation at promoters in vivo. We also performed a Gene Ontology (GO) term analysis in promoter upstream regions. For G4s, Z-DNA motifs, and MRs, we found multiple terms associated with developmental processes, such as pattern specification process (GO: 0007389), embryonic organ development (GO: 0048568), and positive regulation of neuron differentiation (GO: 0045666) (Figure S6A). As these analyses suggest that some non-B DNA motifs could control tissue-specific gene expression, we used TissueEnrich to calculate the enrichment of tissue-specific genes and found sets of tissue-specific genes where a set of neuronal-specific genes were enriched for genes containing G4, MR, DR, and STR at their upstream promoter regions (Figure S6B). Altogether, these results demonstrate that promoters are enriched for non-B DNA motifs relative to other regulatory elements and to other genic compartments and that some non-B DNA motifs are more likely to occur at developmental and neuronal genes. Therefore, the excess of genetic variants at non-B DNA motifs identified earlier could have broad implications on gene regulation expression levels across tissues and developmental stages.

The enrichment of non-B DNA motifs at promoters and the excess of eQTLs localized within certain non-B DNA motifs prompted us to investigate their functional impact on gene transcription utilizing MPRAs. We first analyzed two lentiMPRA datasets generated by our group as part of the ENCODE consortium,⁴⁷ where a total of 14,625 and 7,346 candidate promoter sequences were examined in both orientations in K562 and HepG2 cell lines. We identified non-B DNA motifs across the lentiMPRA tested sequences (STAR Methods) and examined their association with gene expression. We found that sequences with G4 and Z-DNA motifs showed significantly increased expression levels in both cell lines (Figures 3A and 3B; t tests, Bonferroni correction, p < 0.001), while for IRs, DRs, STRs, and MRs we did not observe consistent results (Figure S7A). As there is a known positive correlation between expression and guanine-cytosine (GC) content,⁴⁸ which was also observed in our lentiMPRA datasets (Pearson r = 0.398 and 0.261 in K562 and HepG2, respectively), we constructed a linear model to account for the contribution of GC content toward expression (Figure S7B). Sequences with Z-DNA motifs had substantially elevated expression levels relative to sequences without them, even after controlling for GC content in both cell lines (t tests, Bonferroni correction p < 0.001; Figures 3C and S7C). However, after GC-content correction, G4s were not associated with increased expression, and in HepG2, they were instead significantly associated with reduced expression levels (Figures 3C and S7C). Similar results were obtained after removing outliers from the linear model (absolute Ζ score >2.5). Also, G4s on the template strand were associated with reduced expression relative to non-template strands in both cell lines, but the difference reached statistical significance only in the HepG2 cell line (Figure S7D). For the other non-B DNA motifs, we could not find consistent effects in both cell lines, suggesting that nucleotide composition contributed to the observed effects before GC-content correction.

Finally, we identified transcription factor binding sites (TFBSs) across the MPRA sequences using the JASPAR vertebrate non-redundant list of transcription factor motifs.⁴⁹ We compared the contribution of non-B DNA motifs relative to TFBSs toward expression levels, both before and after GC-content correction. We found that G4 and Z-DNA motifs had similar contributions to known TFBSs, such as EGR1, YY1, and SP9, resulting in increased expression levels relative to sequences without them (Figure S8). However, only Z-DNA motifs had comparable effects when we accounted for GC content (Figures 3D and 3E), and the results were consistent between HepG2 and K562 lentiMPRAs.

To further validate our findings, we analyzed lentiMPRA results from a library that characterized the effect of 3,623 de novo promoter mutations that were identified in the Simons Simplex Collection.⁵⁰ This library tested both alleles, centered around the variant, totaling 7,246 sequences along with 150 positive and 150 negative controls for their effect on promoter activity in neural progenitor cells (NPCs) (Figures S9A–S9C). This library had 1,234 sequences harboring one or more non-B DNA motifs (Figure S9D). We observed that sequences harboring G4, DR, and Z-DNA motifs displayed a significantly higher expression than sequences without them (t tests, Bonferroni corrected p values, G4s, DRs, and Z-DNA p < 0.001), whereas sequences with IRs, MRs, and STRs did not show a significant association (p > 0.05) (Figure 4A).

Similar to the analysis of the ENCODE MPRA libraries, we observed a significant contribution of the GC content toward the effects on expression of certain non-B DNA motifs. After constructing a linear model to adjust for GC content, we observed that G4 motifs are associated with decreased expression, while only Z-DNA sequences remained associated with higher expression (Figure 4B), consistent with previous results. In this case, removing outliers maintained a positive association with G4s and gene expression. We also observed a substantial difference in the expression dependent on the orientation of G4 motifs, with G4s on the template strand having lower expression than those on the non-template strand before and after GC-content adjustment (Figures 4C and S9E; Mann-Whitney U, p < 0.001). The primary sequence comprising consecutive G-runs that are interspersed by loop elements can form G4 structures (Figure 1). The association between G-runs and gene expression was further investigated, finding that consecutive G-runs result in decreased expression when accounting for their GC-content contribution (Figure S9F). Furthermore, we found that the length of the Z-DNA motif was positively associated with the expression levels (Kruskal-Wallis H test, p < 0.001; Figure 4D).

Similar to the previous MPRAs, we identified TFBSs across the MPRA sequences and compared the contribution of non-B DNA motifs relative to TFBSs toward expression levels before and after GC-content correction. We found that G4 and Z-DNA motifs had comparable contributions to TFBSs toward increasing expression levels with increases of 1.27- and 1.51-fold over sequences without them (Figure S10A). However, when we accounted for GC content, the effect of non-B DNA motifs was not comparable to the best TFBS motifs (Figure S10B). Therefore, we find substantial differences in the results in NPCs relative to HepG2 and K562 cell lines, with a lower contribution of Z-DNA motifs in NPCs, which might be due to the selection of loci that were not necessarily proximal to the TSS or due to the lower number of Z-DNA-containing sequences, with only 311 sequences having them.

To validate if the G4s we observed in this NPC lentiMPRA form these structures, we selected ten candidate promoter-proximal sequences with the lowest and highest expression among sequences with G4s (Table S1) and performed multiple spectroscopic assays to characterize their structures (Figures 4E and 4F), as G4 structures possess distinct spectroscopic features.^51,52 We first used circular dichroism spectroscopy measurements of the G4-containing DNA oligonucleotides, in the presence of lithium ions (non-G4 stabilizing) or potassium ions (G4 stabilizing), to examine the formation potential of DNA G4s, which indicated that our candidate sequences can fold into G4 structures (Figures 4E, 4F, S11A, and S11B). In addition, we conducted UV melting and found a hypochromic shift at 295 nm for the potassium-ion condition, which supported the formation of the G4 structure, with a melting temperature above physiological temperature (Figures 4E, 4F, S11A, and S11B).

To confirm the results from the circular dichroism and UV-melting experiments, we used fluorescent-based arrays, including N-methyl mesoporphyrin IX (NMM)-ligand-enhanced fluorescence and intrinsic fluorescence experiments (Figures 4G, 4H,S12, and S12B). In the absence of NMM ligand, no fluorescence was observed at ∼610 nm. Upon NMM addition, weak fluorescence was observed under Li⁺, which was substantially enhanced when substituted with K⁺, supporting the formation of G4 that allows recognition of NMM and enhances its fluorescence (Figure 4G). Similarly, the intrinsic fluorescence of G4s was increased when replacing Li⁺ with K⁺, highlighting the formation of DNA G4s (Figure 4H). Corroborating our results, we observed increased fluorescence intensity under conditions that promote G4 formation for all candidates. We also carried out two positive G4 controls and a negative B-DNA control to verify our findings above (Figure S13). Combined, these results validate that these sequences form G4 structures in vitro.

To directly test the effect of non-B DNA structures on promoter activity, we generated an MPRA library that introduces various non-B DNA perturbations to ten disease-associated genes. This set of genes included cancer oncogenes (CMYC, CKIT, BCL2, KRAS) and genes associated with different cancer types (ADAM12, ALOX5, SRSF6, VEGF12) as well as FMR1, associated with fragile X syndrome (OMIM: 300624), and SNX12, which is associated with neurodegenerative diseases (Table S2). As our MPRA-tested sequences are 200 bp in length, we first validated whether our selected 200 bp sequences could drive promoter activity using luciferase assays in K562, MCF7, IMR90, and HEK293T cells, finding the majority to be active in most cell lines (Figures S14A and S14B).

Following validation of these 200 bp sequences, we next generated an MPRA library that included the following manipulations: (1) disruption of existing non-B DNA motifs and (2) introduction of different non-B DNA motifs with varied biophysical properties, including spacer- and arm-length changes in IRs, DRs, and MRs, orientation and loop length in G4s, and length in Z-DNA motifs. lentiMPRAs and subsequent computational analyses were carried out as previously described.⁵³ Briefly, oligonucleotides were synthesized and cloned into a lentiviral MPRA promoter vector (Figure 5A; Table S2), and lentivirus libraries were generated. Libraries were used to infect both K562 and HEK293T cells for 3 days, to allow non-integrating lentivirus to degenerate, and DNA and RNA barcodes were sequenced. Since previous work in our lab showed that lower basal activity can have a significant effect on MPRA results,⁵⁴ these two cell lines were chosen as almost all the selected promoters showed ≥2-fold activity compared with empty vectors (except for CKIT in HEK293T). All experiments were done in triplicate, and computational analyses were carried out using MPRAflow⁵³ and MPRAnalyze.⁵⁵ We observed a strong correlation between all three replicates (Pearson r ≥ 0.9 in all cases; Figure S15A) and between the two cell lines (Pearson r = 0.87; Figure S15B).

The promoters in our MPRAs showed variable expression, with the highest levels observed for SRSF6 and the lowest for ADAM12 (Figure 5B). We investigated the contribution of each non-B DNA motif toward expression in both cell lines across the promoters, adjusting across genes using Z score normalization. Specifically, for each gene we calculated the Z score of each sequence, which was calculated by subtracting the expression levels of that sequence from the mean across all sequences of that gene and dividing by the standard deviation. In concordance with our previous MPRA analyses, we observed that sequences with Z-DNA and G4 motifs had significantly higher expression (Figures 5C and 5D). Interestingly, while we did not observe consistent results in our previous MPRA analyses for MRs, DRs, and IRs, here, we observed significantly higher expression levels when MRs and DRs were present, whereas for IRs we found significantly lower expression (Figure 5C). For STRs, we did not find consistent patterns in the two cell lines. The above results across non-B DNA motifs did not change when we accounted for GC content; however, this was most likely due to our experimental design having only a small number of loci targeted, which, as a result, had a narrow and uninformative GC-content range.

For G4s, we introduced a single, two, or three mutations in one, two, three, or every G-run at the original G4 genomic sites. We compared the mutated sequences with the original sequence and found that sequences with disruptions in the G-runs did not display significant expression differences from the original sequences (Figure S15C). We designed MPRA sequences with scrambled Z-DNA motifs or with disruptions of purines to pyrimidines in the alternating purine-pyrimidine tract, which served as Z-DNA controls. We found that there was a statistically significant reduction in expression following the disruption of Z-DNA motifs (Figure 5D), supporting the notion that they are activating sequences. We also observed that non-template G4s had higher expression than those at the template strand in both cell lines and both before and after GC-content correction (Mann-Whitney U, Bonferroni corrected; Figure 5E), consistent with our earlier results. For Z-DNA, longer motifs resulted in higher expression (Figure 5F). These results suggest that the non-B DNA motifs and their biophysical properties contribute to expression across promoter templates.

Our study has multiple limitations. First, the examination of the regulatory roles of non-B DNA motifs through MPRA experiments did not investigate how molecules that stabilize their formation affect the conclusions reached. Secondly, the MPRA results are based on specific cell lines, and it would be of interest to examine which of these findings can be generalized across cell types and which effects are cell-type specific. We also cannot exclude the influence of the experimental design in our findings. Furthermore, additional experiments and mechanistic work are required to further our understanding, including biophysical and molecular experiments. Lastly, future work would be needed to resolve the relevance of mutations at non-B DNA motifs in the development and progression of human diseases. The aforementioned limitations could be of high interest for future work.

Cell culturing was performed for HEK293T (RRID CVCL_0063), K562 (RRID CVCL_0004), MCF-7 (RRID: CVCL_0031) and IMR-90 (RRID: CVCL_0347) cell lines. Human HEK293T embryonic kidney cells were cultured in Dulbecco’s modified Eagle’s medium (DMEM, Sigma) supplemented with 10% FBS and 2mmol/ L L-glutamine. Human K562 erythroleukemia cells were cultured in Iscove’s modified Dulbecco’s medium (IMDM, Sigma-Aldrich) supplemented with 10% FBS. Human MCF-7 breast cancer cells were cultured in Eagle’s minimal essential medium (MEM, Sigma-Aldrich) supplemented with 10% FBS, 10μg/ml insulin, 1mM sodium pyruvate and 0.1 mM non-essential amino acids. Human IMR-90 fibroblasts were cultured in MEM supplemented with 10% FBS and 0.1 mM non-essential amino acids. Neural progenitor cells were differentiated from H1 hESCs following the dual-Smad inhibition protocol as described in ⁹¹. All cell lines were grown at 37°C and 5% CO₂.

I.G.-S., M.H., and N.A. conceived the study. I.G.-S., G.E.P., V.A., A.M., and J.Z. wrote the code and performed the analyses. I.G.-S., G.E.P., J.Z., and J.V. generated the visualizations. J.V., O.E., and F.I. performed the MPRA experiments. H.Y.W. and M.I.U. performed the circular dichroism (CD) titration, UV melting, and fluorescence assays, and H.Y.W., M.I.U., and C.K.K. analyzed and interpreted the spectroscopic data. J.-Y.A. and S.J.S worked on the NPC MPRA design. M.H. and N.A. supervised the research. I.G.-S., M.H., and N.A. wrote the manuscript with input from all authors.

The authors declare no competing interests.

The MPRA data for the NPC cell line targeting autism-related loci and the MPRA data for the non-B DNA associated loci in HEK-293T and K562 cell lines are deposited in NCBI BioProject with accession number PRJNA763774. The MPRA data for HEPG2 and K562 cell lines (Figure 3) have been deposited in the ENCODE portal with IDs ENCSR463IRX and ENCSR460LZI.

All original code and data tables to perform the analyses can be found on the GitHub page (https://github.com/IliasGeoSo/High_Throughput_MPRA_Non_B_DNA) and are publicly available.

Any additional information required to reanalyze the data reported in this paper is available from the lead contact upon request.