Circulating Serum MicroRNAs as Potential Diagnostic Biomarkers of Posttraumatic Stress Disorder: A Pilot Study

Posttraumatic stress disorder (PTSD) is a psychiatric disorder that can develop upon exposure to a traumatic event. While most people are able to recover promptly, others are at increased risk of developing PTSD. However, the exact underlying biological mechanisms of differential susceptibility are unknown. Identifying biomarkers of PTSD could assist in its diagnosis and facilitate treatment planning. Here, we identified serum microRNAs (miRNAs) of subjects that underwent a traumatic event and aimed to assess their potential to serve as diagnostic biomarkers of PTSD. Next-generation sequencing was performed to examine circulating miRNA profiles of 24 members belonging to the Dutch military cohort Prospective Research in Stress-Related Military Operations (PRISMO). Three groups were selected: “susceptible” subjects who developed PTSD after combat exposure, “resilient” subjects without PTSD, and nonexposed control subjects (N = 8 per group). Differential expression analysis revealed 22 differentially expressed miRNAs in PTSD subjects compared to controls and 1 in PTSD subjects compared to resilient individuals (after multiple testing correction and a log2 fold-change cutoff of ≥|1|). Weighted Gene Coexpression Network Analysis (WGCNA) identified a module of coexpressed miRNAs which could distinguish between the three groups. In addition, receiver operating characteristic curve analyses suggest that the miRNAs with the highest module memberships could have a strong diagnostic accuracy as reflected by high areas under the curves. Overall, the results of our pilot study suggest that serum miRNAs could potentially serve as diagnostic biomarkers of PTSD, both individually or grouped within a cluster of coexpressed miRNAs. Larger studies are now needed to validate and build upon these preliminary findings.

Background

Posttraumatic stress disorder (PTSD) is a psychiatric disorder that can develop upon exposure to a life-threatening traumatic event, i.e., an event capable of producing intense feelings of fear, helplessness, and horror (Association, 2013). Symptoms associated with PTSD include re-experiencing of the traumatic event, avoidance behavior, overall negative mood, and hyperarousal (Association, 2013). The economic burden associated with PTSD is substantial, and patients with PTSD are at increased risk of committing suicide and having familial issues such as marital problems (Fontana and Rosenheck, 1994; Tarrier and Gregg, 2004; Nock and Kessler, 2006; Ferry et al., 2015). Although ∼60% of individuals within Western Europe will one day be exposed to a traumatic event, only ∼6% of these individuals develop PTSD while others show a positive psychological adaptation process denoted as resilience (Kalisch et al., 2017; Koenen et al., 2017). However, some populations such as military soldiers are at elevated risk for trauma exposure, making PTSD a relatively common chronic disorder in the combat Veteran population (Thomas et al., 2017). Currently, a variety of treatment options exist for PTSD, without one being clearly superior to another (Yehuda et al., 2014). Moreover, pharmacological treatment options for PTSD are at best moderately effective and only work for a subset of patients (Richter-Levin et al., 2018). Therefore, increasing efforts are being made to unravel the biological underpinnings of PTSD in order to develop more efficient therapeutic strategies. It is now becoming clear that epigenetic mechanisms are involved in the lasting behavioral and molecular effects of trauma exposure (Schmidt et al., 2011; Snijders et al., 2018a).

Epigenetics refers to a variety of processes that are triggered by environmental factors and cause lasting but reversible alterations in gene expression (Goldberg et al., 2007). Among epigenetic mechanisms, noncoding RNA molecules such as microRNAs (miRNAs) are involved in the posttranscriptional regulation of gene expression by binding to specific messenger RNAs (Peschansky and Wahlestedt, 2014). Several miRNAs have been found implicated in PTSD, shedding much needed light on the underlying pathophysiological underpinnings of this disorder (Wingo et al., 2015; Bam et al., 2016a, Bam et al., 2016b; Martin et al., 2017). Such findings emphasize the notion that expression profiles of miRNAs could one day serve as relatively easily accessible biomarkers or be embedded within a network of several relevant biological processes that together could more accurately reflect the complexity of PTSD. For those individuals who have difficulties recognizing or properly describing their symptoms, identifying such markers could be of use in clinical contexts in order to objectively confirm the presence of the disorder and establish appropriate treatment plans when needed (Lehrner and Yehuda, 2014). Using these markers could be equally relevant during postdeployment medical screenings since military service members may have secondary reasons to not fully disclose their symptoms (Yehuda et al., 2013).

Here, we aimed to identify serum miRNAs that could one day serve as diagnostic biomarkers of PTSD. We further aimed to gain insights in the coexpression patterns of these miRNAs, their predicted gene targets and underlying biological pathways, along with their diagnostic accuracy. We hypothesized that specific miRNAs are differentially expressed between subjects with PTSD, trauma-exposed healthy individuals (referred to as “resilient” subjects in this paper), and nonexposed healthy controls. For this, we performed next-generation sequencing (NGS) on serum samples of 24 military members belonging to a Dutch military cohort, and we compared miRNA profiles between the three groups. Our findings suggest that miRNAs could potentially serve as biomarkers of PTSD, both individually or grouped within a cluster of coexpressed miRNAs. Larger studies are now needed in order to further validate and build upon these preliminary findings.

Materials and Methods

Participants

A subset of military personnel (24 males) was selected from the larger Prospective Research in Stress-Related Military Operations (PRISMO) study, a prospective cohort of Dutch military members deployed to Afghanistan for 4 months (Reijnen et al., 2015; Eekhout et al., 2016). Based on the level of combat exposure during deployment and the severity of postdeployment PTSD symptoms, three subgroups were identified: 1) susceptible individuals, i.e., trauma-exposed subjects with deployment-related PTSD symptoms at 6 months follow-up; 2) resilient individuals, i.e., trauma-exposed soldiers with no PTSD diagnosis at follow-up; and 3) controls, i.e., deployed, but nonexposed and mentally healthy military members. Blood samples were collected at the Utrecht University Medical Center at 6 months postdeployment. Trauma exposure was assessed using a 19-item deployment experiences checklist (van Zuiden et al., 2011). The severity of PTSD symptoms was established using the 22-item Self-Rating Inventory for PTSD (SRIP) (Hovens et al., 2002). Information on smoking and alcohol was collected using self-report measures. This study was approved by the ethical committee of University Medical Center Utrecht (01-333/0) and conducted in accordance with the Declaration of Helsinki. All participants gave written informed consent.

RNA Isolation

Total RNA was isolated from 300μl human serum using the mirVana PARIS kit (Ambion) according to the manufacturer’s instructions. Briefly, the samples were incubated with an equal volume of denaturing solution, acid-phenol/chloroform was added, and the samples spun for 5min at 10,000×g. The aqueous phase was recovered and passed through a filter which was washed three times with the provided wash solutions. Final RNA was eluted in 100µl nuclease-free water. The concentrations and quality of the recovered RNA were measured using the Agilent Bioanalyzer 2100 (Agilent Technologies, Inc., CA, USA). All eluates were stored at −80°C until further use.

Small RNA Library Preparation and Next-Generation Sequencing

Barcoded libraries (N = 24, 8 per group) were prepared with an input of 25ng total RNA using the Illumina Small RNA TruSeq kit (Illumina, CA, USA). Briefly, 3′ and 5′ RNA adapters were added, the samples were reverse transcribed and amplified using 11 PCR cycles. All samples were processed in parallel and received a unique barcode. The complementary DNA constructs were gel purified and concentrated by ethanol precipitation. The quality control was performed using Agilent's 2100 Bioanalyzer with a High-Sensitivity DNA Chip. The 24 samples were pooled (N = 8 per group) and sequenced in duplicate using the Illumina HiSeq 2000 DNA sequence platform according to the manufacturer’s protocol (GEO accession: GSE137624).

Small RNA Sequencing Data Analysis

Quality control of the raw sequences was done using FastQC (v. 0.11.3), and reads were preprocessed and mapped to the latest release of miRBase (v. 21) (Baras et al., 2015) utilizing miRge with default settings (Kozomara and Griffiths-Jones, 2014). In order to compensate for bias introduced by very low abundant sequences, only those miRNAs with an average of 50 counts (or more) across samples were considered for further analyses.

Differential Expression Analysis

Data normalization and differential expression analysis was conducted with the DESeq2 package in R (v. 3.5.2) (Love et al., 2014) thereby correcting for age, alcohol use, and smoking status. Resulting p-values were controlled by the false discovery rate (FDR) at 5% (Benjamini and Hochberg, 1995).

Weighted Gene Coexpression Network Construction and Module Detection

The identified miRNAs were used to construct coexpression networks using the Weighted Gene Coexpression Network Analysis (WGCNA) R package (Langfelder and Horvath, 2008). Normalized miRNA data was used as input. An adjacency matrix was generated by calculating Pearson’s correlations between all miRNAs. Next, topological overlap between miRNAs was calculated using a power of 9. We performed 200 rounds of bootstrapping in order to construct a network that is robust to outliers. The cutreeDynamic function in the dynamicTreeCut R package was then used to identify coexpression modules of positively correlated miRNAs with high topological overlap. Modules with at least 30 miRNAs were assigned a color. Modules with highly correlated eigengenes were merged using the mergeCloseModules function in R. Pearson correlations between module eigengenes, age, smoking status, and alcohol were calculated. Welch’s t-tests were performed in order to detect differences between module eigengenes of the control subjects and the trauma-exposed individuals. One-way ANOVAs were performed to detect differences between the three groups. When significant, the post-hoc Tukey HSD test was used to detect pairwise group differences.

Target Gene Pathway and Enrichment Analyses

The experimentally validated miRNA–target interactions database miRTarBase 6.0 (Chou et al., 2018) was used to identify gene targets of miRNAs. In order to narrow down the amount of target genes for further analyses, one-sided Fisher tests (with FDR multiple correction) were performed to evaluate whether the amount of miRNAs targeting a specific gene was significantly higher than expected by chance. Those genes were then analyzed for enriched Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways and Gene Ontology terms (GO terms) using the online Database for Annotation, Visualization and Integrated Discovery (DAVID) v6.8 (Huang da et al., 2009; Huang da et al., 2009).

Statistical Analyses

To detect differences in age, number of previous deployments, cigarette smoking, alcohol use, trauma exposure scores, and SRIP scores between the groups, the Welch ANOVA with Games–Howell post-hoc test was applied. Since data on alcohol use at the 6 months follow-up time point was not available for all subjects, predeployment values were used instead. For each individual, smoking status was estimated based on their unique methylation patterns in 183CpGs, as previously described (Zeilinger et al., 2013). Finally, the classification accuracy of specific miRNAs was determined by calculating the area under the receiver operating characteristic (ROC) curve (AUC) in R.

Results

Demographic Characteristics

A total of 24 subjects were included in the present study, of which 8 developed PTSD following deployment, 8 were resilient, and 8 were nonexposed controls (Supplementary Table 1). Based on the sequencing results, four subjects were excluded due to having a distinctively lower amount of reads causing great variation in expression data between samples. The three groups did not differ in terms of age, number of previous deployments, smoking status, and alcohol use (Table 1). On average, subjects with PTSD and resilient individuals were exposed to a similar amount of traumatic events, which was significantly more than the nonexposed controls [F(2, 8.8) = 54.67, p < 0.001. Games–Howell post-hoc showed p < 0.001 for PTSD versus control, and resilient versus control]. Finally, resilient and control subjects had similar postdeployment PTSD scores as measured by the SRIP, which were significantly lower than the average score of the PTSD group [F(2, 11.15) = 25.23, p < 0.001. Games–Howell post-hoc showed p < 0.001 for PTSD versus resilient, and PTSD versus control].

TABLE 1

Table 1 Demographic characteristics of the 20 subjects remaining after outlier exclusion.

miRNAs Sequencing and Differential Expression Analysis

Small RNA sequencing yielded an average of 9.5 million unfiltered sequencing reads across all samples. After adaptor trimming and size selection, 1.9 million high-quality reads remained, which were aligned to miRNA sequences from miRBase (release 21). As mentioned earlier, principal component analysis revealed the presence of four outliers, which were excluded from further analysis. The count data were then filtered for miRNAs that showed an average of 50 reads or more across all samples. This resulted in the identification of 306 different miRNAs. Differential expression analysis in DESeq2 revealed that a total of 123 miRNAs showed differential expression between PTSD cases and nonexposed controls, while 4 were downregulated in PTSD cases compared to resilient individuals (Supplementary Table 2). Selecting those miRNAs with a log2 fold-change (FC) value≥|1.0| and FDR adjusted p < 0.05 revealed that one miRNA, miR-1246, was downregulated in PTSD subjects compared to resilient subjects and 22 were differentially expressed between PTSD subjects and nonexposed controls (Table 2, Figure 1). Of these, 4 were downregulated and 18 were upregulated. We used the Venn tool to identify those differentially expressed miRNAs that are specific for PTSD only (Figure 2). Two miRNAs were identified at the intersection of the blue and yellow circles, i.e., miR-4454 and miR-210-3p. Both miRNAs were significantly downregulated in PTSD subjects compared to resilient subjects and controls and not differentially expressed between resilient subjects and controls, suggesting that these could be more specific to PTSD (Supplementary Table 2). However, both miRNAs had log2 FC values of −0.61 and −0.54, which does not pass our threshold of ≥|1.0|.

TABLE 2

Table 2 Differentially expressed microRNAs (miRNAs) between posttraumatic stress disorder (PTSD) cases versus controls and PTSD cases versus resilient individuals with a log2 fold-change value≥|1.0| and FDR adjusted p < 0.05.

FIGURE 1

Figure 1 Volcano plots of differentially expressed microRNAs (miRNAs) between posttraumatic stress disorder (PTSD) cases and controls (A) and PTSD cases and resilient subjects (B). Black dots represent nonsignificantly differentially expressed miRNAs, red dots represent significant miRNAs with a log2 FC<|1|, orange dots represent nonsignificant miRNAs with a log2 FC≥|1|, and green dots represent significantly differentially expressed miRNAs with a log2 FC≥|1|. Significance is declared when adjusted p<0.05.

FIGURE 2

Figure 2 Venn diagram showing overlapping microRNAs (miRNAs). Of particular interest here are the two miRNAs at the intersection of the blue and yellow circle (but not the green), i.e., miR-4454 and miR-210-3p. Obtained using https://bioinfogp.cnb.csic.es/tools/venny/index.html

Weighted Gene Coexpression Network Analysis

WGCNA was applied using the 306 identified miRNAs in order to detect clusters of coexpressed miRNAs. Based on the sample dendrogram, one outlier was removed from further analyses (Supplementary Figure 1). We identified three miRNA modules (Figure 3). The turquoise, blue, and brown modules each had 84, 79, and 65 miRNAs, respectively. None of the modules were associated with the potential covariates age, smoking status, or alcohol use (Supplementary Figure 2). Within each module, the module eigengenes were significantly different between trauma-exposed individuals and nonexposed controls for the turquoise and blue modules (p = 2.67×10⁻⁰⁴, p = 2.51×10⁻⁰⁶, respectively) but not for the brown module (p = 0.196; Figure 3A). When stratifying the trauma-exposed individuals into PTSD subjects and resilient subjects, the individual eigengenes of the blue module were significantly different between PTSD subjects and resilient individuals (p = 1.46×10⁻⁰³; Figure 3B), which was not the case for the other modules. We therefore focused on the blue module for further analyses.

FIGURE 3

Figure 3 Significant modules of coexpressed miRNAs. The y-axis displays the module eigengene values. Groups are stratified by trauma exposure (A), and trauma exposure and PTSD status (B). Significances were detected using Welch tests (A) or one-way ANOVAs (B) *p < 0.05, ***p < 0.001.

Out of the 79 miRNAs belonging to this module, 67 were differentially expressed between PTSD subjects and controls (Table 1), including miR-138-5p, the hub miRNA (Table 3). In order to evaluate the diagnostic accuracy of some of these miRNAs, we performed ROC analysis for those miRNAs with the highest absolute module memberships. The five most contributing miRNAs, i.e., miR-221-3p, miR-335-5p, miR-138-5p, miR-222-3p, and miR-146-5p (Table 3), could perfectly distinguish PTSD subjects and controls (AUC of 1 for all miRNAs; Supplementary Figure 3 A.1 and A.2 for miR-221-3p). These miRNAs could equally well differentiate PTSD subjects from resilient subjects, except for miR-221-3p and miR-222-3p (AUC of 0.95 and 0.98, respectively). When obtaining ROC curves using miRNA expression levels adjusted for confounders (i.e., age, smoking, and alcohol use), all miRNAs could still distinguish PTSD subjects from controls (Supplementary Figure 3 B.1 and B.2 for miR-221-3p). However, differentiating PTSD from resilience was less accurate as reflected by AUCs of 0.625, 0.775, 0.725, 0.675, and 0.775 for miR-221-3p, miR-335-5p, miR-138-5p, miR-222-3p, and miR-146-5p, respectively (Supplementary Figure 3 B.1 and B.2 for miR-221-3p).

TABLE 3

Table 3 MicroRNAs (miRNAs) belonging to the Weighted Gene Coexpression Network Analysis (WGCNA) blue module.

Target Gene Pathway and GO Enrichment Analyses

Validated gene targets of the 79 miRNAs comprised within the blue module were obtained from the online database miRTarBase (Chou et al., 2018). In order to narrow down this extensive set of target genes (N = 9270), Fisher tests were performed to select only those genes that were targeted by significantly more miRNAs than expected by chance. This revealed a set of 146 genes, which were considered for pathway and enrichment analyses (Supplementary Table 3). After FDR adjustment, 15 significantly enriched KEGG pathways were identified of which most were cancer-related (Table 4). GO enrichment analyses of these target genes further identified eight significant biological processes, five molecular functions, and six cellular components (Table 5). The most enriched GO terms were related to apoptotic processes, protein binding, and intracellular compartments, respectively (Table 5).

TABLE 4

Table 4 Significant Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways.

TABLE 5

Table 5 Significant gene ontology (GO) terms enriched for a subset of target genes (N = 218) of the coexpressed microRNAs (miRNAs) from the blue module (N = 70).

Discussion

In this study, we aimed to identify the diagnostic biomarker potential of circulating miRNAs for PTSD using serum samples from Dutch military subjects. We further aimed to gain insights in the coexpression patterns of these miRNAs, their predicted gene targets, and underlying biological pathways. Our preliminary findings suggest that 1) certain miRNAs could potentially serve as individual biomarkers of susceptibility, and 2) the coexpression of a specific set of miRNAs could accurately distinguish between subjects with PTSD, resilient individuals, and nonexposed controls. Such markers could be useful in clinical settings for accurate diagnosis and treatment planning, which is especially relevant for individuals who have that have difficulties associating their symptoms to a traumatic event, are unable to describe their symptoms, or are unwilling to fully disclose them (Yehuda et al., 2013).

Differential expression analysis identified 1 differentially expressed miRNA between subjects with PTSD and resilient individuals and 22 between subjects with PTSD and nonexposed controls (after multiple testing correction and a log2 FC cutoff of ≥|1|). Of these, miR-138-5p was significantly overexpressed in subjects with PTSD as compared to controls, and WGCNA revealed that this was the hub miRNA of the blue module. Serum levels of this miRNA were previously found altered in a rat model of restraint stress (Balakathiresan et al., 2014), while hippocampal miR-138-5p levels were associated with the formation of fear memories in mice (Li et al., 2018). Another miRNA, miR-1246, was the only significant miRNA that was downregulated in PTSD cases compared to resilient subjects and had a log2 FC>|1|. This miRNA was previously found downregulated in peripheral blood mononuclear cells of war veterans suffering from PTSD as compared to healthy nontrauma-exposed controls (Bam et al., 2016b). Such findings suggest that these miRNAs could be implicated in PTSD and potentially aid in diagnosing this disorder.

The three modules of coexpressed miRNAs identified by WGCNA revealed that most of the detected miRNAs could be clustered based on similarities in their expression patterns. The blue module contained 79 miRNAs which could significantly differentiate between trauma-exposed individuals and nonexposed controls. Interestingly, within the trauma-exposed individuals, the expression profiles of these miRNAs were significantly different between individuals with and without PTSD. This highlights the importance of including and studying not only non-trauma exposed controls but also trauma-exposed healthy individuals in order to disentangle PTSD effects from trauma-related effects. Moreover, 67 of the miRNAs of the blue module, including its hub miRNA, were significantly differentially expressed between PTSD cases and controls, which enhances the notion that these miRNAs could be relevant for PTSD.

Of the five miRNAs with the highest module membership, we calculated the AUCs to assess their diagnostic accuracy (Grund and Sabin, 2010). In order to determine the biomarker potential of these miRNAs, i.e., their potential to reflect PTSD regardless of any other confounding condition, we used uncorrected miRNA expression values. Interestingly, the results suggest that these miRNAs could almost perfectly distinguish PTSD subjects from resilient individuals and controls. However, these results were not reflected by the DESeq2 analyses in which the expression levels of these miRNAs were not different between PTSD cases and resilient individuals. Part of this discrepancy can most likely be attributed to confounding effects, as DESeq2 analyses were corrected for age, alcohol, and smoking status. When obtaining the ROC curves using confounder-adjusted miRNA expression values, the AUCs more accurately corresponded to the DESeq2 results. Although these results suggest that the expression of our selected miRNAs fluctuates with confounders, they mostly strengthen the need of replication in larger cohorts. This will further be valuable in determining whether these miRNAs could be specific for PTSD only as opposed to trauma more broadly.

Enrichment of GO terms indicated that target genes of the coexpressed miRNAs in the blue module are enriched in several KEGG pathways of which most were cancer-related. This suggests that these miRNAs could be implicated in cancer pathways that are also involved in signaling cascades possibly related to PTSD. The target genes were also involved in several biological processes of which most were involved in apoptotic processes. Previous studies found reduced level of apoptotic markers in the serum of subjects with PTSD (Mkrtchian et al., 2013) and abnormal apoptosis in specific brain regions of animals undergoing single prolonged stress as a model for PTSD (Han et al., 2013; Li et al., 2013; Jia et al., 2018). These findings indicate a potential apoptosis dysfunction that could contribute to the inflammation pattern frequently observed within PTSD (Mkrtchian et al., 2013). Furthermore, the involvement of the identified genes in cellular responses after mechanical stimuli could indicate the need to correct for traumatic brain injuries, which are not uncommon among military members. Unfortunately, this information was not available for the present study. Finally, enriched molecular function GO terms suggest their involvement in the binding of proteins and RNA, while the significant cellular component GO terms show involvement of intracellular compartments such as the cytosol and the nucleoplasm.

Of note, the present paper refers to trauma-exposed healthy individuals as being “resilient” in order to create a clear differentiation between trauma-exposed healthy subjects and nonexposed control subjects. However, we do acknowledge and emphasize that resilience is more than just the reverse side of PTSD or the absence of symptomatology (Kalisch et al., 2017; Snijders et al., 2018b). Instead, resilience is an active and dynamic process that needs to remain separated from the multifaceted and complex nature of PTSD. This complexity further suggests that identifying one true and valid biomarker of susceptibility is likely not realistic. We therefore urge future studies to combine findings such as the ones presented in this paper with several other biological networks and phenotypic profiles in order to develop a cross-dimensional, global understanding of PTSD.

The main strength of this study lies in the inclusion of three different groups, i.e., PTSD subjects, resilient subjects, and nonexposed healthy controls, which allows us to disentangle PTSD- from trauma-related effects. However, the study is mainly limited by its relatively small sample size consisting of male subjects only. Given the existence of female- and male-biased miRNAs as recently reported by Cui, Yang et al. (2018) (Cui et al., 2018), these findings may not be applicable to the female population. This study population may also differ from other cohorts such as civilians in terms of demographics, psychological characteristics, and type of experienced trauma, which limits the extrapolation potential. Next, one could question the validity of self-report PTSD measures and whether the observed markers are specific to PTSD since certain comorbidities such as (history of) traumatic brain injuries were not available and thus not accounted for.

In conclusion, this paper presents preliminary evidence for using specific miRNAs as diagnostic biomarkers of PTSD, either individually or grouped within coexpressed clusters. Identifying reliable biomarkers of PTSD is essential for accurate diagnosis and treatment planning. We therefore encourage future studies to build upon these findings by aiming to replicate these in larger cohorts and thus pave the way for functional studies to gain insights into the precise roles of these miRNAs in stress susceptibility.

Data Availability Statement

The Data is available at GEO, accession: GSE137624.

Ethics Statement

The studies involving human participants were reviewed and approved by The ethical committee of University Medical Center Utrecht (01-333/0). The patients/participants provided their written informed consent to participate in this study.

Author Contributions

LN, BR, GK, DH, JKl, and MK participated in the conception and the design of the study. MB, CV, EV, and EG recruited the participants and collected the blood samples. CS and BM performed the experiments. CS, JKr, EP, and LE performed the bioinformatics analysis. CS wrote the manuscript. LN, BR, and MK provided the necessary fundings. All authors critically read, commented, provided scientific content and approved the final manuscript.

Funding

This work has been funded by the European Union’s Horizon 2020 research and innovation program under the Marie Sklodowska-Curie grant agreement N0. 707362 (LN) and by a VIDI award number 91718336 from the Netherlands Scientific Organization (BR).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2019.01042/full#supplementary-material

Supplementary Figure S1 | Sample dendrogram to detect outliers. Clustering was based on miRNA expression data. Sample names were relabeled as 1-20.

Supplementary Figure S2 | Correlations of modules detected by WGCNA and the following potential covariates: age, smoking status and alcohol use. P-values are presented between brackets. ME: module eigengene.

Supplementary Figure S3 | Receiver operating characteristic (ROC) curve for the miRNA with the highest module membership in the blue module, i.e. miR-221-3p. The graphs represent PTSD vs control without confounders (A.1) or with confounders (B.1), and PTSD vs resilient without confounders (A.2) or with confounders (B.2).

References

American Psychiatric Association. (2013). Diagnostic and Statistical Manual of Mental Disorders (5th ed.). Arlington, VA: American Psychiatric Association.

Google Scholar

Balakathiresan, N. S., Chandran, R., Bhomia, M., Jia, M., Li, H., Maheshwari, R. K., et al. (2014). Serum and amygdala microRNA signatures of posttraumatic stress: fear correlation and biomarker potential. J. Psychiatric Res. 57, 65–73. doi: 10.1016/j.jpsychires.2014.05.020

CrossRef Full Text | Google Scholar

Bam, M., Yang, X., Zumbrun, EE., Zhong, Y., Zhou, J., Ginsberg, J. P., et al. (2016).Dysregulated immune system networks in war veterans with PTSD is an outcome of altered miRNA expression and DNA methylation. Sci. Rep. 6, 31209. doi: 10.1038/srep31209

PubMed Abstract | CrossRef Full Text | Google Scholar

Bam, M., Yang, X., Zhou, J., Ginsberg, J. P., Leyden, Q., Nagarkatti, P. S., et al. (2016). Evidence for epigenetic regulation of pro-inflammatory cytokines, interleukin-12 and interferon gamma, in peripheral blood mononuclear cells from PTSD patients. J. Neuroimmune Pharmacol. 11 (1), 168–181. doi: 10.1007/s11481-015-9643-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Baras, A. S., Mitchell, C. J., Myers, J. R., Gupta, S., Weng, L. C., Ashton, J. M., et al. (2015). miRge - A Multiplexed Method of Processing Small RNA-Seq Data to Determine MicroRNA Entropy. PLoS One 10p (11), e0143066. doi: 10.1371/journal.pone.0143066

PubMed Abstract | CrossRef Full Text | Google Scholar

Benjamini, Y., Hochberg, Y. (1995). Controlling the false discovery rate: a practical and powerful approach to multiple testing. J. R. Stat. Soc. 57, 289–300. doi: 10.1111/j.2517-6161.1995.tb02031.x

CrossRef Full Text | Google Scholar

Chou, C. H., Shrestha, S., Yang, C. D., Chang, N. W., Lin, Y. L., Liao, K. W., et al. (2018). miRTarBase update 2018: a resource for experimentally validated microRNA-target interactions. Nucleic Acids Res. 46, D296–D302. doi: 10.1093/nar/gkx1067

PubMed Abstract | CrossRef Full Text | Google Scholar

Cui, C., Yang, W., Shi, J., Zhou, Y., Yang, J., Cui, Q., et al. (2018). Identification and analysis of human sex-biased MicroRNAs. Genomics Proteomics Bioinf. 16 (3), 200–211. doi: 10.1016/j.gpb.2018.03.004

CrossRef Full Text | Google Scholar

Eekhout, I., Reijnen, A., Vermetten, E., Geuze, E., (2016). Post-traumatic stress symptoms 5 years after military deployment to Afghanistan: an observational cohort study. Lancet Psychiatry 3 (1), 58–64. doi: 10.1016/S2215-0366(15)00368-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Ferry, F. R., Brady, S. E., Bunting, B. P., Murphy, S. D., Bolton, D., O’Neill, S. M., et al. (2015). The economic burden of PTSD in Northern Ireland. J Trauma Stress 28 (3), 191–197. doi: 10.1002/jts.22008

PubMed Abstract | CrossRef Full Text | Google Scholar

Fontana, A., Rosenheck, R. (1994). Traumatic war stressors and psychiatric symptoms among World War II, Korean, and Vietnam War veterans. Psychology Aging 9 (1), 27–33. doi: 10.1037/0882-7974.9.1.27

PubMed Abstract | CrossRef Full Text | Google Scholar

Goldberg, A. D., Allis, C. D., Bernstein, E. (2007). Epigenetics: a landscape takes shape. Cell 128 (4), 635–638. doi: 10.1016/j.cell.2007.02.006

PubMed Abstract | CrossRef Full Text | Google Scholar

Grund, B., Sabin, C. (2010). Analysis of biomarker data: logs, odds ratios, and receiver operating characteristic curves. Curr. Opin. HIV AIDS 5 (6), 473–479. doi: 10.1097/COH.0b013e32833ed742

PubMed Abstract | CrossRef Full Text | Google Scholar

Han, F., Yan, S., Shi, Y. (2013). Single-prolonged stress induces endoplasmic reticulum-dependent apoptosis in the hippocampus in a rat model of post-traumatic stress disorder. PLoS One 8 (7), e69340. doi: 10.1371/journal.pone.0069340

PubMed Abstract | CrossRef Full Text | Google Scholar

Hovens, J. E., Bramsen, I., van der Ploeg, H. M. (2002). Self-rating inventory for posttraumatic stress disorder: review of the psychometric properties of a new brief Dutch screening instrument. Percept. Motor Skills 94 (3), 996–1008. doi: 10.2466/pms.2002.94.3.996

CrossRef Full Text | Google Scholar

Huang da, W., Sherman, B. T., Lempicki, R. A. (2009). Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources. Nat. Protoc. 4 (1), 44–57. doi: 10.1038/nprot.2008.211

PubMed Abstract | CrossRef Full Text | Google Scholar

Huang da, W., Sherman, B. T., Lempicki, R. A. (2009). Bioinformatics enrichment tools: paths toward the comprehensive functional analysis of large gene lists. Nucleic Acids Res. 37 (1), 1–13. doi: 10.1093/nar/gkn923

PubMed Abstract | CrossRef Full Text | Google Scholar

Jia, Y., Han, Y., Wang, X., Han, F. (2018). Role of apoptosis in the post-traumatic stress disorder model-single prolonged stressed rats. Psychoneuroendocrinology 95, 97–105. doi: 10.1016/j.psyneuen.2018.05.015

PubMed Abstract | CrossRef Full Text | Google Scholar

Kalisch, R., Baker, D. G., Basten, U., Boks, M. P., Bonanno, G. A., Brummelman, E., et al. (2017). The resilience framework as a strategy to combat stress-related disorders. Nat. Hum. Behav. 1 (11), 784–790. doi: 10.1038/s41562-017-0200-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Koenen, K. C., Ratanatharathorn, A., Ng, L., McLaughlin, K. A., Bromet, E. J., Stein, D. J., et al. (2017). Posttraumatic stress disorder in the world mental health surveys. Psychol. Med. 47 (13), 2260–2274. doi: 10.1017/S0033291717000708

PubMed Abstract | CrossRef Full Text | Google Scholar

Kozomara, A., Griffiths-Jones, S. (2014). miRBase: annotating high confidence microRNAs using deep sequencing data. Nucleic Acids Res. 42, D68–D73. doi: 10.1093/nar/gkt1181

PubMed Abstract | CrossRef Full Text | Google Scholar

Langfelder, P., Horvath, S. (2008). WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics 9, 559. doi: 10.1186/1471-2105-9-559

PubMed Abstract | CrossRef Full Text | Google Scholar

Lehrner, A., Yehuda, R. (2014). Biomarkers of PTSD: military applications and considerations. Eur. J. Psychotraumatol. 5, 10.3402/ejpt.v5.23797. doi: 10.3402/ejpt.v5.23797

CrossRef Full Text | Google Scholar

Li, Y., Han, F., Shi, Y. (2013). Increased neuronal apoptosis in medial prefrontal cortex is accompanied with changes of Bcl-2 and bax in a rat model of post-traumatic stress disorder. J. Mol. Neurosci. 51 (1), 127–137. doi: 10.1007/s12031-013-9965-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, D. W., Liu, J. Z., Li, S. C., Yang, J. B., Sun, H. H., Wang, A. H., et al. (2018).Fear conditioning downregulates miR-138 expression in the hippocampus to facilitate the formation of fear memory. Neuroreport 29 (16), 1418–1424. doi: 10.1097/WNR.0000000000001129

PubMed Abstract | CrossRef Full Text | Google Scholar

Love, M. I., Huber, W., Anders, S. (2014). Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 15 (12), 550. doi: 10.1186/s13059-014-0550-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Martin, C. G., Kim, H., Yun, S., Livingston, W., Fetta, J., Mysliwiec, V., et al. (2017). Circulating miRNA associated with posttraumatic stress disorder in a cohort of military combat veterans. Psychiatry Res. 251, 261–265. doi: 10.1016/j.psychres.2017.01.081

PubMed Abstract | CrossRef Full Text | Google Scholar

Mkrtchyan, G. M., Boyadzhyan, A. S., Avetyan, D. G., Sukiasyan, S. G., (2014). The involvement of abnormal apoptosis in the disturbance of synaptic plasticity in posttraumatic stress disorder. Neurosci. Behav. Physiol. 44 (4), 26–29.

Google Scholar

Nock, M. K., Kessler, R. C. (2006). Prevalence of and risk factors for suicide attempts versus suicide gestures: analysis of the national comorbidity survey. J. Abnorm. Psychol. 115 (3), 616–623. doi: 10.1037/0021-843X.115.3.616

PubMed Abstract | CrossRef Full Text | Google Scholar

Peschansky, V. J., Wahlestedt, C. (2014). Non-coding RNAs as direct and indirect modulators of epigenetic regulation. Epigenetics 9 (1), 3–12. doi: 10.4161/epi.27473

PubMed Abstract | CrossRef Full Text | Google Scholar

Reijnen, A., Rademaker, A. R., Vermetten, E., Geuze, E. (2015). Prevalence of mental health symptoms in Dutch military personnel returning from deployment to Afghanistan: a 2-year longitudinal analysis. Eur. Psychiatry 30 (2), 341–346. doi: 10.1016/j.eurpsy.2014.05.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Richter-Levin, G., Stork, O., Schmidt, M. V. (2018). Animal models of PTSD: a challenge to be met. Mol. Psychiatry 24, 1135–1156. doi: 10.1038/s41380-018-0272-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Schmidt, U., Holsboer, F., Rein, T. (2011). Epigenetic aspects of posttraumatic stress disorder. Dis. Markers 30 (2-3), 77–87. doi: 10.1155/2011/343616

PubMed Abstract | CrossRef Full Text | Google Scholar

Snijders, C., de Nijs, L., Baker, D. G., Hauger, R. L., van den Hove, D., Kenis, G., et al. (2018a). MicroRNAs in post-traumatic stress disorder. Curr. Top Behav. Neurosci. 38, 23–46. doi: 10.1007/7854_2017_32

PubMed Abstract | CrossRef Full Text | Google Scholar

Snijders, C., Pries, L. K., Sgammeglia, N., Al Jowf, G., Youssef, N. A., de Nijs, L., et al. (2018b). Resilience against traumatic stress: current developments and future directions. Front. Psychiatry 9, 676. doi: 10.3389/fpsyt.2018.00676

PubMed Abstract | CrossRef Full Text | Google Scholar

Tarrier, N., Gregg, L. (2004). Suicide risk in civilian PTSD patients–predictors of suicidal ideation, planning and attempts. Soc. Psychiatry Psychiatr. Epidemiol. 39 (8), 655–661. doi: 10.1007/s00127-004-0799-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Thomas, M. M., Harpaz-Rotem, I., Tsai, J., Southwick, S. M., Pietrzak, R.H. (2017). Mental and physical health conditions in us combat veterans: results from the national health and resilience in veterans study. Prim. Care Companion CNS Disord. 19 (3). doi: 10.4088/PCC.17m02118

CrossRef Full Text | Google Scholar

van Zuiden, M., Geuze, E., Willemen, H. L., Vermetten, E., Maas, M., Heijnen, C. J., et al. (2011). Pre-existing high glucocorticoid receptor number predicting development of posttraumatic stress symptoms after military deployment. Am. J. Psychiatry 168 (1), 89–96. doi: 10.1176/appi.ajp.2010.10050706

PubMed Abstract | CrossRef Full Text | Google Scholar

Wingo, A. P., Almli, L. M., Stevens, J. S., Klengel, T., Uddin, M., Li, Y., et al. (2015). DICER1 and microRNA regulation in post-traumatic stress disorder with comorbid depression. Nat. Commun. 6, 10106. doi: 10.1038/ncomms10106

PubMed Abstract | CrossRef Full Text | Google Scholar

Yehuda, R., Neylan, T. C., Flory, J. D., McFarlane, A. C. (2013). The use of biomarkers in the military: from theory to practice. Psychoneuroendocrinology 38 (9), 1912–1922. doi: 10.1016/j.psyneuen.2013.06.009

PubMed Abstract | CrossRef Full Text | Google Scholar

Yehuda, R., Vermetten, E., McFarlane, A. C., Lehrner, A. (2014). PTSD in the military: special considerations for understanding prevalence, pathophysiology and treatment following deployment. Eur. J. Psychotraumatol. 5. doi: 10.3402/ejpt.v5.25322

CrossRef Full Text | Google Scholar

Zeilinger, S., Kühnel, B., Klopp, N., Baurecht, H., Kleinschmidt, A., Gieger, C., et al. (2013). Tobacco smoking leads to extensive genome-wide changes in DNA methylation. PLoS One 8 (5), e63812. doi: 10.1371/journal.pone.0063812

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhou, J., Nagarkatti, P., Zhong, Y., Ginsberg, J. P., Singh, N. P., Zhang, J., et al. (2014). Dysregulation in microRNA expression is associated with alterations in immune functions in combat veterans with post-traumatic stress disorder. PLoS One 9 (4), e94075. doi: 10.1371/journal.pone.0094075

PubMed Abstract | CrossRef Full Text | Google Scholar