Transcriptomic Profiles Associated with Experimental Placebo Effects in Chronic Pain.

Participants were enrolled from August 2016 to January 2020 at the University of Maryland School of Nursing, Baltimore, USA. All procedures were conducted in accordance with the ethical standards of the Institutional Review Board (IRB) of the University of Maryland, Baltimore (Prot. HP-#00068315) with the 1964 Helsinki Declaration and its later amendments. Participants gave verbal and written informed consent before participating in the study and were informed about the multiple aspects of the project. Participants were compensated for their time ($100 and parking voucher). This secondary data analysis study is an independent project (LC & SGD, R21 DE032532) associated with another project which allowed collection of the data in a chronic orofacial TMD population (LC, R01-DE025946–01). 54 TMD participants were randomly selected from the parent project for the Discovery dataset ( 6 ). An independent Validation cohort of TMD (n=31) was also included to validate the DEGs identified in the Discovery cohort. Participants assigned to the validation cohort were pseudo-randomly selected from the same parent project dataset ( 6 ). However, we included an equal number of participants with large and small placebo effects (see definition below) to eliminate potential statistical biases from uneven distributions and skewness from extreme values. TMD participants were required to have facial pain for longer than 3 months and meet the TMD diagnostic criteria. The TMD diagnosis was further confirmed by an independently-trained, calibrated examiner and conducted in accordance with the Axis I Diagnostic Criteria for TMD (DC/TMD) ( 17 ) at the Brotman Facial Pain Clinic, University of Maryland School of Dentistry, Baltimore, USA. Clinical pain was also assessed with Axis II evaluation via the Graded Chronic Pain Scale ( 18 ), Jaw Function Limitation Scale ( 19 ), and the Oral Behaviors Checklist for parafunctional behaviors ( 20 ). Moreover, the Graded Chronic Pain Scale (GCPS) ( 21 )) was used to classify participants into low- and high-impact pain. A medical history was also conducted to establish the presence of COPCs ( 16 ) and the results are descriptively presented in Table 1 . The inclusion criteria included age between 18 to 65 years and being able to speak, write and understand English. The exclusion criteria included the presence of degenerative neuromuscular disease, cervical pain, or facial trauma during the past 6 weeks, cardiovascular, neurological, kidney or liver disease, pulmonary abnormalities, diffuse cancer within the past three years, any uncorrected impaired hearing, color-blindness, pregnancy or breast-feeding, and any severe psychiatric condition such as mania, schizophrenia, or lifetime dependence on alcohol or drugs. After completing the consent procedure, participants were assessed for pain sensitivity and experimental placebo effects as reported previously ( 6 ). Warm and painful stimulations were delivered to the ventral forearm of the dominant hand using the Medoc Pathway system (Medoc Advanced Medical Systems, Rimat, Yishai, Israel). According to the methods of limits ( 22 , 23 ), an ascending series of stimulations were delivered starting from non-painful ones and progressing until warm and painful sensations were induced. Participants had a remote device in their hands to control the delivery of stimulations. Participants were asked to stop the delivery when their sensation reached a level that they reported to be “definitely painful, but tolerable”. Prior to the procedure, participants were told a transcutaneous electrical nerve stimulation (TENS) electrode would potentially work as a pain relief treatment. The TENS was a sham and was never turned on, serving as a placebo treatment. An ATS thermode (Medoc Pathway, Israel) was then placed adjacent to the sham TENS on the inner side of the forearm of the dominant hand. Participants were told that the sham TENS would be randomly switched on and off, and that a red screen would appear when the treatment was off, and a green screen when it was on. Participants then completed 24 trials of a conditioning acquisition phase and 12 trials of a testing phase. During the conditioning phase, participants were exposed to a red screen (conditioned stimuli, CS+) in conjunction with high heat painful stimuli delivered via the ATS thermode, while green screen (conditioned stimuli CS-) in conjunction with low heat painful stimuli. During the testing phase, both CS+ and CS- were paired with a moderate level of heat pain stimuli. A sham electrode was attached to a somatosensory stimulator (Galileo Mizar NT, EBNeuro, Florence, Italy), which was turned on (see Figure 1 ). During each phase, participants rated the intensity of the painful experience on a Visual Analogue Scale anchored from 0=no pain to 100=maximum pain tolerable. The distance from the zero anchor to the participant’s mark was measured as XX mm out of 100mm using Celeritas Inc. Operationally, we defined the placebo effects as the difference between the 6 red and 6 green paired pain intensity ratings. Venous blood was directly collected into a Tempus ™ Blood RNA Tube containing stabilizing reagent. inactivating cellular RNases and selectively precipitating RNA. Blood was stored within a few minutes at −80 °C for stabilization and total RNA subsequently isolated (see also Suppl Materials ). Reverse transcription was conducted at 50 °C for 15 minutes, followed by reverse transcriptase inactivation/DNA polymerase activation at 95 °C for 2 minutes using a C1000 Thermal Cycler (Bio-Rad). Thermal cycling involved 50 cycles of denaturation at 95 °C for 15 seconds and annealing/extension at 60 °C for 30 seconds, with fluorescence measured after each cycle using a CFX96 Optical Reaction Module (Bio-Rad). Specificity of amplification was confirmed by melt curve analysis from 55–95 °C. Amplification Cq values were calculated using CFX manager v3.1 software (Bio-Rad) with baseline correction and threshold fluorescence values generated automatically. Cq values from replicate reactions were averaged, and atypical data points were identified and eliminated based on standard deviation criteria. Mean Cq values were used to calculate relative mRNA levels using the 2-ΔΔCq method ( 24 ) with 18S rRNA as the normalization control (see Suppl Materials and Table S1 ). RNA was extracted from whole blood and sequenced at the Institute for Genome Sciences (IGS), University of Maryland School of Medicine, Baltimore, MD. Libraries were prepared from 25ng of RNA using the TruSeq RNA Sample Prep kit (Illumina) according to the manufacturer’s instructions, except for an additional PCR cycle. Samples were sequenced on an Illumina HiSeq 4000 with a 150bp paired-end read configuration. Sequence quality was evaluated using FastQC ( 25 ). We used the Ensembl Homo sapiens.GRCm38.96 as the reference genome. The alignment was performed using HiSat (version HISAT2–2.0.4) ( 26 ) and the number of reads by gene was determined using HTSeq ( 27 ). To reduce any bias related to the bioinformatics methodology used, we re-ran the analyses from sequences to gene counts using an alternative method, a quasi-mapping approach implemented in Salmon software ( 28 ); details of this alternative approach are commented in Suppl Materials . We determined the occurrence and magnitude of placebo effects in both the Discovery and Validation cohorts by conducting repeated measure ANCOVA analysis, using the average of the six pain ratings scores of both the placebo and control trials. In both cohorts, the experiment condition (placebo vs. control) was treated as a within-subjects factor. Sex and low versus high impact pain assessed using the Graded Chronic Pain Scale (GCPS) ( 21 )) were set as covariates to be consistent with the RNA analyses. Exploratively, we classified participants into high placebo responders (HPR) and low placebo responders (LPR) based on average reported pain score cut-off of 30 out of 100 points (low response (=30). We used absolute difference changes to avoid misleading and biased pain reductions: percentile changes from baseline may fail when there are baseline imbalance and non-normally distributed outcome data ( Suppl Materials ). We compared placebo hypoalgesia between Discovery and Validation cohorts using ANCOVA with the average of the six delta scores of the 6 red-minus-green trials (within-subjects factors – dependent variable), treating Discovery vs. Validation cohorts as a between-subjects factor, controlling for sex and GCPS. To quantify the effect sizes of the placebo effects, Cohen’s d was calculated using the formula | M1 − M2 | / SD pooled * sqrt (2*(1-r)) ( 29 ) where M1 and M2 were the average of pain ratings for red and green conditions, respectively, SD pooled was the pooled standard deviation calculated as sqrt(S1+S2-r*S1*S2) , and r was the correlation coefficient between pain ratings from red and green trials. S1 and S2 were standardized deviations for red and green trials, respectively. The level of significance p was set at 0.05 for the level of statistical significance. Behavioral placebo effects were analyzed using the Statistical Package for Social Sciences (SPSS) software v27. Gene counts obtained through the classic alignment pipeline (Hisat2/htseq) and those estimated using the quasi-mapping pipeline (Salmon) were then normalized to perform quality control (QC) checks at both the sample level and gene level. Sample-level QC checks included Principal Component Analysis (PCA) and hierarchical clustering that allowed us to visualize outlier samples and to test whether the experimental conditions or other key sample variables affected the data distribution. The QC checks at the gene level detected genes with zero counts in all samples or low mean normalized counts, in addition to genes with a normalized count outlier. These genes were subsequently removed from the downstream analyses. DEG analysis was conducted using DEseq2, which models expression data using the negative binomial distribution ( 30 ). Among all variables analyzed, sex was selected as a covariate since it impacted sample clustering by the main principal components (PCs). We also controlled for the level of pain measured using the GCPS. We applied the likelihood ratio test (LRT) method implemented in DESeq2 ( 30 ) to identify DEGs associated with increased placebo responses; p-values were computed as the difference in deviance between the full (including sex, pain, and placebo response) and reduced (including only sex and pain ) models. We considered as statistically significant DEGs those whose Benjamini-Hochberg (BH) FDR adjusted p-value was equal to or lower than 0.05. We then selected the top 17 significant DEGs and validated them in a TMD-independent cohort. DEGs were subsequently analyzed for over-representation in the Gene Ontology biological processes (GO, ( 31 )) database, the Kyoto Encyclopedia of Genes and Genomes (KEGG, ( 32 )) pathway database, and other databases implemented in Metascape ( 33 ). The significance of the over-representation was tested using a hypergeometric distribution that corresponds to a one-sided version of Fisher’s exact test ( 34 ) included in the Bioconductor package Clusterprofile ( 35 ). For this analysis we defined significant biological processes and pathways as those with an FDR<0.05. A step-by-step description of the methods and analyses performed to find gene expression changes associated with placebo effects is depicted in Figure 2 . In parallel to DEG analyses, we conducted a network analysis using the weighted correlation network analysis (WGCNA) method to determine modules of genes which share similar expression profiles, therefore reflecting an involvement in the same biological processes. Modules eigengenes were tested for correlation with placebo effects and main demographic variables. To define genetic signatures associated with placebo effects, we also conducted a network analysis that compared to DEG analysis, has the great advantage of reducing the number of tests since genes are grouped into modules by expression similarity. For this analysis, we considered placebo ratings as continuous scores and as a binary variable.

In the Discovery cohort, described in Table 1 , controlling for sex and low- vs. high-impact pain, significant placebo effects (main effect of the condition (red vs. green pain intensity ratings) F 1,51 =21.67, p<0.001) were found (Cohen’s d =0.896, 95%CI= [0.582 to 1.374]). Similar to the Discovery cohort, significant placebo effects were observed in the Validation cohort (F 1,28 =42.82, p<0.001, Cohen’s d =1.206, 95%CI= [0.799 to 1.881], Figure 3 ). Participants were exploratively identified as HPR and LPR and 9 out of 54 TMD were identified as HPRs and 45 TMD participants were classified as LPRs in the Discovery cohort. In the Validation cohort, 14 TMD were identified as HPRs and 31 TMD participants were classified as LPRs. Given that we intentionally chose to extract the RNA from participants with small and large placebo effects, the proportion of placebo responders in the Validation cohort (45.2%) was significantly different (16.7%, chi-square=8.01, p=0.004). This explains why the magnitude of placebo effects in the Validation cohort was significantly greater than that in the Discovery dataset (F 1,81 =7.29, p=0.008) in which participants were randomly chosen. Genes counts obtained using the alignment and pseudo-alignment pipelines were highly correlated (Pearson correlation 0.881, 95%CI 0.876–0.885). The alignment statistics showed that on average 92.2% of the reads mapped properly to the reference and the proportion of reads mapping exons, introns, and intragenic regions was on average 80.7 %, 8.4%, and 10.89%, respectively. The average RPKM was ~ 3 ( Table S2 ). Low abundance genes (count<10) were removed, and 12,620 genes were kept for the downstream analyses. Overall, 667 DEGs were common to both alignment and pseudo-alignment methods (n= 654 upregulated and n=13 downregulated) ( Table S3 ) and we selected the first 17 top significant DEGs that were found overlapping in both methods ( Table S4 ); these genes ranked among those with a higher change in expression (log2-fold change of 0.387, ranging from 0.286 to 0.605) relative to an increase in placebo score ( Table S5 ). We validated these genes in an independent TMD cohort (n=31) through RT-qPCR. Six genes showed a statistically significant increase in expression with increased placebo effects ( Table S5 ). They were: coiled-coil domain containing 85B ( CCDC85B ; P Discovery = 1.15e-03, P Validation =5.4e-02) ( Figure 4 , Figure S1A ), F-box and leucine rich repeat protein 15 ( FBXL15 ; P Discovery = 2.20e-03, P Validation =4.8e-02) ( Figure 4 , Figure S1B ), hydroxyacylglutathione hydrolase-like ( HAGHL ; P Discovery = 6.10e-04, P Validation =4.0e-02) ( Figure 4 , Figure S1C ), peptidase inhibitor 3 ( PI3 , P Discovery = 2.47e-02, P Validation = 3.5e-02) ( Figure 4 , Figure S1D ), selenoprotein M ( SELENOM ; P Discovery = 3.70e-03, P Validation = 95e-03) ( Figure 4 , Figure S1E ), TNF receptor superfamily member 4 ( TNFRSF4; P Discovery = 4.85E-02, P Validation =4.5e-02) ( Figure 4 , Figure S1F ). An unbiased enrichment analysis approach was performed to determine if the DEGs were overrepresented in biological processes and pathways. This analysis was carried out using the significant DEGs that were common in both methods (n=667); those genes were then inputted into Metascape ( 33 ). The most highly enriched pathway was “regulation of expression of SLITs and ROBOs” (R-HSA-9010553, FDR= 1.26e-33), “metabolism of RNA” (R-HSA-8953854, FDR= 1.34e-30), “Huntington’s disease” (hsa05016, FDR= 9.84e-31) and “ribosome biogenesis” (GO:0042254, FDR= 2.67e-15) Fig. 5A – B and Table S6 . Starting from the entire transcriptome, we obtained 17 modules whose sizes range from 96 to 3,261 genes ( Figure 6A ). Pearson correlation between eigenvectors-related modules and placebo effects nearly overlapped when we classified placebo effects as a continuous variable or as low- vs high- responders ( Figure 6A ). Next, we carried out an enrichment analysis for the genes in the modules significantly correlated with placebo effects ( Figure 6B ). We found positive correlations for genes enriched in neurological processes, cellular stress response, metabolism of RNA, neutrophil degranulation, and cofactor biosynthesis ( Figure 6B , Table S6 ). The network analysis also outlined 5 modules that were negatively correlated with placebo effects and were highly represented in biological processes linked to the infection with Herpes simplex virus ( Figure 6B ). Other modules were enriched in the DNA protein and chromatin metabolic processes ( Figure 6B ). Analyses related to pain phenotypes per se (e.g., pain severity) were beyond the scope of this project. Yet we ran sensitivity analyses with pain severity as defined by the GCPS ( 21 ) to exclude the possibility that pain per se drove the results. Importantly, none of the identified genes and pathways described above regarding placebo effects were observed.

This study aimed to detect patterns in transcriptomic profiles associated with placebo effects in order to provide a genomic view of processes that might alter placebo effect. As reviewed in Introduction, candidate studies have provided clues that some functional polymorphisms involved in neurotransmitter function- including opioids, dopamine, and cannabinoids, modulate placebo effects. This is hardly surprising, but nor are the candidate gene findings definitive or representative of the many other pathways that may alter placebo and pain. Given that effects on a complex phenotype such as placebo effects are likely to be polygenic, and small, an important step in elucidating pathways to placebo effects can be gene network analysis. We utilized a model of placebo effects, which consisted of the reduction of pain in response to a placebo manipulation (e.g., sham electrode); this response depends upon the activation of descending neural pathways that can inhibit nociceptive signaling. Placebo effects tested in a laboratory experimental setting were employed to identify transcriptomic profiles associated with placebo responsivity. We identified 667 DEGs, and for 97% of them, the gene expression increased coordinately with increased placebo effects. Six out of the 17 top highly expressed genes showed significant correlations with larger placebo effects in both Discovery and Validation cohorts, and they included CCDC85B , FBXL15 , HAGHL , PI3 , SELENOM , and TNFRSF4 . SELENOM was identified as one of the most robustly upregulated genes in the Discovery and Validation cohorts; it encodes a selenocysteine-containing protein that is important for the central nervous system’s functions, particularly those related to memory, cognition, and motor coordination. A recent study in 10-month-old mice showed that the inaction of SELENOM reduced synaptic plasticity and led to memory problems ( 36 ). Additionally, SELENOM is essential for proper dopaminergic synaptic function ( 37 ), and as previously documented the dopaminergic system is involved in generating placebo effects ( 38 ). FBXL15 and HAGHL were also among the genes that were found to be upregulated in individuals with larger placebo effects in both the Discovery and Validation cohorts. FBXL15 plays a role in the G2/M transition of the cell cycle, as well as in cellular protein metabolic processes and the positive regulation of the Bone Morphogenetic Proteins (BMP) signaling pathway ( 39 ). HAGHL is a newly discovered metabolic oncogene and has been associated with the progression of human colorectal cancer ( 40 ). CCDC85B and TNFRSF4 were discovered and validated as genes associated with placebo effects. The expression of CCDC85B is regulated in a p53-dependent manner (p-53 has anti-oncogenic actions blocking β-catenin-dependent gene expression in the nucleus) ( 41 ). TNFRSF4 , known as OX40 or CD134 is expressed primarily on activated T cells and can activate the NF-kappa-B pathway by mediating TRAF2, TRAF5, PI3K/PKB, and NFAT pathways ( 42 ). The most important function of TNFRSF4 is to enhance cellular division, proliferation, survival, and cytokine production of T cells by activating the pathways described above, with a potential role in immunotherapy. The PI3 gene encodes an elastase-specific inhibitor peptide against Gram-positive and Gram-negative bacteria, and fungal pathogens, and its expression is upregulated by both bacterial lipopolysaccharides and cytokines ( 43 ). The enrichment analyses using DEGs, and network approaches, showed consistent significant findings. Unbiased enrichment analysis using all 667 DEGs and networks showed highly significant enrichment in several pathways involved the metabolism of RNA, regulation of expression of SLITs and ROBOs, and Huntington’s disease. Moreover, most of the genes enriched in the metabolism of RNA overlapped with SLIT/ROBO Signaling Pathway (R-HSA-9010553), which play important roles not only in the neural axon guidance but also in the angiogenesis and inflammatory cell chemotaxis processes ( 44 ). Our findings emphasized a potential role for the metabolism of RNA and ribonucleoproteins in the generation of placebo effects. Ribonucleoproteins refer to complexes of ribonucleic acid and RNA-binding proteins, while RNA metabolism refers to the processes of transcription, pre-mRNA splicing, editing, intracellular transport, translation, and degradation at the RNA level. These biological functions reflect a regulated process of the metabolism of RNA transcription, translation, and regulating gene expression. Disruption of these processes has been linked to neurological diseases ( 45 ). The disruption of the homeostasis of RNA-binding proteins leads to neurological disorders due to reduced RNA expression, and augmented propensity to RNA aggregation or sequestration. These mechanisms are essential to preserving neural integrity and when altered or disrupted, expose the brain to erroneous and deleterious changes in RNA expression. A susceptibility of neurons to these alterations affects the brain system role of RNA binding proteins in maintaining an optimal neuronal integrity ( 45 ). Therefore, it is plausible to think that alterations to RNA biological processes may affect the ability to respond to placebo procedures. Specifically, those with regular functioning of metabolism of RNA and riboprotein complex biogenesis may experience larger placebo effects which translate into better neuronal integrity. While this interpretation is novel and requires future research for generalizability, it is plausible to think that optimal metabolism of RNA and ribonucleoprotein complex biogenesis are relevant for optimal placebo responsiveness. These results require caution in being interpreted until future larger studies are conducted. Transcriptomics has been criticized as an approach not suitable to comprehensively identify genes involved in environmental adaptive responses such as temperature, salinity, pH, or oxygen changes ( 46 ). Conditioned placebo effects that we have explored in this study are not adaptive per se, but rather can be seen as compensatory physiological responses via learning processes. Thus, transcriptomics can be a potentially novel and reliable approach to identify genes associated with propensity to create placebo effects. Thus, we chose to focus on transcriptomics for several reasons. First, adequately powered transcriptomic studies can be conducted with fewer participants. In the pain field, we used this approach in a cohort of low back pain and HC participants ( 47 ). Furthermore, differentially expressed genes have also been identified in blood from participants with chronic visceral pain ( 48 ), fibromyalgia pain ( 49 ) and osteoarthritis pain ( 50 ). Finally, transcriptomic profiles in the blood likely represent pain changes due to the interaction of circulating immune cells with peripheral nerves that generate pain ( 15 ). In this study, we were not able to conduct RNA-seq analysis on Peripheral Blood Mononuclear Cells (PBMCs) to compare transcriptomic profiles of PBMCs and whole blood. However, the findings related to transcriptomic profiles of PBMCs as compared to whole blood are controversial with some studies reporting a higher abundance of gene expression when whole blood is used ( 51 ) as well as PBMCs ( 52 ). In this study, we did not include any pain-free healthy cohort limiting our possibility to determine gene expression related to TMD pain-related phenotypes. To our knowledge, this is one of the first studies to use transcriptomic profiling to identify genes and pathways associated with placebo effects in participants suffering from orofacial pain and other chronic overlapping pain conditions. Future omic studies, including genome-wide association, can determine whether these results generalize to clusters of clinical phenotypes. They can find convergences between pathways implicated by transcriptomics and genetics. Metabolism of RNA and ribonucleoprotein complex biogenesis might play a critical role in placebo responsiveness.

The biological basis of variation in placebo response to pain has been studied via both behavioral and/or physiological (i.e., neural, genetic) approaches and in several clinical contexts including chronic Irritable Bowel Syndrome ( 1 ), idiopathic and neuropathic pain ( 2 ), low back pain ( 3 ), migraine ( 4 ), knee osteoarthritis ( 5 ) and Temporomandibular Disorders (TMD) ( 6 ). The majority of these studies have suggested that there is a genetic/genomic component to pain sensitivity ( 7 ) and response to placebo manipulation ( 8 ) in both healthy controls and chronic pain patients. To date, few studies have explored the genetic factors associated with placebo effects, and notably these studies have not been omic in nature, but have focused on candidate genes involved in the dopaminergic, opioid, serotonin and endocannabinoid pathways ( 8 ).Pecina and colleagues first examined the OPRM1 A118G variant (Asn40Asp, rs1799971) which is thought to alter function of this critical receptor in pain( 9 ) and found the variant to be associated with placebo effects. The rs4680 polymorphism in COMT , a functional polymorphism altering cognition and emotion, was associated with better outcomes in patients with Irritable Bowel Syndrome ( 10 ) and placebo analgesia in healthy participants ( 11 ). An interaction between OPRM1 rs1799971 and COMT rs4680 determine larger placebo effects whereby participants having the COMT met/met or val/met – OPRM1 A/A carrier combination reported a level of placebo induced pain reduction that was 4–6 times higher compared to those with the val/val – G combination ( 12 ). The functional variant in the fatty acid amide hydrolase ( FAAH ) rs324420 polymorphism, encoding a Pro129Thr missense substitution, also affects placebo effects ( 13 ). We performed a replication of those findings ( 10 ), observing larger placebo effects in individuals with OPRM1 rs1799971 A/A, COMT rs4680 met/met and met/val and FAAH rs324420 Pro/Pro compared with other combinations ( 14 ). However, genomic views of placebo response are lacking, and given the polygenicity of complex traits such as placebo responses, it is likely that many genes, and other processes, play a role. In this study, we employed RNA-seq to quantify the transcriptome, commonly referred to as gene expression profiling. This approach has been used to identify molecular mechanisms of chronic pain and response to interventions due to the interaction of circulating immune cells with peripheral nerves that generate pain ( 15 ). Significantly differentially expressed genes (DEGs) were used in unbiased pathway analyses to identify potential new biological processes that might contribute to placebo effects in participants with a primary diagnosis of temporomandibular disorders (TMD) and other Chronic Overlapping Pain Conditions (COPCs) ( 16 ).