Genomic landscape of CpG rich elements in human.

BACKGROUND
The studies on CpG islands (CGI) and Alu elements functions, evolution, and distribution in the genome started since the discovery in nineteen eighties (1981, 1986, correspondingly). Their highly skewed genome wide distribution implies the non-random retrotransposition pattern. Besides CGIs in gene promoters, CGIs clusters were observed in the homeobox gene regions and in the macrosatellites, but the whole picture of their distribution specifics was not grasped. Attempts to identify any causative features upon their (genome wide) distribution, such as the DNA context mediated preferred insertion sites of Alu repeats, have been made to ascribe their clusters location.


METHODS
Recent emergence of high resolution 3D map of human genome allowed segregating the genome into the large scale chromatin domains of naturally observable nuclear subcompartments, or Topologically Associated Domains (TADs), designated by spatial chromatin distribution. We utilized the chromatin map to elucidate relations between large scale chromatin state and CpG rich elements landscape. In the course of analysis it was confirmed that genes, Alu and CGI clusters maintain obvious, albeit different in strength, preference for open chromatin. For the first time it was clearly shown that the clusters density of the Alu and CGIs monotonically depend on the chromatin accessibility rate. In particular, the highest density of these elements is found in A1 euchromatin regions characterized by a high density of small length genes replicating in the early S-phase. It implies that these elements mediate (CGIs) or are a side element (Alus) of chromatin accessibility.


RESULTS
We elucidated that both methylated and non-methylated CGIs display the affinity to chromatin accessibility. As a part of comparative genomics section, we elucidated that the dog's genome non-canonical structure, outstanding in mammals for its high CGIs abundance compared to gene number, is explained by the presence of dense tandem CGI extended hotspots (500 kb on average) in subtelomeric and pericentromeric regions with highly skewed CG content, and not by CGIs global distribution pattern shift.


CONCLUSIONS
The study underlines the close association of CG-rich elements distribution with the newly introduced large scale chromatin state map, proposing a refined standpoint on interrelation of aforementioned genome elements and the chromatin state. To our expertise, the TAD-associated partition model employed in the study is likely the most substantial one regarding CpG rich clusters distribution among the whole genome chromatin/isochores maps available.


Background
The total number of CpG sites in human genome comprises around 28.3 mln instances [1]. That is less than 1% of genome compared with 4.4% expected given 42% GC content due to C-> T mutation shift following the frequent 5′ cytosine deamination in CG doublet [2]. Aside from randomly distributed CpG sites in mammalians, there are two major expansions of them in primates, namely: Alu retrotransposons [3] and CpG cluster units, CGIs [2]. The share of Alu CpGs was estimated up to 30% as early as in 1993 [4] and now is confirmed to be more than 25% [1], while CGIs account for only 2% of CpG content in human [5].
With that, CGIs proved to represent a highly specific marker for open chromatin [2]. In particular, unmethylated CGI and genes association is reported as early as in 1987 [6].
The average CpG content of CGI is 18%, and Alu is about 3.6% (still higher than 1% genome average), which allows considering both of them CpG-rich elements. Three major Alu families differ significantly in CpG content: from 2.5% on average in AluJ, to 3.3% in the most abundant AluS family, and up 6% in the young AluY sequences. Notably, CGIs and Alu complement in methylation pattern: 70% of CGIs are completely unmethylated, while around 70% of Alus are methylated in all tissues/cell types.
The CGIs play an overwhelming role within promoters of about 70% genes. One of current views of the promoter CpG enrichment implications is the alteration of DNA conformation to non B-DNA (Z-DNA) structure and thus reduced nucleosome affinity [5,[7][8][9][10][11][12]. It was shown that the length of CGI and nucleosome depletion rate significantly correlate: long CGIs are nucleosome free, while small CGIs can attract/position nucleosomes [10]. DXZ4 macrosatellite spanning 45-120 kb [13] and other CGIs in intergenic regions and gene deserts display basal transcription in a majority of cell lines featuring the accessible (open) non-compacted chromatin [13]. This was also confirmed by DHS and CGI colocation, which was also evolutionary conserved [14].
Intriguingly, CpG-rich Alu families tend to distribute their CpG sites at the distance of 31-32 bp, which make them prone to nucleosome binding in mammals [11,15]. It is worth mentioning that CG dinucleotide signal is the only one among all 16 dinucleotides which manifests the significance detected by autocorrelation analysis of nucleosome positioning in mammals [15]. Also, the CGI mediated nucleosome depletion doesn't essentially depend on methylation state [11], but that way they become inaccessible to transcription factors binding and protected from spurious RNA PolII complexes assembly. Besides, there were reports of nucleosome dips around polyadenilation sites (PAS) [16,17].
Analysis of causes of CGIs conservation phenomena during evolution revealed that they evolve in three distinct regimes: a) hypodeaminated (non-methylated promoters); b) methylated with stable CpG content due to Biased Gene Conversion; c) CG-rich exons. There is a small fraction of pseudo-CGIs, which arise due to spurious clustering [18].
The largest CGIs class represents hypomethylated 5′ located CGI promoter regions responsible for PolII based transcription initiation in 70% of genes. Still, there are smaller (comprising around 25% of CGIs total), but nevertheless essential CGI classes within the vicinity of genes that are differentially or constitutively methylated and are presumably involved in tissue and temporal specificity of gene expression regulation. The functional implication of such CGIs is supported by non-random synonymous substitutions inferred from CGI-CDS overlapping instances [19,20].
Based on previous studies, the number of CGIs and protein coding gene numbers are approximately equal in the majority of mammals and correlate chromosome wise [21][22][23]. But there is certainly not one to one correspondence of the CGIs and genes: a range of CGI clusters can be observed in gene deserts, e.g. macrosatellites D4Z4 [24,25], DXZ4 [13] and others. Backwards, there is plenty of tissue specific gene clusters lack of CGIs (e.g. olfactory receptors). For HOX genes, there are multiple CGIs along each gene [20]. Thus, the issue of genes and CGI relation renders further elaboration.
Advances in whole genome epigenetic marks mapping by novel Chip-Seq [26] and Chia-PET [27] techniques have culminated into the elucidation of 3D chromatin conformation. To date it was recognized that the chromatin domains are the ultimate units of chromosome organization [28]. While Chip-seq experiments feaure high resolution maps of 200-400 bp spans usually corresponding to transcription factors and chromatin remodeling binding sites along with the more extent histone modifications areas, the chromatin 3D conformation map maintains large scale domain size of several hundred kbs [28][29][30]. They were named as Topologically Associated Domains (TADs). The TADs preserve essentially cell-type invariant architecture [31], though, since TADs maintain hierarchical structure within, the cell type specific chromatin conformations also take place on a minor scale [29].
Large scale chromatin domains were primarily segregated into Active domains (A type) of early replication timing chromatin, and lamina associated domains (LAD) (B type) conferring heterochromatic domains with late replication [30,32]. The replication timing analysis also identified the timing transition regions (TTR) as a separate domain type [32]. At the same year HiC conformation capture analysis study [31] elaborated A and B domains into 6 classess (A1, A2, B1, B2, B3, B4) based on analysis of several epigenetic markers profiles within observable chromatin contact domains related to the similar nuclear subcompartments [31]. The attempts to adequately segregate the domains into the chromatin state variants accounting for multiple factors are ongoing in chromatin conformation studies [28]. Notable obstacles herein are that certain domains change its state in the course of embryonic development and cell speciation [33]. Currently, 3D chromatin conformation research reached the point of making possible to model chromatin architecture based on primary data of classification of loci into chromatin types and a catalog of the positions of chromatin loops [34,35].
Concerning our study of CpG rich elements and TADS relation, the differential enrichment of A and B TADs with SINEs, such as Alu in human and B1 in mouse, has been reported previously [36]. The extended TADs classes repertoire employed in our study allowed more elaborated quantification of the feature. To our knowledge, no straightforward TADs mediated CGIs content analysis has been made to date.

Genome data sources
We downloaded data set from UCSC genome browser database (genome.ucsc.edu; [37]). We used human genome version hg19; mouse genome version mm10; dog genome version canFam3. For CGIs annotation we used tables named Cpgislandext in all cases. Genes were downloaded from the same genome versions, using refGene tables for human and mouse, and ensGene table for the dog annotation. Transposon locations for human were retrieved from RMSK table (hg19). Genes were defined by the distinct transcription start site.

Correlation significance
Correlation analysis was performed with Pearson correlation coefficient, its significance was ascertained by using Student t test in the form: where r -Pearson correlation coefficient, nsample size (number of bins). Df = n-1.

ANOVA analysis
ANOVA analysis of TADs subcompartment classes for elements density has been performed using XLStat software (www.XLStat.com). To test the deviation significance of the elements density between the TAD subcompartment classes we used Tukey HSD test.

DNA methylation datasets
Methylation profiles for 63 cell lines and one primary liver cell in ENCODE were downloaded from UCSC genome browser HAIB Methyl450 track [38]. The score of the methylation value associated with each CpG site was defined as the beta value (1) (Illumina's Bead Studio software with the Methylation Module v3.2) multiplied by 1000.
where I kmethylation intensity value on the particular CpG site, the sum is genome wide intensity (approximately 450,000 CpGs) CGI methylation score was calculated as the average of inner CpG sites methylation scores presuming CGI methylation homogeneity [39].

DNase Hypersensitive Sites (DHS) set
DHSs are ENCODE elements represented by short (100-200 bp) DNA fragments. They indicate open or accessible chromatin where DNA is not tightly wrapped within a nucleosome, leaving the sequence accessible to DNA-binding proteins [40]. DHSs largely correspond to transcription factor binding sites, chromatin remodelers and other DNA binding proteins sites on DNA. In a database, they are supplied with scores corresponding to intensity rate specific for each cell line considered.
We used 2.8 mln DHS sites compiled in [14] from 112 human samples representing 72 cell types, to characterize 100 kb genome bins by DHS density.

Distribution features of the elements considered
We used 100 kb bins for the global domain wide analysis of human genome. The distribution basic statistics of Alus, DHSs, CGIs, and genes in 100 kb bins are presented in the Table 1.
We ascertained that Alu, transcripts and CGIs feature exponential distribution (see distributions in Additional file 1: Figure S2), which implies significant number of dense clusters [18,41,42], while L1, DHS maintain binomial distributions ( [43]; Additional file 1: Figure S3). To perform DHS vs other elements comparison we split overall 29,381 100 kb euchromatic bins into 292 clusters ordered by overall DHS number spanning from 0 (26 bins) to 341 (1 bin) DHS entries. The number of DHS per bin was treated as a rough scale chromatin accessibility signature attributable for 100 kb genome segments, and was used as a grouping factor for other DNA elements comparison (Alu, transcripts, CGI distributions).

Nuclear subcompartments map of genome (Topologically associated domains, TADs)
In 2014 Cell report [31] a high resolution 3D map of human chromatin in a range of cell lines was reconstructed. In the course of analysis the authors delineated 6 major classes of nuclear subcompartments (TADs), which partition whole genome contact domains volume by chromatin state analysis. Six classes comprise two euchromatic (A1, A2) and 4 heterochromatic (B1, B2, B3, B4) subcompartments [31]. Each class is characterized with specific histone modification profiles [31]. Since these are largescale domains (median length 185 kb; 450 kb on average; [31]) they are compatible with our 100 kb resolution. The overall state of the nuclear subcompartments is presented in Table 2. We excluded B4 class from consideration since it comprised only 25 domains specific for chromosome 19 [31], as well as NA domains.

A1, A2 -open chromatin, B1-B4 -heterochromatin
According to description in [31], "A" euchromatin segment features enrichment in open chromatin histone marks: H3K36me3, H3K79me2, H3K27ac, and H3K4me1. A1 chromatin state completes replication at early S phase, while A2 proceeds replicating up to the mid -S phase. A2 is enriched with H3K9me3 and contains longer genes. Subcompartment B1 is enriched with H3K27me3 and depleted of H3K36me3 marks, representing facultative heterochromatin. Subcompartment B2 includes the majority of pericentromeric heterochromatin and is enriched at the nuclear lamina and at Nucleolar Associated Domains (NAD). Subcompartment B3 is enriched at the nuclear lamina and is depleted at NADs, thus corresponding to constitutive heterochromatin. B4 comprises a range of marks representing highly ambiguous chromatin pattern. It features strong enrichment for both activating chromatin marks, such as H3K36me3, H3K4me3, H3K27ac, and heterochromatin-associated marks, such as H3K9me3 and H4K20me3. B4 contains 130 of the 278 KRAB-ZNF genes in the genome which is highly non-random [31]. More details on the subject could be ascertained from the original work [31].

CGIs association with genes
The chromosome wise dependence of CGIs and genes was reported previously for a range of species genomes, including dog genome featured with highly skewed CGI to gene ratio due to abundance of CGIs [22]. We replicated these results for three mammalian species (Additional file 1: Figure S1; P < 1.3E-25 for human, P < 3.5E-30 for mouse and P < 1.7E-31 for dog). This correlation could apparently be expected as 60-70% of CGIs overlap promoters [44] but we aimed to assess it more explicitly.
We first approached the genes-CGI relations by considering 100 kb bins genome wide. The vast number of genome 100 kb bins lack both CGIs and genes. To assess this, we calculated 2x2 contingency tables for human, mouse, and dog, correspondingly (Table 3). One may see that non-randomness of joint gene/CGI deserts as well as their co-occurrence is highly significant ( Table 1). In particular, the concordant classes for human (no genes and CGIs and at least 1 gene and at least 1 CGI) comprise 76% of bins leaving only 24% of discordant bins (Table 3). This joint distribution leaves no doubt of overall interrelation of CGIs and genes locations.
We further extended 2x2 tables analysis and built up the distributions of CGIs and genes based on their density across 3 species. The CGI vs gene numbers density per 100 kb distributions were highly correlated with approximately the same rate as chromosome wise (Additional file 1: Figure S2; P < 1.6E-31; df = 20 for human; P < 1.4E-23; df = 15 for mouse and P < 3.6E-30; df = 18 for dog). Notably, CGIs and genes correlate also in dog genome irrespective of nearly 2-fold excess of CGIs number over genes number in this species [22]. Thus, the co-variation of CGIs and genes densities across 100 kb segments in mammalian genomes is significant on a coarse grain even given the multiple (but not predominant) spurious independent clusters (Additional file 1: Figure S2).
Further we proceeded with human genome only in elucidating other features in elements distribution due to the lack of appropriate data for other species.

Open chromatin and elements density
Next we assessed if CGIs vs genes co-clustering may be mediated by chromatin accessibility factor. To check that, we built up the joint distributions of DHSs density against CGIs and genes (Additional file 1: Figure S3, S4, S5). The regression plots of the elements against 292 bins of DHS density rate (see materials and methods) revealed that the correlation of the DHS and genes number per DHS bin is r = 0.92; P < 1.4E-120 (Additional file 1: Figure S3a). The correlation between DHS and CGIs was the highest due to the inherent features of CGIs: r = 0.97; P < 1.1E-188 (Additional file 1: Figure S5a). So, the chromatin accessibility is an ultimate factor for clustering of CGIs and genes (Additional file 1: Figure S7). To confirm it in another way, we employed large scale whole genome chromatin state segmentation map from [31] to assess the distribution of CGIs and genes across different chromatin states.
As can be seen from Fig. 1, the zero and one gene and CGI bins consist largely of B2-B3 heterochromatin, and also of A2 euchromatin type. As long as gene density increases, the major gene and CGI clusters containing class confines to A1 open chromatin (Fig. 1).
We compared the distribution of CGIs and genes across 3 largest classes of chromatin (Fig. 2). We observed the gene dense open chromatin of A1 type preference for both CGIs and genes, while A2 and B3 chromatin classes distributions decay rapidly with increasing elements density. High correlation between CGIs and genes are observed (P < 1E-9 for all three cases).

Methylation state and chromatin accessibility
Based on average methylation state computed across 63 HAIB Methyl cell lines data, we found that approximately 20% of CGI (4057) exceeded empirically chosen threshold: average methylation score >400. They were assigned as "(hyper)methylated" states, and "hypomethylated" ones otherwise. Due to the U-shaped methylation score distribution across the data we unintentionally included a range of differentially methylated CGIs into "hypermethylated" sample, but they are not abundant.
Chromosome wise analysis revealed high correlation of genes number with both unmethylated and methylated CGI clusters (Additional file 1: Figure S11). Next we plotted the total, hypomethylated and hypermethylated samples against DHS densities (Additional file 1: Figure S5). To assess the significance of correlation between the CGIs number and DHS densities, we performed regression analysis of the corresponding plots (Additional file 1: Figure S6). We assessed the correlation significance as r = 0.97, P < 1.1E-188 for total CGIs set; r = 0.96, P < 1E-169 for hypomethylated CGIs and r = 0.95, P < 1.7E-154 for hypermethylated CGIs sample. We also observe twice less DHS overall density in hypermethylated set compared to hypomethylated ones probably due to the observed 4-fold smaller size of hypermethylated vs hypomethylated CGIs sets.
These results indicate that there is no straightforward correlation between CGI methylation and binding sites repulsion, which is compliant to current observations in [45,46]. It is plausible that though not all transcription factors are blocked from binding to CGIs by methylation, still the chromatin status could be altered [45].

Alu sequences and chromatin content
While the CGIs and genes exhibit highly specific exponential distribution, Alu sequences are a lot more abundant and, to our knowledge, no global highly specific discriminating factor for their clustering was reported. In particular, 42% of Alu sequences reside within gene loci, which is very close to random insertion pattern since genes occupy about 40% of genome [42]. Chromosome wide Alu distribution closely follows chromosome length, also implying random nature of Alu distribution independent of genes and CGI clusters, at least chromosome wide. Still, there are Alu clusters that were reported abundant genome wide and were linked to various factors including recombination rate and others [41].
We analyzed DHS and Alu distributions (Additional file 1: Figure S3b, S4b, S8a). It underlies that, in contrast with CGIs/genes, the noticeable number of Alu clusters reside in DHSpoor regions that, in accordance with previous randomly based chromosome wise entities, point to a largely stochastic nature of Alu insertions and/or probably other factors besides the open chromatin state, involved in Alu insertion preference and cluster expansions.
Nevertheless, Additional file 1: Figure S8a underscores the overall trend of Alu clusters to open chromatin, which is consistent, but with small statistical significance due to a big standard variation and small slope of regression line. Our interconnection link of this TE with CGIs and genes is thus rather speculative, based on the specific trends of particular subfamilies vs the chromatin state.
When we applied the chromatin breakdown procedure for Alu subfamilies, we observed that AluY clusters  Figure S10). We also employed L1 retrotransposons analysis to assess chromatin distribution. We found that L1 retrotransposons are randomly distributed in the genome except for highly DHS dense regions (Additional file 1: Figures S3c and S4c), and, consequently, gene dense A1 segments (Additional file 1: Figures S8b and S9b). Considering long L1 clusters, they distinctly reside in constitutive heterochromatin lack of DHS (Additional file 1: Figure S8c). Thus, Alus and L1 retrotransposons partially complement each other both in mode of distribution and in chromatin preference.

ANOVA analysis of the genes, CGIs, and retrotransposons distribution across TADs
ANOVA analysis was employed to complement the regression analysis results. We considered elements variation across 4064 TADs (Table 2) to see the overall elements density.
From Fig. 4 and ANOVA tests we may make 3 essential conclusions concerning TADs content of 4 elements: Based on the Tukey HSD test, no dramatic differences were observed between chromatin states for L1 distribution (Fig. 4c), confirming its random retrotransposition and broad variation in density in all classes. Still, when considering long L1, which comprise only 10% of the total pool of L1 genome instances, the situation changes sharply in favor of B2/B3 preference (Additional file 1: Figure S8c). A1 and A2 TADs significantly differ in all cases except (d), which makes the point of quite distinct patterns of euchromatin in human genome given both of them maintain open chromatin signature in histone and other marks [31]. Notably, while B4 maintains reported repressive histone marks [31], and in our case it comprises significant number of L1 transposable elements [47], still, many of its features are also similar to gene dense A1 chromatin (Fig. 4a-c), pointing at its reported duality [31,47].

CG rich elements and clustering properties
We selected gene density distribution bins to plot CGIs and Alu clusters distribution against (Fig. 5).

Discussion
The study aims at gaining the insight on the distribution specifics of several DNA elements, specifically genes, CGIs and Alus, from the large-scale chromatin landscape standpoint. We found that chromatin accessibility rate is the major factor of joint clustering of these elements. Heterochromatic genomic segments spanning nearly a half of the human genome are void of both CGIs and genes, making their distribution highly skewed genome wide. This keeps valid for others species such as mouse and dog (Table 3, Additional file 1: Figures S1 and S2) and might be extended on other mammals. With the help of coarse genome partition by 100 kb non-overlapping segments we showed that CGIs extremely closely associate with the genes. It was shown that this association is consistent and independent of the proportion of the CGIs vs genes number, which is approximately equal in the mammals majority. In this regard, for the dog genome, characterized by an unusual twofold excess of CGIs instances over the genes number, it was found that an excess can be explained by the relatively few instances of CGI extended clusters. At the same time the proportion and localization of genes and the vast majority of CGIs in the dog genome are coordinated as in other mammals (Additional file 1: Figure S2c).
Indeed, when we mapped the dense CGI clusters onto the dog genome, we clearly observed their highly non random grouping at subtelomeric and pericentromeric regions sparse of genes (Additional file 1: Figure S12). These specific telomeric and centromeric superclusters spanning 500 kb on average were observed at virtually all dog chromosomes (subtelomeric regions in particular) and comprise at least 30% of the annotated CGIs (16,000 instances by rough approximation). The subtelomeric CGI clusters have been observed in cat, horse, and bovine, as well as occasionally in human and have been reported in previous studies [21], but not with the same expansion rate and abundance as observed in dog. The apparent function of these CGI clusters is not currently elucidated, but may be related to the repeated nature of the DNA, characteristic for these regions.
Employment of the expanded chromatin spectrum using chromatin signatures inferred from topological associated domains [31] revealed that the distributions of CGIs and genes are highly similar chromatin wise both on the fine, DHS-mediated resolution (Additional file 1: Figures S3-S6) as well as on the coarse-grained topological chromatin domains ( Figs. 1 and 2). This makes us suggest that CGIs are the inherent elements of the genes irrespective of their methylation state and location. The independent tandem expansions of CGIs and gene families do not affect the total trend due to their minority.
ANOVA analysis of the elements content in 6 TAD subcompartments (Fig. 4) corroborated the inferences derived from regression analysis. Additionally, TAD classes interrelations could be observed for the distinct elements content underlining the specific features of TAD classes.
It was elucidated that the genes, Alu and CGIs clusters density is monotonically dependent on the chromatin state derived by DHS densities (Additional file 1: Figures S7 and 8). Considering the TADs mediated chromatin partition, the highest density of these elements is found in A1 euchromatin (Additional file 1: Figures S3-S6 and S10). These regions are characterized by a high density of small length genes, and replicating in the early S-phase (G1, S1 phases) [31]. Thus, these elements mediate or are a side element of chromatin accessibility. It is likely that Alus use open chromatin, but CGIs create it.
Using L1 transposable elements as a background we showed by ANOVA (Fig. 4) and distribution analysis (Additional file 1: Figures S8 and S9) that L1 instances are distributed independently of chromatin state, while subset of long L1 sequences are distinctly resided in heterochromatin.
In this context we'd like to draw attention to chromosome 19 due to the extraordinary CG-rich elements density, both for CGIs (Additional file 1: Figures S1a and S11) and Alus (98 Alus per 100 kb; next highest is chr17 with 70 per 100 kb, and genome average is 40 Alus per 100 kb; P < 0.004 by ESD test). Given it is more than twofold dense in gene number than any other chromosome [48], it thus implicitly confirms the close coordination of genes and CG-rich elements. Importantly, chromosome 19 also is the most hypermethylated one chromosome wide (Additional file 1: Figure S11b), while hypomethylated at fetal stage [49] (Additional file 1: Table S1), which may imply specific mode of evolution of this particular chromosome by gene duplications [48,50], and the abundance of development genes [47]. In particular, it was proposed that the random promoter methylation within the newly emerged paralogous gene pair in the course of embryonic reprogramming stage may highly increase the chance for both genes to keep the functionality upon duplication [50].

Conclusions
One of the possible reasons of chromosome 19 observed evolutionary mode may be that it comprises more than a hundred of KRAB-ZNF genes organized in clusters and represents "defense" system against hypomethylated transposable elements in embryogenesis. It is responsible for identification of non-methylated CpG elements located at retrotransposons. Upon recognition they recruit TET protein, which, in turn, recruits heterochromatin modifiers to repress them, and subsequently change methylation status [51,52]. These particular clusters reside in chromosome specific B4-heterochromatin type nuclear compartments [31] which demonstrate even higher than A1 abundance with Alu clusters (Fig. 4). It may imply B4 accessibility in early embryogenesis, and thus a high fixation rate. The latter hypothesis is also supported by B4 non-random abundance with other transposable elements clusters, such as L1M3a [49] and others.
Coherent with CpG abundance, chromosome 19 contains single A1 euchromatin class, and virtually no constitutive heterochromatin B3 and euchromatin A2 domains, manifested as the largest chromatin classes occupying 50% of the genome ( Table 2). The A1 "style" of gene arrangement on chromosome 19 underlies the fact that the genes located within this chromosome maintain the average length shorter than the genome wide average with a high statistical significance.
Notably, A1 chromatin state relates to early replicated domains [31] and, hence, comprises multiple origin of replication sites (ORI) that fire shortly after G phase. As long as there could be some analogy with the insertion pattern of P-elements in D.melanogaster, where ORI sites are shown to be the hotspots of P-element insertions [52], we speculate that Alu retrotransposons may display the same affinity to ORI sites attributable to open A1 chromatin, and, hence populate the ORI-related chromatin.
Large scale features of open chromatin may affect the gene context at a lower resolution level. Besides well known CGIs preference to 5′ gene regions, we observed reported previously high density of CGIs in 3′ gene regions [17,22]. Significant portion of Alus reside both in promoter regions [53] and in 3′ UTR elements [54]. Given both transcription start and end sites maintain specific epigenetic markers [26], we put forward the hypothesis that both of these elements form the gene locus architecture accommodating gene expression maintenance. In particular, CGIs and Alu may cooperate in promoter regions given the recent observation that highly expressed housekeeping genes maintain their CpGpromoter flanked with hypermethylated "shells", possibly accommodating more accurate RNA PolII positioning [55] and involving methylated Alu retrotransposons as CGI "scaffolding" factor [53].
The non-random distribution of the elements considered mediated by chromatin state may be practically utilized by applications in the association studies for refining their statistical significance in location/clustering by incorporating prior knowledge on the particular chromatin state in a model.
The evolutionary implications of the elements within chromatin state context may confer that L1 and CGIs clusters could lead to the chromatin state setting in the course of expansion, while Alu elements seem rather opportunistic in location preference due to nonautonomous nature.
Finally, the observed CGIs expansion dynamics may impact the issue of the genome CpG content maintenance. CpG methylation rate was recently proposed to be the driving force of genome size evolution within a tradeoff in methyltransferases methylation rate vs genome size [56]. As an example, D. melanogaster and C. elegance minimal genomes phenomena could be explained by the absence of methyltransferases, and, as a consequence, CpG content equaling to the expected one. Thus, it looks like CGI duplication mechanism is of evolutionary importance for maintaining sufficient CpG number. In this regard, some subtle mechanisms apparently 'compel' the dog genome to amplify CG-rich tandems within subtelomeric/pericentromeric regions chromosome wide, even given the significantly reduced gene repertoire and corresponding CGI promoters [22], this way aligning CpG (or just GC) content with that of other mammals.

Declaration
This article has been published as part of BMC Evolutionary Biology Vol 17 Suppl 1, 2017: Selected articles from BGRS\SB-2016: evolutionary biology. The full contents of the supplement are available online at https://bmcevolbiol. biomedcentral.com/articles/supplements/volume-17-supplement-1.

Funding
Data compilation was supported by ICG SB RAS budget project 0324-2015-0003 for IC. YO acknowledges the support of the software development by the Russian Science Foundation (http://rscf.ru/en/), project 14-24-00123. VB received support for transcriptome data analysis in mouse and rat from the Russian Science Foundation, project 14-14-00269. The Russian Science Foundation has funded publication of this article. No funding bodies had any role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Availability of data and materials
The datasets used and analyzed during the current study available from the corresponding author on reasonable request.