Gene expression correlates of advanced epigenetic age and psychopathology in postmortem cortical tissue.

Psychiatric stress has been associated with accelerated epigenetic aging (i.e., when estimates of cellular age based on DNA methylation exceed chronological age) in both blood and brain tissue. Little is known about the downstream biological effects of accelerated epigenetic age on gene expression. In this study we examined associations between DNA methylation-derived estimates of cellular age that range from decelerated to accelerated relative to chronological age ("DNAm age residuals") and transcriptome-wide gene expression. This was examined using tissue from three post-mortem cortical regions (ventromedial and dorsolateral prefrontal cortex and motor cortex, n = 97) from the VA National PTSD Brain Bank. In addition, we examined how posttraumatic stress disorder (PTSD) and alcohol-use disorders (AUD) moderated the association between DNAm age residuals and gene expression. Transcriptome-wide results across brain regions, psychiatric diagnoses, and cohorts (full sample and male and female subsets) revealed experiment-wide differential expression of 11 genes in association with PTSD or AUD in interaction with DNAm age residuals. This included the inflammation-related genes IL1B, RCOR2, and GCNT1. Candidate gene class analyses and gene network enrichment analyses further supported differential expression of inflammation/immune gene networks as well as glucocorticoid, circadian, and oxidative stress-related genes. Gene co-expression network modules suggested enrichment of myelination related processes and oligodendrocyte enrichment in association with DNAm age residuals in the presence of psychopathology. Collectively, results suggest that psychiatric stress accentuates the association between advanced epigenetic age and expression of inflammation genes in the brain. This highlights the role of inflammatory processes in the pathophysiology of accelerated cellular aging and suggests that inflammatory pathways may link accelerated cellular aging to premature disease onset and neurodegeneration, particularly in stressed populations. This suggests that anti-inflammatory interventions may be an important direction to pursue in evaluating ways to prevent or delay cellular aging and increase resilience to diseases of aging.

Multiple forms of psychopathology and life adversity, including posttraumatic stress disorder (PTSD), depression, alcohol-use disorders (AUD), and trauma exposure have been associated with advanced epigenetic age relative to chronological age (i.e., “accelerated aging; ” Han et al., 2018; Jovanovic et al., 2017; Marini et al., 2020; Wolf et al., 2016; Wolf et al., 2018a, Wolf et al., 2018b; 2019). This may, in turn, contribute to the association between psychiatric stress and premature onset of disease and early death (Gradus et al., 2015; Scott et al., 2013; Trivedi et al., 2020). Most studies evaluating stress-related advanced epigenetic age have generated estimates of DNA methylation (DNAm) age from blood samples. Only a few have examined associations between psychopathology and advanced epigenetic age in brain tissue (Han et al., 2018; Wolf et al., 2021), which is critical for understanding how psychiatric conditions relate to neurobiological aging.

One possible consequence of changes in DNA methylation at a given locus in the genome is alterations in gene expression in the corresponding region. However, little is known about how advanced epigenetic age, which is based on a weighted summary index representing DNAm loci from across the epigenome, affects gene expression. To date, only one study has addressed this question: Levine et al. (2018) examined an epigenetic index developed to predict disease and death (rather than to predict chronological age) in association with transcriptome-wide gene expression in peripheral monocytes. The age-adjusted index showed similar patterns of association with gene expression as did chronological age. Upregulated genes were enriched for pro-inflammatory signaling pathways and downregulated genes were enriched for those involved in transcription and DNA damage and repair. How advanced epigenetic age might alter gene transcription in the brain, and how this relationship might be affected by psychopathology is unknown.

The primary aim of this study was to examine the association between epigenetic age (relative to chronological age) and gene expression in three cortical regions using postmortem tissue--dorsolateral prefrontal cortex (dlPFC), ventromedial prefrontal cortex (vmPFC), and motor cortex--and to then further evaluate how psychiatric diagnoses of PTSD and AUD, which have previously shown associations with advanced epigenetic age (Wolf et al., 2018a, Wolf et al., 2018b; 2019), altered the association between epigenetic age estimates and gene expression. These brain regions were selected on the basis of evidence of their relevance to PTSD and stress-related disorders (Fonzo et al., 2017; Hayes et al., 2012) and/or aging (Salat et al., 2004, 2005). To operationalize epigenetic age relative to chronological age we used a commonly employed metric in which epigenetic age estimates per the Horvath (2013) algorithm were regressed on chronological age and the residuals from the equation (“DNAm age residuals”) were saved to form a dimensional index ranging from negative values (underestimated epigenetic age relative to chronological age, i.e., slowed cellular age) to positive values (overestimated epigenetic age relative to chronological age, i.e., advanced cellular age). We examined this algorithm as it is the only DNAm age algorithm developed as a multi-tissue predictor, including validation in brain tissue. We conducted both unbiased and hypothesis-driven analyses. In unbiased analyses we examined the main effects of DNAm age residuals and their interaction with PTSD and AUD on transcriptome-wide expression and tested for enrichment of associated biological networks. Our hypothesis-driven analyses were motivated by previous studies that have shown aging to be associated with increased inflammation (Frasca and Blomberg, 2016), oxidative stress (Finkel and Holbrook, 2000), and stress responding (e.g., glucocorticoid responding; Sapolsky et al., 1987), as well as decreased immune efficiency (Pawelec, 2006), and alterations in circadian rhythms (Lananna and Musiek, 2020). Similarly, advanced DNAm age in blood has been associated with biomarkers of glucocorticoid responding (Jovanovic et al., 2017; Zannas et al., 2015), sleep disruption (Carroll et al., 2017), and inflammation, immune, and metabolic dysregulation (Quach et al., 2017), providing further support for the potential effects of advanced DNAm age on gene expression in these biological systems. We hypothesized that expression and network model results would reveal that advanced DNAm age, and advanced cellular aging X PTSD/AUD, would be associated with differential expression of genes relating to these biological processes and systems. Of note, the main effects of psychiatric disease on gene expression in this brain bank have been reported elsewhere (Logue et al., 2021) and thus are not reiterated in this study.

Method

Participants and procedure

We obtained post-mortem left hemisphere brain tissue from 117 donors to the VA National PTSD Brain Bank. The tissue and accompanying clinical data was originally acquired from the Lieber Institute for Brain Development at Johns Hopkins University (Mighdoll et al., 2018). The brain bank included PTSD cases, depressed cases, and age-matched controls; as detailed below, psychiatric comorbidity among cases was assessed and analyzed, rather than treated as an exclusionary criterion. Of the sample of N = 117, a total of 97 had available RNA sequence data, Horvath DNAm age estimates, and genome-wide genotype data for summarizing ancestry that passed all quality control metrics and were the focus of this investigation (see Fig. S1 for sample size flow chart). Sample sizes differed slightly by brain region (Table 1, Fig. S1) as a result of RNA quality control metrics for each region. As shown in Table 1, the full sample comprised 55 men (56.70%) and the mean age at death was 42.38 years (SD: 11.08). Race is listed in Table 1. The full sample included 42 PTSD cases (43.30%), 30 AUD cases (30.93%) and 24 controls (24.74%) who did not meet criteria for any mental health diagnosis. These diagnostic groups were not mutually exclusive (Table 1; Fig. S2), consistent with the high rates of psychiatric comorbidity in the broader PTSD population (Kessler et al., 1995). Each donor was coded as 0 or 1 on each psychiatric diagnosis. We did not evaluate depression independently from PTSD or AUD as over 90% of the PTSD and AUD groups also met criteria for depression (but we did address the added contribution of depression and other comorbidity in follow-up sensitivity analyses; see below).

Table 1

Sample Demographics.

Characteristic	Total	Control	PTSD	AUD
Sample size; n (%)	97	24 (24.74)	42 (43.30)	30 (30.93)
Age at Death; Mean (SD)	42.38 (11.08)	46.53 (9.97)	40.93 (11.34)	43.44 (10.28)
Sex; n (%)
Female	42 (43.30)	8 (33.33)	23 (54.76)	11 (36.67)
Male	55 (56.70)	16 (66.67)	19 (45.24)	19 (63.33)
Race; n (%)
AA	20 (20.62)	8 (33.33)	8 (19.05)	4 (13.33)
CAUC	76 (78.35)	15 (62.50)	34 (80.95)	26 (86.67)
Multi-Racial	1 (1.03)	1 (4.17)	0 (−)	0 (−)
PMI, hours; Mean (SD)	28.55 (8.14)	29.56 (7.02)	28.89 (9.11)	27.67 (7.70)
Horvath DNAm Age Residuals†; Mean (SD)
dlPFC (n = 96)	0.43 (4.37)	0.39 (5.45)	1.27 (3.88)	0.30 (3.46)
vmPFC (n = 89)	0.65 (3.89)	1.74 (4.07)	0.66 (3.93)	0.48 (3.27)
Motor cortex (n = 92)	0.15 (4.47)	0.62 (4.84)	0.12 (4.35)	1.18 (3.64)
Manner of death; n (%)
Suicide	13 (13.40)	0 (−)	9 (2.14)	7 (23.33)
Alcohol or drug use	47 (48.45)	0 (−)	28 (66.67)	17 (56.67)
Comorbidity; n (%)
Comorbid MDD	NA	0 (−)	38 (90.48)	28 (93.33)
Comorbid PTSD	NA	0 (−)	NA	14 (46.67)
Comorbid AUD	NA	0 (−)	14 (33.33)	NA
Smoking Status; n (%)
Positive	61 (62.89)	7 (29.17)	32 (76.19)	24 (80.00)

Note. SD = standard deviation; AA = African-American; CAUC = Caucasian; PMI = post-mortem interval; DNAm = DNA methylation; dlPFC = dorsolateral prefrontal cortex; PTSD = post-traumatic stress disorder; AUD = Alcohol use disorder.

†The Horvath DNAm age residuals were computed on the methylation data from the sample with DNAm data (n = 116 for dlPFC and vmPFC and = 114 for motor cortex), which is larger than the subset with RNA data, thus the mean residuals for the subset with RNA data is not 0. T-tests revealed that mean DNAm age residuals did not differ significantly between the control versus the PTSD or AUD groups in any brain region. RNA data that passed quality control metrics differed slightly by region, as indicated in the table and in Fig. S1. As a result, the female/male breakdown by brain region was as follows: dlPFC = 41 female/55 male; vmPFC = 39 female/50 male; and motor = 39 female/53 male.

Psychiatric diagnoses were based on next-of-kin interviews and review of medical records; the interviews included the MINI International Neuropsychiatric Interview 6.0, the PTSD checklist for DSM-5 (adapted for postmortem studies), and the Lieber Psychological Autopsy Interview, which was conducted by mental health clinicians where possible (Mighdoll et al., 2018). PTSD diagnoses were reviewed by at least two board-certified psychiatrists who made confidence ratings in the diagnosis on a 1–5 scale (a 3 or better was required to be considered a PTSD case; Morrison et al., 2019). AUD diagnoses were determined via medical record review in conjunction with next-of-kin interviews; diagnoses were reviewed by two board-certified psychiatrists to independently confirm each diagnostic determination. Recent (time-of-death) smoking status was reported on during the next-of-kin interview. The neuropathology and toxicology evaluations were conducted by board-certified neuropathologists and cause and manner of death were determined by Maryland state medical examiners. Samples were excluded if there was evidence of history of severe traumatic brain injury, neuritic pathology, or neurodegenerative disease (Morrison et al., 2019). Left dlPFC or Brodmann Area (BA) 9/46, and left vmPFC (BA 12/32) were acquired from the level of the genu of the corpus callosum, and left motor cortex (BA 4) from the level of the superior central sulcus.

DNA and DNAm

DNA extraction was conducted from all three regions (see Supplementary Materials for details), with genotypes determined from motor cortex samples. Genotypes were assayed on the Illumina HumanOmni2.5-8 array and DNAm from each region was assessed with the Illumina Infinium MethylationEPIC array. We calculated the Horvath (2013) DNAm age estimates, a multi-tissue age predictor that included brain tissue in the development of the algorithm, using the 335 probes on the EPIC array that overlap those in the algorithm which was originally developed for the Illumina 450 K BeadChip. Horvath DNAm age estimates correlated with age at death in each region at r = 0.91 to 0.93. As reported in Wolf et al. (2021), DNAm age estimates across brain regions were highly correlated (r = 0.92 to 0.93, ps < 0.001). Chronological age was regressed out from each DNAm age estimate in the sample of n = 116 with DNAm data and the unstandardized residuals were saved (“DNAm age residuals”) to index slowed to advanced epigenetic age relative to chronological age. As per Wolf et al. (2021), DNAm age residuals were moderately correlated with each other across brain regions (r = 0.50 to 0.51, ps < 0.001). Ancestral variation was estimated via principal components (PC) analysis of 100,000 common polymorphisms, with the first three PCs retained as covariates in analyses.

RNA extraction and sequencing

RNA from each of the three brain regions was extracted from 25 mg of tissue using Qiagen RNeasy Fibrous Tissue Minikit. We obtained RNA integrity (RIN) values for a subset of the data, which were found to be acceptable to proceed with library preparation. Illumina TruSeq Stranded total RNA kit with globin depletion was used. A Hiseq 2500 which produced paired-end 75bp reads was used for library sequencing. The Hiseq was performed in two different manners in order to avoid empty lanes: the “high output” mode (flow cells run over eight lanes that contain unique library pools) and “rapid” mode (single cell over two lanes). Trimmomatic (Bolger et al., 2014) was used to eliminate adapters and remove short or low-quality reads in conjunction with aligning the results with the hg38 human reference genome via STAR (Dobin et al., 2013) and Kallisto (Bray et al., 2016a, 2016b) for transcriptome quantification. To further evaluate quality control, the aligned reads were also examined using FastQC, RSeQC, and MultiQC (Ewels et al., 2016). Data were eliminated if there was evidence of less than 50% uniquely mapped reads. Samples were also eliminated if there was evidence that they were outliers in the PCs of regularized log transformed (rLog) expression values. This was accomplished by collapsing Kallisto transcript abundance estimates to the gene level via tximport Bioconductor package and rLog values as estimated in DEseq2 (Love et al., 2014). The PCs were estimated from the rLog values and a threshold of 6 SDs from the mean on the first 10 PCs was used as the cut-point for outlier identification and removal.

For additional quality control metrics, we confirmed expression of X and Y chromosome genes (Shi et al., 2016) against self-reported sex and compared genotypes from the sequenced RNA data (identified via the GATK HaplotypeCaller) against genotype calls from the Illumina HumanOmni 2.5–8 beadchips (i.e., to ensure the data were all correctly aligned to each other). Cell type proportion estimates (weights) were generated using BrainInABlender (Hagenauer et al., 2018) and included as covariates.

Data analyses

We conducted unbiased, transcriptome-wide linear regression analyses in the full sample. Genes with more than one read count in at least 30 subjects were included in analyses. The first set of analyses evaluated the main effects of DNAm age residuals, as predictors of gene expression across the transcriptome, controlling for age, sex, postmortem interval (PMI), cell type estimates (astrocytes, endothelial cells, microglia, mural, neurons, oligodendrocytes and red blood cells), top 3 ancestry PCs, top 3 quality surrogate variables (qSVs), and sequencing-run ID (see Supplementary Materials for details). These models were evaluated for each brain region separately. These analyses were then repeated in sex-stratified cohorts, limiting the transcripts evaluated to those shown to be expressed in each sex (genes with more than one read count in at least 16 male subjects and 14 female subjects were included in corresponding sex-specific analyses, see Supplementary Materials). After examining the main effects of DNAm age residuals, we added each diagnosis (in separate models) and their interaction with DNAm age residuals to the models to examine differential associations between DNAm age residuals and gene expression as a function of psychopathology (i.e., moderation). Analyses were conducted using DESeq2 package in R (Love et al., 2014). All subjects with available data were included in each analysis.

We followed a two-step approach to multiple-testing correction. First, we used an FDR-corrected threshold in each individual analysis (i.e., for a given brain region and cohort) to determine which genes were differentially expressed after correction across the transcriptome. Second, we generated experiment-wide corrected p-values for the main effects of DNAm age residuals which were adjusted across brain region and cohort (i.e., 3 regions X 3 cohorts = 9 sets of analyses of all expressed transcripts) using FDR. Similarly, for the interaction effects, we corrected across all the tests (i.e., 3 regions X 2 diagnoses X 3 cohorts = 18 sets of analyses for all transcripts). Only results that were FDR significant at the experiment-wide level were considered further. To address concerns about the possible confounding effects of psychiatric comorbidity (e.g., the effects of AUD and depression in PTSD analyses), we conducted secondary follow-up linear regressions for genes with multiple-testing corrected differential expression as a function of PTSD, AUD, or their interaction with DNAm age residuals. For example, for all genes that evidenced a corrected significant effect for DNAm age residuals X PTSD, we added the AUD variable and DNAm age residuals X AUD to the model to determine if doing so altered the significance of the DNAm age residuals X PTSD effect. The same approach was followed to examine the potential confounding effects of major depression. Additional follow-up analyses further controlled for smoking status and its interaction with DNAm age residuals.

After the transcriptome-wide analyses, we examined the significance of the main effects of DNAm age residuals and their interactive effects with psychopathology in the five classes of genes hypothesized to be associated with advanced epigenetic age: inflammation, immune, oxidative stress, glucocorticoid, and circadian gene classes. The results in each gene class were extracted from the transcriptome-wide analyses and corrected across the number of genes in each gene class using an FDR correction, rather than across the transcriptome. The genes which were FDR corrected significant within their class were considered significant candidate genes. The genes that were evaluated in each class are listed in Table S1 and were selected in concert with curated gene lists developed by ThermoFisher Scientific (Applied Biosystems TaqMan Gene Expression assays).

We next conducted a gene ontology (GO) overrepresentation analysis of the 200 top differentially expressed genes (in association with the main effect of DNAm age residuals and in association with each interaction term) from each model using GOseq (Young et al., 2010). This allowed us to test for enrichment of particular biological pathways in association with DNAm age residuals and their interaction with PTSD or AUD. We also examined the top 200 up- and down-regulated genes from each model. P-value thresholds were corrected for the number of GO terms examined.

Finally, we conducted a weighted gene co-expression network analysis (WGCNA) to examine gene networks (in the full and sex-stratified samples, using gene lists specific to each sex and brain region). We examined the associations between these gene expression network modules (represented as the first PC of the corresponding module, reflecting the co-expression of each network module) and the main effects of DNAm age residuals and the interactive effects of DNAm age residuals X PTSD or AUD. Gene networks were constructed from the gene correlation matrix computed from regularized log transformed expression values of a set of highly expressed genes which were selected based on the genes with at least 10 read counts in 30 subjects for the full cohort, and the genes with at least 10 read counts in 16 male subjects and 14 female subjects for corresponding cohorts. Co-expression network modules were labeled arbitrarily with a color and were regressed on our main and interaction terms. To interpret each network module, we entered the genes in the co-expression network into DAVID (Sherman and Lempicki, 2009) to examine overrepresentation compared to a background gene set which includes all genes highly expressed enough to be included in the WGCNA analysis. The p-values for DAVID were corrected for the number of pathways examined using the Benjamini method. In addition, we examined if the associated co-expression networks (PCs) were enriched for particular cell type markers with a hypergeometric test using the expressed cell-type marker genes as the background gene set. P-values were corrected for the number of cell type markers tested within each module using FDR.

Results

Transcriptome-wide analyses

Q-Q plots for the main effects of DNAm age residuals did not yield evidence of p-value inflation (Table S2). There were no experiment-wide corrected significant main effects of DNAm age residuals on expression across the cohorts or regions. The addition of the interaction terms (DNAm age residuals X PTSD or DNAm age residuals X AUD), yielded 11 experiment-wide significant interaction effects on expression (Table 2; smallest p = 9.372E-12, smallest pcor = 5.795E-06). Four of these effects were in interaction with PTSD (SNORA73B, COL6A3, GCNT1, and GPRIN3) and seven were in interaction with AUD (ADGRG6, IL1B, NUTM2A-AS1, CES3, ADAMTS18, LINC00643, RCOR2). Effects were distributed across the three regions and across the full cohort and sex-stratified models. Interaction plots for five genes (GCNT1, GPRIN3, IL1B, CES3, RCOR2) of particular interest due to their associations with age-related processes (Franceschi et al., 2000; Alvarez-López et al., 2014; Nolz and Harty, 2014; Karadurmus et al., 2019; Dominguez et al., 2014) are shown in Fig. 1 (see Fig. S3 for the remaining interaction plots). As shown in Fig. 1, advanced epigenetic age (i.e., increasing DNAm age residuals) was associated with increased expression of GCNT1 (in the full cohort) and GPRIN3 (in the men) in dlPFC among those with PTSD. IL1B motor cortex expression among the women was also positively associated with DNAm age residuals among those with (but not without) AUD. In contrast, motor cortex expression of CES3 (in the full cohort) and RCOR2 (in the women) was negatively related to DNAm age residuals as a function of AUD. Fig. 2 shows the results for the same five genes across all analyses, not limited to the experiment-wide region, cohort, or diagnosis in order to evaluate the consistency of effect across models. Of note, while expression of IL1B was upregulated in motor cortex as a function of DNAm age residuals X AUD among the women in the experiment-wide results, the same gene was nominally significantly downregulated among women in vmPFC as a function of DNAm age residuals X AUD and was also nominally downregulated in dlPFC, vmPFC, and motor cortex as a function of DNAm age residuals X PTSD (Fig. 2). All experiment-wide corrected significant effects remained nominally significant in follow-up analyses which additionally controlled for the main effects of various comorbid psychiatric diagnoses (PTSD, AUD, depression, bipolar) and their interaction with DNAm age residuals (Table S3), suggesting that comorbid psychiatric diagnoses did not confound the reported effects. Likewise, experiment-wide effects remained significant with the addition of smoking status and smoking X DNAm age residuals in the model (Table S3). Within-analysis (as opposed to across all analyses) corrected significant effects are listed in Table S4. Overlapping nominally significant main and interaction effects across brain regions are shown in Fig. S4. There were a small number of at least nominally significant effects that overlapped all three brain regions in each model. Notably, IL1B was at least nominally associated with the PTSD X DNAm age residuals interaction term across all three brain regions.

Table 2

Experiment-wide Significantly Expressed Genes in Association with DNAm age Residuals X Diagnosis Interaction Terms.

Gene	Type	Region	Cohort	Moderator	Log2FC	p	pcor-experiment
SNORA73B	snoRNA	vmPFC	Female	PTSD	0.365	9.372E-12	5.795E-06
COL6A3	Protein coding	dlPFC	Male	PTSD	−0.280	1.224E-10	3.782E-05
ADGRG6	Protein coding	dlPFC	Male	AUD	−0.306	6.815E-10	1.405E-04
IL1B	Protein coding	motor cortex	Female	AUD	0.677	1.495E-08	0.002
NUTM2A-AS1	lncRNA	dlPFC	Male	AUD	0.304	1.008E-07	0.012
GCNT1	Protein coding	dlPFC	Full	PTSD	0.106	1.314E-07	0.014
GPRIN3	Protein coding	dlPFC	Male	PTSD	0.114	1.661E-07	0.015
CES3	Protein coding	motor cortex	Full	AUD	−0.150	3.661E-07	0.027
ADAMTS18	Protein coding	vmPFC	Female	AUD	−1.495	3.901E-07	0.027
LINC00643	lncRNA	dlPFC	Male	AUD	−0.188	5.739E-07	0.035
RCOR2	Protein coding	motor cortex	Female	AUD	−0.145	7.323E-07	0.041

Note. Log2FC = log2 fold change; pcor-experiment = experiment-wide FDR adjusted p-value; dlPFC = dorsolateral prefrontal cortex; vmPFC = ventromedial prefrontal cortex; PTSD = posttraumatic stress disorder; AUD = alcohol use disorder; snoRNA = small nucleolar RNA; lncRNA = long non-coding RNA.

Fig. 1

The figure shows the association between DNAm age residuals (X-axis) and gene expression (Y-axis, as regularized log-transformed expression values) as a function of PTSD or AUD diagnosis (the moderator) for five of the eleven genes that achieved experiment-wide significance across all transcriptome-wide analyses (across brain regions, diagnoses, and cohorts). Interaction plots for the remaining differentially expressed genes are shown in Fig. S3. PTSD = posttraumatic stress disorder; AUD = alcohol-use disorder.

Fig. 2

The figure shows the results for the same five experiment-wide significant genes of interest across all analyses, not limited to the experiment-wide region, cohort, or diagnosis. This shows the pattern of results across models for these genes. The cohort is listed along the X-axis and corresponding log2 fold change in expression along the Y axis. PTSD = posttraumatic stress disorder; AUD = alcohol-use disorder. *nominal significance (p < .05). **experiment-wide significance (pcor-experiment < .05).

Candidate gene classes

We next examined the main effects of DNAm age residuals from each region on expression of the genes in the ThermoFisher Scientific curated gene lists for each hypothesized class. There was one corrected significant main effect of DNAm age residuals from the candidate class analyses: expression of AQP1 was significantly associated with DNAm age residuals in motor cortex in the full cohort after correction for 80 genes in the glucocorticoid candidate class (p = 1.722E-04, pcor = 0.014; Table 3, Part A). There were 17 additional corrected significant interaction effects of DNAm age residuals X PTSD or AUD on expression of genes in the candidate class lists (Table 3, Part B; several genes were significant in more than one cohort or region). Specifically, four genes from the inflammation class (IL1B, IGF1, VCAM1, and IL-6), three genes from the glucocorticoid class (VLDLR, BMPER, and ZFP36), three genes from the circadian class (CARTPT, HTR7, PRKACB), three from the oxidative stress class (PTGS2, ALB, PDLIM1), and one from the immune class (CCL19) evidenced differential associations between DNAm age residuals and expression in the presence of PTSD and/or AUD.

Table 3

Significantly Expressed Genes from Candidate Class Analyses.

A. Main Effects of DNAm Age Residuals
Gene	Region	Cohort	Log2FC	p	pcor-class	Gene class	N
AQP1	Motor cortex	Full	-0.107	1.722E-04	0.014	Glucocorticoid	80

Note. Log2FC = log2 fold change; pcor-class = gene class-wide FDR adjusted p-value; N = number of genes included in each gene class; dlPFC = dorsolateral prefrontal cortex; vmPFC = ventromedial prefrontal cortex; PTSD = posttraumatic stress disorder; AUD = alcohol use disorder.

Gene ontology (GO) term overrepresentation

Table 4 shows the top 5 GO terms associated with the main effects of DNAm age residuals and the interaction terms (for each region and cohort; the genes associated with each GO term are listed in Tables S5A–S5C for the top 200 differentially expressed genes and the top 200 up- and down-regulated genes). Each ID in Table 4 was at least nominally significantly associated with each main effect and the interaction terms. Immune-system pathways were among the top enriched for the main effects of DNAm age residuals in the dlPFC and vmPCF among the full cohort, but the results did not withstand correction for multiple testing (smallest pcor = .21). In the female dlPFC model, there was significant enrichment of protein-folding related pathways (smallest pcor = 1.667E-04) in association with the main effects of DNAm age residuals. In the male vmPFC model, there was significant enrichment of cellular development pathways associated with DNAm age residuals (smallest pcor = 7.978E-05). With respect to associations with the PTSD interaction terms, there was significant enrichment of pathways involved in organ development in the dlPFC male model (smallest pcor = 4.601E-06). Immune-response pathways were nominally associated with the PTSD interaction term in the motor cortex in the male cohort (smallest pcor = 0.095) and these associations were driven by the top 200 downregulated genes (smallest pcor = 0.001; Supplementary Table S5C). With respect to associations with AUD interaction terms, there was significant enrichment that was localized to membrane structures in the dlPFC male model (smallest pcor = 3.611E-05) and in inflammation pathways in the female motor model (smallest pcor = .033). This latter association was also likely driven by the top 200 downregulated genes in the female motor model as inflammation pathways were significantly associated with the AUD interaction term in the downregulated genes analysis (smallest pcor = 2.071E-05; Supplementary Table S5C). Also of note, the main effects of DNAm age residuals and the AUD interaction term were associated with enrichment of the adenylate cyclase-activating adrenergic receptor signaling pathway (i.e., stress responding) in the full sample motor cortex and female vmPFC models, though neither term yielded a corrected-significant p-value.

Table 4

Gene Ontology (GO) Overrepresentation Analysis of the top 200 Differentially Expressed Genes.

	Main effect of DNAm age residuals		PTSD interaction		AUD interaction
	GO term	pcor	GO term	pcor	GO term	pcor
dlPFC Full Cohort	pos reg of lymphocyte mediated immunity	0.213	metencephalon dev	0.229	digestion	1
	pos reg of adaptive immune response based on somatic recombination of immune receptors built from immunoglobulin superfamily domains	0.213	sprouting angiogenesis	0.342	lipid digestion	1
	pos reg of adaptive immune response	0.213	cerebellum dev	0.345	osteoblast dev	1
	reg of lymphocyte mediated immunity	0.213	retinol metabolic process	0.345	flagellated sperm motility	1
	reg of adaptive immune response based on somatic recombination of immune receptors built from immunoglobulin superfamily domains	0.213	dendrite dev	0.345	sperm motility	1
dlPFC Male Cohort	single-stranded DNA binding	1	multicellular organism dev	4.601E-06	intrinsic component of membrane	3.611E-05
	DNA repair	1	anatomical structure dev	4.601E-06	integral component of membrane	3.611E-05
	heterocycle catabolic process	1	system dev	6.273E-06	membrane part	3.611E-05
	DNA ligation	1	animal organ dev	1.173E-05	cell projection	3.731E-04
	pos reg of sterol transport	1	dev process	1.241E-05	nervous system dev	3.731E-04
dlPFC Female Cohort	response to topologically incorrect protein	1.667E-04	cellular response to endogenous stimulus	0.213	mitogen-activated protein kinase kinase binding	0.645
	response to unfolded protein	0.002	cell-substrate adhesion	0.213	I band	0.645
	cellular response to topologically incorrect protein	0.013	sprouting angiogenesis	0.213	primary lysosome	0.645
	neg reg of RNA metabolic process	0.040	cell-matrix adhesion	0.267	azurophil granule	0.645
	neg reg of transcription by RNA polymerase II	0.040	response to endogenous stimulus	0.505	protein kinase binding	0.645
vmPFC Full Cohort	innate immune response in mucosa	0.846	purinergic receptor signaling pathway	1	receptor ligand activity	1
	cellular response to glucose stimulus	0.846	sperm midpiece	1	intramolecular oxidoreductase activity	1
	cellular response to hexose stimulus	0.846	sperm flagellum	1	phosphatidylinositol phospholipase C activity	1
	cellular response to monosaccharide stimulus	0.846	purinergic receptor activity	1	response to retinoic acid	1
	response to glucose	0.846	Hsp70 protein binding	1	receptor regulator activity	1
vmPFC Male Cohort	dev process	7.978E-05	endothelium dev	0.439	reg of transport	0.151
	anatomical structure dev	1.574E-04	endothelial cell differentiation	0.439	reg of secretion	0.151
	multicellular organism dev	0.001	lipopolysaccharide-mediated signaling pathway	0.611	reg of signaling receptor activity	0.151
	cellular dev process	0.001	signaling pattern recognition receptor activity	0.611	receptor regulator activity	0.151
	system dev	0.004	determination of heart left/right asymmetry	0.611	pos reg of secretion	0.151
vmPFC Female Cohort	autophagosome	1	RNA processing	0.113	system process	0.280
	anion:anion antiporter activity	1	nucleolus	0.113	anterior/posterior axis specification	0.280
	meiotic cell cycle	1	membrane-enclosed lumen	0.113	monovalent inorganic cation transmembrane transporter activity	0.298
	reg of meiotic nuclear division	1	organelle lumen	0.113	structural molecule activity conferring elasticity	0.298
	C4-dicarboxylate transport	1	intracellular organelle lumen	0.113	adenylate cyclase-activating adrenergic receptor signaling pathway	0.298
Motor Full Cohort	intrinsic component of plasma membrane	0.439	response to interleukin-1	0.186	cytoplasmic part	0.195
	anchored component of membrane	0.439	DNA damage response, detection of DNA damage	0.186	reg of I-kappaB kinase/NF-kappaB signaling	0.195
	sodium channel regulator activity	0.439	cellular response to interleukin-1	0.186	syntaxin-1 binding	0.211
	adenylate cyclase-activating adrenergic receptor signaling pathway	0.439	translesion synthesis	0.237	I-kappaB kinase/NF-kappaB signaling	0.217
	integral component of plasma membrane	0.439	nucleotide-excision repair, DNA gap filling	0.237	pos reg of I-kappaB kinase/NF-kappaB signaling	0.217
Motor Male Cohort	protein targeting to membrane	0.136	immune response	0.095	reg of epidermis dev	0.100
	cotranslational protein targeting to membrane	0.136	exocytosis	0.095	reg of epidermal cell differentiation	0.100
	nuclear-transcribed mRNA catabolic process, nonsense-mediated decay	0.216	endocytic vesicle	0.095	pos reg of epidermis dev	0.142
	nuclear-transcribed mRNA catabolic process	0.216	myeloid leukocyte mediated immunity	0.095	neg reg of cell proliferation	0.147
	protein localization to endoplasmic reticulum	0.216	immune system process	0.099	neg reg of cellular component movement	0.147
Motor Female Cohort	carbohydrate metabolic process	0.012	transcytosis	1	response to interferon-gamma	0.033
	reg of synapse organization	0.179	apoptotic cell clearance	1	cytokine-mediated signaling pathway	0.033
	reg of synapse structure or activity	0.179	macromolecule transmembrane transporter activity	1	cellular response to interferon-gamma	0.034
	protein kinase C signaling	0.563	pos reg of endothelial cell apoptotic process	1	cellular response to cytokine stimulus	0.034
	reg of carbohydrate metabolic process	0.563	response to oxygen-containing compound	1	reg of epithelial cell apoptotic process	0.034

Note. The top 5 GO terms for each analysis are listed along with the FDR corrected p-values. All of the top 5 GO terms were nominally significant at the p < .05 level. Associated GO term IDs are listed in Tables S5a–S5c. Pos = positive; neg = negative; dev = development(al), reg = regulation.

Association with cell type markers

We examined our main and interaction terms in association with the cell type markers in each region and cohort. We found that in the motor cortex, the PTSD X DNAm age residuals interaction term evidenced a corrected significant positive association with oligodendrocytes (β = .056, p = .003, pcor = .034; Table S6). No other effects were significantly associated with cell type markers after correction for multiple testing.

Weighted gene co-expression network analysis

WGCNA analyses identified networks in each brain region. Several of these networks were significantly associated with the main effects of DNAm age residuals and/or the interaction terms (Table 5). Enrichment analyses were used to indicate the function and location of the expression networks (Table 5, Table S7). Several co-expression networks were localized to basic cellular components, such as membrane structure and cellular junctions. In addition, two networks were enriched for myelination pathways and showed corresponding cell type marker enrichment for oligodendrocytes (Table S8). Specifically, the motor cortex pink network was associated with the DNAm age residuals X AUD interaction and the vmPFC turquoise network (generated in the full cohort) was associated with DNAm age residuals X PTSD among women in sex-stratified analyses; both showed enrichment for the myelination biological process. Gene membership in the networks is listed in Table S9.

Table 5

Network Modules Significantly Associated with DNAm Age Residuals and Interaction Terms.

Module	Description	Cell type marker	Na	Main effect or Moderator & Cohort	Beta	SE	p	n (%)b
dlPFC-Full-orangered4	Membrane/cell junction/dendrite	Ex	45	PTSD/	−0.009/	0.003/	0.001/	20 (44)/
				PTSD♂︎	−0.009	0.003	0.010	8 (18)♂︎
dlPFC-Full-royalblue	Membrane/chemical synaptic transmission/cell junction	Ex/In	96	AUD♂︎/	−0.018/	0.005/	0.002/	63 (66)♂︎/
				PTSD♂︎/	−0.013/	0.004/	0.004/	48 (50)♂︎/
				PTSD	−0.007	0.003	0.010	28 (29)
dlPFC-Male-darkseagreen4	Post-synaptic & plasma membrane	Ex/In	109	AUD♂︎/	−0.026/	0.009/	0.004/	61 (56)♂︎/
				PTSD♂︎	−0.020	0.007	0.006	47 (43)♂︎
motor-Full-lightcyan1	cell junction/membrane/cardiac conduction	Ex	210	AUD/	−0.008/	0.003/	0.005/	77 (37)/
				Main effect♂︎	−0.004	0.002	0.027	49 (23)♂︎
dlPFC-Full-darkorange2	Calcium ion binding/membrane	Ex	83	PTSD♂︎/	0.016/	0.005/	0.005/	47 (57)♂︎/
				PTSD	0.009	0.003	0.006	47 (57)
dlPFC-Male-coral1	Membrane/presynaptic active zone	Ex	108	PTSD♂︎	0.021	0.007	0.006	48 (44)♂︎
dlPFC-Male-darkred	Plasma membrane	In	93	AUD♂︎	−0.025	0.009	0.007	38 (41)♂︎
dlPFC-Full-violet	Nucleus/transcription	NA	58	PTSD	0.010	0.004	0.008	16 (28)
dlPFC-Male-lightcyan1	Ion transport	NA	103	AUD♂︎	−0.015	0.005	0.011	33 (32)♂︎
motor-Full-thistle2	NA	NA	25	AUD♂︎	−0.003	0.001	0.011	2 (8)♂︎
motor-Full-plum1	NA	Ex	34	Main effect	−0.005	0.002	0.011	18 (53)
vmPFC-Full-yellowgreen	NA	Ex	40	Main effect♂︎/	0.007/	0.003/	0.013/	18 (45)♂︎/
				Main effect	0.005	0.002	0.031	13 (33)
dlPFC-Male-darkgreen	RNA binding	NA	787	AUD♂︎	−0.020	0.008	0.014	235 (30)♂︎
motor-Full-mediumpurple3	NA	NA	31	Main effect	0.005	0.002	0.014	13 (42)/
motor-Full-darkorange2	Membrane	Oligo	114	AUD	0.008	0.003	0.015	42 (37)
motor-Full-violet	Ubiquitin-protein transferase activity	NA	281	AUD	0.010	0.004	0.015	124 (44)
vmPFC-Full-darkred	Translation, RNA processing, ribosomes	Ex	135	Main effect	−0.004	0.002	0.015	49 (36)
motor-Full-pink	Myelination	Oligo	299	AUD	0.002	0.001	0.016	22 (7)
vmPFC-Male-orangered4	NA	Ex	48	Main effect♂︎	0.010	0.004	0.018	22 (46)♂︎
dlPFC-Female-yellow	Transcription/dendrite	Ex	1616	PTSD	0.012	0.005	0.020	242 (15)
dlPFC-Full-pink	Ligand-gated ion channel, plasma membrane, nicotine response	Ex/In	210	PTSD/	−0.008/	0.003/	0.024/	49 (23)/
				PTSD♂︎/	−0.013/	0.005/	0.027/	49 (23)♂︎/
				AUD♂︎	−0.015	0.007	0.037	54 (26)♂︎
dlPFC-Full-cyan	Transcription/zinc ion binding/dendrite	NA	1403	PTSD♀︎	0.008	0.003	0.027	213 (15)♀︎
vmPFC-Full-turquoise	Plasma membrane/myelination/actin filament binding	Oligo	1039	PTSD♀︎	−0.007	0.003	0.027	47 (5)♀︎
dlPFC-Full-darkgreen	NA	NA	91	AUD♂︎/	−0.016/	0.007/	0.030/	28 (31)♂︎/
				Main effect♀︎	0.006	0.003	0.032	18 (20)♀︎
dlPFC-Male-saddlebrown	Membrane/paranode region of axon	Oligo	72	PTSD♂︎	−0.003	0.001	0.031	3 (4)♂︎
motor-Full-plum2	Neurofilament/axon	NA	26	Main effect	0.005	0.002	0.031	13 (50)
motor-Full-black	Nucleic acid binding/transcription	NA	344	AUD	0.010	0.005	0.034	91 (26)
dlPFC-Male-darkgrey	Cis-Golgi network/nucleic acid binding/nucleoplasm	NA	249	AUD♂︎	−0.018	0.008	0.034	51 (20)♂︎
vmPFC-Full-darkorange	NA	Astro	71	Main effect♂︎	0.003	0.002	0.036	16 (23)♂︎
dlPFC-Full-brown4	Protein targeting to membrane/mRNA catabolic process	Ex	354	PTSD♀︎	−0.008	0.003	0.038	23 (7)♀︎
dlPFC-Female-lightgreen	NA	Oligo/OPC	539	PTSD♀︎	0.009	0.004	0.040	51 (9)♀︎
dlPFC-Full-darkred	Mitochondrion/proteasome complex/cytosol	Ex	1353	AUD♂︎	−0.007	0.003	0.040	352 (26)♂︎
vmPFC-Full-greenyellow	NA	Astro/OPC	272	Main effect	0.002	0.001	0.043	37 (14)
vmPFC-Full-saddlebrown	NA	Oligo	62	PTSD♀︎	−0.017	0.008	0.047	7 (11)♀︎
vmPFC-Full-bisque4	NA	NA	29	Main effect	−0.001	0.0003	0.048	2 (7)

Note. Beta = estimated regression coefficient; SE = standard error; dlPFC = dorsolateral prefrontal cortex; vmPFC = ventromedial prefrontal cortex; PTSD = posttraumatic stress disorder; AUD = alcohol use disorder; ♂︎ = male-specific network association analysis; ♀︎ = female-specific network association analysis; Astro = astrocytes; Endo = endothelial cells; Ex = excitatory neurons; In = inhibitory neurons; Oligo = oligodendrocytes; OPC = oligodendrocyte progenitor cells; Per = pericytes.

aN = number of genes included in each network module; bn(%) = number and percent of nominally significant genes from the transcriptome-wide analysis examined in the same region, cohort and with the same independent variable as in corresponding network analysis.

Discussion

Advanced DNAm age is a biomarker for increased risk for premature disease onset, including neuropathology, such as Alzheimer's disease (Levine et al., 2015) and reduced morphological integrity of the brain (Proskovec et al., 2020; Wolf et al., 2019). Despite the direct role of DNAm on gene expression, no study to date has evaluated the expression correlates of slowed or advanced DNAm age in neural tissue. In this study we found that the relationship between DNAm age residuals and the expression of 11 genes was dependent on psychiatric disease status. That is, advanced methylation age and psychopathology exerted synergistic effects on gene expression, with particular enrichment of inflammation-related genes and pathways.

Inflammation and accelerated aging

The results of several analyses converged to suggest differential expression of inflammation genes in association with both the main effects of DNAm age residuals and the interaction between the residuals and PTSD and AUD. Specifically, IL1B, which has been associated with PTSD (Passos et al., 2015), alcohol use (Szabo and Lippai, 2014), and aging (a.k.a. “inflammaging; ” Franceschi et al., 2000), evidenced an experiment-wide significant positive association (i.e., upregulation) with DNAm age residuals X AUD in the motor cortex among women. Its expression was also at least nominally associated with DNAm age residuals X PTSD across all three brain regions in the full cohort, although it was downregulated in those analyses. IL1B is an inflammatory cytokine produced and secreted primarily by microglia and astrocytes. It is considered a “master regulator” of neuroinflammation due to its hierarchical role in signaling the expression of other inflammatory cytokines that are neurotoxic and contribute to neuroinflammation (Basu et al., 2004) and neurodegeneration (Liu and Chan, 2014). Preclinical research suggests that exposure to severe stress can lead to time-dependent increases in IL1B immunoreactivity and mRNA expression within the dentate gyrus of the dorsal hippocampus (Jones et al., 2015). IL1B has been shown to be significantly upregulated in PTSD cases compared to controls (Guardado et al., 2016; Logue et al., 2021). In this study, we found that psychiatric disease seems to amplify the effects of accelerated DNAm age on expression of IL1B.

Another gene which emerged as an experiment-wide significant effect was RCOR2, which was downregulated in motor cortex among women with advanced DNAm age and AUD. RCOR2 (repressor element 1-silencing transcription [REST] factor corepressor-2) suppresses inflammation, including IL-6 (Hanzu et al., 2013). REST silences genes involved in apoptosis and is protective against Alzheimer's-related neurodegeneration and neural oxidative stress (Lu et al., 2014). A mouse study of accelerated aging and Alzheimer's-like neurodegeneration found that RCOR2 expression was downregulated in the cortex and hippocampus of mice exhibiting an accelerated aging phenotype; in astrocytes, this downregulation preceded onset of age-related degenerative phenotypes (Alvarez-López et al., 2014). Downregulation of RCOR2 in this study may signal loss of protection against neuroinflammation, oxidative stress, and neurodegeneration. The gene may be part of an important pathway to pursue in understanding the mechanisms of early disease onset in which both neurotoxic effects of inflammation and loss of neuroprotection simultaneously amplify risk for neurodegeneration.

Experiment-wide significant results also revealed that expression of GCNT1 was increased in dlPFC as a function of DNAm age residuals among men and women with PTSD. GCNT1 encodes an enzyme that is critical for core 2 O-glycan synthesis, a component of lymphocyte trafficking and memory T-cell responses to inflammation in the endothelium (Nolz and Harty, 2014; Perkey et al., 2020). High-endothelial venules (HEV) are tissues that allow for lymphocytes to be recirculated through blood vessels, including through endothelial cells (HECs). The GCNT1 enzyme that is expressed by HECs directs lymphocyte distribution to allow for effective response to infection (Veerman et al., 2019). Both PTSD and aging are associated with decreased T-cell mediated immune responding, including increases in the number of end-stage memory cells relative to naïve cells (Aiello et al., 2016) and altered expression of immune gene networks (Breen et al., 2018; Mehta et al., 2018). There is preclinical evidence that peripheral inflammation, as indexed by endothelial cells, results in increased macrophage and microglial activity in the brain (i.e., neuroinflammation; Wohleb et al., 2014). Given the role of endothelial cells in regulating the blood-brain-barrier, it is possible that altered expression of GCNT1 in neural tissue could be a reflection of increased peripheral inflammation crossing the blood-brain-barrier via endothelial cell trafficking of lymphocytes among those with increased cellular aging and PTSD. GCNT1 may be a novel target for better understanding and treating the pathophysiology of accelerated aging in stressed populations.

Results of candidate gene class analyses converged with those from the unbiased analyses to reveal significant associations with three additional inflammation genes (IGF1, VCAM1, IL6), all of which have previously been associated with psychiatric stress or PTSD in protein expression studies (Passos et al., 2015; Santi et al., 2018; Sumner et al., 2018). In addition, CCL19, which is part of the adaptive immune response, was downregulated with increasing DNAm age residuals in those with AUD in the male dlPFC. Overrepresentation analyses in our study showed nominally significant enrichment of immune system GO terms in association with the main effects of DNAm age residuals in the dlPFC and vmPFC. Similarly, there was corrected significant enrichment of down-regulated genes in immune pathways in association with the PTSD interaction term in the male motor cortex and in association with the AUD interaction term in the female motor cortex. Related to this, our analysis of cell type markers found that DNAm age residuals were positively associated with oligodendrocyte markers in motor cortex among those with PTSD. Oligodendrocytes are primarily responsible for the production of myelin, which is critical for axon conduction and highly sensitive to aging and diseases of aging, including neuroinflammation and Alzheimer's disease (Cai and Xiao, 2016; Nasrabady et al., 2018). Collectively, this pattern of results, across variables, analytic approaches, brain regions, and cohorts strongly suggests that accelerated epigenetic aging in the context of psychopathology yields alterations in the expression of inflammatory genes.

Given the role of inflammation in numerous diseases of aging (Furman et al., 2019) and neurodegeneration (Newcombe et al., 2018), these results suggest that accelerated epigenetic age and psychopathology exert individual and synergistic effects on expression of inflammation genes, which in turn, may serve as a mechanism for early onset of disease and health decline. Future research could evaluate the impact of anti-inflammatory treatments on epigenetic aging, particularly among those with psychopathology. This might include pharmacological approaches to counter inflammation. For example, a large trial of canakinumab, which inhibits IL1B, among individuals who had previously experienced myocardial infarction found that the drug reduced risk of subsequent nonfatal myocardial infarction and stroke relative to placebo (Ridker et al., 2017). Caloric restriction is also known to dramatically slow the aging process (Sinclair, 2005) and has immediate beneficial effects on inflammation and immunity (Calder et al., 2011; Sierra Rojas et al., 2016; Wu et al., 2019). There is preclinical evidence that caloric restriction may slow age-related methylation changes (Maegawa et al., 2017). Likewise, a recent study found that non-steroidal anti-inflammatory (NSAIDs) medications as well as calcium channel blockers were cross-sectionally associated with decreased Horvath age residuals, though NSAIDs appeared to accelerate aging in longitudinal analyses (Kho et al., 2021). Collectively, this suggests that anti-inflammatory medications and anti-inflammatory diets are worthy of further investigation to test their ability to alter epigenetic aging, gene expression, and disease onset, particularly in stressed populations who are already at risk for increased inflammation (Marsland et al., 2017; Michopoulos et al., 2017).

Stress responding, oxidative stress, and circadian effects

Candidate class analyses also revealed effects for stress response, oxidative stress, and circadian genes in association with the main effects of DNAm age residuals and the interaction terms. These three biological systems and processes, in concert with immune and inflammation responses, influence each other, thus these effects are unlikely to be independent of one other. Circadian processes, for example, have known coordinating homeostatic effects on immune, inflammatory, stress response, and oxidative stress processes (Buxton et al., 2012; Irwin et al., 2015; Rijo-Ferreira and Takahashi, 2019; Trivedi et al., 2017). Circadian genes synchronize multiple daily metabolic and neural functions via regularized and dynamic gene expression (Mazzoccoli et al., 2012). The automated sequenced patterns of gene transcription that coordinate daily biological processes can be disrupted with advancing age (Terzibasi-Tozzini et al., 2017), leading to reduced neurogenesis (Malik et al., 2015). Poor sleep, a common transdiagnostic feature of psychopathology, alters the rate at which the circadian clock oscillates (Wells et al., 2017), and changes the expression of core clock genes (Cedernaes et al., 2015). This has downstream effects on metabolic and neurodegenerative health outcomes (Rijo-Ferreira and Takahashi, 2019). One possibility is that an increased pace of cellular aging in methylation is reflected in shorter circadian gene expression cycles. Efforts to delay age-related disease onset and the pace of cellular aging via therapeutics that influence chronobiology may need to consider circadian pharmacodynamics to maximize efficacy (Dallmann et al., 2014).

Genes implicated in age-related diseases

The transcriptome-wide analyses also yielded evidence of involvement of two additional genes of interest. Specifically, in the full sample, expression of CES3 was decreased in motor cortex as a function of DNAm age residuals among individuals with AUD. Because of its role in fatty acyl and cholesterol ester metabolism, CES3 is frequently implicated in age-related diseases characterized by a surplus of fatty acids, such as obesity and diabetes (Dominguez et al., 2014), steatohepatitis (Lian et al., 2016; Matsubara et al., 2012), and atherosclerosis (Wang et al., 2012). CES3 is involved in the detoxification of xenobiotics and in drug metabolism (Sanghani et al., 2009); its expression in liver cells has been shown to change depending on level of ethanol exposure (Bardag-Gorce et al., 2006). Its role in accelerated aging could help to explain early onset of metabolic diseases among those with psychiatric stress. Finally, GPRIN3 (G protein-regulated inducer of neurite growth) has been associated with dopamine receptor activation and knock out of this gene in preclinical research is related to increased anxiety (Mototani et al., 2018) and proclivity for substance use (Karadurmus et al., 2019). Its effect on dopaminergic receptors carries downstream implications for disorder of aging, including Parkinson's disease (Karadurmus et al., 2019). These are potentially additional novel contributors to accelerated aging and may hold the key to new treatments to slow the biological aging process.

Study limitations

Results should be interpreted in light of a number of study limitations. First, given the nature of the tissue, this was a small cohort and statistical power was therefore limited. Second, as these are cross-sectional data, we cannot gain leverage on the direction of association or clearly differentiate risks versus consequences of accelerated cellular aging and psychiatric disease. Third, we did not have peripheral biomarkers from these samples so we could not evaluate the consistency of results across the peripheral and central nervous systems. We also did not have access to a second cohort to test for replication of these effects.

Conclusions

This is the first study to evaluate the gene expression correlates of accelerated epigenetic age in brain tissue as a function of PTSD and AUD. Results from unbiased, hypothesis-driven, and enrichment analyses converged on the association between accelerated epigenetic age, alone and interaction with psychopathology, on differential expression of inflammatory and immune system related genes. Effects for IL1B were particularly robust across analytic approaches. Results, in concert with prior research, suggest the importance of evaluating anti-inflammatory interventions in future studies aimed at slowing the pace of cellular aging and increasing resilience to diseases of aging. This could contribute to meaningful extensions in lifespan, healthspan, and functionality, yielding both personal and societal benefits.

Funding

This work was supported by the grant number 1R21AG061367 to EJW, and grant number 1I01BX003477 to MWL. This work was also supported by a Presidential Early Career Award for Scientists and Engineers (PECASE 2013A) to EJW, as administered by U.S. Department of Veterans Affairs Office of Research and Development and by the National Center for PTSD. FGM's and SEH's contribution to this work was supported by award number 5T32MH019836-16. Genotype and methylation data were generated with the support of resources at the Pharmacogenomics Analysis Laboratory (Research Service, Central Arkansas Veterans Healthcare System, Little Rock, Arkansas), a core research laboratory funded by the Cooperative Studies Program, Research and Development, Department of Veterans Affairs. The contents of this manuscript do not represent the views of the U.S. Department of Veterans Affairs, the National Institutes of Health, or the United States Government.

CRediT authorship contribution statement

Erika J. Wolf: Conceptualization, Methodology, Formal analysis, Writing – original draft, Writing – review & editing, Project administration, Funding acquisition. Xiang Zhao: Methodology, Formal analysis, Data curation, Writing – original draft, Writing – review & editing, Visualization. Sage E. Hawn: Writing – original draft, Writing – review & editing. Filomene G. Morrison: Investigation, Writing – review & editing. Zhenwei Zhou: Methodology, Formal analysis, Data curation. Dana Fein-Schaffer: Writing – review & editing. Bertrand Huber: Resources, Data curation, Supervision, Project administration, Writing – review & editing. Mark W. Miller: Resources, Writing – original draft, Writing – review & editing, Conceptualization, Project administration, Funding acquisition. Mark W. Logue: Conceptualization, Resources, Writing – original draft, Writing – review & editing, Supervision, Project administration.

Declaration of competing interest

Dr. Wolf owns stock in Illumina, Inc. All other named authors report no financial or other conflicts of interest in relationship to the contents of this article. Filomene G. Morrison's contribution to this work was completed as a post-doctoral fellow at Boston University School of Medicine and the National Center for PTSD. Dr. Morrison is currently an employee of BlackThorn Therapeutics.