The study was approved by the institutional review board of the University of Illinois at Chicago. Tissues were obtained from the Quebec Suicide Brain Bank. The activities of the Quebec Suicide Brain Bank are approved by the Douglas Hospital McGill University IRB. Family members/informants signed written informed consents. For identification and dissection of neuroanatomical regions, we relied on experienced histopathologists using reference neuroanatomical maps [23]. Gyri and sulci were used to landmark specific frontal cortical areas. More precisely, the frontal lobes were sliced into 1–1.5 mm thick coronal sections. The dorsomedial PFC (Brodmann's area 9) was taken just dorsal to the frontopolar area including the most polar portion of the superior and partly the middle superior gyrus between the superior and intermediate frontal sulci. In the sections of the dissected cortical area, the gray and white matters were separated and the gray matter was used in this study. All the tissues were screened for evidence of neuropathology by experienced neuropathologists. pH of the brain was measured in cerebellum in all cases as described by Harrison et al. [24].

Psychiatric diagnoses of subjects were made by means of the psychological autopsy method. This technique, which has been validated for Axis I and II diagnoses [25], [26] consists in part of selecting a family member who is best acquainted with the deceased to serve as an informant and undergo the interview process. Informants included primarily mother, father, siblings, and significant other friends, or relatives. In our previous studies, the type of informant made no significant difference in rating the specific disorders identified [27]. Families are recruited at the Montreal Morgue and are interviewed, on average, 4 months after the death. The activities of the Quebec Suicide Brain Bank are approved by the Douglas Hospital McGill University IRB. Family members/informants sign written informed consents. Psychiatric diagnoses are obtained using the SCID [28] interview for axis I DSM-IV diagnoses, as well as the SCID II. Family members are requested to give written permission for clinical records to be obtained from prior mental health treatment providers. An attempt is made to collect all available records on each case, and then the appropriate data are extracted from the records and collected. Information collected through SCID I and II interviews and from the Coroner's notes and medical records is used by interviewers to write a case history for each subject. These case histories are then reviewed by a clinical panel, in order to reach a consensus DSM-IV diagnosis for each subject. All normal controls are verified as being free of mental illness by the same diagnostic process. Inter-rater reliability is routinely checked and periodic interviewer retraining is performed to avoid drifting. Toxicology screening for alcohol, a comprehensive battery for illicit drug use, and screening for antidepressant or psychoactive drugs taken prior to death were performed in blood/urine of each subject. Subjects were negative for all the toxicology screens. Characteristics of normal controls and suicide subjects are provided in Table 1. There were no significant differences in age, PMI, or pH of the brain between depressed suicide and normal control subjects.

Total RNA was isolated in samples of prefrontal cortex (Brodmann Area 9) using a modified protocol designed to optimize recovery of small RNAs [29]. Total RNA was isolated with Trizol reagent (Invitrogen Life Technologies, Carlsbad, CA, USA) according to manufacturer's directions with a few modifications to maximize yield of small RNAs. Glycoblue 20 µg (Ambion) was added to the RNA precipitation step which was allowed to proceed overnight at −20°C. The RNA pellet was spun down at 20,000 g for 25 min at 4°C; rinsed with 80% ethanol in DEPC-treated water (Invitrogen Life Technologies, Carlsbad, CA, USA); resuspended and treated with RNAsecure (Ambion); and treated with DNase I using DNA-free (Ambion). RNA was treated with DNAse I and checked for purity by OD 260/280 ratio (NanoDrop 1000 spectrophotometer, Thermoscientific, Wilmington, DE, USA). The quality of total RNA was verified for all subjects by running on agarose gels. Quality of total RNA was also checked by determining RNA integrity number (RIN) and only samples showing RIN>8.5 were included.

Reverse transcription was performed with the TaqMan® miRNA Reverse Transcription kit (ABI: Applied Biosystems) and the multiplex RT for TaqMan® miRNA Assays (ABI) that consists of eight predefined reverse transcription primer pools (ABI) following the manufacturer's protocol. For each RT pool, 100 ng of total RNA was used and the product was diluted 1∶62.5 and 55 µl mixed with 55 µl of TaqMan® Universal PCR Master Mix, No AmpErase® UNG. 100 µl of each mix was dispensed in the appropriate well in the TaqMan® Human miRNA Array v1.0 (TLDA, ABI) and run to 50 cycles as per manufacturer's protocol on a ABI 7900HT. An equal number of control and depressed suicide samples were processed and assayed on the same day, by an individual who was unaware of group identity. A sample processed without reverse transcription showed no detectable miRNA values.

Using 4 samples run on duplicate plates to monitor inter-plate reliability, we observed that Ct values above 35 were markedly less reliable, and so Ct = 35 was set as the threshold of detectability. Four data points were removed as outliers (defined as measurements that were at least 3 SD different than all other measurements of the same miRNA; in all cases, there was no more than one outlier observed for a single miRNA); however, removal of outliers did not affect the results reported here. Ct values were normalized using U6 RNA, whose mean value did not vary across groups and showed low variability. In the first 15 samples in each group, exogenous synthetic mir-122a was spiked into samples at the time of RNA extraction as a further control to check for possible differences in sample yield, integrity or loading; results were similar whether normalizing to the endogenous (U6) or exogenous (mir-122a) RNA. Note that Ct values cannot be compared directly across different miRNAs, but when comparing the same miRNA in different samples, a higher Ct value indicates lower abundance, and a difference in Ct values of 1 indicates a two-fold difference in abundance.

Tissues were homogenized on ice with a glass-Teflon homogenizer in lysis buffer (50 mM Tris–hydrogen chloride, pH 7.5, 100 mM sodium chloride, 1% Triton X-100, 50 mM sodium fluoride, 10 mM sodium phosphate, 5 mM EDTA, 1 µg/mL aprotinin, 1 µg/mL leupeptin, 1 mM sodium orthovanadate, and 1 mM 4-[2-aminoethyl] benzenesulfonyl fluoride hydrochloride). The homogenates were centrifuged at 14,000 g for 10 minutes. The protein content in the supernatant was determined by the Bradford method (Bio-Rad, Hercules, CA).

Samples containing equal amount of protein (30 µg) were resolved onto 10% (wt/v) sodium dodecyl sulfate (SDS)–polyacrylamide gel and blotted on an enhanced chemiluminescence (ECL) membrane (Amersham) as described earlier (30). Membranes were incubated with human-specific DNMT3b (LifeSpan Bio, Seattle, WA), BCL2 (Abcam, San Francisco, CA), or VEGFA (LifeSpan Bio, Seattle, WA) antibody overnight at 4°C. The dilution for each antibody was as follows: DNMT3b, 1∶500; BCL2, 1∶1000; VEGFA 1∶1000. The ECL membranes were then incubated with horseradish peroxidase–linked secondary antibody (anti-mouse or anti-rabbit IgG, 1∶1000) for 5 hours at room temperature and exposed to ECL autoradiography film. Membranes were stripped and reprobed with endogenous protein β-actin antibody (Sigma Chemical Co, St. Louis, MO) as described earlier [30]. The optical densities of the bands were quantified using the Loats Image Analysis System (Loats Associates, Inc, Westminster, MD). Ratios of optical densities of DNMT3b, BCL2, or VEGFA to β-actin were calculated. DNMT3b, BCL2, and VEGFA migrated to 96, 26, and 22 kDa, respectively. Immunolabeling of β-actin was not significantly different in depressed subjects (1.1±0.2 AU) compared with normal controls (1.2±0.2 AU).

The specificity of DNMT3b, BCL2, and VEGFA antisera were checked by using a 100-fold excess of blocking peptide (relative to the molarity of the antisera) corresponding to the epitope used to generate DNMT3b, BCL2, and VEGFA. In addition, to validate our data, we initially determined the immunolabeling of DNMT3b, BCL2, and VEGFA in the PFC of depressed suicide and control subjects using 5 different concentrations of protein (5–50 µg). It was observed that the optical density of the band increased linearly with increased concentration of protein and that the curve shifted toward the right when a decrease in immunolabeling was observed (data not shown).

The non-parametric Wilcoxon sign-rank test, 2-tailed was utilized when making lists of miRNAs whose mean expression levels differed significantly across treatment groups, since this test is appropriate for all miRNAs whether or not they follow a normal distribution [31]. The Bonferroni correction of statistical significance values was not appropriate in this study because it assumes that the expression of the vast majority of genes is independent of each other. This situation does not apply in the case of miRNAs, which form extensive cross-correlation networks –in some cases due to co-transcription from the same primary gene transcripts [29], [32], and in some cases due to shared factors that co-regulate their expression within a given tissue or even within a disease condition (see Results). Instead, we performed SAM analysis (Significance Analysis of Microarrays software, version 3.08, Stanford University, http://www-stat.stanford.edu/~tibs/SAM/),which estimated statistical significance by subjecting the data to multiple random permutations; parameters were set for two-class unpaired analysis, Wilcoxon statistic and 5000 permutations.

To measure possible correlations between miRNAs and age, postmortem interval and pH, the latter factors were converted to log2 scores to be comparable with Ct values, and all miRNAs were tested for correlations with each factor across all 35 samples. Given the total sample size (n = 35), a correlation of r = 0.333 is significantly different from r = 0 at the p = 0.05 significance level. The protein levels of DNMT3b, BCL2, and VEGFA were compared between depressed suicide and normal controls using independent sample t-test. The correlations between protein levels of DNMT3b, BCL2, and VEGFA and corresponding miRNAs were calculated using Pearson product-moment correlation analysis.

Despite active research in the miRNA field, it is currently not possible to predict accurately which individual brain-expressed mRNAs are most likely to be targeted by a given miRNA [33]. However, we attempted to glean further insights by combining evidence from 11 miRNA target prediction servers using the miRecords meta-server [34]. For each down-regulated miRNA, the meta-server was queried to identify targets that are predicted by at least 3 different prediction servers, and the top 30 ranked targets were tabulated (this is a small fraction of the total number of predictions). We then asked whether certain targets of potential interest were common to many different miRNAs in the down-regulated set of 37 miRNAs listed in Table 2.

Total RNA (1 µg) was reverse transcribed using 1.5 µl of a 10 µM mix solution of random hexanucleotides (experiment 1) or short specific RT primers (experiment 2), 10 mM dNTP mix, 0.1 M DTT, 40 units RNaseOUT and 200 units SuperScript™ III. Denaturation was performed at 90°C for 5 min. RT reaction was carried out by 5 min at 38°C, 30 min at 42°C, 30 min at 55°C, in a final reaction volume of 20 µl. Realtime PCR reaction was carried out in a final volume of 20 µl, containing 5 µl of using cDNA (diluted 1∶10), 20 µM each of specific primers, and 1× SYBR green PCR Master Mix (Applied Biosystems). To determine the linear range and sensitivity of each assay, a standard curve was generated using a serial 10-fold dilution of pooled cDNA derived from all subjects. All experiments are performed in duplicate and include a non-RT and no template control. Reactions were run on Mx3005p (Strategene, La Jolla, CA, USA) with cycling conditions: 95°C for 10 min, 40 cycles 95°C for 15 sec and 60°C for 1 min and followed by a melting curve cycle; 95°C for 1 min, 55°C for 30 sec, 95°C for 30 sec. Primer specificity was verified by a single peak in the melting curve analysis and a single band of the right size on gel electrophoresis for the all genes analyzed. Short RT primers: Dicer (TM 40.2 C) TTGGTGGACCAA; Drosha (TM 39 C) CGTCCAAATAACTG; DGCR8 (TM 43.6 C) ACGTCCACGGT; U6 (TM 36.2 C) TATGGAACGCTT; U44 (TM 37.2 C) CAGTTAGAGC TAATT; PRI-155 (TM 41.9 C) ACAGCCTACAGC; PRI-494 (TM 37.7 C) TCAAAACAGA ACTG. qPCR Primers: DICER (154 bp) Forward: CACATCAATAGATACTGTGCT, Reverse: TTGGTGGACCAACAATGGAGG. DROSHA (64 bp) Forward: AAGCGTTAATA GGAGCTGTTTACT, Reverse: CGTCCAAATAACTGCTTGGCT. DGCR8 (92 bp) Forward: GCTGAGGA AAGGGAGGAG, Reverse: ACGTCCACGGTGCACAG. U6 (95 bp) Forward: CTCGCTTCGGCAGCACA, Reverse: AACGCTTCACGAATTTGCGT. U44 (48 bp) Forward: TGATAGCAAATGCTGACTGA, Reverse: CAGTTAGAGCTA ATTAAGACC. Pri-155 (135 bp) Forward: CCAGCTTTATAACCGCATGT, Reverse: ACAGCCTACAGCAAGC CTTC. Pri-494 (134 bp) Forward: TGACTTCCCAA AAGCAACCT, Reverse: TCAAAACA GAACTGCCGAATC.

The multiplex RT-PCR plate assayed 367 miRNAs in parallel, of which 196 miRNAs exhibited mean expression values at least 2-fold greater than the threshold for detectability (i.e., mean Ct value ≤34 before normalization) in the control group. This set of 196 miRNAs was analyzed further. The complete set of normalized Ct values is included as Supporting Information (File S1).

First, miRNA expression values were normalized relative to the small RNA U6; this appears to be a suitable normalizer since U6 showed no significant mean difference across groups and exhibited low, equal variance in both groups. Twenty-one miRNAs were significantly down-regulated (none up-regulated) at p = 0.05. This is a large and consistent response; note that one would expect to observe only 5 miRNAs down-regulated and 5 miRNAs up-regulated by chance alone, at p = 0.05. Similarly, 24 miRNAs were down-regulated but none up-regulated by 30% or more (i.e., delta Ct>0.5). The results were even more striking when SAM analysis (Significance Analysis of Microarrays, Stanford University) was carried out: 100 miRNAs were identified as significantly decreased at p = 0.05 or better, with a false-discovery rate of 2.45% (i.e., fewer than 3 miRNAs would be expected to be significant by chance alone).

In large part, this down-regulation is due to a general, global down-regulation of miRNAs levels in the depressed suicide group. In almost every case, the mean abundance of a given miRNA in the depressed group was lower than its mean abundance in controls (17% decrease on average). This finding was highly significant (p = 4.5×10⁻⁴³). Furthermore, this global down-regulation was also apparent when the data were normalized instead to the geometric mean of U6, U44 and U48 (p = 8.3×10⁻²¹). Finally, to make sure that the apparent finding was not due to an up-regulation of small RNAs in the depressed group, we examined the 15 samples in each group that had mir-122a spiked-into samples as an exogenous control (mir-122a levels showed no difference across groups). We still observed a highly significant global downregulation of miRNA levels in the depressed group (p = 1×10⁻⁶).

Because the depressed group might be expected to be heterogeneous in terms of underlying etiologies or other variables, we hypothesized that the measured miRNA values in the depressed group might an overall increased extent of heterogeneity relative to the control group. However, by plotting the relation between the standard deviation and the mean expression value across all expressed miRNAs, we found that the miRNA expression in the depressed suicide group was significantly less variable than in the control group (Fig. 1).

Using data normalized to U6, twenty-one miRNAs were significantly decreased (none were increased) when tested individually at p = 0.05 or better (Table 2). The majority of these were also significant when the data were normalized instead to the geometric mean of U6, U44 and U48 (data not shown). To examine whether confounding variables had any effect on the expression of individual miRNA, we performed correlational analysis using Pearson product-moment analysis. We did not find any effect of age, PMI, gender or pH of the brain on the expression of miRNAs (data not shown).

The 21 down-regulated miRNAs include many miRNAs that have been implicated in cellular growth and differentiation. The set comprises both very abundant and low-abundance miRNAs, and those with high and low synaptic enrichment as studied in mouse homologues [29]; the miRNAs arise from both intronic and intergenic loci, and are located on numerous chromosomes. An additional 16 miRNAs were decreased by 30% or more that failed to meet the criterion for significance (Table 2). This may be due, at least in part, to their much lower expression – only 3 of 16 miRNAs in the non-significant set were abundantly expressed (i.e., mean Ct value<30), vs. 18 of 21 miRNAs in the significant set. Nevertheless, because many of the non-significant miRNAs shared chromosomal loci, 5′-seeds, and predicted targets with the significant miRNAs, they were also included in the analysis.

Almost half of the combined down-regulated miRNAs are encoded at chromosomal loci near another miRNA on the list (Table 2) and are possibly transcribed by the same primary miRNA gene transcripts [32], [35]: a) mir-142-5p and 142-3p; b) mir-494, 376a*, 496, and 369-3p; c) mir-23b, 27b and 24-1*; d) mir-34b* and 34c; and e) mir-17* and 20a. In addition, three pairs of miRNAs are encoded at distances greater than 100 kb but still lie within the same chromosomal region: a) mir-424 and 20b at Xq26.2-3, 377 kb apart; b) mir-142 and 301a at 17q22, 820 kb apart; and c) mir-324-5p and 497 at 17p13.1, 205 kb apart. This suggests that at least some of the down-regulated miRNA expression is due to decreased transcription.

Many of the down-regulated miRNAs also shared 5′-seed sequences (particularly bases 2–8, Table 2) that are involved in target recognition [36]. For example, identical seed sequences are shared by a) mir-20a and 20b; b) mir-301a and 130a; and c) mir-424 and 497. As well, a 6-mer nucleotide motif is shared by mir-34a, 34b* and 34c, and strikingly, a 5-mer motif (AGUGC) within the 5′-seed is shared by 5 of the affected miRNAs (mir-148b, 301a, 130a, 20a, 20b) that is predicted to bind Alu sequences within the 3′-UTR region of target mRNAs [37]. This suggests that the down-regulated miRNAs should exhibit extensive overlap among their mRNA targets.

At the time that the analysis was carried out, 42 experimentally validated targets were currently known for miRNAs in the down-regulated set, according to curated data compiled at the miRecords server [38], http://mirecords.umn.edu/miRecords/, though it should be noted that validation experiments were generally conducted in non-neural systems. Many of these are transcription factors (e.g., E2F1, E2F6, BACH1, SP1, HOXA5, RUNX1) and other nuclear proteins, but include transmembrane and signaling proteins as well. Intriguingly, 4 different down-regulated miRNAs all target the same validated growth factor, VEGFA (mir-20b, 20a, 34a, 34b*), a molecule implicated in depression in both humans and animal models (see Discussion). Other validated targets include BCL2 (mir-34a), DNMT3B (mir-148b), and MYCN (mir-101, 34a).

Among predicted targets, estrogen receptor alpha, ESR1, was predicted to be targeted by 3 different down-regulated miRNAs (mir-148b, 301a, 496). Others targeted by 3 or more affected miRNAs include ubiquitin ligases (UBE2D1 and UBE2W), signal transduction mediators (CAMK2G, AKAP1), the splicing factor NOVA1 that regulates brain-specific alternative splicing; the GABA-A receptor sub-unit GABRA4; calcium channel CACNA1C; and brain-active transcription factors including SMAD5, MITF, BACH2, MYCN, and ARID4A. Several of these predicted targets interact with validated targets; for example, ARIA4A binds E2F1; SMAD5 binds RUNX1; and estradiol treatment decreases E2F1 levels in prefrontal cortex [39]. BACH2 transcription factor binding sites have been identified upstream of many brain-expressed miRNAs [40]. Retinoblastoma binding protein 1 (ARIA4A) is of interest because it recruits histone deacetylases and regulates gene expression via chromatin-based silencing.

We measured protein levels of several selected target proteins implicated in depression (DMNT3b, VEGFA, and BCL2) in all individuals by Western blotting (normalizing values using beta-acting immunoreactivity). Whereas mean levels of DMNT3b was strongly up-regulated (2.4-fold) in the depressed suicide group (p = 1×10⁻¹¹), VEGFA showed no significant change (1.03-fold, p = 0.4) and BCL2 was strongly down-regulated (0.57 of control, p = 2.4×10⁻¹⁴) (Fig. 2). Changes in mean expression levels may reflect a variety of influences (e.g. transcription factor activity and differential turnover), as well as possible exerting regulatory effects by miRNAs listed in Table 2 and others. When we calculated the correlations between miRNA abundance and protein abundance across individual samples within the control group and (separately) within the depressed suicide group, surprisingly, DNMT3b levels showed an extremely strong positive correlation with mir-148b across individuals (r = 0.91 in controls, r = 0.94 in the depressed suicide group). Similarly, BCL2 was strongly and positively correlated with mir-34a (r = 0.92 in controls, r = 0.82 in the depressed suicide group). These findings strongly suggest that these miRNAs are likely to be actively co-regulated with their targets, as often observed in brain [7]. In contrast, VEGFA showed a mix of smaller positive and inverse correlations: (mir-20a: r = 0.08 in controls and r = 0.26 in depressed; mir-20b: r = 0.34 in controls and r = 0.24 in depressed; mir-34a: r = 0.34 in controls and r = −0.39 in depressed; mir-34b-5p: r = −0.363 in controls and r = 0.08 in depressed). Of these, mir-34a is probably the most interesting finding, since the correlation was positive in controls but inverse in the depressed suicide group. This may reflect a reorganization of miRNA-target networks (see below).

Quantile rank analysis is a non-parametric technique that identifies miRNAs which change their expression significantly relative to the other miRNAs in the population. This is an important parameter since different miRNAs often compete for binding to RISC complexes, and therefore changes in quantile rank may be as (or more) relevant for miRNA functional effects as absolute changes in miRNA abundance across groups. The miRNAs in each group were rank-ordered according to their mean expression levels (i.e., were assigned quantile ranks), and we then looked for miRNAs that show significant changes in quantile rank across groups. Overall, changes in quantile rank were relatively small and normally distributed (mean change = 0, SD = 3.6). Although six miRNAs showed a relative increase in quantile ranks relative to other miRNAs in the population (Table 3), this finding is not different from what would be expected by chance alone (expected number at p<0.05 or better = 5.3, observed = 6). In contrast, ten miRNAs showed a relative decrease in quantile ranks in excess of what is expected by chance (expected number at p<0.05 or better = 5.3, observed = 10) (Table 3). Most of these miRNAs were also listed in Table 2, and may be regarded as especially likely to exert functional effects in the depressed suicide group.

Another complementary method of analyzing the miRNA expression data is to identify pairs of miRNAs that are co-regulated in their expression, up or down, across individuals within a single group. The pairwise correlations for all miRNAs were computed separately within the control group and within the depression group. We were particularly interested in identifying miRNA pairs that showed no significant correlation in the control group, yet did exhibit a significant positive correlation in the depressed group. This phenomenon allows one to detect coordinated changes in the underlying driving forces for expression (e.g., due to changes in transcription factors or miRNA processing) in the depressed suicide group.

In this analysis, the list of miRNAs was first filtered to include only those miRNAs that had robust expression (i.e., whose Ct value <35 in all individuals in both control and depressed groups). The Ct values for each miRNA were normalized to the global mean Ct value of the same individual (averaged over all miRNAs), in order to remove correlations due to inter-individual differences in overall miRNA content. Next, all miRNAs were examined pair-wise and the Pearson correlation coefficient was computed for each group separately. (Note that for n = 17 subjects per group, only correlations of r = 0.5 or greater are significantly different from 0 at p = 0.05.) Finally, we identified pairs of miRNAs that satisfied the following criteria: 1. They exhibited a significant positive correlation in the depressed group (r-dep>0.5), but no significant correlation (positive or negative) in the control group (−0.5 <r-control≤0.2). 2. The correlation coefficients in the control group were much lower than in the depressed suicide group (r-dep minus r-control >0.75).

This procedure revealed a set of 29 miRNAs, none of whom were pairwise correlated in the normal control group, but which formed a very extensive inter-connected network in the depressed group (Fig. 3a). Several of the miRNAs (let-7b, mir-132, 181b, 338-3p, 486-5p, and 650) were “hubs” correlated with four to nine other miRNAs in the network. Only two of the miRNAs were significantly down-regulated in the depressed group (mir-142-5p and 146a, neither of them hubs). Of the network miRNAs, mir-181b was reported to be upregulated in schizophrenia [14], and mir-132 and 134 are known to be CREB-regulated and involved in neurite growth and dendritic spine growth, respectively [41], [42]. Thus, correlation analysis revealed the existence of factors that drive co-expression of a group of miRNAs across individuals specifically in the depressed group (Fig. 3a). A similar analysis also revealed a smaller set of miRNAs that were pairwise correlated in the control group but NOT in the depressed suicide group (Fig. 3b).

Three genes involved in miRNA biogenesis (dicer, drosha and DGCR8) were measured in the entire cohort using RT-PCR, using the same total RNA employed for the miRNA plate measurements (Table 4). No significant difference was observed between the control and depressed groups for any of these mRNAs (Table 4). Two of the significantly down-regulated miRNAs were also examined, that are encoded by single primary transcripts at intergenic loci. Pri-miR-155 levels showed a 19% decrease in the depressed group, but the finding was not statistically significant (p = 0.30; Table 4). Pri-miR-494 showed a 18% increase in the depressed group, but this did not achieve significance either (p = 0.088).