Epigenetic supersimilarity in MZ twins

Rather than being predominantly determined by genetics, interindividual variation in DNA methylation at MEs is determined, at least in part, stochastically [6] and influenced by the nutritional milieu of the preimplantation embryo [10,11,12]. We therefore expected that, at MEs, methylation concordance within MZ twin pairs would be greater than that of unrelated individuals, but comparable to that within DZ twin pairs. To test this, we analyzed a genome-scale DNA methylation data set from Grundberg et al. [9], who used the HM450 array to assess methylation in adipose tissue from adult female twins of European-descent (97 MZ twin pairs and 162 DZ twin pairs). As did Grundberg et al., we discarded low-quality probes potentially affected by single nucleotide polymorphisms (SNPs) and, of the remaining 344,303 probes, focused our analysis on the 10% (34,405) with the highest interindividual variance (hereafter referred to as the top 10%).

Within regions previously identified as candidate or bona fide MEs [13, 14], we assessed twin–twin methylation concordance inversely by probe-specific mean square error (MSE) of β values. MSE assesses the deviation of a twin pair from the line of identity, providing a direct measure of discordance. Contrary to our expectation, MZ twin concordance in putative ME regions was between 2.5- and 16.5-fold higher than that of DZ twins (Fig. 1). This suggested that establishment of DNA methylation at these regions is under genetic control. To test this, we examined the probe-specific narrow-sense heritability (h2) estimates (based on the ACE method [15]) from Grundberg et al. [9]; h2 is the proportion of phenotypic variation in a population that is attributable to genetic variation [16]. Strikingly, 1058 probes (3% of total) showed h2 estimates > 1 (Fig. 2a). Most of the probes within the candidate MEs featured in Fig. 1 were among them (Fig. 2a), indicating that these superordinate h2 values are not simply a result of sampling error.

Fig. 1 Methylation in MZ twin pairs is highly concordant at candidate metastable epialleles (MEs). Each plot shows probe-specific β values for 97 MZ (blue, twin 2 > twin 1) and 162 DZ (red, twin 2 < twin 1) twin pairs, at loci previously identified as bona fide or candidate MEs [13, 14]. Insets show locus-average mean square errors (MSE) across all the MZ and DZ twins. MSE is much lower in MZ compared to DZ twins. a VTRNA2-1, 15 probes, 10.6-fold lower MSE (in MZ vs. DZ). b DUSP22, 11 probes, 16.5-fold lower. c PAX8, eight probes, 2.5-fold lower. d CYP2E1, three probes, 10.8-fold lower. e SFT2D3, four probes, 3.1-fold lower. f CFD, one probe, 6.6-fold lower Full size image

Fig. 2 Some HM450 probes exhibit epigenetic supersimilarity (ESS). a Distribution of probe-specific narrow-sense heritability (h2) estimates from [9]. (Shown are data on 24,839 probes; 9566 probes with h2 < 0.001 were excluded for clarity.) Of the probes, 1058 show h2 > 1, including most of the probes illustrated in Fig. 1 (red box plot). b Normalized DZ MSE vs. MZ MSE for the 34,405 probes (top 10%) from Grundberg et al. [9]. Histograms (right and top) show distribution; red curves show best normal fit. Normalized DZ MSE (mean ± standard deviation = 0.76 ± 0.13) is normally distributed, but normalized MZ MSE (0.63 ± 0.23) is skewed left (P = 7.0 × 10–66). Probes with h2 > 1 are shown in blue. Probes to the left of the green line (y = 2x) are classified as ESS. c Associations between probe-level mQTL and heritability estimates (both from Grundberg et al. [9]). Among the 9708 probes that are both in the top 10% of interindividual variance and positive for mQTL (top panel) mean heritability is 0.64 (gray vertical line) and positively associated with the strength of mQTL. Among ESS probes positive for mQTL (middle panel), mean heritability is 0.90 and not associated with mQTL. Mean heritability of ESS probes negative for mQTL (0.99, bottom panel) is similar to that of mQTL-positive ESS probes. d Model to explain ESS in MZ twins. Numbers on the dice represent different methylation states at a specific locus. If de novo methylation occurs after embryo cleavage (top), each MZ embryo undergoes independent establishment. If de novo methylation occurs prior to embryo cleavage (bottom), both MZ embryos inherit the same methylation state. e Consistent with this model, bisulfite pyrosequencing in three tissues of 17 cadavers indicates that ESS probes also show systemic interindividual variation. Two examples are shown, OR2L13 and HLA-DQB2 Full size image

Clearly, h2 values > 1 are difficult to interpret. To better understand this observation, we calculated MSE for all 34,405 top 10% probes [9]. To elucidate the extent to which DZ and MZ twins are more similar than pairs of unrelated individuals, probe-specific MSEs were normalized relative to randomized pairs (RZ), simulating pairwise MSE within the general population. DZ/RZ MSE and MZ/RZ MSE were generally < 1, as expected (Fig. 2b). Genetic influences on CpG methylation generally occur when the local sequence context in cis (i.e., a haplotype) affects establishment of methylation [17]. Given that DZ twins are identical by descent at 50% of haplotypes [18] and MZ twins at 100% of haplotypes, a model based on genetic determination predicts that the mean normalized DZ MSE should be no more than twice the mean normalized MZ MSE. Hence, for probes to the left of the green line (y = 2x) in Fig. 2b, MZ twin pairs show greater-than-expected similarity in DNA methylation. We refer to this phenomenon as “epigenetic supersimilarity” (ESS). According to the central limit theorem, assuming that probe-specific methylation is determined by many unobserved (genetic) factors, the mean intra-pair errors should be normally distributed. Indeed, normalized DZ MSE are, but normalized MZ MSE are skewed to the left (P = 7.0 × 10–66) (Fig. 2b). Each probe with DZ/MZ MSE > 2 (corresponding to those left of the green line in Fig. 2b) is > 5 standard deviations (sd) away from the expected normal mean (P < 0.0001) (Additional file 1: Figure S1a), well beyond the range of sampling error. Most of the probes for which Grundberg et al. estimated h2 > 1 are characterized as ESS (Fig. 2b). Our initial validation studies found that many ESS probes with interindividual β range < 0.4 in the Grundberg et al. data set [9] are essentially unmethylated in several human primary tissues. We therefore refined the selection criteria to MSE DZ/MZ > 2 and an interindividual β range > 0.4, identifying 1580 probes (4.6% of the 34,405) as ESS (Additional file 1: Figure S1b and Additional file 2: Table S1). Across all probes with β range > 0.4, normalized DZ MSE remained normally distributed, but normalized MZ MSE were shifted even further to the left (Additional file 1: Figures S1b, c).

To further test whether the superordinate heritability estimates of Grundberg et al. might somehow result from the genetic identity of MZ twin pairs, we analyzed their data on methylation quantitative trait loci (mQTL), i.e., sequence variants correlated with methylation at specific CpG sites [19]. Grundberg et al. [9] combined their genotyping and HM450 data on 603 adipose tissue samples and applied a conservative significance threshold (P < 1.2 × 10–9), identifying 9708 mQTL probes within the top 10% of interindividual variance. Among these, as expected, the strength of the mQTL association was positively associated with heritability (Fig. 2c, top). There was no such association across ESS probes (Fig. 2c, middle). If the superordinate heritability associated with ESS results from the genetic identity of MZ twins, the mean heritability of ESS probes with mQTL should be higher than that of those without mQTL. This was not the case (Fig. 2c, middle and bottom). This analysis, using mQTL data from the same samples in which we identified ESS, provides strong evidence that ESS is not simply a consequence of the isogenicity of MZ twins.

Testing a model for ESS

During MZ twinning, if de novo DNA methylation at a particular locus occurs prior to embryo cleavage, both twins will inherit the same epigenotype at the locus simply because of developmental timing, rather than as a consequence of their genetic identity [14]. This provides a potential explanation for ESS (Fig. 2d). If correct, methylation at ESS loci must be established in the cleavage-stage embryo. If the epigenetic state is maintained during subsequent cellular differentiation, these loci should show systemic interindividual variation in DNA methylation.

To test this, we selected 13 ESS regions and assessed systemic interindividual variation (SIV) by bisulfite pyrosequencing in liver, kidney, and brain of cadaver tissues [13]. Methylation tended to be correlated in these tissues derived from the different embryonic germ layers (Fig. 2e; Additional file 1: Figure S2). Overall, 9 (69%) of the 13 loci showed evidence of SIV (Additional file 2: Table S2). For a broader evaluation, we analyzed a previously published data set from Lokk et al. [20], who profiled multiple tissues from each of several individuals using the HM450 platform. From each of four individuals, we considered data for all SNP-free and high-quality HM450 probes for tissues representing the three embryonic germ layer lineages: gall bladder (endodermal), abdominal aorta (mesodermal), and sciatic nerve (ectodermal). Individual- and tissue-specific methylation were estimated as the average across the three tissues and the four individuals, respectively, and variation was quantified as the range of these averages (Fig. 3a). Though most probes showed little of either (Fig. 3b, histograms), tissue-specific was generally greater than interindividual variation (Fig. 3b). To focus on robust SIV we restricted our analysis to probes with interindividual variation that was at least 0.2 delta β and three times greater than tissue-specific variation (Fig. 3b, shaded region). These cutoffs identified 1042 probes with evidence of SIV (Additional file 2: Table S1). Bisulfite pyrosequencing in cadaver tissues (Fig. 3c; Additional file 1: Figure S3) confirmed SIV at 8 (67%) of 12 regions evaluated (Additional file 2: Table S3).

Fig. 3 Epigenetically supersimilar (ESS) probes are enriched for systemic interindividual variation (SIV). a Analytical strategy applied to data of Lokk et al. [20] on abdominal aorta, gall bladder, and sciatic nerve from each of four individuals. Interindividual and tissue-specific variation were quantified as the range of the individuals’ mean beta values (μ 1 , μ 2 , μ 3 , μ 4 ) and the tissues’ mean beta values (μ aa , μ gb , μ sn ), respectively. b Tissue-specific vs. interindividual variation for 344,151 probes. Histograms (top and side) indicate the density distribution. Green lines illustrate cutoffs used to identify 1042 SIV probes (shaded region, interindividual > 0.2 and tissue-specific < 1/3 of interindividual variation). c Examples of bisulfite pyrosequencing data confirming systemic interindividual variation in selected loci: PF4 and LDHC. d The 1580 probes with evidence of ESS are 6.3-fold enriched for SIV (P < 10–10, chi-squared test) Full size image

Perfect overlap between ESS and SIV probe sets was not anticipated for two reasons. First, as they survey only four individuals, the Lokk et al. data cannot capture all interindividual variation. Second, epigenetic states established prior to gastrulation may not be maintained in all differentiated lineages (i.e., early embryonic establishment is necessary but not sufficient for SIV). Nonetheless, relative to the 5388 non-ESS probes with interindividual range > 0.4, the 1580 ESS probes were 6.3-fold enriched for SIV (P < 10–10, chi-squared test; Fig. 3d), supporting our model for the developmental basis of ESS (Fig. 2d).

ESS and SIV sites share genomic and epigenomic features, and are enriched for MEs

ESS appears to be a marker for individual-specific epigenetic states that are established in the cleavage-stage embryo. Such states could be established under genetic influence, or stochastically; only the latter are consistent with epigenetic metastability [6]. The mQTL data of Grundberg et al. (Fig. 2c) demonstrate that ESS is not generally associated with genetic effects. To test this more generally we evaluated additional data sets in which the HM450 platform was used to assess mQTL in at least 100 individuals [21]. Volkov et al. [22] profiled SNPs and DNA methylation in adipose tissue of 119 men and identified 15,208 CpG sites with significant cis-mQTL. Shi et al. [23] assessed mQTL in histologically normal lung tissue from 210 individuals and reported estimates of the proportion of methylation variance explained by neighboring SNPs (which we refer to as β SNP ). We considered probes with β SNP > 0.33 as exhibiting substantial mQTL. Of the 34,304 probes Shi et al. identified with statistically significant cis-mQTL, only 4306 (12.6%) showed substantial mQTL (Additional file 1: Figure S4). Although both the Grundberg et al. [9] and Volkov et al. [22] data were based on adipose tissue, less than half of the mQTL probes identified by either were identified in both (Additional file 1: Figure S5). Conversely, most of the Shi et al. substantial mQTL probes were also identified by the other two studies (Additional file 1: Figure S5). Moreover, > 80% of the probes Shi et al. reported as substantial mQTL in lung also exhibited significant mQTL in independent studies of breast and kidney [23]. For these reasons, we focus our subsequent analyses on the Shi et al. substantial mQTL probe set. (Nonetheless, we have included data on all three mQTL lists in our annotation of ESS and SIV probes in Additional file 2: Table S1.)

Relative to the probe sets from which they were drawn, those with evidence of either ESS (Fig. 4a) or SIV (Fig. 4b) were enriched for substantial mQTL (15- and 24-fold, respectively, both P < 10–10, chi-squared test). Substantial mQTL affected 25.1% of ESS and 17.9% of SIV probes (Additional file 2: Table S1). We consider ESS probes without evidence of substantial mQTL to be candidate MEs. Likewise, since our SIV analysis is analogous to previous ME screens [13, 14], SIV probes without evidence of substantial mQTL are also candidate MEs. Indeed, most of the HM450 probes identified as MEs in a previous screen that employed genome-wide bisulfite sequencing [13] overlap with mQTL-filtered ESS or SIV hits (Fig. 4c). After excluding those with substantial mQTL, ESS probes remained 5.6-fold enriched for SIV (P < 10–10; Fig. 4c), indicating that the common epigenetic behavior of these probes sets is not due to genetic effects. Importantly, most ESS probes do not show substantial mQTL (Fig. 4a), further evidence that ESS is not simply a consequence of the genetic identity of MZ twins. To directly assess the influence of local sequence on interindividual variation at ESS loci we validated several top hits, performing genotyping and methylation analysis by pyrosequencing in peripheral blood DNA of 64 Gambian children [13]. Each genotyping assay targeted a nearby common SNP with minor allele frequency > 25%. Two regions negative for mQTL (CYP2E1 and DUSP22) and two with some mQTL-positive probes (SPATC1L and ZFP57) showed substantial interindividual variation even among individuals of the same genotype (Fig. 5a–d). These regions show strong linkage disequilibrium, indicating that SNP genotype is generally an indicator of haplotype. Notably, the SNP we genotyped at ZFP57, rs3129057, was recently reported to be the strongest index SNP in phase with haplotype-dependent allele-specific methylation (Hap-ASM) in the region [17]. Significant mQTL was detected for ESS CpGs at DUSP22, SPATC1L, and ZFP57 (Fig. 5b–d). At these same loci, however, interindividual variance of methylation was associated with haplotype, providing the novel insight that the local sequence context can influence epigenetic metastability. The pyrosequencing results were further validated by clonal bisulfite sequencing for selected individuals at SPATC1L and ZFP57 (Fig. 5e and f), confirming that even in regions of substantial mQTL, individuals with the same local sequence context can exhibit dramatic interindividual variation in DNA methylation.

Fig. 4 Regions of epigenetic supersimilarity (ESS) and systemic interindividual variation (SIV) share genomic and epigenomic features. a Normalized DZ MSE vs. MZ MSE for the 6968 probes with range > 0.4, of which 489 (red) show substantial mQTL. Inset: ESS probes are 15-fold enriched for substantial mQTL (P < 10–10, chi-squared test). b Tissue-specific vs. interindividual variation at 344,151 probes, of which 2702 (red) are substantial mQTL. Inset: SIV probes are 24-fold enriched for substantial mQTL (P < 10–10, chi-squared test). c After filtering out substantial mQTL, ESS and SIV hits overlap more than two-thirds of probes at previously identified MEs [13]. d Relative to all probes in the top 10% of interindividual variance, ESS and SIV probe sets are enriched for CpG islands (both comparisons P < 10–10, chi-squared test). e Gene set enrichment analysis shows that both ESS and SIV probes are enriched for genes expressed in cancer (P = 4.7 × 10–8 and 4.8 × 10–9, respectively). Each row represents a different type of cancer in The Cancer Genome Atlas [24] (key to abbreviations in Additional file 2: Table S5). f Association of probe sets with epigenomic feature annotations derived from 111 reference epigenomes [25]. ESS and SIV probes are enriched for active promoters (TssA) and underrepresented at enhancers (Enh) (all four comparisons P < 10–10) Full size image

Fig. 5 Interactions between DNA methylation and local sequence context at some top ESS regions. a–d Average methylation vs. SNP genotype at ESS regions within CYP2E1, DUSP22, SPATC1L, and ZFP57. In each panel, gene diagram (top) shows location of ESS region where methylation analysis was performed (asterisk) relative to that of a SNP that was genotyped in 64 Gambian children. Grid summarizes normalized linkage disequilibrium (D’) across these ~3-kb regions in a Gambian population in Western Gambia (GWD, 1000 Genomes Project [80]). With the exception of G/G individuals at rs3129057 (ZFP57), there is substantial interindividual variation in average methylation within each genotype class. At CYP2E1 (a), average methylation is not associated with SNP genotype (P = 0.31). At DUSP22, SPATC1L, and ZFP57 (b–d) average methylation is associated with genotype (P = 0.002, 0.02, and 0.0001, respectively). At these same loci, interindividual variance differs between the two homozygous genotypes; i.e., C/C vs. T/T at DUSP22 (P = 0.02), G/G vs. A/A at SPATC1L (P = 0.04), and G/G vs. A/A at ZFP57 (P = 1.9 × 10–6). e, f Clonal bisulfite sequencing data at two homozygous individuals at each of SPATC1L and ZFP57, respectively, confirm dramatic interindividual variation in DNA methylation in the absence of local sequence variation. Black, empty, and gray circles represent methylated, unmethylated, and indeterminate CpG sites, respectively. Vertical red line indicates the position of the SNP Full size image

Relative to negative control probes with interindividual variation comparable to ESS probes but no evidence of ESS or SIV (Additional file 1: Figure S6 and Additional file 2: Table S4), ESS and SIV probes were 3.6- and 5.0-fold enriched for CpG islands, respectively (Fig. 4d; P < 10–10 for both comparisons). Likewise, ESS and SIV probes were 3.3- and 2.4-fold enriched in subtelomeric regions (<2 Mb from chromosome ends; Additional file 1: Figure S7a; P < 10–10 for both comparisons). Since subtelomeric regions are rich in genetic variation, we tested whether the subtelomeric enrichment might be due to mQTL. However, similar enrichments were found in the ESS and SIV probe subsets not associated with substantial mQTL (Additional file 1: Figure S7b). The ESS and SIV gene lists each included six genomically imprinted genes, no different from what is expected by chance; imprinted loci among these two classes are ANO1, GNAS, GRB10, NAP1L5, NLRP, and VTRNA2-1 (ESS) and DLGAP2, KCNQ1OT1, NAP1L5, NLRP2, and VTRNA2-1 (SIV) (http://www.geneimprint.com/). Gene set enrichment analysis (GSEA) using data from The Cancer Genome Atlas [24] showed that, relative to negative controls, both ESS and SIV probes are more likely to be annotated to genes expressed in a wide range of tumors (Fig. 4e). Across 111 reference epigenomes encompassing a wide range of cell lines and primary tissues [25], both probe sets were enriched for active promoters and depleted for enhancers (Fig. 4f). ESS and SIV CpGs were identified independently but exhibit highly overlapping genomic and epigenomic features, indicating that they share similar fundamental biological properties.

Periconceptional environment affects establishment of methylation at ESS and SIV CpGs

Mouse [10, 11, 26] and human studies [13, 14, 27, 28] have shown that establishment of DNA methylation at MEs is sensitive to periconceptional environment. Previous studies tested this using a “natural experiment” exploiting seasonal variation in maternal nutritional status in The Gambia [12]. Here, we analyzed an independent set of 128 blood samples collected from 2-year-old Gambian participants in the Early Nutrition and Immune Development (ENID) trial [29] who were conceived at the peak of either the rainy or the dry season [30]. Based on the notion that MEs are largely free of genetic influence, we set out to test whether ESS and SIV probes without substantial mQTL show season-of-conception effects. To our surprise, we found comparable and highly significant enrichments for season-of-conception effects in ESS and SIV probe sets regardless of substantial mQTL (Additional file 2: Tables S17 and S18). We therefore examined the unfiltered probe sets and found that both ESS and SIV probes, but not negative control probes, were significantly enriched for season-of-conception effects (P = 3.3 × 10–9 and 1.4 × 10–23, respectively; Fig. 6a). Consistent with previous studies of candidate MEs in independent cohorts [12,13,14], being conceived during the rainy season generally resulted in higher levels of DNA methylation at both ESS and SIV loci (Fig. 6b). These findings further support the conjecture that, regardless of mQTL, many ESS and SIV probes are MEs. This season of conception effect also provides independent support for our model (Fig. 2d) that ESS arises due to early embryonic establishment of epigenotype.

Fig. 6 Sites of epigenetic supersimilarity (ESS) and systemic interindividual variation (SIV) are enriched for effects of periconceptional environment on DNA methylation. a Relative to all probes on the HM450 array, mQTL-filtered ESS and SIV probes (but not negative control probes) are highly enriched for significant (FDR < 10%) associations with season of conception in rural Gambia. b Heat map of average effect of season of conception at loci that show a significant seasonal difference in methylation (FDR < 10%). At both ESS and SIV probes, as in previous studies of MEs in independent cohorts [13, 14], children conceived in the rainy season have higher methylation Full size image

Prospective associations between DNA methylation in blood and later cancer

Although the probe signature of ESS was identified from a study of adult twins, it appears to be a consequence of methylation establishment in the early embryo and hence must be stable from embryonic development to adulthood. Since ESS is associated with genes expressed in tumors (Fig. 4e), we asked whether interindividual variation in DNA methylation at ESS loci predicts risk of later cancer in adults. To test this, we examined data from the Melbourne Collaborative Cohort Study (MCCS), which enrolled 41,514 healthy adult volunteers between 1990 and 1994 [31]. Peripheral blood samples and information on health-related behaviors were collected at enrollment, and incident cases of cancer were ascertained prospectively by linkage to the Victorian Cancer Registry, which receives mandatory notification of all new cancer cases in Victoria, Australia. The systemic nature of interindividual variation at ESS probes enabled us to use DNA methylation in peripheral blood as an indicator of methylation in various tissues. A control was matched to each incident case on sex, country of birth, and age at enrollment, using density sampling. Using the Illumina HM450 platform, DNA methylation at baseline was assessed on 3464 case-control pairs overall in studies of seven types of cancer [32] (Additional file 2: Table S6); average time from sample collection to diagnosis was 9.2 ± 4.9 years (mean ± sd).

Regardless of potential genetic influences, our data indicate that interindividual epigenetic variation at ESS probes occurs systemically and is stable over time. We therefore evaluated ESS probes without regard to mQTL. Combined effects across multiple CpGs (i.e., differentially methylated regions) are more likely to demonstrate long-term stability and affect gene expression [33]. Hence, rather than analyze individual probes, we focused on the 198 clusters of multiple ESS probes separated by no more than 500 bp (523 CpGs total; Additional file 2: Table S7). Analysis of expression [34] and methylation data [9] from Grundberg et al. showed that at many ESS clusters, average methylation in adipose tissue is associated with gene expression, not only in adipose tissue but also in skin and lymphoblastoid cell lines (Additional file 1: Figure S8 and Additional file 2: Table S8). These results provide evidence that methylation at ESS loci in one tissue yields information about epigenetic regulation in additional tissues. To test for probe-specific associations between peripheral blood DNA methylation at baseline and risk of specific cancer diagnosis we performed conditional logistic regression, adjusting for estimated leukocyte composition (using the Houseman algorithm [35]) and other covariates. Statistical significance of associations at the cluster level were then evaluated by permutation testing (tabulated results in Additional file 2: Table S11). Relative to negative control clusters (Additional file 2: Table S12), the 198 ESS clusters were four-fold enriched for associations with later cancer (P = 1.5 × 10–5). To minimize multiple testing, we focused on the ten ESS clusters with the largest number of CpG sites. Remarkably, at seven of these, peripheral blood DNA methylation at baseline was significantly associated with later cancer (Fig. 7; Additional file 2: Table S13); three (SPATC1L, VTRNA2-1, and DUSP22) were significantly associated with multiple types of cancer. Elevated methylation in a cluster of six CpG sites at SPATC1L was associated with reduced risk of colorectal and prostate cancer (Fig. 7b, f), and elevated methylation in a cluster of 12 CpGs at VTRNA2-1 was associated with higher risk of lung cancer and mature B-cell neoplasm (Fig. 7d, e). Interestingly, elevated methylation in a cluster of eight CpG sites at DUSP22 was associated with increased risk of mature B-cell neoplasm (Fig. 6e) yet reduced risk of urothelial cell carcinoma (Fig. 7g). The 154 negative control clusters showed few and relatively weak associations with later cancer (Additional file 1: Figure S9). Results are also shown for the 128 ESS clusters that included no probes with substantial mQTL (Additional file 1: Figure S10 and Additional file 2: Table S14).