DNA methylation signature on phosphatidylethanol, not on self-reported alcohol consumption, predicts hazardous alcohol consumption in two distinct populations

The process of diagnosing hazardous alcohol drinking (HAD) is based on self-reported data and is thereby vulnerable to bias. There has been an interest in developing epigenetic biomarkers for HAD that might complement clinical assessment. Because alcohol consumption has been previously linked to DNA methylation (DNAm), we aimed to select DNAm signatures in blood to predict HAD from two demographically and clinically distinct populations (Ntotal = 1,549). We first separately conducted an epigenome-wide association study (EWAS) for phosphatidylethanol (PEth), an objective measure of alcohol consumption, and for self-reported alcohol consumption in Cohort 1. We identified 83 PEth-associated CpGs, including 23 CpGs previously associated with alcohol consumption or alcohol use disorder. In contrast, no CpG reached epigenome-wide significance on self-reported alcohol consumption. Using a machine learning approach, two CpG subsets from EWAS on PEth and on self-reported alcohol consumption from Cohort 1 were separately tested for the prediction of HAD in Cohort 2. We found that a subset of 143 CpGs selected from the EWAS on PEth showed an excellent prediction of HAD with the area under the receiver operating characteristic curve (AUC) of 89.4% in training set and 73.9% in validation set of Cohort 2. However, CpGs preselected from the EWAS on self-reported alcohol consumption showed a poor prediction of HAD with AUC 75.2% in training set and 57.6% in validation set. Our results demonstrate that an objective measure for alcohol consumption is a more informative phenotype than self-reported data for revealing epigenetic mechanisms. The PEth-associated DNAm signature in blood could serve as a robust biomarker for alcohol consumption.


Introduction
Hazardous alcohol drinking (HAD) is detrimental to health and is highly correlated with medical comorbidities and psychiatric diseases [1,2]. Diagnosing HAD is challenging due to a lack of stable and objective measures for chronic heavy alcohol consumption [3]. Phosphatidylethanol (PEth) is a lipid metabolite of ethanol formed from phosphatidylcholine in erythrocytes and has been proposed as a biomarker for alcohol consumption. Compared with selfreported data, PEth reliably detects ethanol levels up to 21 days after the last drink [4], and the PEth level is highly correlated with alcohol consumption [5]. However, the clinical applicability of PEth is limited because its half-life is~4-7 days and its window of detection is considered to be 21 days [6]. Thus, other longer-term biomarkers for alcohol consumption are needed to inform clinical practice.
Epigenetic signatures have emerged as attractive biomarkers for complex diseases such as cancers and neurodegenerative diseases [7]. Epigenetic markers may reflect environmental exposures, including alcohol consumption. Among these epigenetic markers, DNA methylation (DNAm) biomarkers are particularly attractive because they are relatively stable and capture an early stage of pathophysiological changes [8,9]. A recent longitudinal study on DNAm showed that most DNA methylome changes occurred 80-90 days before clinically detectable glucose elevation [10], suggesting that DNAm is involved in an early stage of diabetes. Finally, epigenetic modifications can be reliably detected in noninvasive fluids and biospecimens [11]. Thus, the utility of epigenetic alterations has motivated the biomarker research field to develop epigenetic signatures derived from easily accessible cells for clinical use [12][13][14].
DNAm markers are emerging as diagnostic biomarkers in many areas of medicine and are applied to predict complex diseases [15]. For example, DNAm markers on the promoters of several genes, including BMP3, NDRG4, and SPEPT9, in blood or stool samples have been approved by the Food and Drug Administration as biomarkers for colorectal cancer screening [16]. DNAm markers also distinguish smokers and nonsmokers [17,18]. However, we do not yet have validated DNAm biomarkers for the diagnosis of HAD.
DNAm is directly altered by HAD in the following manner. HAD often causes folate and vitamin B deficiency, resulting in the reduction of S-adenosylmethionine (SAM). DNAm is modulated by DNA methyltransferase which transfers a methyl group from SAM to the 5-position of cytosine in the context of cytosine-phosphate-guanine (CpG) dinucleotide. Reduced methyl transfer reaction cofactors (folate and vitamin B) reduce methyltransferase activity that may lead to alteration in DNAm. Recent studies have shown that alcohol consumption modifies DNAm [19] in animals and in the human epigenome from blood, liver, and saliva cells [17,[20][21][22][23][24]. As a result, DNAm in peripheral cells can serve as a robust biomarker for HAD.
Epigenome-wide association studies (EWAS) have identified hundreds of DNAm CpGs in blood for alcohol consumption [25][26][27][28], alcohol use disorder (AUD) [29,30], stress-related alcohol consumption [31], and fetal alcohol syndrome [32][33][34][35]. A large number of CpGs in the blood have recently been reported to have associations with dietary folate and alcohol intake [36]. CpGs have been associated with alcohol consumption in different cell types, ethnic groups, and phenotypic assessments [28,29,37]. More than a dozen CpGs for alcohol phenotypes have been replicated. For example, cg11376147 on SLC43A1 has been linked to alcohol consumption and HAD diagnosis in several studies [17,28,29]. Thus, DNAm in blood has been proposed as a diagnostic and prognostic biomarker of alcohol consumption for clinical use [38]. For this purpose, a previous study identified a panel of 144 CpGs as biomarkers for alcohol consumption [29]. However, these CpGs have not been validated in independent studies.
Differentially methylated CpG sites have also been associated with differential gene expression for alcohol exposure in both animals and humans. Alcohol exposure is associated with hypomethylation in the promoter of the proprotein convertase subtilisin/kexin type 9 (PCSK9) gene [30] that is also correlated with PCSK9 expression for heavy alcohol consumption in humans and mice. Because PCSK9 is well known to regulate low-density lipoprotein cholesterol, DNAm alteration and dysfunction of PCSK9 is thought to be a mechanism for alcoholrelated abnormalities in lipid metabolism. Most recently, Gatta et al. [31] reported the hypermethylation of DNA 5methylcytosine at the promoter regions of NR3C1 (Nuclear Receptor Subfamily 3 Group C Member 1), the glucocorticoid receptor, that was correlated with the reduction of mRNA expression of NR3C1 in human brains with AUD. The expression levels of several stressresponsive genes within the NR3C1 gene network were also decreased in brain samples from individuals with AUD. This evidence further supports the feasibility of DNAm biomarkers for HAD that may have both clinical utilities and help elucidate underlying pathophysiologic mechanisms of heavy alcohol consumption.
One of the limitations of previous EWAS is that alcohol consumption was assessed by self-report, which may lead to inaccurate assessment and introduce bias [29,39,40]. A selfreported phenotype may, in part, explain the discrepancy of EWAS findings on alcohol consumption or alcohol userelated phenotypes observed in previous studies. Objective measures such as PEth may improve the association signals for alcohol consumption in EWAS because PEth-associated DNAm markers are more proximal to the biological changes and pathological processes underlying HAD.
In this study, we hypothesized that the DNAm signatures associated with PEth would be a more robust predictor of HAD than self-reported drinking data. We conducted a 2stage study with the goal of identifying DNAm CpGs for alcohol consumption and linking the CpG features to HAD (N total = 1,549). We first identified CpGs associated with PEth in Cohort 1. Then, the informative CpGs were selected to predict HAD by using elastic net regularization (ENR) in a demographically and clinically independent sample (Cohort 2). We also compared the findings of DNAm markers for PEth with those for self-reported alcohol consumption. The analytical strategy is presented in Fig. 1.

Sample descriptions
Cohort 1 (N = 1,047): The DNA samples in Cohort 1 were from the Veterans Aging Cohort Study (VACS) [41]. Data were obtained from the patients after they provided written consent; data were collected via telephone interviews, focus groups, and full access to the national Veterans Affairs electronic medical record system. All subjects in this subset of the VACS cohort were men.
Samples in Cohort 1 were divided into a discovery set (N = 580) and a replication set (N = 467) for EWAS. A majority of discovery samples were HIV-positive (~85.34%), and all replication samples were HIV-positive.
Cohort 2 (N = 502): We recruited 502 HIV-negative healthy community volunteers who responded to advertisements placed either online or in local newspapers and at a community center in New Haven, CT [42]. Phenotypic assessment including alcohol consumption was obtained through the in-person interview. To reduce confounding effects, we excluded subjects who met the Diagnostic and Statistical Manual of Mental Disorders, 4th Edition (DSM-IVTR) (American Psychiatric Association, 1994) criteria for substance dependence on any drug or alcohol other than nicotine. Subjects with head injury or those who used prescribed medications for any psychiatric or medical disorders were also excluded.
All phenotypic data in Cohort 1 and Cohort 2 were obtained in the same time window as the blood draws for DNAm profiling. Genomic DNA was extracted from whole blood using a standard method [43]. The study was approved by the committee of the Human Research Subject Protection at Yale University and the IRB committee of the Connecticut Veteran Healthcare System.

Phosphatidylethanol (PEth) measurement
In this study, PEth was only measured in Cohort 1 using dried blood spot samples derived from frozen peripheral blood mononuclear cells stored at −80°C [5]. We analyzed the PEth levels from dry blood spots at the U.S. Drug Testing Laboratory (Des Plaines, IL) via LC-MS/MS, as described in previous studies [44][45][46]. The LC-MS/MS method has a high capacity and is costeffective and clinically reliable [46]. PEth can be detected at concentrations as low as 2 ng/ml. A study showed that the PEth value is linearly related to alcohol consumption [47]. In forensics, 20 ng/ml of PEth was used as a cutoff to detect harmful alcohol use [48]. The sensitivity of PEth has been reported to be 99% [47], with several studies showing the assay to have perfect specificity, including in the presence of liver disease and hypertension. We previously reported that PEth was highly correlated with the Alcohol Use Disorders Identification Test-Consumption (AUDIT-C, first three items of AUDIT) score from electronic records [49].

Definition of hazardous alcohol drinking (HAD)
Alcohol consumption was measured by both PEth and AUDIT-C in Cohort 1 and was only measured by AUDIT in Cohort 2 (Supplementary Table S1). In Cohort 1, HAD was defined as PEth level ≥20 based on previous studies and non-HAD was defined as PEth <20 [48]. HAD was corresponding to AUDIT-C score ≥4 and non-HAD was corresponding to AUDIT-C score <4 for men [5]. In the discovery set of Cohort 1, there were 166 HADs and 414 non-HADs. In the replication set of Cohort 1, there were 135 HADs and 332 non-HADs. In Cohort 2, alcohol consumption was assessed by a full scale of 10-item AUDIT with a total score of 40. HAD was defined as AUDIT ≥8 for men and ≥7 for women based on previous studies. Non-HAD was defined as AUDIT <8 for men and <7 for women [50]. There were 148 HADs and 354 non-HADs. Cohort 2 was divided into a training set (N = 402) and a testing set (N = 100), with an 80-20 split, for machine learning prediction of HAD. Demographic and clinical variables for HAD versus non-HAD participants in Cohort 1 and Cohort 2 are presented in Table 1.

DNA methylation and data quality control (QC)
In Cohort 1, DNAm for the discovery sample was profiled by using the Illumina Infinium HumanMethylation450 Beadchip (Illumina HM450K) (San Diego, CA, USA). DNAm for the replication sample was assessed by using the Illumina Infinium MethylationEPIC Beadchip (Illumina EPIC) (San Diego, CA, USA). In Cohort 2, DNAm was measured by using Illumina HM450K. All samples in Cohorts 1 and 2 were processed at the Yale Center for Genomic Analysis [14]. After QC (details in Supplementary Information), in Cohort 1, a total of 437,722 CpGs from 450K array remained in the discovery sample and 846,604 CpGs from EPIC array remained for the replication sample. A total of 48.26% common CpGs (408,583) were analyzed in meta-analyses. In Cohort 2, we applied the same QC criteria. A total of 437,722 CpGs remained for analysis.

Discovery and replication EWAS in Cohort 1
EWAS were separately performed to test the association of each CpG methylation with PEth and AUDIT-C score in the discovery and replication samples. To adjust for significant global confounding factors, we followed a comprehensive analysis pipeline developed by Lehne et al. [51]. Since previous studies have shown that a large number of CpGs were significantly associated with age [52], smoking status [12], race [53], HIV status, and HIV-1 VL [14], these variables were adjusted in the models. The cell proportions of six cell types were also adjusted in the models [54]. The log10 of viral load (log 10 VL) and ART adherence were adjusted in the replication sample. In addition, a recent study reported by Jiao et al. [55] demonstrated that sample position affected the measurement of DNAm in Illumina methylation arrays and may introduce biases and increase batch effects. Thus, we adjusted positional effects in the models to further address confounding effects. Epigenomewide significance was set at a Benjamini-Hochberg false discovery rate (FDR) < 0.05 in the discovery sample. Significance in the replication sample was set at p < 0:05 number of CpGs being tested . Analytical models present as the following:

First generalized linear model
For discovery,

Principal component analysis (PCA) of intermediary residuals
We then performed a PCA on the resulting regression residuals. The top five principal components (PCs) on the residuals (PC1 residuals ,…,PC5 residuals ) were adjusted in the final model.
A final generalized linear model for identifying differential methylation We performed a final generalized linear regression analysis for each CpG predicting the β as a function of the In(PEth) adjusted for confounders and the top five residual PCs derived from the model above.
The same models were also used for EWAS on AUDIT-C score in discovery and replication samples, where the independent variable ln(PEth) was replaced by the AUDIT-C. To evaluate whether the residual DNAm was adjusted for confounding effect in the above model, we tested the correlations between the top 30 PCs and position, batch, age, race, smoking status, WBC, and six cell-type proportions in before and after QC, respectively.

Meta-analysis of EWAS in Cohort 1
An EWAS meta-analysis was conducted by combining the findings from the discovery and the replication stages. For each CpG, we obtained effect size estimates and p values from the two samples and weighted the effect size estimates by their estimated standard errors. Then, the summary statistics of the two samples were combined using a samplesize weighted meta-analysis using the METAL program [56]. Epigenome-wide significance was set at an FDR < 0.05.
As a comparison of meta-EWAS findings, we conducted a single EWAS in a total of 1,047 samples combining the discovery sample and replication samples together. The batch effect and positional effect were removed by using ComBat [57]. The analytical models and covariables were the same as described above.

Polygenic methylation score (PGMS)
We constructed a PGMS for each individual as a weighted sum of the CpG β values using the effect size estimated from the EWAS as weights [13]. In detail, the PEth-related CpGs identified in the meta-analysis were chosen to construct the PGMS. Then, the PGMS was applied to establish a prediction model for HAD in Cohort 2.
M i : the PGMS of individual i; a j : the estimated coefficient for CpG probe j; β ij : the methylation β value for individual i at CpG probe j.
Adjusted R 2 and incremental adjusted R 2 We used the adjusted R 2 to estimate the phenotypic variances explained by the DNAm. The adjusted R 2 represented the percentage of variation explained by only the independent variables that affected the dependent variable. Here, the adjusted R 2 was the proportion of the variance of the PEth values, AUDIT-C scores, or AUDIT scores that were explained by the individual CpG or the linear combination of CpGs.
We applied the incremental adjusted R 2 (incremental R 2 ) as one of the parameters for feature selection as described below. The incremental R 2 was used to determine whether a new predictor increases the predictive ability above and beyond that provided by an existing model. It was calculated for each selected CpG or the linear combination of selected CpGs.

Feature selection using elastic net regularization (ENR)
CpG features were separately preselected from the EWAS results on PEth and on AUDIT-C in Cohort 1. The selected features were used to evaluate the prediction of HAD. Using the ENR method, we performed 10-fold cross-validation for feature selection in the training sample of Cohort 2. Here, we randomly selected 80% of the samples in Cohort 2 and cross-validated them to obtain the values for the ENR tuning parameters. The following steps were taken to select the CpG features and to evaluate their performance.
Step 1. Preselection CpGs. Because DNAm of CpGs under the epigenome-wide significance threshold may collectively account for phenotype variation and may improve prediction of a phenotype, we preselected PEth-associated CpGs with a meta p < 1E-04 from the meta-EWAS in Cohort 1 for both PEth and AUDIT-C. The preselected CpGs were used to establish the predictive model in the training set of Cohort 2.
Step 2. Importance ranking CpGs. In the training set of Cohort 2, we performed an ENR for feature selection among the preselected CpGs. We extracted the coefficients for the model with the lambda value corresponding to the minimum mean cross-validated error. This procedure was repeated N times. We excluded the CpGs with the percentage of zero coefficients larger than 95%. All selected CpGs were ranked according to the summation of the absolute value of the N coefficients.
Step 3. Model building by ENR in the training set. CpG features were selected based on the area under the receiver operating characteristic curve (AUC) and the incremental R 2 for different numbers of CpG sets. The model with the best performance was determined, and the optimal values of the parameters in the best model were found by performing cross-validation in ENR.
Step 4. Model performance testing in the testing set. The performance of the CpG features selected from the training set was evaluated in the testing set using AUC, balanced accuracy, and incremental R 2 .
A sensitivity test using different cutoffs of p values was performed to select the model with the best performance. Different sets of CpGs with p value <1E−06, <1E−05, <1E −4, <1E−3 were selected for feature selection in the training sample and model evaluation in the testing sample. The CpG set with the best performance was determined in the final model.
All analyses were performed using R software (https://www.r-project.org/). ENR was performed using the function "cv.glmnet" in the "glmnet" package.

Biological interpretation of the prediction model
Gene enrichment analysis was performed using the CpGs from the final prediction model to understand the underlying biological significance. We applied the web-accessible, gene annotation term-based Database for Annotation, Visualization and Integrated Discovery (DAVID) for gene enrichment analysis (http://david.niaid.nih.gov) [58]. The expanded DAVID Knowledgebase integrates almost all major and wellknown public bioinformatics resources [59]. A significant pathway was set as an FDR < 0.05.

Results
EWAS identifies new DNA methylation CpGs for PEth but not for self-reported alcohol consumption Two analyses of EWAS on PEth values and on AUDIT-C scores were separately conducted in Cohort 1. Phenotypically, as expected, PEth level and AUDIT-C score were highly correlated (r = 0.45, p < 2.00E−16) ( Supplementary  Fig. S1a). Compared with the non-HAD group, the HAD group had a greater AUDIT-C score (p = 5.42E−132) and a higher level of PEth (p = 3.47E−33) ( Supplementary Fig.  S1b). Demographic and clinical variables are presented in Table 1.

Discovery EWAS on PEth and on AUDIT-C
Prior to data QC, we found 10 PCs out of 30 PCs in DNAm was significantly correlated with position and batch effect, 4 PCs correlated with WBC, 2 PCs correlated with CD8T, 1 PC correlated with CD4T, and 2 PCs correlated with monocyte (p < 0:05 30 ¼ 1:67E À 03) ( Supplementary Fig.  S2a). After adjusted confounding effects in the model, residual methylation showed no correlations with position, batch, age, race, smoking status, WBC, or six cell-type proportions ( Supplementary Fig. S2b), suggesting that the EWAS findings below are unlikely contributed by nonspecific variables in the cohort.
As expected, the analysis of the EWAS on AUDIT-C scores revealed no CpG reaching epigenome-wide significance in the replication sample. Only one of nine CpGs associated with PEth were associated with AUDIT-C score (p < 0:05 9 ¼ 5:56E À 03) (cg11376147 in SLC43A1 with p = 2.74E−04) and showed the same direction as the association of PEth.

Meta-analysis
A meta-analysis revealed 83 epigenome-wide significant CpGs on PEth (FDR = 4.94E−06~4.97E−02) ( Table 2 and Fig. 2a). Of note, despite removing batch effects and position effects, a single EWAS conducted in the combining the discovery and replication samples revealed a greater number of λ than meta-EWAS (1.442 for the EWAS for combining samples and 1.130 for meta-EWAS) (Supplementary Fig. S5), suggesting that meta-EWAS was less likely inflated and biased than the single EWAS.
Interestingly, even with an increased sample size in the meta-analysis, we found no epigenome-wide significant CpG site of the meta-EWAS on AUDIT-C scores (Fig. 2b).
We further tested the correlation between the β values of the 12 CpGs with Bonferroni significance and PEth. All 12 CpGs were significantly correlated with PEth levels with p < 0:05 12 ¼ 4:17E À 03 after the model was adjusted for confounding factors (Fig. 2c), 4 of the 12 CpGs were positively correlated with PEth, and the remaining 8 CpGs were negatively correlated with PEth.

PEth-associated CpG sites improve the prediction of HAD in Cohort 1
Because PEth itself was highly correlated with AUDIT-C scores and differed significantly between the HAD and the non-HAD groups, we were interested in whether PEthassociated CpG DNAm improved the prediction of HAD compared with the prediction of HAD using PEth alone. We found that the AUC was 74.2% for PEth alone, 76.8% with the 12 Bonferroni significant CpGs and PEth, and 87.2% with the 83 epigenome-wide significant CpGs and PEth ( Supplementary Fig. S6). Thus, DNAm features improved the prediction of hazardous alcohol consumption compared with PEth alone in the same cohort.

PGMS derived from 83 PEth-associated CpGs is correlated with alcohol consumption in an independent sample
To be consistent with the analysis in Cohort 1, we performed an EWAS on AUDIT-C score in Cohort 2. We found no epigenome-wide significant CpG for AUDIT-C. An EWAS for a full scale of AUDIT score also revealed no significant CpG.
A PGMS constructed from the 83 PEth-associated CpGs was highly correlated with the self-reported 10-item AUDIT score in Cohort 2 (r = 0.40, p = 5.47E−19). The incremental R 2 of the association between the PGMS corresponding to 83 PEth-related CpGs and the 10-item AUDIT score was 0.0976, which implied that the PGMS explained 9.8% of the variance of the full AUDIT score in an independent population (Supplementary Fig. S7a).
We further tested whether the PGMS derived from the PEth-associated CpGs was separately correlated with selfreported alcohol consumption (AUDIT-C, first three items of AUDIT) and self-reported problem alcohol drinking behaviors (AUDIT-P, item 4-10 of full AUDIT). We found that the PGMS was significantly correlated with AUDIT-C score (r = 0.36, p = 3.36E−15) ( Supplementary Fig. S7b) and AUDIT-P score (r = 0.34, p = 1.29E−10) ( Supplementary  Fig. S7c). The correlation of the PGMS was slightly stronger with the AUDIT-C score than with the AUDIT-P score.

PEth-associated DNA methylation CpG sites predict HAD in Cohort 2
We found no statistically significant difference in the characteristics between the training set and the testing set in The CpGs that are significant using African Ancestry sample in Liu et al. [29]. Table S3). As shown in Supplementary Fig. S8, a sensitive test showed that the best performance model was a panel of CpGs preselected at p < 1E −04 assessed by AUC and incremental R 2 . Of note, although a larger cutoff value, e.g., p < 1E−03, showed a greater incremental R 2 , the AUC was less than the CpG set at the cutoff of p < 1E−04, which may be due to the increased background noise with a larger number of preselected CpGs at p < 1E−03. Therefore, the panel of CpGs with p < 1E−04 from the meta-EWAS in Cohort 1 were preselected for feature selection. A total of 259 CpGs were preselected to build a predictive model in the training set of Cohort 2. All 259 CpGs were ranked according to the summation of the absolute value of the N coefficients. As shown in Fig. 3a, a panel of 143 CpGs (Supplementary Table S4) showed the greatest AUC with 89.4% and the highest incremental R 2 with 19.3% in the training set. Therefore, a model of 143 CpGs was validated in the testing set.

Cohort 2 (Supplementary
In the testing set, we found that the 143 CpGs model showed an AUC of 73.9%, a balanced accuracy of 62.3%, and an incremental R 2 of 5.9% (Fig. 3b). The results show that the 143 selected PEth-associated CpGs enabled the good prediction of HAD. Notably, the panel of 143 CpGs included 44 epigenome-wide significant CpGs for meta-EWAS on PEth in Cohort 1.
Using the same approach for the analysis of feature selection of AUDIT-C-associated CpGs from Cohort 1 to predict HAD in Cohort 2, a panel of 18 CpGs were selected from 54 CpGs with p < 1E−04. In the training set, the AUC was 70.2%, and the incremental R 2 was 2.2%. In the testing set, the AUC was 57.6% (46.1-69.1%), and the incremental R 2 was 1.1%.

Discussion
Using samples from two distinct populations, we have demonstrated that an objective phenotype, PEth, is a robust phenotype for identifying DNAm in blood associated with HAD and that PEth-associated CpGs are predictive of HAD. We revealed 83 CpGs associated with PEth, while none of the CpGs were associated with self-reported alcohol consumption. A PGMS derived from the 83 CpGs explained 9.8% of the variance of alcohol consumption in a demographically and clinically independent sample. We further showed that the 83 CpGs combined with PEth improved 13% of AUC of predicting HAD compared with the AUC of predicting HAD by PEth alone. Importantly, we identified a panel of 143 CpGs that were relevant to PEth levels in a mostly HIV-positive sample and that predicted self-reported HAD in an HIV-negative sample. The 143 CpGs included five CpGs that were previously included in the DNAm biomarker panel for prediction of alcohol consumption and five CpGs were Bonferroni significant associated with alcohol consumption in an African Ancestry sample in Liu et al. (Supplementary Fig. S10) [29]. Interestingly, a panel of CpGs related to self-reported AUDIT-C score showed poor predictive performance for HAD. Together, these findings suggest that PEth-associated DNAm features, but not DNAm for self-reported alcohol consumption, is a robust biomarker in predicting hazardous alcohol consumption that may have potential clinical utility. Emerging evidence suggests that a set of epigenetic modification markers across different tissues is more stable and reproducible than we previously expected [60]. In this study, we replicated 24 CpGs that had previously reported associations with alcohol consumption or alcohol use disorders (Supplementary Table S4). For example, three promoter CpGs, cg19731612 on NSD1 [28,29], cg03523740 on TXLNA [28,29], cg18121224 on NSD1 [28], and cg00407659 on ANXA6 [28] that were associated with alcohol consumption in previous studies were also significantly associated with PEth in our study. In addition, we revealed multiple new PEth-associated CpGs that are located on the genes involved in tyrosine autophosphorylation, catalyzed phosphorylation of histones H3 and H2B (DYRK2) and the serine/threonine p21-activating kinases (PAK1), sequence-specific serine/arginine splicing factor (TRA2B) functions, and extracellular matrix protein (FBLN2). These results suggest that alcohol consumption alters DNAm on the genes involved in the cellular process and epigenetic programming. One intriguing question is whether the significant CpGs for HAD detected from the blood sample is relevant to methylation alteration by alcohol consumption in the brain. Several studies have reported associations of methylation of CpGs with alcohol Fig. 3 Feature selection using elastic net regularization (ENR) for hazardous alcohol drinking (HAD). a The area under the receiver operating characteristic curve (AUC) and the incremental adjusted R 2 (incremental R 2 ) of the selected CpGs using the ENR method. A set of CpGs associated with the natural logarithm of phosphatidylethanol (In (PEth)) in cohort 1 (p < 1.00E−04) was preselected for ENR analysis in training samples of Cohort 2. Incremental R 2 denotes the difference in adjusted R 2 between the model with the predicted variable and the model without the predicted variable. b The ROC curve for HAD prediction in the testing set of Cohort 2 using the 143 ENR-selected CpGs from the training samples. consumption in the human postmortem brain samples. Several of those CpGs showed nominal significance in our blood samples. For example, cg18362496 in H19 that was previously reported hypermethylation in AUD brain samples [61] showed a positive association with PEth in our blood sample (p = 0.03). A stress-related gene, KCNK6, that was previously associated with AUD in brain [31] was nominal significant in the same direction in our study (p = 0.04). However, the majority of significant CpGs for alcohol consumption differs between brain and blood samples. The discrepancy is not unexpected considering distinct methylome architectures between the brain and peripheral tissues. Although the findings do not elucidate the etiology of alcohol drinking behavior in brain, our results suggest a peripheral mechanism of how alcohol consumption changes the epigenome, which may lead to medical disorders. Given the inaccessibility of brain tissues in living humans, biomarkers from peripheral cells could be of benefit in the clinical care of HAD patients.
The 83 PEth-associated CpGs identified in a mostly HIVpositive population collectively explained 9.8% of the variance of HAD in an HIV-negative population, suggesting the stability of the DNAm effect of alcohol exposure. Notably, the 9.8% effect size of the PGMS in our study is comparable with the previously reported 12-13.8% effect size of a PGMS in a study with a larger sample size (N = 13,317) than our study [29]. We further showed that PGMS was not only significantly associated with recent alcohol consumption (AUDIT-C) (r = 0.36, p = 3.36E−15) but was strongly associated with the problematic consequences of alcohol use (AUDIT-P) (r = 0.34, p = 1.29E−10), further indicating that DNAm is a relatively stable marker for the long-term effects of alcohol consumption.
The reproducible CpGs suggest a robust, consistent epigenetic response to alcohol consumption that can serve as biomarkers for clinical use. Using a machine learning approach, we identified a set of 143 CpGs that enables the distinction of HAD and non-HAD individuals. One of the common challenges for machine learning prediction is model overfitting. We took several steps to address this concern: (1) feature preselection and selection were conducted in two different cohorts; (2) the processes of feature selection and model evaluation were carried out in the same cohort but in different sets without overlapping samples; and (3) we applied a newly developed machine learning ENR method to select features in a combination of 10-fold cross-validation. Compared with two traditional penalized regression methods, Ridge [62] and the least absolute shrinkage and selection operator (LASSO) [63], ENR has the advantage of selecting informative features without compromising predictive accuracy and has been shown to outperform both the Ridge and LASSO methods [64]. With these strengths of the analytical approach, we showed that a panel of 143 CpGs performed fairly well in the testing sample set.
Compared with the findings from the largest DNAm biomarker study for alcohol consumption up to date by Liu et al. [29], we found that a small proportion of CpGs is indeed overlapping between two studies. Despite many differences between ours and Liu's studies, e.g., sample size, phenotype assessment, DNAm profiling array, and analytical strategy, nine epigenome-wide significant CpGs are identical between ours and Liu et al.'s studies in African Americans. These nine overlapped CpGs are located in five genes: SLC43A1, FBLN2, HNRNPA1, CAND2, and GAS5. Five biomarker-CpGs are overlapping between the two studies. The overlapped CpGs are located on SLC7A11, DYRK2, TRA2B, NCOA2, and GPR133. Enrichment analysis suggests the overlapped CpGs are not discovered by chance across the two studies (T χ 2 = 2400 and p = 0 for epigenome-wide significant CpGs; T χ 2 = 486.4 and p = 0 for biomarker-CpGs). Therefore, the overlapped CpGs across two very different studies further underscore the stability and reproducibility of DNAm as a biomarker for alcohol consumption.
Several limitations should be considered in interpreting the current findings. (1) There was a lack of power to detect sex-specific associations between CpGs and HAD. It is well known that HAD in men and women is epidemiologically and mechanistically different. The individuals in Cohort 1 were all men and~50% of the individuals in Cohort 2 were women. These samples are insufficient to seek sex-specific DNAm markers. (2) The DNAm signatures were identified from whole blood samples that lacked cell-type-specific profiles. Future analyses using cell-type-specific CpGs may improve prediction performance. (3) The 143 CpGs in the DNAm signature were preselected from an HIV-positive sample, while the prediction model was built and validated in an HIV-negative sample. We expect to improve the predictive efficiency in a relatively homogenous sample in future studies. (4) Other psychiatric disorders such as depression are common in HAD, which might have confounded the findings. Validation of the prediction panel on other alcohol use-related phenotypes, e.g., alcohol use disorder, and address other psychiatric disorders are necessary to confidently claim the predictive performance and accuracy for clinical use.
In summary, to the best of our knowledge, this is the first study to demonstrate that PEth is a robust phenotype for detecting subtle DNAm changes associated with alcohol consumption compared with self-reported alcohol use data. PEth-associated DNAm markers predicted HAD with a good accuracy. These findings suggest that DNAm signatures may have clinical utility as biomarkers for alcohol consumption, and further development and testing of these biomarkers are warranted.

Data availability
Demographic variables, clinical variables, and methylation status for the VACS samples were submitted to the GEO dataset (GSE117861) and are available to the public. All codes for analysis are also available upon a request to the corresponding author.
Publisher's note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons. org/licenses/by/4.0/.