VTRNA2-1: Genetic Variation, Heritable Methylation and Disease Association

VTRNA2-1 is a metastable epiallele with accumulating evidence that methylation at this region is heritable, modifiable and associated with disease including risk and progression of cancer. This study investigated the influence of genetic variation and other factors such as age and adult lifestyle on blood DNA methylation in this region. We first sequenced the VTRNA2-1 gene region in multiple-case breast cancer families in which VTRNA2-1 methylation was identified as heritable and associated with breast cancer risk. Methylation quantitative trait loci (mQTL) were investigated using a prospective cohort study (4500 participants with genotyping and methylation data). The cis-mQTL analysis (334 variants ± 50 kb of the most heritable CpG site) identified 43 variants associated with VTRNA2-1 methylation (p < 1.5 × 10−4); however, these explained little of the methylation variation (R2 < 0.5% for each of these variants). No genetic variants elsewhere in the genome were found to strongly influence VTRNA2-1 methylation. SNP-based heritability estimates were consistent with the mQTL findings (h2 = 0, 95%CI: −0.14 to 0.14). We found no evidence that age, sex, country of birth, smoking, body mass index, alcohol consumption or diet influenced blood DNA methylation at VTRNA2-1. Genetic factors and adult lifestyle play a minimal role in explaining methylation variability at the heritable VTRNA2-1 cluster.


Introduction
Mendelian-like inheritance of germline DNA methylation can be due to cisor transacting genetic factors known as methylation Quantitative Trait Loci (mQTL) or epimutations (heritable change in gene activity that is not associated with a DNA mutation but rather with gain or loss of DNA methylation or other heritable modification of chromatin). Both can mimic germline pathogenic variants in their effect on gene function and disease association and discriminating between the two possibilities (mQTL or epimutation) in specific genomic regions and disease context is often challenging.
We previously made a genome-wide assessment of heritable methylation using a family design [1]; probes were ranked by a methylation-heritability metric and 24 of the 1000 most heritable CpGs were identified to be associated with breast cancer risk in these families. Several CpGs within VTRNA2-1 were among those that appeared to be most heritable and most strongly associated with breast cancer risk, including five within the gene promoter [1]. Evidence that DNA methylation can be transmitted from parent to offspring in the absence of a genetic explanation is scarce and controversial [2][3][4]. It is therefore important to assess whether any genetic variation may influence DNA methylation at the VTRNA2-1 region, i.e., cisor trans-acting mQTL.
VTRNA2-1 has been demonstrated to adopt two structurally and functionally distinct RNA conformations, one which strongly inhibits protein kinase R (PKR) and downstream eukaryotic translation initiation factor 2 subunit α (eIF2α) phosphorylation, and one which acts as a pseudo-inhibitor of PKR when competing with other double-stranded (ds) RNA molecules [5]. The~2 kb region overlapping the VTRNA2-1 locus has been reported to be polymorphic or atypically imprinted, and somatically acquires DNA methylation of the maternal allele in the majority of cases [6,7]. PKR is an interferon-induced kinase consisting of 551 amino acids that acts as an intracellular stress sensor, primarily associated with viral infection. dsRNA produced by viral replication binds to and activates PKR, causing dimerisation and subsequent phosphorylation of its substrate, eIF2α. Phosphorylation of eIF2α converts eIF2 from a substrate to an inhibitor of its GDP-GTP exchange factor eIF2B, inhibiting mRNA translation and arresting global protein synthesis [8,9]. In addition to eIF2α phosphorylation, PKR can activate the nuclear factor kappa-light-chain-enhancer of the activated B cells (NF-κB) signalling pathway, which is known to play an oncogenic role in tumourigenesis [10]. PKR is reported to act on the NF-κB pathway by inducing phosphorylation of IκBα, which interestingly requires the expression of VTRNA2-1 to occur [9].
There is evidence that VTRNA2-1 may act as a tumour suppressor and is a metastable epiallele [11,12]. The study by van Baak et al. used data from the Melbourne Collaborative Cohort Study (MCCS) to assess the association between methylation at metastable epialleles and cancer risk and concluded that methylation at VTRNA2-1 was potentially associated with risk of lung cancer and B-cell lymphoma [13]. We also showed that VTRNA2-1 promoter methylation was associated with prostate cancer risk, and these associations appeared stronger for aggressive disease [14].
Consistent with a familial aggregation of VTRNA2-1 methylation not being due to genetic factors are observations that DNA methylation at this locus is sensitive to the perinatal environment, including factors such as season of conception [13,15] or maternal folic acid supplementation during pregnancy [16]. That VTRNA2-1 methylation at birth was found to be associated with childhood overweight/obesity [17] may signal another plausible link with cancer risk in adulthood, in addition to the tumour suppressor role [15].
DNA methylation at VTRNA2-1 has therefore been hypothesised to show the following pattern [7]: (1) The paternally-inherited allele seems to always be unmethylated, as observed in several studies [1,14], and (2) the maternally-inherited allele is methylated in~75% and unmethylated in~25% of individuals. Methylation at this locus is thought to be influenced by the aforementioned pre-/perinatal environmental factors, as well as by genetic variants, the latter possibly due the role of CTCF (transcription factor, CCCTC-binding) in imprinting via the influence of rs2346018. A recently published study of genome-wide mQTLs in 27,750 European participants [18] revealed relatively weak associations between genetic variants and methylation at cg26328633 (CpG site identified as part of the strong heritable cluster in our Australian families), but no data were available for neighbouring CpGs, and no apparent association was found for other single nucleotide polymorphisms (SNP) of interest within this region, including rs2346018. No mQTLs were found for any of the most heritable VTRNA2-1 DNA methylation marks in a previous genome-wide assessment that included~2000 participants [19].
In this study, our aims were three-fold: first, to sequence the VTRNA2-1 region to assess the presence of rare genetic variation at this locus; second, to conduct a genomewide assessment of mQTLs and SNP-based methylation heritability in the VTRNA2-1 region (previously identified heritable marks); third, to assess whether any genetic variants associated with DNA methylation in this region contribute to the previously observed associations with breast cancer risk.

Results
This study used data from (1) 179 participants of multiple-case breast cancer family studies to assess the presence of genetic variants in the VTRNA2-1 heritable region, (2) 4500 participants in a prospective study to assess mQTLs and SNP-based heritability (genome-wide and focusing on cis-variants), and (3) 2141 participants in breast cancer family-based studies to adjust the VTRNA2-1 results of our previous publication [1] for a nearby SNP.
The characteristics of participants in the prospective cohort (Melbourne Collaborative Cohort Study) and multiple-case family (Australian Breast Cancer Family Registry and Kathleen Cuningham Foundation Consortium for research into Familial Breast Cancer) studies (see Methods Section) included in this study are shown in Table 1. In the prospective cohort study, the majority were male, aged between 50 and 70 and never or former smokers. The distribution of methylation beta values at the VTRNA2-1 region is shown in Figure 1b,c. Nine genetic variants were identified via sequencing (Figure 1a) in members of multiplecase breast cancer families. Of these, one was rare (identified in only one participant) and was excluded from further analysis. The other eight variants were used to estimate carrier probabilities in members of multiple-case breast cancer families (see Methods Section; 4.4.4. Associations with breast cancer risk). All eight variants were available after genotype imputation in the prospective cohort study (see Methods Section; 4.2. Genetic and methylation data) and therefore included in the mQTL analysis. None of these variants were found to directly overlap with the most heritable VTRNA2-1 methylation site (cg06536614).  Methylation values at the five most heritable methylation marks (cg06536614, cg26328633, cg25340688, cg26896946 and cg00124993) were highly correlated (in the MCCS all r ≥ 0.88); we therefore focused on cg06536614, which was found to be the most heritable methylation mark [1]. The percentage methylation at cg06536614 was lower than 50% for 1946 (43%) of participants and lower than 60% for 4414 (98%) of them; 23% had less than 30% methylation ( Figure 2). The mQTL analysis was performed in 4500 participants in the prospective, populationbased study (MCCS) for a total of 10,484,498 genetic variants; a cis-mQTL analysis was then carried out by focusing on genetic variants within 50 kb of the most heritable methylation mark (cg06536614). Genome-wide, we found no evidence that any included genetic variant was associated with VTRNA2-1 methylation (all p > 5 × 10 −9 ), with similar results obtained for the M-value (logit transformation of beta value) or RINT (rank-based inverse normal transformation, which was applied previously in the context of DNA methylation analyses and provides a Gaussian methylation distribution, which is not always the case for the Mvalues) transformation (Supplementary Table S1 and S2, respectively, showing the 100 CpGs with smallest p-values). Results from the cis-mQTL analysis, i.e., genetic variants within 50 kb of the CpG, are shown in Supplementary Table S3 (M-values) and Supplementary  Table S4 (RINT-values). There were 43 variants with p < 1.5 × 10 −4 , indicating evidence of genetic influences on VTRNA2-1 methylation (Table 2 and Figure 3). The strongest evidence of association was observed for rs2190622 (p = 5 × 10 −6 ) ( Figure 2), although this association was not substantially stronger than for other neighbouring SNPs. A significant association was also observed for rs2346018 (p = 8 × 10 −5 ) (Figure 2), which was previously reported to modify methylation in this region [7]. Consistent with this, associations appeared stronger for the variants located closer to the CpG of interest (Supplementary Table S3 and Figure 3). Associations were qualitatively similar but appeared somewhat stronger for the RINT-transformed values (strongest hit: rs2190622, p = 8 × 10 −8 ; rs2346018: p = 8 × 10 −6 and Supplementary Table S4 and Figure 3). These variants appeared to explain very little of variation in VTRNA2-1 methylation ( Figure 2) with a variance explained ranging from 0.33% to 0.47%. The findings were similar when restricting the analyses to MCCS participants selected as controls (not shown).   Table 2. All variant names, positions, and quantitative results are shown in Supplementary Table S3 and S4.
Results from SNP-based heritability analyses (i.e., taking into account > 1M variants) are shown in Table 3. These were consistent with the mQTL analyses in showing minimal influence of SNPs on methylation in the region: h 2 = 0, 95%CI: −0.14 to 0.14. These results were virtually the same when the RINT transformation was used instead of M-values (Table 3) or restricting the analyses to MCCS participants selected as controls (not shown). Non-null or high heritability was observed for methylation sites distant from the heritable VTRNA2-1 region (>2-8 kb, Table 3).

Discussion
Our study provides further evidence that DNA methylation at VTRNA2-1 is minimally influenced by genetic factors, and thus, the Mendelian-like inheritance of germline DNA methylation at this locus is likely to be via a true epimutation mechanism rather than via a mQTL. Therefore, the "missing heritability" (approximately the difference between familybased and SNP-based heritability) appears to be substantial, which confirms the findings of Joo et al. [1]. Genetic variants are therefore unlikely to fully explain any associations between methylation at this locus and disease risk, including breast cancer as we found previously in the context of a multiple-case breast cancer family study.
Some genetic variants at this locus were statistically significantly associated with blood DNA methylation (carriers of the minor allele being less likely to show hypomethylation), but the effect sizes and variance explained were small (variance explained 0.4% to 0.5% for the strongest individual associations). Rs2346018 [6,7], which was previously implicated with methylation at this locus, showed a significant association with methylation in our data, thereby confirming it may exert a small influence. Although our study does not allow disentangling which of the identified SNPs might causally affect DNA methylation, it should be noted that (1) rs2346018 was one of the SNPs with smallest p-value among the 334 tested, and (2) many SNPs that appeared most strongly associated with methylation did not have a clear functional interpretation, for example rs2190622, the strongest observed association, is located in an intergenic region and has unknown regulatory function (Supplementary Table S5).
In their genome-wide assessment of imprinting in the methylome, Zink et al. concluded that VTRNA2-1 is an example of a region with polymorphic imprinted methylation unrelated to SNP genotypes. In our study, a large number of participants (45%) had a percentage of methylation between 50% and 60%. Although this might appear to be inconsistent with maternal imprinting, this could indicate the limitation of the HM450 assay to measure DNA methylation with sufficient precision or DNA methylation at the paternal allele accumulated over the lifetime. However, technical validation using Pyro-Mark (Pyrosequencing technical validation) produced similar methylation values, and we did not find evidence that VTRNA2-1 methylation was influenced by factors for which our data shows widespread methylation changes, such as age, sex, country of birth, or other factors shown to strongly affect DNA methylation such as tobacco smoking, alcohol consumption or body mass index [20][21][22]. In this study, we did not have information on early-life factors. Although several factors in utero and early in life were shown to modulate DNA methylation at VTRNA2-1, none of the findings presented for season of conception [15], maternal folate during pregnancy [16], or gestational famine exposure [23] appear to fully explain methylation variation in this region; the mechanisms by which VTRNA2-1 methylation is inherited therefore appear to be essentially non-genetic and only partially explained by the factors studied in the literature so far.
Previous studies of mQTLs and heritability have only provided a partial assessment of genetic influences on VTRNA2-1 methylation. The study by Gaunt et al. [24] also did not report SNP-based heritability but assessed mQTLs. Various potential trans-mQTLs were identified but evidence of associations appeared weak (all p > 3 × 10 −8 , http:// www.mqtldb.org/search.htm, accessed on 1 December 2020). The large meta-analysis by Min et al. only reported significant associations for cg26328633 based on 27,750 European participants; 12 SNPs located close to VTRNA2-1 had a p-value ranging from 10 −22 to 10 −26 (http://mqtldb.godmc.org.uk/, accessed on 1 December 2020), but the effect sizes were virtually the same as those obtained in our study (beta~0.1) so that the variance explained was likely similarly small. Results for rs2346018 or other variants of interest were not provided. In the study by McRae et al. [25], the family-based heritability (peripheral blood leukocytes) of our five most heritable marks associated with breast cancer was 0.50. The study by van Dongen et al. [26] found very high twin-based heritability (~0.97) for the five VTRNA2-1 CpGs but SNP-based heritability was reported as "NA" (due to convergence problems using GCTA, which might mean those values were in fact equal to zero). The p value cut-off we used for detecting mQTLs genome-wide was conservative (strict Bonferroni correction, p = 5 × 10 −9 ), but using other cut-offs commonly used in genome-wide association studies, such as p = 5 × 10 −8 , would have resulted in the same conclusion (M-values: one intergenic variant in chromosome 1, p = 3 × 10 −8 and explaining 0.7% of methylation variability; RINT values: no variant with p < 5 × 10 −8 ). It is possible that some true associations were not detected, but these would likely have a weak influence on DNA methylation. Other studies used various cut-offs to declare statistically significant mQTLs, e.g., p = 10 −14 in Gaunt et al. [24], p = 10 −11 in McRae et al. [19], and p = 10 −8 in Min et al. [18], but these studies carried out more tests because they investigated >400,000 CpGs.
Although findings are difficult to compare across studies, they appear to be consistent with ours in showing weak influences of genetic variants on VTRNA2-1 methylation. Consistent with this, the associations of DNA methylation with breast cancer risk after adjustment for rs2346018 were only slightly attenuated. It should also be noted that none of rs2346018 or other sequenced genetic variants were found to be associated with breast cancer risk in the largest genome-wide association studies to date [27,28].
Although all methylation measures in this study were made on blood samples, it is worth noting that VTRNA2-1 methylation has been implicated as playing a role in, or being influenced during, carcinogenesis. Fort et al. [29] sought to find direct association between VTRNA2-1 transcript levels and methylation of its promoter in prostate tumour samples. Average VTRNA2-1 promoter methylation was found to be substantially increased in both low-grade and metastatic tumour tissue compared with normal prostate tissue. Additionally, average VTRNA2-1 promoter methylation appeared to correlate with Gleason score, clinical T-value and biochemical relapse [29]. The levels of VTRNA2-1 transcript were found to be inversely correlated with average promoter methylation. The relationship between VTRNA2-1 and cancer growth appears to be tissue specific, with several studies suggesting a tumour suppressive role, e.g., cholangiocarcinoma [30], oesophageal carcinoma [31], small cell lung cancer [32], gastric cancer [33] and acute myeloid leukaemia [34], and some suggesting an oncogenic role, e.g., in endometrial cancer [35] and thyroid cancer [36].
We conclude that the genetic and non-genetic factors we investigated play a minimal role in explaining variation in blood DNA methylation at VTRNA2-1, so these are unlikely to play a strong role in observed associations between VTRNA2-1 methylation and disease risk. The mechanism of inheritance of DNA methylation in this region remains to be elucidated.

Prospective Cohort Study
The Melbourne Collaborative Cohort Study (MCCS) is a community-based study that recruited 41,513 participants in 1990-1994 [37]. Several nested case-control studies were conducted to assess associations between DNA methylation in blood and the risk of eight types of cancer. Incident cases were matched to controls on age, sex, country of birth and sample type (buffy coats/dried blood spots/peripheral blood mononuclear cells) using incidence density sampling [37,38]. We also used questionnaire-collected data on smoking and alcohol consumption [20][21][22], measures of body mass index and derived a healthy eating index using a validated 121-item food frequency questionnaire [39,40].

Multiple-Case Breast Cancer Families
A total of 210 individuals from 25 multi-generational multiple-case breast cancer families, including 20 from Kathleen Cuningham Foundation Consortium for research into Familial Breast Cancer (kConFab) and 5 from the Australian Breast Cancer Family Registry (ABCFR)), were included in this study [41][42][43]. Among these family members, there were 87 breast cancer cases and 123 unaffected relatives.

Genetic and DNA Methylation Data
The VTRNA2-1 region (GRCh37, ch5:g.135414615-135417597) was screened in 179 of 210 individuals from a multiple-case breast cancer family by targeted-sequencing using a custom-designed HaloPlexHS panel (Agilent, Santa Clara, CA, USA). Libraries were prepared from blood-derived DNA according to the manufacturer's instructions and sequenced on a 2 × 150 bp high-output flow cell on the HiSeq3000 (Illumina, San Diego, CA, USA). Paired-end reads were aligned to the human reference genome GRCh37 using BWA-mem 0.7.17 [44]. Adapter sequences were removed, and unique molecular indices were marked for downstream read-deduplication using the Agilent Genomics NextGen Toolkit (Agilent, Santa Clara, CA, USA). Target coverage was calculated using bedtools v.2.27.1 [45], and variant calling was performed using VarDict v.1.7 [46]. Genetic variants filtering was performed as described previously [47]. Variant annotation was performed on variants with a read depth ≥30× and a variant allele frequency (VAF) ≥0.15, using VarSeq VSClinical v2.2 (Golden Helix Inc., Bozeman, MT, USA). Finally, an additional 23 participants in ABCFR/kConFab had genetic measures made using OncoArray, using the same method as for the MCCS for genotype imputation (see next paragraph).
Genome-wide genotyping was conducted on blood DNA samples from 12,584 MCCS participants using the Infinium OncoArray-500K BeadChip (Illumina, San Diego, CA, USA) [37,48]. Following previous standardised protocols [28], we imputed autosomal genotypes using the Michigan imputation server [49] and IMPUTE version 2 [50] with the 1000 Genomes Project dataset (phase 3) as the reference panel. The genotype probabilities from imputation were used to hard-call (uncertainty < 0.1) the genotypes for variants with an imputation info score > 0.3. For the current analysis, we included 4748 participants for whom DNA methylation data was also available. We then retained the hard-called variants with minor allele frequency > 0.001, missing genotype rate < 0.2 and Hardy-Weinberg equilibrium p-value > 10 −6 . Furthermore, to avoid bias due to confounding by shared environment among close relatives, participants were removed based on relatedness, i.e., excluding one participant randomly selected from any pair with a genetic relationship ≥ 0.05 (4th-degree or closer relationship) [51,52]. This procedure also removed duplicated methylation samples (genetic relationship = 1) [38]. After these quality control steps, 4500 paired genetic-methylation samples were retained (including 2228 cancer cases and 2272 controls) and 10,484,498 genetic variants (including 9,551,474 SNPs) for the analysis.
For all samples, DNA methylation was measured using the HumanMethylation450 (HM450) BeadChip (Illumina, San Diego, CA, USA) using methods described previously [1,38,53]. We used methylation M-values as their distribution is usually closer to Gaussian than methylation beta-values [54]. As a sensitivity analysis, we performed a more direct normalization of beta-values using rank-based inverse normal transformation (RINT) which was applied previously in the context of DNA methylation analyses and provides a Gaussian methylation distribution, which is not the case for the M-values [55,56].

Technical Validation of Methylation Measures Using Pyrosequencing
Pyrosequencing (PSQ) conducted on the PyroMark Q48 (Qiagen, Hilden, Germany) was used to validate the methylation measures made on the HM450 assay. DNA was extracted as described previously [1] and bisulfite converted using Zymo (Zymo Research, Irvine, CA, USA). Forward, reverse and sequencing PCR primers were designed using PyroMark Assay Design Software version 2.0. (Forward primer sequence: 5Bios 5 -G/GGAGGAATTGAGAGTTTTTTTAGGATA-3 ; Reverse primer sequence: 5 -CCTTCAA-AATAACACCAACTTATATTATCA-3 ; Sequencing primer sequence: 5 -ACATAAAAAAA-TCAATAAACACC-3 ) to target cg04481923 and were synthesised by Integrated DNA Technologies (Coralville, IA, USA). The EpiTect PCR control DNA set (Qiagen, Hilden, Germany), which includes a completely methylated and a completely unmethylated bisulfiteconverted control DNAs, was used to generate a standard curve. The EpiTect control DNAs were mixed in known ratios (0, 0.25, 0.50, 0.75 and 1) and run with ten test samples, along with a non-converted DNA sample and a no-template control. The PyroMark PCR cycling protocol was as follows: denaturation for 15 min at 95 • C, then 45 cycles of 30 s at 94 • C, 30 s at 56 C, and 30 s at 72 • C, then final extension for 10 min at 72 • C. For each sample, the raw percentage of methylation was determined at cg04481923 (VTRNA2-1) and calibrated using the EpiTect control DNAs standard curve.

Assessment of mQTLs
After QC of methylation and genetic data, 4500 participants (including 2228 cancer cases and 2272 controls) and 10,484,498 genetic variants (including 9551,474 SNPs) were available for the analysis. We first removed factors that may confound DNA methylation values using linear mixed models with methylation M-values (or RINT-values) as the outcome and as covariates: age, sex, sample type, white blood cell proportions (estimated using the Houseman algorithm [57,58]) and 20 genetic principal components to account for population structure/ancestry as fixed effects; and as random effects: study, plate and slide of the assay. Our sample therefore included both cancer cases and controls. The inclusion of cancer cases may bias mQTL associations because of collider bias [59]; collider bias is usually considered to be small [60,61], and this may be particularly true in our setting because no strong associations of individual methylation markers with cancer risk were observed, but we nevertheless assessed consistency of associations in controls by analysing them separately.
As slightly over 10M genetic variants were tested, we used the Bonferroni correction for multiple testing and considered associations with a p-value less than 0.05/10 −7 = 5 × 10 −9 to be potential true signals. Further, because cis acting genetic variants are considerably more likely than trans acting variants to influence DNA methylation, we considered all SNPs within 50 kb pairs of the methylation sites analysed. A total of 334 variants were identified, so we corrected the cis-mQTL analyses for multiple testing using the Bonferroni cut-off p = 1.5 × 10 −4 (0.05/334).

SNP-Based Heritability
The univariate genome-based restricted maximum likelihood (GREML) method [62,63] was used to estimate the SNP-based heritability of methylation values in the sample of 4500 participants and a subsample of 2272 controls, respectively. The M-values (or RINTvalues) after removing confounding effects were used as phenotypes in these analyses. We used only 1050,921 HapMap3 SNPs as they have been shown reliable and robust to bias in estimating SNP-based heritability and genetic correlations [64][65][66]. A genetic relationship matrix based on these SNPs was created and implemented in GREML. The heritability analyses were performed using the software GCTA [62].

Association of Non-Genetic Factors with VTRNA2-1 Methylation
We used mixed linear regressions similar to our previous publications [20,21] to assess the association of age, sex, BMI, smoking, alcohol consumption with methylation at VTRNA2-1. This analysis was undertaken using the same set of 4500 participants used for the genetic analyses, as well as separately in MCCS controls only.

Associations with Breast Cancer Risk
Cox proportional hazards survival analysis was used to test for associations between variants associated with VTRNA2-1 methylation and breast cancer risk using all participants from the 25 multi-generational multiple-case breast cancer families (n = 2141). This analysis was based on the phenotype and relationships data of these 2141 participants, and the methylation and genetic data on 202 of them. Unobserved methylation and SNP data were replaced by estimated carrier probabilities using the methods presented in [1]. As the families in this study were ascertained because they each contained multiple breast cancer cases, and no adjustment for this ascertainment criterion was made, hazard ratio estimates are biased, but since the ascertainment criterion has no effect on the test statistic under the null hypothesis, the p values for association with breast cancer are valid. These p values were based on the likelihood ratio test, not the Wald test, so variances for the hazard ratios were not needed and hence were not estimated using either standard maximum likelihood or robust variance. The same models as in Joo et al. [1] were performed, with additional adjustment for genetic variant carrier probabilities at rs2346018, which was one of the strongest mQTLs in this study (see Results), and previously reported to influence VTRNA2-1 promoter methylation [7]). Similar results were obtained after adjustment for any of the other eight variants for which carrier probabilities could be estimated (data not shown).