The effect of listening to music on human transcriptome

Ethics statement

The Ethics Committee of Helsinki University Central Hospital approved this study. Written informed consent was obtained from all the participants.

Participants and phenotypes

A total of 48 individuals (aged 18–73; mean 42.5) participated in the study (Table S1). All the 48 participants were characterized for musical aptitude and music education. Musical aptitude was measured using three tests: the auditory structuring ability test (Karma music test) (Karma, 2007), and the Seashore tests for pitch and for time (Seashore, Lewis & Saetveit, 1960). Details of the tests have been described in Oikkonen et al. (2014). A combined score of the three tests (COMB score; range 0–150) was used to define musical aptitude. COMB scores were classified as either high or low based on the upper and lower quartiles of score distribution. Data on music education were collected using a questionnaire. The self-reported levels of music education (referred to as edu classes 1–4), received through studies or degrees from music schools, institutes, or universities, were classified as follows: 1 represents no music education; 2 represents music education of less than two years; 3 represents music education of more than ten years; and 4 represents a professional musician. The approximate time of systematic music education and training was 21.42 years on average, in participants of edu classes 3–4. Table 1 shows the sample’s phenotype characteristics.

Exposure to music

To our knowledge, no previous studies have systematically studied the effect of listening to music on genome-wide transcriptional profiles of humans. We have previously shown that music-listening habits vary a lot among listeners (Ukkola-Vuoti et al., 2011). To start with, we chose to study the effect of classical music, Wolfgang Amadeus Mozart’s Violin Concerto No. 3 in G major, K.216 because it is relatively familiar in the western culture. As the human brain perceives complex sounds in a millisecond-level time frame (Wang et al., 2009; Kayser, Logothetis & Panzeri, 2010), we expected that the 20 min-listening session (duration that the concerto lasts) will induce an effect of music on human transcriptome. In studies on the effect of pain in humans, very short durations of pain induction (8 to 150 s) have been used (Hubbard et al., 2011). The participants were unaware of the type of music that was intended for the listening session. Peripheral blood samples were collected from all the 48 participants just before and after the listening session. From here on, participants who listened to music are referred to as listeners throughout the text.

Control study

The same 48 participants were invited to a control study. Of these, 15 participants could attend. The participants were advised to avoid listening to music and hard exercise the day before the control study. The control study was performed in a ‘music-free’ environment, where the participants had an opportunity to converse, read a magazine, or take a walk outside (no exercise) during the session. Peripheral blood samples were collected from the participants just before and after 20 min in the control session (the same duration as in the listening session). From here on, participants of the control session are referred to as controls throughout the text.

Genome-wide expression profiling

For this procedure, 2 × 2.5 ml samples of peripheral blood were drawn into PAXgene blood RNA tubes (PreAnalytiX GmbH, Hombrechtikon, Switzerland) as per the kit instructions, in both of the sessions. Total RNA was isolated using the PAXgene blood miRNA Kit (PreAnalytiX GmbH, Hombrechtikon, Switzerland) as per the kit manual. Purified RNA samples were measured for concentration and purity on the NanoDrop 1000 v.3.7 (Thermo Fisher Scientific, Waltham, Massachusetts, USA). Globin mRNA was depleted from our samples using Ambion’s Human GLOBINclear™ kit (Applied Biosystems, Carlsbad, California, USA) as per the kit insert. The samples were measured on the NanoDrop 1000 to determine the sample concentration and purity and for integrity on the 2100 Bioanalyzer (Agilent Technologies, Waldbronn, Germany) before being diluted to 50 ng/µl using RNase-free water. A total of 2 µg of RNA was assayed on the Illumina HumanHT-12 v4 bead array (Illumina Inc., San Diego, California, USA), which targets more than 47,000 probes. The gene expression profiling assays for the listening and control sessions were conducted in two separate batches. To account for the batch effect corrections, six samples from the listening session were assayed together with the control session samples. Intensity data were exported through Bead Studio software. The data reported in this article have been deposited in the Gene Expression Omnibus database, www.ncbi.nlm.nih.gov/geo (accession no. GSE48624).

Data preprocessing

The Lumi package was used to read and preprocess the signal intensity data. Specifically, pre-processing included background correction, variance stabilizing transformation, and quantile normalization. Data from both the listening (n = 48) and control (n = 15) sessions were normalized separately. Five samples from the control session were excluded from further analyses owing to data quality. In addition, we used the ComBat method (Johnson, Li & Rabinovic, 2007) to adjust for batch effects and determine if this correction affected the pre-post fold-changes across conditions. However, we did not find significant differences between the fold-changes of corrected and uncorrected data over time across conditions. Therefore, we chose to proceed with the uncorrected data owing to the strengths of our analysis methods as described below. After normalization, duplicate and un-annotated probes were excluded using the genefilter package (R package version 1.40.0). Before extracting the expression values from the normalized data, Illumina’s detection p-values (threshold: 0.01) were used to filter out probes with low intensities corresponding to the background signal. Finally, only those probes that were expressed in at least half of all of the arrays (listening and control sessions) were chosen for the study.

Differential expression analysis

The choice of an appropriate statistical test for the identification of differentially expressed genes depends upon several aspects of the data including the underlying distribution, homogeneity, and the sample size. As the statistical tests for normality are known to be sensitive to sample size, we used a normal Q-Q plot to get a glimpse of the distribution of the data. For this, we randomly visualized the distribution of transcriptional responses of control samples (n = 10) for several transcripts (Fig. S1) using normal Q-Q plots. We observed that the data appeared to deviate from normality in several instances. As the central limit theorem does not always hold true for small sample sizes, a cautious approach here would be to employ a non-parametric test (better being safe than sorry). Non-parametric tests are deemed to be appropriate analysis tools when the distribution of data is difficult to characterize, because they make less stringent distributional assumptions. Therefore, we chose to use the rank product (Breitling et al., 2004) non-parametric test statistic, which is relatively powerful especially for small sample sizes, and when the data is heterogeneous and does not meet normality. Rank product-based methods outperformed several other methods including empirical Bayes statistic (limma) and SAM, when the sample size is small and when the data is non-homogeneous (Jeffery, Higgins & Culhane, 2006). Also, a comparison of eight gene ranking methods using Microarray Quality Control datasets (golden-standard) has demonstrated the high sensitivity and specificity of the rank product method (Kadota & Shimizu, 2011).

To identify the differentially expressed genes, we compared the magnitude of pre-post changes in gene expression across conditions using the rank product method implemented in the RankProd Bioconductor package. Based on the estimated percentage of false predictions (pfp), RankProd employs a non-parametric statistic to identify genes that are consistently ranked high among the most up- or down-regulated genes in replicate experiments. Instead of analyzing the actual expression value, this method utilizes the ranks of genes in each sample. The strength of rank product method allows us to compare and combine the datasets of the listening and control studies. After the identification of differentially expressed genes using a pfp of 0.05 in RankProd, we selected only those genes that exceeded an effect-size cut-off (>1.2-fold change over time across conditions, and at least a pre-post change of 10% in gene expression in the listening session). Here, we would like to point out a couple of aspects of these selection criteria. First, pfp employed by RankProd is equivalent to the standard false discovery rate (FDR). Second, there exists a widespread misconception that only two-fold changes are significant (Hoheisel, 2006) and that false notion is based on the very initial publications of microarray studies, which used a two-fold change criteria for a particular group of experiments because of biological relevance. Fold-change thresholds are completely arbitrary and in the majority of the cases they depend upon the underlying biological question. For example, studies that investigated the effect of gene-environment interactions (socio-environmental effect (Cole et al., 2007), yogic meditation effect (Black et al., 2012)) used unorthodox fold-change thresholds.

We further performed successive functional annotation analyses using GeneTrail (http://genetrail.bioinf.uni-sb.de/) and IPA (Ingenuity^® Systems, www.ingenuity.com). The gene ontological overrepresentation analysis in GeneTrail uses a hypergeometric distribution test along with a conservative multiple testing correction method (FDR < 0.05) to assess whether genes belonging to certain functional categories are enriched in the dataset. In addition to the standard gene ontology analyses, we performed upstream transcription regulator analysis, which essentially predicts all the upstream transcription regulators (transcription factors, receptors, cytokines, microRNA, and kinases) that could have possibly mediated the observed differential expression. Based on the overlap between known targets of a transcription regulator and the set of differentially expressed genes, an overlap p-value is computed using Fisher’s exact test (p < 0.01). Further, the activation states of the predicted transcription regulators are also inferred using an activation Z-score, which is based on literature-derived knowledge on the direction of regulation (either activating or inhibiting).

Transcriptional responses after a control session

First we assessed the homogeneity of transcriptional responses in the control session. For this, we analyzed whether the transcriptional responses of participants belonging to edu classes 3–4 (n = 4) and edu classes 1–2 (n = 6) differed in the control session. Only 4 genes were found to be differentially expressed, suggesting no major differences in the transcriptional responses between groups to a ‘non-musical activity.’ Therefore, we used the expression data of the whole ‘music-free’ control group as a reference for the comparative analyses.

We used Spearman’s rank-based correlation to check for each gene, whether the transcriptional responses correlated with the age or sex of the participants in either the listening session (n = 48) or the control sessions (n = 10). After multiple testing corrections, we found no significant effects of age or sex on the transcriptional responses.

Transcriptional response after listening to music

Based on neuroscientific studies, the brains of musicians and non-musicians differ structurally and functionally (Elbert et al., 1995; Gaser & Schlaug, 2003). This led us to ask whether the transcriptional responses of musically experienced participants would differ from those of musically inexperienced participants when listening to music. Therefore, we compared the transcriptional responses of listening to music separately for musically experienced and inexperienced participants vs controls. Comparisons were made based on two subphenotypes of musical experience: musical aptitude and music education.

First, we compared the magnitude of pre-post fold-changes in the genome-wide transcriptional profiles of listeners of edu classes 3–4 (n = 29) and controls (n = 10). Using RankProd non-parametric statistics and stringent selection criteria, we identified 45 differentially expressed genes (27 up-regulated and 18 down-regulated). Next, we compared the genome-wide transcriptional profiles of listeners with high COMB scores (n = 12) and controls (n = 10). Similar statistical analysis identified 97 differentially expressed genes (75 up-regulated and 22 down-regulated). The differentially expressed genes from both the comparisons are listed in Table S2, and a comparison of the pre-post changes in both conditions is shown in Fig. 1.