Responder and nonresponder patients exhibit different peripheral transcriptional signatures during major depressive episode

To date, it remains impossible to guarantee that short-term treatment given to a patient suffering from a major depressive episode (MDE) will improve long-term efficacy. Objective biological measurements and biomarkers that could help in predicting the clinical evolution of MDE are still warranted. To better understand the reason nearly half of MDE patients respond poorly to current antidepressive treatments, we examined the gene expression profile of peripheral blood samples collected from 16 severe MDE patients and 13 matched controls. Using a naturalistic and longitudinal design, we ascertained mRNA and microRNA (miRNA) expression at baseline, 2 and 8 weeks later. On a genome-wide scale, we detected transcripts with roles in various biological processes as significantly dysregulated between MDE patients and controls, notably those involved in nucleotide binding and chromatin assembly. We also established putative interactions between dysregulated mRNAs and miRNAs that may contribute to MDE physiopathology. We selected a set of mRNA candidates for quantitative reverse transcriptase PCR (RT-qPCR) to validate that the transcriptional signatures observed in responders is different from nonresponders. Furthermore, we identified a combination of four mRNAs (PPT1, TNF, IL1B and HIST1H1E) that could be predictive of treatment response. Altogether, these results highlight the importance of studies investigating the tight relationship between peripheral transcriptional changes and the dynamic clinical progression of MDE patients to provide biomarkers of MDE evolution and prognosis.


Design setting and subjects
The study design was a naturalistic, prospective, longitudinal and comparative study with assessments of major depressive episode (MDE) patients and healthy controls at baseline (week 0), 2 and 8 weeks after inclusion. Sixteen patients who met the Diagnostic and Statistical Manual of Mental Disorders, fourth edition, text revision (DSM-IV-TR) criteria for major depressive disorder (MDD) participated in the study. 1 Inclusion criteria were: (i) treated or untreated MDE (ii) 17-item Hamilton Rating Scale for Depression (HDRS) score ≥ 20 corresponding to severe or very severe MDE. 2 Patients with bipolar disorder, schizophrenia, personality disorder, substance use disorder, medical problems, abnormal body temperature or laboratory tests (especially blood cell count), and neurological exploration were excluded.
Patients were clinically evaluated using an adapted standardized procedure (SCID-I/P) 3 .
Socio-demographic information including age, gender, marital status and educational level was collected. The severity of MDE was evaluated with a French version of 17-item HDRS at each evaluation. 4 Presence of recurrent MDD, episode duration, familial history of MDE and other severe psychiatric disorders, suicide attempts history, and tobacco smocking were evaluated.
At the end of the 8-week clinical follow-up, patients were classified as responders and nonresponders based on the consensual definition of clinical response corresponding to a minimal reduction of 50 % of the HDRS score from the initial evaluation. 5 The first nine patients, classified as responders, and their matched controls were organized by pairs within cohort A. The other patient/control pairs formed cohort B.
For the control group, age-and sex-matched subjects were evaluated to exclude any patients with a history of psychiatric disorder using the French version of standardized interview validated for health control subjects (SCID-NP). Moreover, all participants were carefully interviewed and had clinical examinations to eliminate any psychiatric or medical conditions. Tobacco smoking was also evaluated.
All experiments on human subjects were conducted in accordance with the latest version of the Declaration of Helsinki. The project was approved by the local ethics committee (Comité de Protection des Personnes, CPP Sud Méditerranée II, Marseille, France, study registered under number 09.025) and written informed consent was obtained after a complete description of the study to the subjects.

Blood mRNA extraction
At each evaluation, 8-10 ml venous blood was collected from fasting patients and matchedcontrols in EDTA tubes between 7:00 a.m. and 9:00 a.m. and processed within two hours.
Peripheral blood mononuclear cells (PBMCs) were isolated from the blood by Ficoll density centrifugation. Total RNAs were extracted from the PBMCs with the mirVana miRNA isolation kit (Ambion, Austin, TX) according to the manufacturer's protocol. RNA concentration was determined using a nanodrop ND-1000 spectrophotometer (NanoDrop Technologies, Wilmington, DE). RNA integrity was assessed on an Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA). To ensure a robust analysis for the following procedures, samples with an RNA integrity number (RIN) inferior to 8 were excluded.

Preparation of samples and micorarray assay
Sample amplification, labeling, and hybridization essentially followed the one-color microarray-based gene expression analysis (low RNA input linear amplification PLUS kit) recommended by Agilent Technologies. In brief, 200 ng of each total RNA sample was reverse transcribed into cDNA using oligo dT-T7 promoter primer. Labeled cRNA was prepared from the cDNA. The reaction was performed in a solution containing dNTP mix, cyanine 3-dCTP, and T7 RNA polymerase, and incubated at 40°C for 2h. It was followed by RNeasy column purification (QIAGEN, Valencia, CA). Dye incorporation and cRNA yield were checked with the NanoDrop ND-1000 Spectrophotometer. 600 ng of cyanine 3-labeled cRNA (specific activity > 6 pmol Cy3/ug cRNA) was fragmented at 60°C for 30 minutes before hybridization, for 17 hours at 65°C in a rotating Agilent hybridization oven, onto Agilent whole human genome oligo microarrays containing 50 599 different oligonucleotide probes (SurePrint G3 Human GE 8x60K, Agilent Technologies, Santa Clara, CA) and covering 27 958 Entrez Gene RNAs and 7 419 lincRNAs. Microarrays were then washed according to manufacturer's instructions. Slides were scanned immediately after washing on the Agilent DNA microarray scanner (G2505C) using one color scan setting for 8x60K array slides. The scanned images were analyzed with Agilent feature extraction software 10.5.1.1 using default parameters (protocol GE1_107_Sep09 and Grid 028004_D_F_20110325) to obtain background subtracted and spatially detrended processed signal intensities. All data were normalized by quantile normalization using limma R/bioconductor package (v.2.16.4).
Only 39 784 oligonucleotide probes with signal intensities detectable (i.e. above background according to the "gIsWellAboveBg" Agilent feature extraction value) in ≥ 70% of samples (from either MDE patients or controls and at either baseline or 8 weeks later) were subsequently analyzed. This filtering step was obtained using the AgiND R package (http://tagc.univ-mrs.fr/tagc/index.php/software/agind, v1.0.5). The microarray data are available from the gene expression omnibus (GEO, http://www.ncbi.nlm.nih.gov/geo/) under the series accession number GSE38206.

Gene and chromosome band enrichment analyses
Lists of genes with significant changes were uploaded on DAVID (database for annotation, visualization and integrated discovery) for identifying statistically relevant signaling pathways 6 using high classification stringency, P-value (p) < 0.05 and false discovery rate (FDR) ≤ 5%. The sequences for the boundaries of the chromosomal regions where the probes to detect differential expression are located were obtained from the Ensembl database (Human reference assembly GRCh37, WTSI/EBI). To determine the chromosomal cytoband enrichment in our lists of genes, we used the ToppFun algorithm and restricted the search to p<0.05. 7 To take into account multiple testings, the FDR was calculated with the PROC MULTTEST statement as previously described. 400  were then normalized based on the expression level of a combination of 2 to 5 reference miRNAs, exhibiting a stable level of expression among all samples, and selected with both the genorm PLUS (Biogazelle, Zulte, Belgium) 8 and DataAssist softwares. In order to normalize the target miRNA with reference miRNAs exhibiting a proximal expression level (i.e. difference between the target and the reference miRNAs < 3 Ct), 5 windows of expression levels (Ct < 21; 21 < Ct < 24; 24 < Ct < 27; 27 < Ct < 30; Ct > 30) were defined in which a specific set of reference miRNAs were selected. Normalized expression levels of each miRNA are also quantified as 2 -Ct , 9 with the DataAssist software (Applied Biosystems, v3.0), relative to the normalized expression level of the same miRNA in a calibrator sample (either the mean of 9 samples from the 9 control subjects when comparing patients to controls, or the week 0 sample of the MDE patient when comparing patient at week 0 to the same patient 8 weeks later).

Target prediction
For each list of differentially expressed miRNAs, the Ingenuity pathway analysis (IPA) software (Ingenuity systems), which relies on 3 popular algorithms (TargetSan, TarBase and miRecords), was queried to identify targets within lists of differentially expressed genes in our microarray analysis. Only highly predicted and experimentally validated targets were considered for further analysis.

Real-time RT-PCR for candidate gene validation
1.6 μg of RNA was reverse transcribed with the High Capacity cDNA archive kit (Applied Biosystems, Foster City, CA, USA). 200 ng of the resulting cDNA was combined with a TaqMan® universal PCR Master Mix (Applied Biosystems) and 48 PCR reactions were simultaneously performed on an ABI PRISM 7900HT thermocycler using tLDA technology according to manufacturer recommendations (Applied Biosystems). All PCR reactions were performed in triplicate. For each tested candidate gene, primer sets and probes were selected using the web portal of the manufacturer (Applied Biosystems) to amplify sequences matching as close as possible the transcripts detected by the Agilent probe in our initial microarray screen. In addition, after a search within our microarray data for probes showing minimal variation of signal intensity among the 36 samples tested from both MDE patients and controls, we selected 3 genes, in addition to very highly expressed gene GAPDH (imposed by the manufacturer), to be used as highly (PAFAH1B1), moderately (ALDOC) and weakly (SV2A) expressed reference genes. After verifying with DataAssist that our 4 selected reference genes were stably expressed among all the samples from MDE patients and controls tested by RT-qPCR, we set-up 4 windows of expression intensities to normalize target gene expression (21<Ct<24; 24<Ct<27; 27<Ct<30; Ct>30) using adequate reference genes alone or in combination (Supplementary Table 17). The expression level of each candidate gene was calculated as 2 -Ct with the DataAssist software (Applied Biosystems, v3.0), in which each candidate gene is quantified relative to the expression of a single, or a combination of two reference genes, to calculate a proximal level (i.e. difference between the target and the reference miRNA < 3 Ct) of expression compared to the target gene, and compared to a calibrator sample (the mean of the samples from the control subjects at week 0 when comparing patients to controls).

Demographic and clinical data
Demographic and clinical variables were compared between patients and matched-controls for cohort A and between responder and nonresponder patients with a Fisher exact test for qualitative variables and a t-test for quantitative variables. All statistical tests were made using the IBM SPSS Statistics v20 software.

Microarray data
To calculate the number of samples required to draw conclusions with enough statistical power for our microarray screen, we determined a minimum sample size (n=9) by using sizepower tool within the R/bioconductor package 10 with the following parameters: FDR = 5 %, number of significant genes (π 0 ) ≥ 0.5%, power of the study = 90%, expected difference between groups = 1.5, and standard deviation = 0.5 (Supplementary Figure 1). Accordingly, samples of the first 9 responder MDE patients and matched controls at the inclusion visit and eight weeks later were assessed in the microarray analysis and constituted 4 different groups: "MDE 0w " for responder patients at the inclusion, "MDE 8w " for responders patients eight weeks later, "C 0w " for matched control subjects at the inclusion and "C 8w " for matched control subjects eight weeks later.

MiRNA study
To select differentially expressed miRNAs (at 0 and 8-week visits), non-parametric unpaired Mann-Whitney tests were used to compare the fold change (FC) between patients and matched-controls with a threshold P-value of 0.05. Non-parametric paired Wilcoxon tests were used to compare the FC between the two evaluations within the MDE group with a threshold P-value of 0.05.

Candidate gene validation
For exploratory purposes and to compare microarray data with RT-qPCR data, parametric unpaired t-tests were used to compare the FC between patients and matched-controls at the first evaluation, and 8 weeks later, within either the whole cohort (FC all ), the intial sub-cohort A (FC A ), or the subsequent sub-cohort B (FC B ).
Since we conducted a prospective study with three visits and have repeated measures for each subject, a mixed linear model method was used for each candidate mRNA to test the differences between responders and nonresponders patients. As potential confounding factors could introduce bias in the observed results, additional characteristics were taken into account and tested as co-variables in the mixed linear model such as sex, age, and treatments received.
To explore differences between responders and nonresponders, we calculated and compared marginal estimated mean of each group.
Finally, to quantitatively explore the predictive value of the best mRNA candidates in the previous analysis, we calculated Spearman's correlations between HDRS evolution (0 to 8-and 2 to 8-week intervals) and FC at inclusion or at 2 weeks, respectively. Moreover, a discriminant function analysis was used with the goal of establishing a predictive score to classify patients in the responder versus nonresponder groups. Sensitivity, specificity, positive/negative predictive value, and the confidence intervals of each selected mRNA candidates, or a combination of the best classifiers, were computed for the combination of selected variables. Receiver-operating characteristic (ROC) curves analysis was used to determine the area under the curve (AUC). This list contains 31 genes differentially overexpressed in MDE at both baseline and at 8-weeks versus control. Changes were assessed by a SAM analysis set with a FDR of 5% and then challenged with a t-test (FC ≥ 1.50; p<0.001). This list contains 39 genes differentially underexpressed in MDE at baseline and at 8-weeks versus control. Changes were assessed by a SAM analysis set with a FDR of 5% and then challenged with a t-test (FC ≥ 1.50; p<0.001). Genes in bold are also significant with a FDR set at 1%.

Supplementary Figure S1: Sample size calculation
To evaluate the number of samples required for the microarray screen, a minimum sample size was determined by using sizepower tool within the R package with the following parameters: FDR = 5 %, power of the study = 90%, expected difference between groups = 1.5, and standard deviation = 0.5. Three curves were drawn according to different values for the percentage of significant genes (1-π 0 ). The dashed line indicate that with a power of 90% and π 0 = 0.095, the minimum sample size is 9.

45
Supplementary Figure S2: Spearman's correlation between HDRS score evolution and mRNA expression level for PPT1 and CELF2 Each circle indicates, for a specified MDE patient (the label attached to the circle refers to the case number in Table 1), the difference of HDRS scores at week 2 and 8 weeks after inclusion as a function of the PPT1 (a) or CELF2 FC (b) for MDE 2w sample, calculated with the 2 -ΔΔCt method.
The calibrator was a mathematical pool of control sample Ct. Normalization was performed using GAPDH and PAFAH1B1 for PPT1, and PAFAH1B1 and ALDOC for CELF2.