Changes in the expression of inflammatory and epigenetic-modulatory genes after an intensive meditation retreat.

Background: Meditation retreats are characterized by intensive or concentrated periods of meditation practice, commonly undertaken in a residential setting. Although research indicates that meditation training can positively influence physical and mental health outcomes, the biological consequences of meditation retreat interventions are relatively understudied. In this study, we examined the influence of a month-long, silent meditation retreat on the expression of genes involved in epigenetic modulation and immune processes. Method: We assessed gene expression changes in experienced meditators attending a month-long Insight meditation retreat (n = 28), as compared to a community control group (n = 34) of experienced practitioners living their everyday lives. Blood samples were collected on day two of the retreat (Time 1) and again 3 weeks later (Time 2). Control participants were also assessed across a 3-week interval, during which they maintained their regular daily routines.
Results: As compared to controls, retreat participants showed differential changes in the expression of several genes involved in chromatin modulation and inflammation. The most substantive finding was downregulation of the TNF pathway in retreat participants, which was not observed in controls. Conclusions: These findings indicate that meditation retreat participation may influence some of the inflammatory mechanisms involved in the development of chronic diseases, and that this style of psychosocial intervention may have therapeutic potential, particularly in experienced practitioners.

ABSTRACT.

Introduction

Decades of research have demonstrated bidirectional relationships between the immune system and psychological well-being, with important implications for health and longevity [[1], [2], [3]]. Accordingly, there is a growing body of research aimed at understanding whether mind-body interventions, including meditation, can influence biological pathways involved in immunity and inflammation [[4], [5], [6], [7]]. Initial evidence suggests that pro-inflammatory genes, which are upregulated in response to social stressors and early life adversity [8,9], may be downregulated by mind-body interventions [10]. Here, we add to this literature by examining changes in the expression of genes involved in inflammation and epigenetic modulation across 21 days of intensive meditation practice in a residential retreat context.

One of the mechanisms by which psychosocial stressors are thought to negatively impact health is through inappropriate or chronic activation of inflammatory pathways [11,12]. Tumor necrosis factor-α (TNF-α) and nuclear factor-κβ (NF-κβ), in particular, have been linked to a variety of psychosocial stressors and poor health outcomes, including neuroinflammation [13] and depression [12,14]. Interestingly, downregulation of the NF-κβ pathway is one of the more consistent molecular outcomes of mind-body interventions [5,10,15], which are hypothesized to improve health outcomes, in part, by promoting adaptive responses to stress [16,17].

Given initial evidence suggesting that meditation training can alter gene expression, there is also growing interest in understanding how meditation influences the epigenetic mechanisms that regulate gene expression. A study comparing one day of mindfulness practice to a day of leisure activities (i.e., resting, reading and watching documentaries), found rapid reductions in several histone deacetylases (HDACs) genes along with reductions in the proinflammatory genes COX2 and RIPK2. Interestingly, the reductions in HDAC2 and RIPK2 were more pronounced in meditators with better cortisol recovery following a psychosocial stressor [7]. These data are consistent with the social signal transduction theory put forth by Black and colleagues, which proposes that mindfulness training mitigates inflammatory gene expression by altering the interpretation of potentially threatening social stimuli [4].

Meditators in this day-long study also showed changes in the acetylation and methylation of histones [7] and increased methylation of other genes involved in chromatin remodeling [18], suggesting that meditation practice may elicit anti-inflammatory effects through a variety of gene regulatory mechanisms. In addition to the methylation changes observed in this day-long study [18], a cross-sectional study comparing long-term meditators to non-mediators also showed 64 differentially methylated sites (DMS) in the experienced meditators and a central role for TNF-α and NF-κβ-regulated networks [19]. Collectively, bioinformatic analyses from these two studies found differentially methylated sites in pathways involved in cell aging, neurotransmission, lipid and glucose metabolism, immunity and inflammation in meditators vs. controls [18,19].

Despite growing evidence that meditation can impact inflammatory gene expression, the existing literature is sparse and limited in a variety of ways. For example, the majority of work has been conducted by a single research group [15] using DNA microarrays, which lack the precision provided by more targeted assays. Here we compliment this work by using highly specific targeted gene probes to quantify gene expression by real-time qPCR, which is one of the most sensitive and accurate methods, broadly applied to validate microarray data.

Additionally, existing studies have examined gene expression changes during shorter-term or less intensive meditation interventions. In contrast, meditation retreats represent an important, but understudied, mode of training, which are typically used to develop meditative expertise through extended periods of concentrated practice [20]. Retreats are intentionally designed to encourage practitioners to use each moment and every activity as an opportunity for continued mindfulness. To reduce the distractions of daily life, retreats are often held in secluded natural environments where meditators follow a rigorous schedule of formal practice under the guidance of experienced teachers, and with the social support of fellow practitioners. Meditators may also adopt the practice of noble silence, which entails temporarily refraining from speaking, communicating, or initiating eye contact with others, except during meetings with teachers or to navigate common spaces. This practice is intended to facilitate quietude and focus on one's inner experience. Together, these conditions afford a unique opportunity for practitioners to observe their mental experience, to cultivate particular qualities of mind, and to experience meditative insights that may have synergistic effects on their overall well-being. Given the profundity of retreats for developing meditative expertise, these interventions may also represent a valuable research tool for understanding psychobiological effects of meditation training.

Here we investigated a month-long meditation retreat intervention to characterize changes in gene expression that occur when experienced meditators engage in high doses of practice in a retreat context. The first aim of this study was to compare gene expression changes observed in experienced practitioners after one day of intensive practice [7] to those observed after 21 days of intensive practice. The genes HDAC2, HDAC3, HDAC9, RIPK2, COX2 and TNFα were chosen for this purpose. For these genes, we expected to see directionally similar patterns to those observed by Kaliman et al. [7].

Our second aim was to expand on these prior findings by examining changes in other genes involved in inflammatory and epigenetic processes. For this aim, we selected 27 additional candidate genes based on 1) existing literature indicating their mechanistic role in the immune system, NF-κB signaling, and chromatin modification; 2) a preliminary analysis using human chromatin modification enzymes and inflammation RT2 Profiler™ PCR arrays (Qiagen); and 3) a database analysis of physical and functional protein-protein interactions (String). This list of candidate genes included several cytokines, chemoattractants, and regulators of chromatin conformation and gene expression (e.g., histone deacetylases, acetyltransferases, and methyltransferases). Although we did not have specific hypotheses for each of these genes, we did expect to see retreat-related reductions in pro-inflammatory genes along with changes in the expression of genes involved in epigenetic processes that might support the suppression of these genes.

Methods

Data presented in this report were collected as part of a larger investigation of intensive meditation training (ClinicalTrials.gov #NCT03056105). The purpose of the broader study was to evaluate the effects of a month-long residential retreat program on psychological well-being and psychobiological outcomes from an ecologically-informed perspective. Candidate genes were selected by the University of Barcelona team (P.K., M.J.A, and M.C.) after samples were collected by the UC Davis team (C.D.S, Q.A.C, A.P.Z., and B.G.K.). Interested readers can find additional details regarding the retreat environment, meditation training, participant recruitment and demographics in our prior report [21].

Participants

To maximize ecological validity, we chose to evaluate a long-standing and annually offered retreat program at a well-established retreat center. Retreat participants (i.e., retreatants, n = 28) were recruited from a pool of individuals pre-enrolled in one of two month-long residential retreats held at Spirit Rock Meditation Center (Woodacre, CA, USA) in February and March of 2013. Given the length of the intervention and the logistical constraints of coordinating with a large-scale retreat center running established retreat programs, it was infeasible to randomly assign participants to groups. Instead, to control for some of the lifestyle factors that might be expected to accompany long-term practice, comparison participants2 (i.e., controls, n = 34) were recruited from the local Spirit Rock community and assessed between May of 2013 and February of 2014. This control group had similar levels of lifetime meditation experience and prior retreat attendance, and did not significantly differ from the retreat group in age, gender, BMI, education, or income (see Ref. [21]).

Applicants who self-reported major medical conditions expected to impact immune function (i.e., cancers, autoimmune diseases, immunodeficiency disorders, or other conditions involving chronic inflammation such as hepatitis) were excluded from the blood collection procedure. Two retreat and two control participants withdrew from the study after the initial assessment. Two additional controls could not attend the second assessment due to scheduling conflicts. Thus, the final sample consisted of 28 retreatants and 34 controls at the first assessment, and 26 retreatants and 30 controls at the second assessment. All participants gave informed consent before taking part in the study and were compensated at a rate of approximately $20 per hour for their time. The study protocol was approved by the University of California, Davis Institutional Review Board.

Meditation retreats

Each retreat was taught by a team of six experienced meditation teachers who were not involved in the study design or implementation. Retreat participants lived at Spirit Rock Meditation Center for the duration of their retreat, and maintained noble silence while practicing formal meditation for up to 10 h per day. Meditation practices included sitting and walking variants of Insight meditation [22], during which practitioners were instructed to direct attention to their bodily sensations (and movements, when walking) and to observe their thoughts and emotions as transient phenomena without grasping or over-identification. Participants were also instructed in practices intended to support the four immeasurables: loving-kindness, compassion, empathetic joy, and equanimity. Together, these practices aim to cultivate beneficial qualities of mind and positive aspirations for oneself and others.

Blood sampling procedure

Fasting blood was collected via antecubital venipuncture. Retreatants gave blood at Spirit Rock between 5:30–6:00 a.m. the morning following their first full day of retreat and again three weeks later (one week before the end of retreat).3 Controls were assessed at the Anubhuti Retreat Center in Novato, CA in cohorts of 4–11 people over the course of the following year. Controls were also assessed at the beginning and end of a 3-week interval, during which they maintained their regular routines. At both assessments, controls gave fasting blood between 9:00–10:00 a.m. after completing a 40-min meditation session.

Quantification of gene expression

Whole blood was collected in Vacutainer Cell Preparation Tubes (CPT, Becton Dickinson) and transported to a field laboratory where peripheral blood mononuclear cells were isolated by density-gradient centrifugation within 30 min of collection. Cell pellets were conserved in RNA Later (Sigma St Louis, MO) at −80 °C and stored at the UC Davis Center for Mind and Brain until they were shipped to the University of Barcelona for RNA extraction in P.K.‘s laboratory. Samples were shipped in an insulated shipping container with ∼55lb of dry ice; they were in transit for four days and arrived frozen with dry ice remaining.

Total RNA was extracted using the mirVana™ RNA Isolation Kit (Applied Biosystems). Yield, purity, and quality of RNA were determined spectrophotometrically (NanoDrop, USA) using Bioanalyzer 2100 (Agilent Technologies). The resultant RNA had 260/280 nm ratios above 1.9 and RNA Integrity Numbers (RIN) higher than 7.5.

Quantitative real time (q-RT) PCR was performed for each sample in duplicate using a Bio-Rad CFX384 real-time PCR system (Bio-Rad) with TaqMan FAM-labeled specific probes (Applied Biosystems). The probes used are listed in Supplementary Table 1. A pre-amplification step (TaqMan PreAmp Master Mix; Applied Biosystems) was performed for genes exhibiting very low levels of expression (i.e., PRMT5).

Q-RT PCR data were analyzed with Bio-Rad CFX Manager (Bio-Rad) using the automatic Ct setting to assign baseline and threshold Ct values. The expression of each gene was normalized to the expression of the TBP reference gene, and then the relative expression was calculated using the 2−ΔΔCt method [23]. Gene expression data were natural log transformed to correct for skew prior to analysis.

Multivariate statistical analyses

First, we examined differences in the overall profile of all 33 genes. Gene expression dissimilarities were compared between groups and assessments using multivariate distance matrix regression (MDMR [24]) with the MDMR package in R. MDMR is a person-centered regression method that allows for the estimation of statistical associations between multivariate outcomes and categorical or continuous predictors based on dissimilarities among sets of data [24]. MDMR is ideally suited for analyzing genes defined by their participation in a biochemical pathway or other a priori grouping [25]. Dissimilarity was quantified using Euclidean distance between sets of genes. We also calculated δ effect size statistics to estimate how much each gene contributed to the multivariate effects. δ represents the decrease in pseudo-R2 when an individual gene is removed from the MDMR model relative to the estimate of pseudo-R2 from the full model [24]. Effect sizes were calculated at Time 1 and Time 2 separately based on the between-groups comparison (Fig. 1).

Fig. 1

A)Dendrograms and heatmaps depicting relations between genes across participants at each time point. The dendrograms drawn to the left show the hierarchical clustering between genes. In the heat maps, the red end of the gradient indicates lower gene expression while blue indicates greater expression. Each column represents the gene expression profile for one participant, with retreat participants plotted to the left of the black line and controls plotted to the right. B) Gene expression heatmaps by Group and Time. Each row represents an individual participant; columns represent individual genes. Here, the most notable differences between groups lie in COX2, IL8 and IL1β, which are lower in the retreat group at Time 1 as shown in the upper left panel. There is also a visible pattern of upregulation in CCR3, CXCR2, and GADD45G across groups at Time 2, shown in the panels on the right. (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)

Univariate statistical analyses

Next, we analyzed changes in individual gene expression via linear mixed-effect models. For each gene, models included the factors of Group, Time, and the Group by Time interaction, as well as age and BMI as covariates. Our primary analyses of interest were the Group by Time interactions (i.e., differential changes between groups over time). We also considered main effects of Group (i.e., group differences across both time points), since our baseline assessment occurred after one full day of retreat and prior work has indicated that gene expression differences can emerge after one full day of meditation practice [7]. Our analyses yielded 66 test statistics in total (2 tests × 33 genes). Resultant p-values were subjected to false discovery rate (FDR) corrections to control for Type I error. The expected proportion of false discoveries was set at .05. Observed (uncorrected) and corrected p-values are reported in Table 1, Table 2.

Table 1

Immune System Genes: Mixed-Model ANOVAs with Observed and Corrected p-values.

Outcome	Gene Name	Predictor	F	pobserved	pcorrected
BCL3	B-cell CLL/lymphoma 3	Age	0.33	0.566
		BMI	0.02	0.885
		Group	3.72	0.059	0.118
		Time	4.76	0.034
		Group by Time	0.44	0.511	0.636
CCR3	C–C motif chemokine receptor 3	Age	1.03	0.314
		BMI	0.84	0.362
		Group	1.02	0.316	0.443
		Time	20.06	<.001
		Group by Time	0.17	0.685	0.779
CCR5	C–C motif chemokine receptor 5	Age	0.02	0.878
		BMI	0.07	0.795
		Group	0.05	0.83	0.869
		Time	3.79	0.056
		Group by Time	0.17	0.683	0.779
CCR7	C–C motif chemokine receptor 7	Age	1.59	0.211
		BMI	2.48	0.12
		Group	4.78	0.033	0.083
		Time	7.14	0.01
		Group by Time	5.48	0.023	0.067
COX2 *	cyclo-oxygenase-2 or prostaglandin-endoperoxide synthase 2 (PTGS2)	Age	1.06	0.308
		BMI	0.91	0.343
		Group	82.45	<.001	<.001
		Time	9.88	0.003
		Group by Time	12.05	0.001	0.006
CXCL10	C-X-C motif chemokine ligand 10 or Interferon gamma-induced protein 10 (IP-10)	Age	3.18	0.079
		BMI	0.38	0.538
		Group	1.16	0.286	0.412
		Time	3.43	0.069
		Group by Time	1.75	0.191	0.301
CXCR2	C-X-C motif chemokine receptor 2	Age	0.23	0.632
		BMI	0.03	0.862
		Group	11.23	0.001	0.008
		Time	20.55	<.001
		Group by Time	0.11	0.745	0.82
FKBP5	FKBP prolyl isomerase 5	Age	0.24	0.629
		BMI	1.46	0.231
		Group	9.13	0.004	0.015
		Time	0.03	0.865
		Group by Time	2.53	0.118	0.205
GADD45G	growth arrest and DNA damage inducible gamma	Age	0.7	0.406
		BMI	0.67	0.417
		Group	3.99	0.05	0.111
		Time	19.52	<.001
		Group by Time	0.12	0.733	0.819
GSTP1	glutathione S-transferase pi 1	Age	0.21	0.65
		BMI	0.08	0.781
		Group	6.62	0.013	0.04
		Time	5.44	0.023
		Group by Time	0.22	0.64	0.768
HIF1A	hypoxia inducible factor 1 subunit alpha	Age	2.25	0.139
		BMI	0.07	0.795
		Group	10.86	0.002	0.009
		Time	7.66	0.008
		Group by Time	4.86	0.032	0.083
IL1β	interleukin 1 beta	Age	0.22	0.639
		BMI	2.22	0.141
		Group	24.04	<.001	<.001
		Time	0.3	0.588
		Group by Time	0	0.96	0.96
IL6R	interleukin 6 receptor	Age	0.14	0.712
		BMI	0.01	0.906
		Group	0.24	0.623	0.761
		Time	6.21	0.016
		Group by Time	5.43	0.023	0.067
IL6ST	interleukin 6 signal transducer	Age	0.11	0.739
		BMI	5.36	0.024
		Group	1.98	0.164	0.271
		Time	2.12	0.151
		Group by Time	2.15	0.148	0.251
IL8	interleukin 8 or C-X-C Motif Chemokine Ligand 8 (CXCL8)	Age	0.83	0.365
		BMI	1.54	0.22
		Group	34.15	<.001	<.001
		Time	4.22	0.044
		Group by Time	9.94	0.003	0.011
IFNγ	interferon gamma	Age	5.68	0.02
		BMI	0.18	0.674
		Group	4.08	0.048	0.109
		Time	0.37	0.544
		Group by Time	4.28	0.043	0.103
NOD2	nucleotide binding oligomerization domain containing 2	Age	1.12	0.294
		BMI	0.52	0.476
		Group	3.38	0.071	0.134
		Time	6.91	0.011
		Group by Time	4.27	0.044	0.103
NR3C1	nuclear receptor subfamily 3 group C member 1	Age	4.6	0.037
		BMI	0.79	0.378
		Group	1.42	0.239	0.359
		Time	3.06	0.087
		Group by Time	0.7	0.408	0.539
RETN	resistin	Age	0	0.947
		BMI	0.01	0.932
		Group	3.6	0.063	0.122
		Time	59.14	<.001
		Group by Time	3.26	0.076	0.136
RIPK2*	receptor interacting serine-threonine kinase 2	Age	0.39	0.536
		BMI	0.13	0.716
		Group	9.61	0.003	0.012
		Time	3.43	0.069
		Group by Time	5.06	0.028	0.078
TNFα*	tumor necrosis factor alpha	Age	6.79	0.011
		BMI	4.02	0.049
		Group	3.77	0.057	0.117
		Time	12.96	0.001
		Group by Time	8.51	0.005	0.018
TNFRSF1A	TNF receptor superfamily member 1A	Age	0.27	0.604
		BMI	0	0.959
		Group	6.92	0.011	0.037
		Time	14.54	<.001
		Group by Time	13.63	0.001	0.004
TNFRSF1B	TNF receptor superfamily member 1B	Age	0.15	0.704
		BMI	0.02	0.903
		Group	0.83	0.367	0.504
		Time	11.13	0.001
		Group by Time	10.6	0.002	0.01

Note: Expression values were natural log transformed prior to analyses. We FDR corrected the p values testing Group and Group by Time effects, with the proportion of false discoveries (Q) across both effects set to .05. Significant effects are in bold. Genes with a priori hypotheses are marked with an asterisk (*).

Table 2

Epigenetic Modulatory Genes: Mixed-Model ANOVAs with Observed and Corrected p-values.

Outcome	Gene Name	Predictor	F	pobserved	pcorrected
AURKC	aurora kinase C	Age	0.04	0.848
		BMI	0.04	0.835
		Group	0.18	0.676	0.779
		Time	0.54	0.465
		Group by Time	1.58	0.213	0.327
DNMT3A*	DNA methyltransferase 3 alpha	Age	14	<.001
		BMI	0.15	0.695
		Group	7.57	0.008	0.027
		Time	0.85	0.361
		Group by Time	3.34	0.073	0.134
HDAC2*	histone deacetylase 2	Age	1.7	0.197
		BMI	1.63	0.206
		Group	1.83	0.181	0.292
		Time	0.79	0.379
		Group by Time	1.15	0.287	0.412
HDAC3*	histone deacetylase 3	Age	0.28	0.599
		BMI	0	0.951
		Group	3.77	0.057	0.117
		Time	0.03	0.862
		Group by Time	0.47	0.497	0.63
HDAC9*	histone deacetylase 9	Age	0.18	0.677
		BMI	1.17	0.283
		Group	15.11	<.001	0.002
		Time	0.14	0.714
		Group by Time	18.53	<.001	0.001
KAT6B	lysine acetyltransferase 6B	Age	0.32	0.576
		BMI	0.12	0.732
		Group	9.88	0.003	0.011
		Time	0.02	0.899
		Group by Time	0.8	0.374	0.504
KAT7	lysine acetyltransferase 7	Age	0.08	0.772
		BMI	0	0.968
		Group	0.09	0.76	0.822
		Time	0.49	0.486
		Group by Time	0.01	0.943	0.957
KMT2G	lysine-specific methyltransferase 2G or SET domain containing 1B (SETDB1)	Age	0.31	0.582
		BMI	3.02	0.087
		Group	23.13	<.001	<.001
		Time	1.17	0.284
		Group by Time	0.48	0.49	0.63
KMT5B	lysine methyltransferase 5B or suppressor of variegation 4–20 homolog 1 (SUV420H1)	Age	0.01	0.915
		BMI	1.6	0.21
		Group	12.13	0.001	0.006
		Time	4.33	0.042
		Group by Time	0.06	0.8	0.851
PRMT5	protein arginine methyltransferase 5	Age	0.06	0.813
		BMI	0.04	0.843
		Group	0.01	0.941	0.957
		Time	0.01	0.921
		Group by Time	15.63	<.001	0.002

Results

Multivariate gene analyses

The retreat and control groups significantly differed in their gene expression profiles at Time 1 (p < .001) and at Time 2 (p < .001) as demonstrated by multivariate distance matrix regression (MDMR) of Euclidean distances between all 33 genes. Both retreat (p < .001) and control (p = .012) participants’ profiles changed over time, but these changes were qualified by a significant Group by Time interaction suggesting that retreat participants changed differently than controls (p = .003). The dendrograms and heat maps depicted in Fig. 1A show that COX2, IL8 and IL1β clustered together, with directionally similar patterns across participants. Effect size measures depicted in Fig. 1B also suggest that Il1β, COX2, and IL8 contributed most to the multivariate differences observed between groups. Interestingly, retreat participants showed a fairly consistent pattern of downregulation in these genes at Time 1, with a subtle but visible increase in expression across retreat participants at Time 2. Controls, on the other hand, showed more variability and upregulation in these genes at Time 1, with three distinct patterns of expression emerging at Time 2 (i.e., upregulation, downregulation, no change). It is also of note that CXCR2, GADD45G, and CCR3 clustered together at both time points, but showed greater upregulation at Time 2 across all participants. Moreover, COX2, IL8 and IL1β expression appears to be inversely related to CXCR2, GADD45G, and CCR3 expression, as these clusters show the greatest distance from one another in the dendrograms. Neither age (p = .347) nor BMI (p = .347) significantly predicted overall gene expression profiles.

Univariate group differences across assessments

Next, we examined individual genes. First, we tested the main effect of Group on each gene. After FDR correction, fourteen analyses remained significant (Table 1, Table 2). The retreat group had lower expression of the immune-system-regulating genes COX2, RIPK2, HIF1A, GSTP1, Il1β, IL8, and TNFRSF1A, but higher expression of CXCR2 and FKBP5. The retreat group also had higher expression of the epigenetic-modulating genes DNMT3A, HDAC9, KAT6B, KMT2G, and KMT5B. Group effects not qualified by an interaction are presented in Fig. 2 (i.e., Il1β, RIPK2, HIF1A, GSTP1, CXCR2, and FKBP5 shown in Fig. 2A; DNMT3A, KMT2G, KMT5B, and KAT6B shown in Fig. 2B). Importantly, the group differences in COX2, IL8, TNFRSF1A, and HDAC9 were qualified by significant Group by Time interactions described below. There were no significant effects of Group or Group by Time interactions in BCL3, CCR3, CCR5, CCR7, CXCL10, GADD45G, IL6R, IL6ST, IFNγ, NOD2, NR3C1, RETN, AURKC, HDAC2, HDAC3, or KAT7.

Fig. 2

Significant Group effects that survived False Discovery Rate correction. Expression values were log transformed. Error bars represent standard errors of the mean. Group differences at Time 1 or Time 2 are indicated: ns = not significant, *p < .05, **p < .01, ***p < .001.

Univariate group differences emerging over time

Our primary analyses resulted in seven Group by Time interaction effects that remained significant after FDR correction (Fig. 3). The groups showed differential changes over time (i.e., interaction effects) in five immune-related genes (COX2, IL8, TNFα, TNFRSF1A, and TNFRSF1B; Table 1) and two genes involved in histone modification (HDAC9 and PRMT5; Table 2). Among these significant interactions, COX2, IL8, HDAC9 and PRMT5 showed group differences at Time 1 with differential trajectories across assessments. TNF genes, on the other hand, showed no group differences at Time 1, but displayed a consistent pattern of retreat-related downregulation from Time 1 to Time 2.

Fig. 3

Significant Group × Time effects that survived False Discovery Rate correction.

Gene expression values were log transformed. Error bars represent standard errors of the mean. Group differences at Time 1 or Time 2 are indicated: ns = not significant, *p < .05, **p < .01, ***p < .001.

Decomposing these effects, the retreat group (M = −1.98, SE = 0.15) showed lower COX2 expression than controls (M = −0.08, SE = 0.13) at Time 1, t (108) = 9.65, p < .001, as well as Time 2, t (110.6) = 5.15, p < .001 (Time 2 Retreat: M = −1.19, SE = 0.15; Control: M = −0.12, SE = 0.14). However, retreatants showed an increase in COX2 across assessments, t (59.1) = −4.5, p < .001, that was not observed in controls, t (60.1) = 0.25, p = .806.

Similarly, the retreat group (M = −1.69, SE = 0.18) showed lower IL8 expression than controls (M = −0.06, SE = 0.16) at Time 1, t (96.7) = 6.80, p < .001, and Time 2, t (101.5) = 3.285, p = .001 (Time 2 Retreat: M = −1.03, SE = 0.18; Control: M = −0.20, SE = 0.17). Like COX2, IL8 expression increased in retreatants over time, t (57.9) = −3.55, p = .001, but not in controls, t (58.4) = 0.81, p = .419.

In contrast to these patterns, retreatants (M = 0.40, SE = 0.05) had higher HDAC9 expression than controls (M = 0.04, SE = 0.04) at Time 1, t (103.1) = −5.78, p < .001. HDAC9 then decreased across assessments in retreatants, t (59.1) = 3.18, p = .002, while it increased in controls, t (59.9) = −2.91, p = .005, resulting in no difference between the retreat (M = 0.23, SE = 0.05) and control (M = 0.18, SE = 0.04) groups at Time 2, t (106.8) = −0.81, p = .417.

For PRMT5, retreatants (M = −0.37, SE = 0.10) had significantly lower expression than controls (M = 0.01, SE = 0.09) at Time 1, t (110.5) = 2.61, p = .01. The opposite pattern emerged at Time 2, with the retreat group having higher expression (M = −0.01, SE = 0.11) than controls (M = −0.38, SE = 0.09), t (110.7) = −2.61, p = .01, due to an increase across assessments in retreatants, t (56.1) = −2.59, p = .012, alongside a decrease in controls, t (57.1) = 3.03, p = .004.

Interestingly, we observed evidence of retreat-related declines in the expression of TNFα, and the TNFα receptor genes, TNFRSF1A, and TNFRSF1B, as demonstrated by the emergence of group differences at Time 2 that were not present at Time 1 (Fig. 3).

There was no difference in TNFα between retreatants (M = −0.05, SE = 0.06) and controls (M = −0.04, SE = 0.06) at Time 1, t (99.2) = 0.11, p = .914. However, retreatants decreased in TNFα over time, t (59.5) = 4.43, p < .001, while controls showed no change, t (60.1) = 0.5, p = .62, resulting in lower TNFα expression in retreatants (M = −0.34, SE = 0.06) than controls (M = −0.07, SE = 0.06) at Time 2, t (103.6) = 3.12, p = .002.

Likewise, retreatants (M = 0.008, SE = 0.05) and controls (M = 0.004, SE = 0.04) showed no difference in TNFRSF1A at Time 1, t (105.3) = −0.05, p = .96. Retreatants, however, decreased across assessments, t (51.7) = 5.10, p < .001, while controls showed no change, t (52.7) = 0.09, p = .933, again resulting in lower expression in retreatants (M = −0.29, SE = 0.05) than controls (M = 0.00, SE = 0.05) at Time 2, t (108.5) = 4.19, p < .001.

This same pattern was observed for TNFRSF1B: there was no difference between retreatants (M = 0.11, SE = 0.06) and controls (M = 0.01, SE = 0.05) at Time 1, t (115.7) = −1.5, p = .138, but retreatants (M = −0.24, SE = 0.06) had significantly lower TNFRSF1B expression than controls (M = −0.01, SE = 0.06) at Time 2, t (116.2) = 2.75, p = .007, as a result of a decrease in retreatants, t (57.8) = 4.48, p < .001, alongside no change in controls, t (59.1) = 0.06, p = .955.

Discussion

Residential retreats are an increasingly popular form of meditation training, yet relatively little is known about how these interventions might alter the molecular underpinnings of health. Here, we addressed this gap by examining changes in the expression of genes directly and indirectly involved in the regulation of inflammation over the course of a month-long meditation retreat. We report a consistent pattern of downregulation in the TNF pathway, as indicated by retreat-related reductions in TNFα, TNFRSF1A, and TNFRSF1B. We also observed differential patterns in several other genes indicative of lower inflammation and greater transcriptional repression in the retreat group.

Downregulation of the TNF pathway

Downregulation of the TNF pathway is indicative of reduced inflammation in the meditation group. Given recent interest in developing anti-TNF therapies to treat depression [14] and inflammatory disorders [26], these results suggest that meditation retreats may be a useful, non-pharmacological tool for modulating the TNF pathway.

Reductions in the TNF pathway are particularly interesting in light of prior findings from this study: As reported elsewhere, we found retreat-related increases in bulk telomere length [21] that were predicted by basal levels of brain derived neurotrophic factor (BDNF). Telomeres are nucleoprotein complexes that flank linear chromosomes, protecting coding DNA from degradation and instability. They play a critical role in initiating cellular senescence, and shorter telomeres have been linked to higher levels of TNF-α and IL-6 [27]. Presumably, the retreat-related downregulation of the TNF pathway reflects reduced inflammation, which may, in turn, have contributed to less telomere degradation and better telomere maintenance. TNF-α also negatively regulates neurogenesis through its TNFRSF1A receptor, and BDNF expression is necessary for neurogenesis [28]. The observed reductions in TNFα and TNFRSF1A may, therefore, suggest an increased capacity for neurogenesis in our retreat participants, despite no group-level change in BDNF.

Our results are also consistent with the finding that one day of meditation practice led to greater methylation, and thus suppression, of TBKBP1 and TNFSF13B, which are involved in the TNF and NF-kB pathways [18]. Moreover, when Epel et al. [29] compared experienced and novice meditators attending a week-long retreat to a vacationing control group residing at the same resort, they found that TNF-α protein levels were maintained in experienced meditators, but increased in the novice and vacation control. Although these findings are not entirely consistent with our observations, they both point to lower TNF-α-mediated inflammation in experienced meditators on retreat.

Other indicators of lower inflammatory burden in retreat participants

Despite increases in IL8 and COX2 across assessments, Il1β, COX2, and IL8 were consistently lower in the retreat group. Effect size measures further suggest that these genes contributed most to the multivariate differences observed between groups (Fig. 1B). Although we cannot determine whether these group differences were the result of early intervention effects or pre-existing group differences, they too are consistent with a lower inflammatory burden in retreatants.

Interleukin-(IL-)1β is a proinflammatory cytokine that activates microglia, regulates the activity of growth factors, and can stimulate immune cells to produce other proinflammatory cytokines [28]. Both TNF-α and IL-1β are implicated in sickness behavior and depressive phenotypes induced by stress [1,13,[30], [31]]. In fact, IL-1β mediates the stress-induced inhibition of neurogenesis in the hippocampus, indicating that it likely counteracts BDNF, which promotes neurogenesis [28,31]. When paired with downregulation of the TNF pathway, lower IL1β expression in our retreat group further suggests a link between retreat experience and lower inflammation.

In the retreat group, consistently lower levels of IL8 were coupled with higher levels of CXCR2, which codes for one of two IL-8 receptors (the other being CXCR1 or IL8RBP). Despite having lower levels of IL8 across time, the retreat group showed increases in both IL8 and CXCR2. IL-8, also known as CXCL8, is a chemoattractant involved in the migration and activation of neutrophils during acute inflammation [32]. IL-8 has been implicated in chronic inflammatory pain states [33] and is upregulated in the anterior cingulate cortex in mice with persistent inflammatory pain [34].

Interestingly, Creswell et al. [35] found lower IL8 and CXCR1 expression in older adults who were lonely, compared to those who were not—but they observed no change in IL8 or its receptor in relation to an 8-week MBSR intervention. Taken together with our findings, these results could indicate that the IL-8 pathway is implicated in loneliness in older adults, and that more intensive meditation interventions are needed to modulate this relationship. On the other hand, Black et al. [36] found IL8 downregulation in PBMCs following an 8-week Kirtan Kriya meditation intervention, which is more consistent with the group difference we observed here.

RIPK2 and HIF1A were also lower in retreatants than controls, and showed similar, but non-significant increases in the retreat group across time. These results are indicative of differences in NF-κβ signaling. Moreover, the group differences we observed in COX2 and RIPK2 after one day of retreat are consistent with reductions observed by Kaliman et al. [7] after one day of practice. It is important to note, however, that we observed an increase in COX2 expression in our retreat group, whereas Kamilan et al. [7] observed a decrease in COX2 after a day of practice. These findings appear to contradict one another, although it is possible that the substantial group difference at Time 1 reflects a rapid effect of one day of practice that then subsided across the course of the retreat. Further research will be needed to determine if these contradictory effects are spurious or reflect differences in the time course of training.

We also observed lower levels of GSTP1 in our retreat group across assessments. GSTP-1 belongs to a family of glutathione S-transferases, which are antioxidant enzymes that play an important role in detoxification. Interestingly, this finding, which may suggest lower antioxidant activity in the retreat group, contradicts a finding by Sharma et al. [37], who report higher GSTP1 expression in regular practitioners of Sudarshan Kriya—a breathing exercise used to induce relaxation.

Finally, we observed higher levels of FKBP5 in our retreat group. FKBP5 modulates glucocorticoid receptor activity in response to stress [38]. This finding contradicts those by Creswell et al. (2012), who found downregulation of FKBP5 after mindfulness-based stress reduction, and by Bishop et al. [[39], [40]] who found increased methylation of the FKBP5 gene in individuals with PTSD who responded to a mindfulness-based stress reduction intervention compared to non-responders. Again, further research will be needed to determine whether these discrepancies are due to differences in the intensity or style of practice, individual differences (e.g., in mental health profiles), or other environmental factors.

Expression of epigenetic modulators

With respect to epigenetic modulators, we found few effects attributable to the retreat intervention. Consistent with Kaliman et al. [7], who observed a decrease in HDAC9 across one day of practice, we observed a decrease in HDAC9 in the retreat group across assessments (Fig. 3). However, in contrast to Kaliman et al. [7], we found no change in HDACs 2 or 3.

The retreat group also had higher expression of the methyltransferases DNMT3A, KMT2G, and KMT5B, and the acetyltransferase KAT6B across both time points, as well as lower expression of the methyltransferase PRMT5 at Time 1, which significantly increased by Time 2. PRMT5 regulates inflammation by recruiting DNMT3A to repress genes downstream [41]. DNMT3A also interacts with HDACs and the repressive histone-methyltransferase KMT2G (also known as SETD1B) to silence gene transcription [[42], [43]]. KMT2G is involved in regulating the activation of genes downstream of NF-κB [44]. KMT5B (also known as SUV420H1) is involved in the methylation of histone 4, and the age-related loss of this methylation is associated with human cancers [45,46]. KAT6B (also known as MORF) can also promote and suppress transcription [47], and is involved in the epigenetic regulation of neurogenesis in the mammalian brain [48]. Although each of these epigenetic regulators is involved in multiple functions, collectively their levels being consistently higher in the retreat group suggests greater transcriptional repression, likely contributing to lower inflammation.

Unfortunately, the logistical constraints of our study design make it difficult to interpret the health-related implications of these findings, because it impossible to determine the source of the group differences observed. On the one hand, these group difference could reflect meditation effects that were initiated during the first full day of retreat and maintained by subsequent retreat practice. However, it is also possible that they reflect pre-existing group differences. They might also reflect seasonal or circadian fluctuations, as the control samples were collected at various times throughout the year and at a different time of day than retreatants. Future research involving randomization to short-term interventions are needed to clarify these possibilities.

Strengths and limitations

One methodological advantage of our study is the ability to examine the effects of a high dose of concentrated practice in experienced meditation practitioners [20]. Most meditation research has focused on the effects of short-term interventions as they are easier to control, and thus often easier to interpret. However, for many people, meditation represents a life-long endeavor to cultivate well-being and to mitigate suffering. Meditation practice can also be effortful or puzzling at first—though a practitioner's experience of, and relation to, meditation is bound to shift over time. It is therefore important to investigate how psychobiological processes differ between short- and long-term interventions, and for practitioners with varying degrees of prior experience. For instance, Kaliman et al. [7] found no change in TNFα expression after one day of practice, but we find significant downregulation after 21 days of practice, suggesting that modulation of the TNF pathway may require more intensive practice.

Another advantage of our study is the use of targeted gene probes. Many of the initial studies examining meditation-related changes in gene expression have relied on genome-wide microarrays. These extensive arrays are valuable exploratory tools to begin to understand which pathways might be influenced by a particular intervention, but require targeted follow-ups using more reliable assays as we have done here.

Unfortunately, studying a well-established and regularly offered retreat also comes with necessary logistical and methodological constraints, such as a limited sample size [20]. Given the length and intensity of the retreat intervention, it was infeasible to randomly assign participants to retreat or control conditions, or to devise a comparable sham intervention. Instead, we aimed to maximize the ecological validity of our participant sample by recruiting a comparison group with an existing interest in meditation, an ongoing practice, and prior retreat experience. Conceivably this allows us to better generalize our results to a population of like-minded individuals as opposed to a population of non-practitioners [49]. However, the lack of randomization and an active control condition preclude causal inference. While our results do indicate retreat-related changes in TNFα and its receptors, true randomization would be necessary to determine whether the retreat caused these changes. It should be further noted that many of the within group changes that we observed were relatively small compared to the between group differences, meaning that the biological significance of these changes is unclear. Our design also lacks a true baseline assessment, which limits our ability to distinguish true group differences from early intervention effects. Finally, we cannot determine the degree to which features of the retreat other than meditation (e.g., vegetarian diet, natural environment) contributed to the changes observed.

A further limitation of this study is our inability to control for potentially confounding seasonal variables. In this case we were unable to assess the retreat and control participants simultaneously, or in a stratified manner. As gene expression can be responsive to environmental factors such as temperature and daylight hours [50], it will be important for future studies to stratify assessments and to devise methods to account for these potentially confounding variables. Repeated studies investigating shorter retreat interventions will be important for overcoming these logistical constraints and bridging gaps in the existing literature.

Conclusions

In sum, the differential patterns observed between our retreat and control groups suggest that engaging in retreat practice may promote beneficial changes in gene expression indicative of lower inflammation. These findings are notable, as chronic, low-grade inflammation is associated with many modern health problems [51]. The data presented here also suggest that the TNF pathway should be further investigated with respect to meditation training. Given the limitations of our study, these effects should be considered preliminary and require replication. Nevertheless, these results set the foundation for future randomized and controlled studies to test whether intensive meditation interventions might contribute to the treatment and prevention of inflammatory conditions.

CRediT author statement

QC: Writing - Original Draft, Investigation, Project administration, Data Curation, Formal analysis, Software, Visualization; GS: Writing - Original Draft, Formal analysis, Visualization; MÁ: Conceptualization, Investigation, Data Curation; Writing - Review & Editing; MC: Conceptualization, Investigation, Data Curation; Writing - Review & Editing; BK: Conceptualization, Investigation, Data Curation, Writing - Review & Editing; AZ: Conceptualization, Investigation, Data Curation, Formal analysis, Software, Visualization; PK: Conceptualization, Funding acquisition, Supervision, Writing - Review & Editing; CS: Conceptualization, Funding acquisition, Investigation, Resources, Supervision, Writing - Review & Editing.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.