{'fraction_non_alphanumeric': 0.03770255093855286, 'fraction_numerical': 0.020985079416011552, 'mean_word_length': 5.0792938652101824, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 2, 'https://': 4, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Non-B DNA: a major contributor to small- and large-scale variation in nucleotide substitution frequencies across the genome

Abstract Approximately 13% of the human genome can fold into non-canonical (non-B) DNA structures (e.g. G-quadruplexes, Z-DNA, etc.), which have been implicated in vital cellular processes. Non-B DNA also hinders replication, increasing errors and facilitating mutagenesis, yet its contribution to genome-wide variation in mutation rates remains unexplored. Here, we conducted a comprehensive analysis of nucleotide substitution frequencies at non-B DNA loci within noncoding, non-repetitive genome regions, their ±2 kb flanking regions, and 1-Megabase windows, using human-orangutan divergence and human single-nucleotide polymorphisms. Functional data analysis at single-base resolution demonstrated that substitution frequencies are usually elevated at non-B DNA, with patterns specific to each non-B DNA type. Mirror, direct and inverted repeats have higher substitution frequencies in spacers than in repeat arms, whereas G-quadruplexes, particularly stable ones, have higher substitution frequencies in loops than in stems. Several non-B DNA types also affect substitution frequencies in their flanking regions. Finally, non-B DNA explains more variation than any other predictor in multiple regression models for diversity or divergence at 1-Megabase scale. Thus, non-B DNA substantially contributes to variation in substitution frequencies at small and large scales. Our results highlight the role of non-B DNA in germline mutagenesis with implications to evolution and genetic diseases.

nomic regions may be at a higher (or at a lower) risk of acquiring mutations important for adaptation and/or genetic diseases (1)(2)(3). In a broad sense, deciphering the causes of regional variation in mutation rates is essential to understanding both evolution and diseases (1,2).
Numerous genomic features contribute to regional variation in mutation rates, but those identified to date cannot account for all such variation. Some features are directly related to DNA sequence and usually act at the scale of single nucleotides, e.g. guanines and cytosines are more mutable than adenines and thymines (4,5). Neighboring nucleotides also have an effect, e.g. methylated cytosines in CpG dinucleotides are 10 times more mutable than other sites because of their spontaneous deamination (6), and several other contexts leading to guanine holes and increased mutagenesis were previously identified (7). Other genomic features--such as recombination rate (8), replication timing (9), chromatin accessibility (10,11), histone modifications (12,13), and Lamina Associated Domains (14)--contribute to regional variation in mutation rates through the variable activity of different enzymatic processes along the genome. These frequently act at larger scales, from several hundreds of kilobases to several megabases (Mbs). The magnitude of regional variation in mutation rates decreases with the increase in the genomic scale considered; most such regional variation in fact occurs at the single-nucleotide scale (1). At the 1-Mb scale, which is considered the natural long-range variation scale for mammalian genomes (15), most regions have mutation rates deviating by ∼2-fold (1). Notably, at this scale, several analyses indicated that the genomic features listed above explain only ∼50% of the regional variation in mutation rates (12,16,17). The correlation in regional variation in mutation rates between human and great apes (18) suggests that the unexplained portion of this variation is not random, and that additional factors remain to be discovered. Non-B DNA may be one such factor.
Certain DNA sequence motifs have the ability to fold (at least over part of their length) into secondary conformations that differ from the canonical right-handed B-DNA helix that has 10 bp per turn (19)(20)(21). Such non-B DNA sequence motifs (henceforth called 'non-B DNA motifs') range in length from a few dozen to a few hundred nucleotides and are non-randomly distributed across the genome (22)(23)(24). A genomic locus harboring a non-B DNA motif is usually referred to as a non-B DNA locus, and the same non-B DNA motif is frequently present at multiple loci in the genome.
Several types of non-B DNA have been identified based on the structures they can form (Figure 1), which in turn depend on their motif sequences (25,26). The Gquadruplex (henceforth called 'G4') structure ( Figure 1A) alternates stems consisting of guanines and loops consisting of unspecified nucleotides, with the canonical sequence of G 3+ N 1−12 G 3+ N 1−12 G 3+ N 1−12 G 3+ (27). Alternating purine and pyrimidine sequences can form Z-DNA ( Figure 1C)--a left-handed double-stranded helix with 12 bases per turn (28). Each mirror, inverted and direct repeat locus contains two repeat arms, usually separated by a non-repetitive spacer. Homopurines and homopyrimidines organized in mirror repeats with or without a ≤100-nucleotide spacer can form H-, or triplex, DNA ( Figure 1B) (29). Inverted sequence repeats with or without a ≤100-nucleotide spacer can form DNA cruciforms ( Figure 1D) (30,31). Direct sequence repeats with or without a ≤10-nucleotide spacer can form slipped-strand structures ( Figure 1E) (32). Finally, Aphased repeats--three or more segments composed of three to nine adenines and/or thymines (A-tracts) and whose centers are separated by 10 bp--can bend, or create a curvature in, the double helix ( Figure 1F) (26,33,34).
In the cell, non-B DNA loci can affect several DNA metabolic processes, and both replication-dependent and independent pathways contribute to elevated mutation rates at such loci. Replicative polymerases encounter many non-B DNA structures that act as natural impediments to DNA synthesis elongation. Specialized DNA polymerases associated with the replication fork, including Pols eta and kappa, perform highly efficient synthesis through non-B DNA structures (reviewed in (56,65,66)), and can take over synthesis from stalled replicative polymerases (67). DNA damage-induced mutagenesis has been associated with non-B structure formation (reviewed in (68)), and damage induced by reactive oxygen species is affected by the inherent structure and local sequence of DNA (reviewed in (69)). Non-B DNA structures are recognized by DNA repair pathways, such as nucleotide excision repair (70)(71)(72)(73), and error-prone repair processing can result in mutations (74)(75)(76). In addition, non-B DNA structures are associated with the formation of double-strand breaks (DSBs), which lead to increased genomic instability (71,(77)(78)(79)(80)(81)(82)(83)(84). Even though non-B DNA loci have the potential to influence mutations by affecting multiple processes in the cell, their effects on mutation rates have not been studied in detail.
While an elevated density of single nucleotide polymorphisms (SNPs) at non-B DNA was previously reported (24), little is known about where mutations concentrate within the non-B DNA motifs, and whether the increase in mutations is limited to the motifs themselves or extends into their flanking sequences. Furthermore, no studies have reported whether the effect of non-B DNA structures on SNP density depends on their stability. In all, genetic variation at non-B DNA and its flanking sequences remains understudied.
We recently demonstrated that, across the human genome, both the speed and accuracy of the modified bacteriophage phi29 polymerase used in Pacific Biosciences (PacBio) sequencing are modulated by the presence of non-B DNA loci and that, for G4 loci, this effect depends on their stability (85). We also showed G4 loci that are most divergent between human and orangutan or are highly diverse among human populations have more pronounced polymerization slowdown and error rates (85). These results suggest that both sequencing errors and germline mutations are elevated at G4 motifs and may have similar mechanisms, i.e. result from polymerization slowdown.
Here, we extend this study, and test a hypothesis that non-B DNA motifs contribute to regional variation in mutation rates. We use the frequency of human SNPs as a measure of intraspecific variation or 'diversity', and the frequency of fixed nucleotide substitutions (FNSs) between human and orangutan reference genomes as a measure of interspecific variation or 'divergence'. These measures proxy germline mutation rates because we minimize selection effects by limiting our analysis to putatively neutrally evolving non-coding regions of the genome. We analyze the influence of non-B DNA on nucleotide substitution frequencies at two different scales: the small scale of single nucleotides and the large scale of 1-Mb genomic windows. Using statistical methods from the functional data analysis domain, which investigate shapes of signals (86,87), we test whether diversity and divergence levels differ between each type of non-B DNA and randomly selected control sequences genome-wide. Furthermore, we test whether the predicted stability of G4 loci affects the levels of diversity and divergence. Moreover, we describe distortions in the nucleotide substitution spectrum surrounding G4s. Lastly, we evaluate whether adding non-B DNA loci to statistical models that include known contributing genomic features (1,12,16,88,89) increases the explained share of regional variation in diversity and divergence. Overall, our study measures the contributions of non-B DNA to germline variation in nucleotide substitution frequencies across the human genome.

Non-coding non-repetitive (NCNR) subgenome
To obtain the NCNR subgenome of the hg19 version of the human genome, we excluded (i) all repeats as annotated with RepeatMasker (90) (rmsk track at the UCSC Genome Browser (91)); (ii) NCBI RefSeq (92) genes (exons and their 1-kb up-and downstream flanking regions, and 5 kb up-stream of the 5 UTRs and downstream of the 3 UTRs); (iii) conserved elements as annotated with phastConsEle-ments100way (93) and (iv) enhancers as annotated with GeneHancerRegElementsDoubleElite (94,95) (Supplementary Figure S1). The hg19 version of the human genome was used because it has a larger number of annotations of genomic features than the more recent hg38 version.

Divergence and diversity nucleotide substitutions datasets
69 329 877 FNSs that occurred between human and orangutan were retrieved from the 100-way Vertebrate Multiz Alignment (96) obtained from the UCSC Genome Browser (91). SNPs from the Simons Genome Diversity Project (97) were acquired from the Seven Bridges Cancer Genomic Cloud (https://cgc.sbgenomics.com/). 44 833 480 SNPs from all individuals in this project were merged into a single BED file. We only considered variants (26 875 194 FNSs and 15 555 617 SNPs) located in the NCNR subgenome in subsequent analyses.

Non-B DNA annotations
G4 loci were annotated in the human genome reference (version hg19) using Quadron software (98). Annotations of the other non-B DNA loci (direct, mirror, and inverted repeats, as well as Z-DNA and A-phased motifs) were downloaded from the non-B Database (non-B DB, https://nonb-abcc. ncifcrf.gov/, as assessed on September 8, 2020) (99).

Small-scale variation in nucleotide substitution frequency
We compared SNP (diversity) and FNS (divergence) frequencies inside non-B DNA loci and their 2-kb up-and downstream flanking regions. We only considered non-B DNA located within the NCNR subgenome (Supplementary Table S1). Non-B DNA loci longer than 100 bp were discarded, to minimize the possibility that longer (and more rare) motifs and their flanking sequences may overlap with multiple shorter (and more common) loci. Even with this filter, some loci did overlap; when such overlaps occurred for loci of the same non-B DNA type, we retained only one locus (and its flanking regions) selected at random (Supplementary Figure S2). Overlaps between non-B DNA loci of different types are rare (23) and do not interfere with statistical analyses, thus we did not filter them out. For each SNP or FNS located in the flanking region of a non-B DNA locus, we computed the distance from the closest annotated motif end. For each SNP or FNS located inside a motif, we computed a position scaled between 0 and 1. This scaling was performed for each section of the non-B locus (stems and loops in G4s, repeat arms and spacers in inverted, direct, and mirror repeats) in order to 'align' all non-B loci of the same type while retaining the information about sections. 180 bins were considered in total, for all sections of the scaled interval to be larger than the ones observed at non-B DNA loci. The sections' maximal sizes were 30 bins for stems and 20 bins for loops in G4s; 50 bins for repeat arms and 80 bins for spacers in inverted, direct, and mirror repeats; 40 bins for A-tracts and 30 bins for spacers in A-phased repeats. Because the resolution in the scaled intervals was higher than in the original one, each SNP or FNS could occupy several scaled positions ( Supplementary Figure S3). For each non-B DNA locus and its flanking regions, a non-overlapping control interval of matching length and chromosome was created using bedtools shuffle (100).

IWTomics
IWTomics (Interval-Wise Testing for Omics data; (101)) was employed to compare SNP and FNS frequencies between non-B DNA and control sequences. For each non-B DNA type and each nucleotide substitution frequency, we performed three separate tests: one for each 2-kb flanking sequence and one for the non-B DNA locus itself. IW-Tomics performs a functional permutation test for the null hypothesis that two groups of curves have the same distribution, versus the alternative hypothesis that their distributions differ (two-sided test). Here, the two groups are represented by non-B DNA and control sets, while each curve consists of the substitution counts measured in contiguous nucleotides (for flanking sequences) or in contiguous bins (for non-B DNA loci). Although substitution counts are discrete measurements and each of these curves might look quite noisy, the substitution frequency difference between the two groups--employed in the test statistics of the test--is a rather smooth curve (e.g. Supplementary Figure  S4, bottom panels). This makes functional methods suitable for our analysis and guarantees that IWTomics has good power and accuracy. For each comparison, we based the test on 1,000 permutations and on the test statistic where f 1 (x), f 2 (x) are the substitution frequencies in the two groups and S is the subinterval where the test is performed. We computed an adjusted P-value curve--controlling the interval-wise error rate--for any scale ranging from the single nucleotide to the entire 2000bp interval (for flanking regions) or from the single bin to the entire 180-bin interval (for non-B DNA loci). An example of a complete IWTomics output is shown in Supplementary Figure S4.

Substitution spectrum at the first flanking position
To account for the potential sequence context bias introduced by the non-B DNA motif annotations, SNP and FNS frequencies were computed in their trinucleotide context (by considering a substituted base with its two flanking bases).
For each trinucleotide, we counted total occurrences surrounding non-B DNA loci and their controls, and identified and counted those harbouring a SNP or an FNS in the central base. Next, we added up substitution counts in trinucleotides for each substitution type (for instance, all NAN→NCN were summed as A→C substitutions), and computed substitution frequencies dividing these by the counts of total corresponding trinucleotides. We tested whether substitution frequencies surrounding non-B DNA loci differed from those of their controls using two-sided Fisher's exact tests. P-values were adjusted for multiple testing using a Bonferroni correction.

Building 1-Mb genomic windows
Our multiple regression analysis utilized 1-Mb consecutive non-overlapping windows obtained partitioning the hg19 human genome. We excluded windows that overlapped with gaps in the hg19 assembly using bedtools (100), resulting in 2511 1-Mb autosomal windows. The coverage of NCNR subgenome in each 1-Mb window was computed with bedtools annotate (100).

Genomic features
For each of our 1-Mb windows, we extracted genomic features from various sources (Supplementary Table S2). GC content was computed directly from the genomic DNA sequence. RNA polymerase II and replication origins were measured as 'coverages', computing the proportion of the window covered by the feature using bedtools (100). Replication timing and recombination rates were measured as weighted averages, i.e. the sum of all annotated intervals' length multiplied by their score. DNA methylation coverages (for CHH, CHG and CPG methylation (102)) were also measured as weighted averages, but using a different method implemented through the Galaxy tool 'Assign Weighted Average Values' (103,104). Specifically, a window partially or totally overlapping multiple feature intervals was assigned the average of the values of the features weighted by the corresponding number of overlapping bases. DNA binding profiles of H2AFZ, H3K9me3, H3K9ac, H3K27ac, H3K36me3, H4K20me1, H3K79me2, H3K4me1, as well as CTCF motifs, were collected from the ENCODE project portal (105). The results of ChIP-seq assays performed on hESC H1 cell lines were downloaded (provided on GitHub) and histone modification features were measured as 'signals'--the average number of reads aligned in each window. DNase Hypersensitive Sites were also collected from the ENCODE project portal (105) and measured as 'coverage', i.e. as the proportion of the window covered by DHS peaks. Mappability (106), CpG islands (107), and Lamina Associated Domains (108) tracks were downloaded from the UCSC Genome Browser (91). The distance of each window to the closest centromere and telomere was computed as described in the code on the GitHub. Finally, the telomeric hexamer TTAGGG was annotated on hg19 using fastaRegexFinder (https://github. com/dariober).

Multiple linear regressions
Multiple regression models were built for both SNP frequency and FNS frequency using the NCNR subgenome (see above). Since these frequencies are not reliable in windows with small NCNR coverage (this ranged from 4.8% to 66.6%, Supplementary Figure S5A), we filtered out all windows where this coverage was <20%. Moreover, since SNP and FNS frequencies reliability increases (their variance decreases) with NCNR coverage, we employed weighted least squares--weighing each retained window with its NCNR coverage. As predictors, we used all genomic features listed above--measured as coverage, count, or signal (Supplementary Table S2)--as well as the coverage of each non-B DNA type (Supplementary Table S1), in the same 1-Mb genomic windows. For G4 loci, instead of simple coverage, we used a weighted coverage, computed weighing each G4 locus in a 1-Mb window with its predicted stability score provided by Quadron (98). Windows with average mappability score <0.8 (Supplementary Figure S5B) were also filtered out. Finally, we identified and filtered out a total of 12 windows that represented strong outliers, which may have been influential points for the regressions (in particular, three windows with no G4s, one with unusually large inverted repeats coverage, one with unusually large mirror repeats coverage, eight windows with very small RNA Pol II coverage, one with very large H2AFZ signal and one with very large H3K79me2 signal). A total of 2203 windows were retained for the regression analysis. Several genomic features, as well as SNP frequency, were log-transformed to obtain more symmetric distributions, while other variables were scaled (see Supplementary Table S3 and the GitHub for details). Telomere hexamer coverage was binarized (0 when the coverage was 0 and 1 when it was >0). In order to avoid strong multicollinearity in the models, we used hierarchical clustering (with 1 -|Spearman's correlation| as dissimilarity and complete linkage) to group predictors (Supplementary Figure S6A). Predictors with Spearman's correlations higher than 0.8 in absolute value were clustered, and only one predictor from each cluster was included in the regression models (Supplementary Figure S6). A total of 15 genomic features and all six non-B DNA predictors were retained through this exercise (using Pearson's correlation produced similar clustering; Supplementary Figure  S6B). Mappability was not included in the clustering, but it was forced to be in the regression models as a predictor, in order to control for it. A two-step procedure was then employed to include the relevant quadratic terms in the models and perform variable selection. First, for each variable (all non-B DNA and genomic features except telomere hexamer, which was binarized), we compared a weighted linear model comprising the variable and mappability as predictors, with a weighted quadratic model comprising the variable, its squared values, and mappability as predictors. The quadratic term was retained if the P-value of the ANOVA test comparing the linear and quadratic models was <10 −10 . Second, we employed a weighted regression model with Elastic Net regularization (109) to select relevant (linear and/or quadratic) terms. The elastic net mixing parameter was set to 0.5 (equal mixture of Lasso and Ridge penalties). The regularization parameter was selected through repeated 10-fold crossvalidation (five repetitions) to minimize mean squared error (we obtained the same results maximizing R-squared).
Mappability was forced to be in the final models, hence it was not included in the Elastic Net penalty. Finally, we refitted the final models comprising the terms selected by the Elastic Net and mappability (Supplementary Table S3), we evaluated their R-squared and we computed the coefficient of partial determination of each selected variable (non-B DNA and genomic feature) as where R 2 f ull is the R-squared of the full model and R 2 red is the R-squared of the reduced model obtained by removing all linear and/or quadratic terms involving the variable.

RESULTS
To study the effects of non-B DNA loci on regional variation in mutation rates, we used sequence-motif-based annotations of Z-DNA loci, and of direct, inverted, mirror, and A-phased repeats (Supplementary Table S1) in the hg19 version of the human genome available in the non-B DNA DataBase (99). G4 loci (Supplementary Table S1) were annotated with Quadron (98)--which, compared with other G4 annotation tools (110), has the advantage of assigning a predicted stability score to each G4 locus by using a machine-learning algorithm trained on the experimental output of G4-seq (111).
To minimize the effects of selection, regional variation in mutation rates was assessed by evaluating the distribution of nucleotide substitutions in putatively neutrally evolving non-coding regions of the genome (see Materials and Methods). From these non-coding regions, we additionally excluded repetitive elements as annotated by RepeatMasker (90) because of difficulties in sequencing, read mapping, and subsequent variant calling, and because of the demonstrated role of some non-B DNA motifs in the life cycle of transposable elements (43). This resulted in ∼900 Mb of the Non-Coding Non-Repetitive (NCNR) subgenome (i.e. approximately one third of the genome). We assumed that these sequences evolve largely neutrally (14,88,89,112) and thus, the levels of diversity and divergence reflect variation in mutation frequencies. Within the NCNR subgenome, we were able to analyze 171 653 direct repeats, 2 017 399 inverted repeats, 213 899 mirror repeats, 146 174 A-phased repeats, 108 279 Z-DNA loci and 178 312 G4 loci (Supplementary Table S1).
Always within the NCNR subgenome, we gathered nucleotide substitutions, the most common type of genetic variants (2), from two independent datasets. First, we considered the frequency of FNSs (the number of fixed differences per site per locus) obtained from whole-genome alignments of human and orangutan genomes (96). Second, we considered the frequency of SNPs (the number of SNPs per site per locus) from the Simons Genome Diversity Project (SGDP) (97). This dataset was generated from 279 human genomes sequenced at a relatively high depth (∼30×). We did not use low-depth sequencing datasets, e.g. the 1000 Genomes Project (113), as they have high probability of false SNPs arising from sequencing errors (85).

Non-B DNA loci affect small-scale variation in nucleotide substitution frequency
We studied small-scale variation in nucleotide substitution frequency in non-B DNA loci and their immediate flanking regions (2 kb upstream and 2 kb downstream), as compared to NCNR subgenome control sequences matching the non-B loci in number, length and chromosome of origin, and chosen at random. Previous studies suggested a mutagenic effect of non-B DNA extending into adjacent hundreds of bases (80,114); we increased this distance to ensure we can detect potential effects even at larger distances. Whenever the flanking regions of non-B DNA loci overlapped, only one of the overlapping loci was chosen (at random) for analysis (Supplementary Table S1; see Materials and Methods). To include non-B DNA loci of different lengths (from 11 to 100 bp) in a unified statistical analysis, we scaled each locus to a common length interval split in 180 bins (Supplementary Figure S3). Such scaling was chosen based on the number and maximum length of different sections within non-B DNA loci (e.g. stems and loops in G4s, and repeat arms and potential spacers in inverted, direct, mirror, and A-phased repeats; see Materials and Methods). After scaling the loci, we aggregated nucleotide substitutions per bin for each non-B DNA type genome-wide, and computed the corresponding fixed and polymorphic nucleotide substitution frequencies dividing the respective per-site nucleotide substitution counts by the number of loci. Differences in polymorphic and fixed nucleotide substitution frequencies between non-B DNA loci and control sequences were evaluated using IW-Tomics (115), a functional data analysis method. IWTomics compares two groups of curves (e.g. the substitution frequencies along non-B DNA loci, or their up-and downstream flanking regions versus the respective control sequences), performs permutation tests on all possible subintervals, and produces adjusted P-value curves for any possible subinterval length (from the single bin to the whole 180 bin interval for non-B loci, and from the single nucleotide to the whole 2000-bp interval for each flanking region). This approach identifies locations of significant differences between the two groups of curves.
G4 loci. We studied variation in nucleotide substitution frequency within G4 loci, separating each locus into its stems (consisting of guanines) and loops (consisting of any nucleotides). In order to more easily parse patterns in stems and loops, we focused on G4s with four stems and three loops. These are the most abundant in the genome; they constitute 49.7% of the G4s we considered in the NCNR subgenome. We found that G4 loci exhibited significantly elevated polymorphic and fixed nucleotide substitution frequencies compared to control sequences (Figure 2A and B). This elevation was evident in both stems and loops, but was particularly striking in loops, which had up to 3.2-fold SNP frequency increase and up to 1.7-fold FNS frequency increase compared to control sequences.
We next divided G4 loci into stable and unstable based on their predicted stability. Following published recommendations (98), we labeled loci with Quadron score >19 as stable (a total of 20 156 loci) and loci with Quadron score ≤19 as unstable (a total of 30 212 loci). We found that SNP and FNS frequencies were significantly higher at stable versus unstable G4 loci, and again this difference was more pronounced in loops (up to 5.1-and 1.8-fold for SNP and FNS frequencies, respectively) than in stems (up to 2.1-and 1.3fold for SNP and FNS frequencies, respectively; Figure 2C and D). To test whether mutations at G4 loci tend to stabilize or destabilize G4 structures, we analyzed the effects of rare variants on G4 stability (we considered SNPs with derived allele frequency ≤0.01, excluding those supported by 1-2 reads only to minimize the effect of sequencing errors). Because genetic drift dominates over natural selection when allele frequencies are low (116), the distribution of rare variants is likely to reflect mutational processes rather than natural selection. Interestingly, we found that, on average, there are as many mutations leading to an increase as there are mutations leading to a decrease in predicted G4 stability. Specifically, we substituted positions in the reference human genome by the corresponding rare variants, creating an alternative human genome, and reran G4 annotation with Quadron (98). A total of 609 013 G4s were present at the same or similar coordinates (requiring ≥90% overlap) in both alternative and reference human genomes. Among these shared G4s, the proportion predicted to change from stable to unstable (201/(201 + 179)) was not significantly different from proportion predicted to change from unstable to stable (179/(201 + 179); P = 0.281, 1-sample test with continuity correction for proportion of G4s with increased stability equal to 50%), based on the Quadron stability threshold of 19 (98).
Other non-B DNA loci. Other types of non-B DNA loci also affected variation in both SNP and FNS frequencies in NCNR (Figures 3-4), with patterns that were specific to the different types of loci. For inverted repeats ( Figures  3A and 4A), SNP and FNS frequencies were both significantly elevated in the spacer (up to 1.3-and 1.2-fold, respectively), but significantly depressed in most repeat arms (up to 1.2-fold for both), as compared to control sequences (the comparisons in the rest of the paragraph are also against control sequences). For direct repeats ( Figures 3B and 4B), SNP and FNS frequencies were both significantly elevated in the spacer (up to 3.0-and 1.8-fold, respectively); however, SNP frequency was significantly elevated (up to 1.4-fold), whereas FNS frequency was significantly depressed (up to 1.5-fold), in the repeat arms. For mirror repeats ( Figures  3C and 4C), both SNP and FNS frequencies were significantly elevated in the spacer (up to 1.4-and 1.2-fold, respectively); however, SNP frequency was significantly elevated (up to 1.3-fold), whereas FNS frequency was significantly depressed (up to 1.2-fold), in the repeat arms. For Z-DNA ( Figures 3D and 4D), both SNP and FNS frequencies were significantly elevated (up to 3.3-and 2.0-fold, respectively). A-phased repeats can be separated into A-tracts, which are required for structure formation, and spacers, which can vary in sequence without affecting structure formation. We focused on those with three A-tracts and two spacers (the most abundant in the genome). Both SNP and FNS frequency in repeat arms were significantly depressed (up to 1.7-and 1.6-fold lower, respectively; Figures 3E and 4E). In spacers, SNP and FNS frequencies were both significantly elevated (up to 1.2-and 1.3-fold, respectively; Figures 3E  and 4E).
We also discovered that spacer length influences variation in SNP frequencies for inverted and mirror repeats (Supplementary Figure S7). Namely, inverted and mirror repeats with <15-bp-long spacers exhibited a more pronounced el-evation in SNP frequency. This result suggests that inverted and mirror repeats with short spacers have a higher probability of forming stable structures as compared to the ones with long spacers.

Nucleotide substitution frequencies at flanking sequences of non-B DNA loci
For most types of non-B DNA loci, the alteration in nucleotide substitution frequencies continued into their flanking regions. Even though we analyzed 2 kb up-and downstream of flanking regions, most of the effect was observed in the first 100 bp, and thus we present them in Figures 2-4 (full 2-kb flanking sequences are presented in Supplementary Figure S8). For all NCNR G4 loci (with four stems and three loops), we measured significantly elevated SNP and FNS frequencies extending into the flanking regions. Contiguous regions of significant elevation extended into the neighboring 3-6 bp, and more scattered regions of smaller, but still significant, elevation existed up to ∼500 bp. However, this flanking region elevation was smaller than that within the loci themselves and was the highest at the immediate first flanking nucleotide position upstream and downstream. Indeed, the first flanking nucleotide position experienced up to 2.2-and 1.5-fold increase in SNP and FNS frequency, respectively, compared to control sequences ( Figure  2A and B). In contrast, the other flanking nucleotides analyzed experienced only up to 1.4-and 1.3-fold increase in SNP and FNS frequency, respectively, compared to control sequences (Figure 2A and B). Importantly, the elevation in SNP and FNS frequency was more pronounced in the flanking regions of stable than unstable G4 loci (Supplementary Figure S9A-D).
In the flanking regions of direct repeats, SNP frequency was significantly elevated (up to 2.5-fold), although gradually decreasing over a distance of ∼450 bp ( Figure 3B), whereas FNS frequency was elevated at the first flanking position only (up to 1.7-fold; Figure 4B), as compared with control sequences. Both SNP and FNS frequencies were significantly elevated (up to 1.4-and 1.1-fold, respectively) in the flanking regions of mirror repeats, mainly in the first few (one to eight) flanking nucleotides ( Figures 3C and 4C). The influence of inverted repeats on SNP and FNS frequency did not extend into their flanking regions (Figures 3A and  4A). The elevated SNP and FNS frequencies at Z-DNA loci extended over a distance of ∼175 and ∼12 bp, respectively ( Figures 3D and 4D), in the surrounding flanking regions. In the flanking regions of A-phased repeats, the decrease in SNP and FNS frequency was concentrated at the first flanking nucleotide (Figures 3E and 4E).
SNP frequencies at the first flanking nucleotide positions (immediately adjacent to non-B motif annotations), both up-and downstream, were particularly high for all types of non-B DNA loci except for inverted repeats (Figures  2-4). In an effort to explain this observation, we asked whether this increase was driven by a particular substitution type (A→T, A→C, etc.), thus skewing the substitution spectrum. To investigate this possibility, we compared the frequencies of different substitutions at the first flanking nucleotide position, separately for each type of non-B DNA.   Only the frequencies of trinucleotides present at the immediate flanking positions of stable G4 loci were compared with those present at control sequences (trinucleotides present only in control sequences were not considered). A correction for the trinucleotide context was applied (see Materials and Methods). Two-sided Fisher's exact test was used to evaluate significant differences, and P-values were adjusted for multiple testing using Bonferroni correction. An asterisk (*) marks significant differences between G4 and control sequences (adjusted P-value < 0.05).
The substitution spectra were indeed different for the first flanking positions of G4 loci as compared to control sequences ( Figure 5; Supplementary Tables S4A and B). Because G4 structures form only on one DNA strand, they might skew the substitution spectrum at the first flanking base in a strand-specific fashion. To account for this, we studied the substitution spectrum at G4 loci annotated on the reference strand only (i.e. restricting G4 motifs to the G-rich strand and not considering C-rich strand motifs; G-quadruplexes annotated on the non-reference strand were analyzed in Supplementary Figure S10, and the results were similar). The requirement of the first and the last bases at G4 loci to be a guanine might lead to a nucleotide substitution spectrum bias at the first flanking positions. For instance, a cytosine in the upstream flanking base followed by a guanine in the annotated G4 can form a CpG dinucleotide, which, if methylated, has elevated rates of CpG→TpG substitutions (117). To correct for this potential bias at the immediate flanking bases of G4 loci, we computed nucleotide substitution frequencies in their trinucleotide context, i.e. the ratio of the number of occurrences of a trinucleotide with the mutated middle base and the total number of occurrences of this trinucleotide. This procedure was performed separately for G4 loci (with four stems and three loops and located in the NCNR subgenome) and control sequences, and effectively corrected for the differences in trinucleotide composition between these two groups of sequences. After applying such corrections, we observed that the T→G substitution frequency was significantly elevated at the first up-and downstream flanking positions of stable G4 loci (odds ratios 14.04 and 8.51, Bonferronicorrected P-values 1.09 × 10 −30 and 1.60 × 10 −13 , respec-tively, Fisher's exact test; Figure 5 and Supplementary Table  S4A). The A→G substitution frequency was also significantly elevated at the first up-and downstream positions of stable G4 loci, albeit to a smaller extent (odds ratios 3.61 and 2.66, Bonferroni-corrected P-values 1.34 × 10 −6 and 3.97 × 10 −3 , respectively, Fisher's exact test). The substitution frequencies at the first flanking positions of unstable G4 loci were not significantly different from those at control sequences (Supplementary Table S4B). These results suggest that the first flanking nucleotides of stable G4s have a tendency to acquire additional guanines, leading to G4 elongation.
Whereas other types of non-B DNA we examined also displayed differences in their nucleotide substitution frequencies at the first flanking positions compared to those of control sequences (Supplementary Table S4C-G), none demonstrated a mutation spectrum as skewed as that observed in G4s. In Z-DNA, all substitution frequencies were elevated except for C→T and G→A (Supplementary Table S4F). All substitution frequencies were evenly (2-3-fold) elevated at both immediate flanking positions of direct repeats (Supplementary Table S4D).

Non-B DNA loci explain a substantial portion of large-scale variation in nucleotide substitution frequencies
To evaluate whether non-B DNA loci contribute to largescale regional variation in nucleotide substitution frequencies across the genome, we studied how divergence, as measured with human-orangutan FNS frequency, and diversity, as measured with human SNP frequency from the SGDP (97), vary with the coverage of each of the six non-B DNA types measured in 1-Mb genomic windows (see Materials and Methods; for G4 loci, we computed a weighted average of per-window coverage, weighting with predicted stability as provided by Quadron scores (98)). Nucleotide substitution frequencies and non-B DNA coverages were measured in the NCNR portion of each 1-Mb window. Non-linear (quadratic) relationships were observed between nucleotide substitution frequencies and coverage of each non-B DNA type studied (Figures 6 and Supplementary Figure S9); however, there were some differences among non-B DNA types. When considering divergence, G4, direct repeat and Z-DNA coverage exhibited an almost flat trend followed by positive associations, Aphased repeat and inverted repeat coverage exhibited negative associations followed by a flat trend, whereas mirror repeat coverage exhibited a negative trend followed by a positive trend, depending on their values ( Figure 6). Similar pat-terns were observed considering diversity (Supplementary Figure S11).
We next examined how non-B DNA contributes to regional variation in nucleotide substitution frequencies in the context of 26 additional genomic features that were previously shown to affect such variation (Supplementary Table S2). To minimize multicollinearity, the overall set of potential predictors was reduced by clustering features based on pairwise Spearman's correlations and retaining only one feature per cluster (Supplementary Figure S6), similar to the procedure used in (118,119). We retained all six non-B DNA types and 15 additional genomic features. Next, we fitted two multiple regressions (separately for divergence and diversity) on this reduced set of predictors and, when needed, their quadratic terms, performing further variable selection with the Elastic Net and measuring the contribution of each retained predictor with a coefficient of partial determina-Nucleic Acids Research, 2021, Vol. 49, No. 3 1509 tion (partial R 2 ; see Materials and Methods). Average mappability score (106) was included as a predictor in all three models to control for its effect.
The resulting models explained 58.0% and 55.6% of variation in divergence and diversity, respectively (Table 1 and  Supplementary Table S3). Importantly, all types of non-B DNA were retained as significant predictors in both models. Among different types of non-B DNA, Z-DNA, direct repeats, and G4 loci were the strongest predictors of divergence (explaining 2.32%, 1.63% and 1.04%, respectively), whereas direct repeats and G4 loci were the strongest predictors of diversity (explaining 8.01% and 1.08% of variation, respectively). Even when considered in the context of other relevant genomic features, as a group, different types of non-B DNA explained as much as 15.94% and 19.95% of variation not explained by the other features in the two regressions, respectively. If we were to exclude non-B DNA, our models would explain 50.0% and 44.5% of variation in divergence and diversity, respectively (note that partial R 2 values usually do not add up because of correlations among predictors). In our full models (Table 1), among other genomic features, distance to telomere, replication timing, H3K9me3, and lamina-associated domains were strong predictors of divergence (explaining 15.3%, 2.08%, 1.72% and 1.35% of variation, respectively), whereas distance to telomere, recombination rate, distance to centromere, H3K9me3, and DNase hypersensitive sites were strong predictors of diversity (explaining 14.2%, 6.12%, 3.45%, 1.52% and 1.08% of variation, respectively). Notably, the contribution of all non-B DNA types taken as a group exceeded that of any other genomic feature in both divergence and diversity regressions. Thus, our results strongly suggest that non-B DNA is a key contributor to large-scale regional variation in nucleotide substitution frequencies.

DISCUSSION
Harnessing publicly available datasets that span different evolutionary times (human-orangutan divergence and human diversity) and sophisticated statistical methods, our genome-wide study illuminates the relationship between non-B DNA loci and variation in nucleotide substitution frequencies. Using IWTomics, we delineated distinct patterns of small-scale variation in and around different types of non-B DNA loci. Using multiple regression analysis, we determined the contribution of non-B DNA to explaining regional variation in nucleotide substitution frequencies at the 1-Mb scale. Our results compellingly argue for a pivotal role of non-B DNA in shaping the mutation dynamics of the genome, which has important biological implications.
On the one hand, our genome-wide study supports a 'unifying hypothesis', whereby non-B DNA plays a role linking genomic architecture and local DNA environment with mutations that underlie human inherited disease (120). Thus, future studies should consider the potential of a locus to form non-B DNA when evaluating candidate genetic variants for human genetic diseases. On the other hand, because of its high mutability, non-B DNA likely represents an inexhaustible source of novel genetic variation in natural populations. In agreement with this, a recent study in stickleback Table 1. Coefficients of partial determination of regional variation in nucleotide substitution and the proportion of variation explained by each model ( fish found that Z-DNA forming in an enhancer of the Pitx1 gene increases the probability of deletions that lead to the loss of pelvic hindfins, and are repeatedly utilized for adaptation to a freshwater environment (3). Based on our findings, we expect that many analogous examples of the use and re-use of non-B DNA in adaptation will be discovered in the near future.

Small-scale variation in nucleotide substitution frequency
We have shown that non-B DNA loci often have nucleotide substitution frequencies significantly different from those characterizing B DNA (usually higher). This corroborates an earlier report of increased SNP density at non-B DNA in the 1000 Genomes Project (24), even though this study might have been biased by a potentially increased Illumina sequencing error rate at non-B DNA that is expected to elevate false positives for genomes sequenced at low depth (85). The present study does not suffer from this limitation as it uses deeply sequenced genomes. The altered nucleotide substitution frequencies at non-B DNA loci we detected here are in agreement with another study from our group, in which we found that PacBio sequencing errors are affected by non-B DNA (85). Thus, non-B DNA loci likely affect both errors in the PacBio sequencer and mutations in germline cells. Other recent studies found non-B DNA loci to be mutation hotspots in cancer genomes (23,64), suggesting that such loci also drive mutagenesis in somatic cells.
G4 loci. Aggregating variants at tens of thousands of G4 loci, we observed an increase in both diversity and divergence, consistent with previous reports of elevated 1000 Genome Project SNP frequency at G4 loci (24) and of G4 motif enrichment at cancer mutations (23). Potential mechanisms underlying heightened mutagenesis at G4 loci include DNA polymerase errors during replication or repair, and/or damage-induced mutations. Importantly, these mechanisms are not mutually exclusive. Due to their G-rich sequence, G4 motifs are subject to oxidative DNA damage (reviewed in (121)), which, if unrepaired, can lead to mutations. Eukaryotic replicative DNA polymerases are inhibited at some G4 structures (reviewed in (122)). Possibly, G4s compromise the high fidelity of replicative DNA polymerases, leading to errors during DNA synthesis. This hypothesis is based on our previously published genome-wide results showing a direct link between polymerization errors during PacBio sequencing and germline mutations at G4 loci (85). Indeed, we found levels of divergence and diversity to be negatively correlated with polymerization speed and accuracy in the sequencer at such loci (85). Alternatively, error-prone polymerases may be engaged at a slowed or stalled fork to ensure complete G4 replication. This mechanism is supported by evidence suggesting that specialized polymerases eta, kappa, and Rev1 are important for G4 stability (122). We demonstrated that stable G4 loci have higher nucleotide substitution frequencies as compared to unstable G4 loci ( Figure 2C and D). Because stability reflects the probability and strength of structure formation, this result suggests that stable G4 structures are a more serious and/or more common obstacle for polymerization progression leading to increased errors, which frequently result in mutations. In support of this claim, we previously found positive relationships between polymerization slowdown and error rates, and between polymerization slowdown and G4 stability, as measured by circular dichroism (85). Potential mechanisms include replication fork restart at G4 motifs by PrimPol, an error prone polymerase (123), or fork degradation and breakage followed by polymerase theta-mediated, error-prone DNA end joining (124). The effect of stability on polymerization accuracy was also found for other types of non-B DNA--Z-DNA and triplex DNA (125).
Stable and unstable G4 loci differed more in divergence comparisons than in diversity comparisons, suggesting that the stability of G4 structures evolves over time. While the classification of individual G4 loci into stable versus unstable likely remains valid for human populations, there is a probability that extrapolating the predicted stability of some human G4 loci to orangutan loci might be incorrect, due to accumulated mutations in the latter. Thus, for instance, some G4 loci annotated as stable in the human genome might be unstable in the orangutan genome. As a consequence, they may be a weaker obstacle for polymerization or accrue less oxidative damage, leading to fewer mutations, and resulting in lower fixed nucleotide substitution frequencies.
Within G4 loci, we found that nucleotide substitution frequencies were more elevated in loops, which can vary without impeding G4 structures, than in guanine stems, which are critical for G4 formation. This suggests that stems and loops differ in mutability and/or selective pressure. Guanines located in both stem tetrads and loop sequences are subject to oxidative damage when folded into quadruplex structures, and the 5 guanine within a tetrad can be a hotspot of damage (126). If G4 stems had low mutability, we would expect guanines to harbor few polymerization errors. However, this was not the case for the polymerase used in PacBio sequencing, which exhibited many polymerase errors at guanines (85). Therefore, the precise mechanisms leading to the observed differences in stem versus loop substitutions remain to be elucidated.
To assess potential selection in stems, we tested whether the site frequency spectra differed between stems and loops in the SGDP SNPs (97). Although the two distributions were similar visually, we observed a significantly higher proportion of SNPs with lower minor allele frequency in stems than in loops (P-value < 2.2e-16, Kolmogorov-Smirnov test; Supplementary Figure S12). Assuming the latter evolve neutrally, this suggests purifying selection acting on the former. Such an observation is consistent with a study in S. cerevisiae, proposing that stems are conserved to safeguard G4 structure formation (127), and contradicts another study, proposing that the underrepresentation of SNPs in G4 guanine tracts has neutral causes (128). Notably, we found evidence of selection acting on G4 stems even though our analyses focused on G4 loci found in the NCNR subgenome, which ought to have a high probability of evolving neutrally. We also found the centers of G4 stems to be more conserved than their edges ( Figure  2), consistent with previous reports of the guanine at the center of stems being the most important nucleotide for the formation of a G4 structure (128)(129)(130). In summary, our results suggest that even though they drive mutagenesis and genome instability interfering with DNA replication (131,132), and even when they are located in the NCNR subgenome, some stable G4 loci might be beneficial for other essential processes and thus preserved in the genome.
A-phased repeats. In contrast to G4 loci, A-phased repeats had depressed overall nucleotide substitution frequencies. This is consistent with slightly lowered rates of polymerization errors observed for A-phased repeats during PacBio sequencing (85). Interestingly, these loci showed a periodic pattern that was in some ways similar to that of G4 loci. In our analysis, nucleotide substitution frequencies were depressed in A-tracts but similar to controls in spacers. Lowered nucleotide substitution frequencies in A-tracts may be explained in part by the previously reported low mutability of adenines (1,4,5); PacBio polymerization error rates at adenines were also reported to be low (85).
Direct, inverted, and mirror repeats. We previously reported elevated PacBio sequencing errors at mirror, direct, and (to a lower extent) inverted repeats (85). Elevated SNP frequencies (1000 Genomes Project, (24)) and cancer somatic mutations (23) were also reported for these loci. Intriguingly, we observed both increases and decreases in nucleotide substitution frequencies at mirror, inverted and direct repeats. For all three repeat types, frequencies were Nucleic Acids Research, 2021, Vol. 49, No. 3 1511 higher in spacers than in repeat arms. This is reminiscent of observations for cancer mutations (23), and consistent with previous experimental studies demonstrating that spacers are more mutable in hairpin structures (84,133). We also found that nucleotide substitution patterns within repeat arms differed among repeat types. (84,133). Non-B DNA loci composed of repeats can pose an obstacle to DNA replicative polymerases, and the observed nucleotide substitution frequencies could reflect processing by low-fidelity, specialized DNA polymerases (65,66,(134)(135)(136) gene conversion and increased DNA damage. Gene conversion is common at cruciform and slipped-strand structures (137), and could explain the decreased nucleotide substitution frequencies in repeat arms of inverted and direct repeats--particularly for fixed nucleotide substitutions, where multiple rounds of gene conversion might have taken place. Oxidative DNA damage was proposed as an important mechanism of increased mutagenesis at mirror repeats (125,138).
Z-DNA. At Z-DNA loci, the elevated diversity and divergence we detected are consistent with previous reports of increased germline and somatic mutations (23,24), but contradict our previous study, in which such loci showed decreased levels of PacBio polymerization errors (85). A potential resolution of this contradiction might lie in the importance of compromised DNA repair (125,138) and of environmental factors (139,140) for Z-DNA mutagenesis. We found that the edges of Z-DNA loci were particularly enriched in nucleotide substitutions ( Figures 3D and 4D), potentially because they coincide with the B-Z junctions (125,138).

Nucleotide substitution frequencies in regions flanking non-B DNA loci
Our results suggest that the effect of most types of non-B DNA loci on diversity and divergence is not restricted to their sequences but extends into their flanking regions--albeit at a smaller magnitude and limited range. We found an elevation in substitution frequencies extending up to ∼500 bp in the flanking regions of G4 loci, which is further than the range of elevated indel frequencies previously reported for humans (23) and C. elegans (80) (150 and 200 bp, respectively). We also found substitution frequencies to be elevated up to ∼300 bp from direct repeats, which is further than the 50 bp previously reported for microsatellites, a similar type of loci (141). While Z-DNA and H-DNA-forming mirror repeats have been reported to increase mutation frequencies in their neighboring regions (125), we found that this effect is limited to at most ∼100 bp and 4 bp, respectively. Several molecular mechanisms have been proposed to explain the influence of non-B DNA loci on nucleotide substitutions in their flanking regions. Double-strand breaks associated with non-B DNA structures and the resulting DNA repair can be mutagenic (124,138). Alternatively, folding into non-B DNA structures may expose downstream sequences to oxidative damage through long-range hole migration, creating hotspots for mutagenesis in the surrounding sequences (125).
We also found a marked elevation in the nucleotide substitution frequencies at the first flanking nucleotides of some types of non-B DNA loci. Diversity and divergence were particularly elevated at the first flanking nucleotides of stable G4 loci, and showed a skewed substitution spectrum favoring T→G and A→G substitutions. This specificity is intriguing, and the mechanism(s) responsible are not clear at this time. Oxidative damage to dGTP precursor pools and subsequent incorporation of 8-oxo-dGTP opposite template A causes T→G substitutions (69). Possibly, the alternative folding of G4 distorts the immediate flanking template bases, and promotes this misincorporation by DNA polymerases during replication. Alternatively, the skewed spectrum may reflect the error specificity of specialized polymerases that are engaged at G4 motifs (122). Among the nuclear DNA polymerases thus far identified for this function, only DNA Pol eta produces T→G and A→G substitutions at measurable frequencies (142). Notably, these substitutions result in the elongation of G4 loci, which should further increase their stability--an observation that is in line with the potential functionality of some G4 loci in the NCNR portion of the human genome. The first flanking nucleotide of Z-DNA loci also showed a skewed substitution spectrum, but this was not the case for other non-B DNA loci considered in our study. These results echo findings of Bacolla et al. (143), who found that SNP frequency and spectrum at the first flanking nucleotides of mononucleotide microsatellites is affected by microsatellite repeat identity and length.

Large-scale variation in nucleotide substitution frequency
Our multiple regression models indicate that non-B DNA loci across the human genome are prominent contributors to regional variation in divergence and diversity measured at the 1-Mb scale. Each type of non-B DNA loci contributed significantly, with particularly strong contributions of direct repeats and G4 loci to variation in diversity, and of direct repeats and Z-DNA loci to variation in divergence (Table 1). Although we analyzed the NCNR genome, many of the studied non-B DNA loci do have a repetitive nature and could elevate mutation rates up to 10 kilobases away from the repetitive tract via repeat-induced mutagenesis (144). When taken together, non-B DNA loci explained ∼20% of the variation in divergence and ∼16% of the variation in diversity--greater than that of any previously reported predictors, which we confirmed to have substantial contributions--e.g. distance to telomeres, distance to centromeres (89,145), recombination rate (89,145,146), replication timing (9,17), histone modification H3K9me3 (12), and DNase hypersensitive sites (11,112,147).
Though comparisons are hindered by differences in methodologies, data sources, and overall sets of predictors considered, including non-B DNA does appear to increase the proportion of regional variation explained in our models as compared to that in prior modeling efforts. We explained 58.0% of variation in human-orangutan divergence, as compared to 52.6% and 52% obtained in prior models for human-chimpanzee (145) and human-macaque (89) divergence, respectively. We explained 55.6% of variation in SNP frequency from SGDP (97), as compared to 35% obtained in prior models of SNP frequency from dbSNP (12).

Limitations of the study and future directions
At any given time, only a subset of annotated non-B DNA loci are expected to fold into non-B DNA structures in vivo (148). The preferential conformation of a motif depends on the precise sequence and stability of B versus non-B DNA (35) and on conditions within the cell. Supercoiling, transcriptional activity, and conditions in the nucleus (e.g. ionic concentrations) can impact transitions between B and non-B DNA conformations at a locus. Due to this transient nature, the effects of non-B DNA structures at individual loci are difficult to study and must be aggregated into genomewide trends. Thus, the results of our study cannot be extrapolated to individual loci. At present, it remains challenging to identify individual loci that form non-B DNA structures at a given time in vivo. However, new methods such as permanganate sequencing (148), kethoxal-assisted single-stranded DNA sequencing (149), and G4-ChIP-seq (36,150) will increase our ability to study the effects of individual non-B loci on nucleotide substitution frequencies in the near future.
While including non-B DNA loci as predictors in our models of large-scale variation in substitution frequencies provided a substantial improvement, we did not evaluate the interaction terms. Including interactions between predictors will complicate the analysis, but may bolster our predictions and further increase the amount of explained variation. For instance, the distribution of non-B DNA is nonuniform along the genome and may correlate with other genomic features used in the models (e.g. recombination rate; note that multicollinearity in our models was controlled by clustering and eliminating some correlated predictors at the outset). Co-occurrence of non-B DNA loci and other features creates an opportunity for joint effects on nucleotide substitution frequencies, which ought to be evaluated in future studies.
Despite our best efforts to restrict our study to neutrally evolving regions of the genome by limiting it to the NCNR subgenome, we captured footprints of purifying selection at G4 loci. Because we are aggregating multiple loci, it may be challenging to distinguish local changes in mutability from selective pressure. However, the possibility of natural selection acting upon some loci previously considered as neutral is an exciting prospect; the potential functions of such loci should be characterized in future work. In fact, natural selection at some non-B DNA loci has been reported previously (24,151)--but the issue remains underexplored to date.
Finally, our study offers only glimpses into the molecular mechanisms underlying our findings. Our results are consistent with DNA polymerase errors being an important mechanism behind increased mutagenesis at non-B DNA loci, and it will therefore be critical to deepen the investigation through further genome-wide studies. Additionally, future work evaluating the role of DNA repair and the effect of single-strandedness on oxidative damage-induced mutagenesis genome-wide will be important for a complete understanding of the mutability of non-B DNA loci and their flanking sequences.

DATA AVAILABILITY
All computational tools used in this study are available as a GitHub repository at the following URL: https://github. com/makovalab-psu/nonB-RegVar.

SUPPLEMENTARY DATA
Supplementary Data are available at NAR Online.