Dysregulated immune system networks in war veterans with PTSD is an outcome of altered miRNA expression and DNA methylation

Post-traumatic stress disorder patients experience chronic systemic inflammation. However, the molecular pathways involved and mechanisms regulating the expression of genes involved in inflammatory pathways in PTSD are reported inadequately. Through RNA sequencing and miRNA microarray, we identified 326 genes and 190 miRNAs that were significantly different in their expression levels in the PBMCs of PTSD patients. Expression pairing of the differentially expressed genes and miRNAs indicated an inverse relationship in their expression. Functional analysis of the differentially expressed genes indicated their involvement in the canonical pathways specific to immune system biology. DNA methylation analysis of differentially expressed genes also showed a gradual trend towards differences between control and PTSD patients, again indicating a possible role of this epigenetic mechanism in PTSD inflammation. Overall, combining data from the three techniques provided a holistic view of several pathways in which the differentially expressed genes were impacted through epigenetic mechanisms, in PTSD. Thus, analysis combining data from RNA-Seq, miRNA array and DNA methylation, can provide key evidence about dysregulated pathways and the controlling mechanism in PTSD. Most importantly, the present study provides further evidence that inflammation in PTSD could be epigenetically regulated.

disorders, and autoimmune diseases 10,11 , all of which have an inflammatory component. Thus, it is important to study the molecular changes occurring in the immune cells and the inflammation manifested in PTSD patients. Peripheral blood mononuclear cells (PBMCs), which include T cells, B cells and monocytes are major players in the peripheral immune system and constitute cells that are both pro-and anti-inflammatory in nature. Additionally, many of the cytokines and chemokines of the pro-and anti-inflammatory milieu are produced by the PBMCs and can give rise to profound changes in the immune response [12][13][14] . The consensus now is that, even though PTSD is a psychiatric disorder, PBMCs play a major role in exacerbating the symptoms 15,16 . However, the knowledge of initiation of inflammation and the canonical pathways dysregulated during PTSD in PBMCs are poorly understood.
Recent use of powerful techniques like Next Generation Sequencing (NGS) and microarrays has led to the identification of differentially expressed mRNAs and miRNAs at the global level. Moreover, combining data from these techniques with various bioinformatics tools for data analysis has made it possible to predict and discern biological pathways that are affected in a diseased state, thereby making it possible to overcome the limitations of single gene-based studies. Furthermore, data from miRNA expression arrays and DNA methylation studies can be used to predict and study the regulation of expression of differentially expressed genes at the epigenetic level. Epigenetic regulation of gene expression includes the influence from processes like DNA methylation, histone modifications and miRNA expression [17][18][19] . Micro-RNAs are small, ~22 nucleotide long, non-coding, regulatory RNAs and are one of the key epigenetic entities that regulate the expression of genes at the post-transcriptional level 20 . The mechanism of gene regulation by miRNAs involve physical interactions with complimentary sequences, typically at the 3′ untranslated regions (UTR) of an mRNA leading to the degradation of the mRNA or inhibition of translation 21,22 . Recently, in mammals, it was reported that most (66%-> 90%) of the mRNA-miRNA interaction leads to destabilization of the mRNA 23 , clearly implying that differential expression of miRNAs can lead to change in the level of transcripts of gene(s). Regulation of numerous immune system genes by miRNAs has already been reported demonstrating that the immune system is tightly regulated by the miRNAs 6,7,[24][25][26] . For example, our lab has shown that elevated expression of interferon gamma (IFNG) in PTSD is regulated by hsa-miR-125a, which is found to be downregulated in PTSD patients 6 . In a related study, we have shown that expression of another pro-inflammatory cytokine, interleukin (IL)-12, is elevated in PTSD and its expression is correlated with downregulation of hsa-miR-193a-5p 7 and many other miRNAs (data not shown). On the other hand, unlike miRNAs, DNA methylation at CpG islands present near the transcription start sites (TSS) can regulate the expression of a gene at the transcriptional level [27][28][29] . Effects brought about by DNA methylation is also reported to influence the expression of several genes including immune system related genes 30,31 . Until now, the majority of the studies intended to identify gene regulators have focused on a single or few genes/miRNAs at a time. Moreover, to our knowledge there is no report on expression paired analysis using RNA-Seq and miRNA array data to correlate differentially expressed genes in PTSD. Therefore, combining RNA-Seq and miRNA array data is a novel approach to simultaneously identify several genes and their biological pathways that are dysregulated and their regulators (miRNAs, in this case) that are possibly causing the differential expression of the genes.
In the present study, we could identify and predict many differentially expressed genes involved in canonical pathways related to the immune system biology, in the PBMCs of PTSD patients. We further provide preliminary evidence at the global level that the differential expression of the genes is possibly an outcome of differential expression of miRNAs and change in DNA methylation level.

Results
Gene expression analysis in PTSD reveals differentially expressed genes. We performed RNA-Seq on RNA samples from PBMCs of PTSD patients ( Table 1 provides the demographics of the controls and the patients included for the RNA-Seq analysis) and, on average, obtained ~60 million reads per sample which resulted in a good sequencing depth upon considering the size of human genome. A total of 40420 mRNA and 11218 non-coding RNA Ids were present in the list after obtaining the RNA-Seq data (Fig. 1a) and 48518 Ids were obtained after initial quality control. The list after quality control included both coding and non-coding RNAs as well as their transcript/splice variants, due to which the number was possibly very high for PBMCs. We selected genes keeping a cutoff value of 2 or more fold change plus a p value < 0.005 and obtained a total of 326 mRNAs (Fig. 1b) and 40 non-coding RNAs which were significantly changed in their expression levels based on RNA-Seq analysis (Supplementary Table S3 provides the list of the differentially expressed Ids). Of the 326 mRNAs, 64 (19.63%) were up-regulated and 262 (80.36%) downregulated (here, the differentially expressed genes included only the protein coding genes and not the non-coding Ids which made a significant proportion of the total Ids obtained after RNA-Seq data analysis). Since miRNAs are the main entities in this study, we looked for possible alterations in the transcript level of primary and pre-miRNAs. None of the significantly differentially expressed non-coding Ids included primary or pre-miRNAs. With regards to mature miRNAs, it is understood that due to their small size (~22 nucleotides), they will not be included during size selection in the RNA-Seq library preparation stage stage. Table 1 provides the demographics of the controls and the patients included for the miRNA microarray analysis. At the time of performing the miRNA array, 847 probes for miRNAs were used (Fig. 1c). To identify the differentially expressed miRNAs, we employed selection criteria as p value less than or equal to 0.05 (significant) and linear fold change of at least ±1.5 or more. Thus, we obtained a total of 190 miRNAs ( Fig. 1d  21.63% of the total miRNAs, however, making it 96.31% of the total miRNAs which were significantly different in their expression in PTSD as compared to control. mRNA functional enrichment reveals dysregulation in immune system pathways in PTSD. We first performed functional enrichment in IPA to identify canonical pathways with the 326 genes that indicated significantly differential expression. Majority of the pathways were related to immune system functioning/biology. To that end, we selected the top 20 canonical pathways for further analysis (Fig. 2a) of which we found Th cell differentiation pathway as the one that was most relevant to our studies (Fig. 2b). In this particular pathway, STAT4, TBX21 and HLA-DQA1 were differentially expressed. STAT4 and TBX21 were up-regulated and HLA-DQA1 was downregulated based on RNA-Seq analysis. The canonical pathways were then categorized on the basis of number of genes present from our dataset. The top canonical pathway revealed was 'agranulocyte adhesion and diapedesis' in which 16 of the genes from our dataset were present ( Table 2, and Supplementary Fig.  S1). Of the 16 genes, 4 were up-and 12 were downregulated. The up-regulated genes included CXCL2, CXCL3, CCL4 and CCL5. To give a clearer picture of the up-or downregulated genes in a specific canonical pathway, the list of genes and the top 20 canonical pathways are provided in Table 2. Thereafter, we used the differentially expressed genes for functional enrichment and gene ontology analysis by analyzing them in Panther pathways and DAVID as well. Based on the Panther pathway analysis (Fig. 2c), 'inflammation mediated by chemokine and cytokine pathway (P00031)' was the top pathway with 10 genes present from our dataset. Gene ontology analysis by DAVID revealed that the most significant (lowest p value, 3.80E-07) ontology was 'immune system process' with 40 genes from our dataset ( Table 3).

Micro-RNA expression analysis indicates a global downregulation in PTSD.
Expression pairing and other miRNA target analysis reveals differentially expressed miRNAs that target genes related to immune system in PTSD. We paired the expression of all the differentially expressed miRNAs (190) and genes (326) by using the IPA Expression pairing tool. From the total 190 differentially expressed miRNAs in the dataset, 44 (23.16%) miRNAs were identified to have targets in the RNA-Seq dataset during the time of analysis in IPA. Of the 326 differentially expressed genes, > 70 (21.48%) were identified to be targets of the 44 paired miRNAs ( Supplementary Fig. S2). A list of the paired genes and miRNAs is provided in supplementary files (Supplementary Table S5) and Fig. 3a shows the miRNAs and the target genes in the IPA-generated interactive network. Functional enrichment in IPA of the genes paired with the 44 miRNAs revealed that most of the canonical pathways were related to the immune system biology (Fig. 3b), similar to that observed in the first analysis (Fig. 2a). It is also worth mentioning that Th cell differentiation pathway (− log p value = 2.86E00) was one of the top canonical pathway obtained after expression pairing. Table 4 provides the list of the genes in all the top pathways shown in Fig. 3b.
As mentioned earlier, because all miRNAs or genes are not covered by one single database or tool, we wanted to see whether there are genes and miRNAs in our dataset that were not picked during Expression pairing in IPA, but are predicted to be target(s) of miRNAs in our dataset. Based on the miRNA target gene information available on www.targetscan.org website, we could identify many more miRNAs (78, list not exhaustive) that target  the genes (top 10 differentially expressed) from our dataset. Table 5 provides names of the top ten differentially expressed genes from our dataset based on the highest fold change values and miRNAs with the highest relevancy scores as provided in www.targetscan.org.
Altered DNA methylation pattern is evident during PTSD and correlates to differential gene expression. We compared the level of DNA methylation at the CpG sites of differentially expressed genes (326) present in our dataset. Based on the Illumina CpG site designations, only 138 genes from our dataset had true CpG sites. As per our set criteria used for this analysis, in the PTSD group, 40 (12.27%) genes had higher DNA methylation average β -values at their corresponding CpG sites when compared to that of controls (Supplementary Table S6) and these 40 genes had lower expression values as per RNA-Seq analysis. On the other hand, only 12 (3.68%) genes had decreased DNA methylation average β -values in PTSD (Supplementary Table S6), and the expression level of those genes was higher in PTSD group. In table 6, only the top 10 up-and downregulated genes with their DNA methylation levels are shown. The DNA methylation values were not significant as per Student's t-test or Wilcoxon test. However, in all the genes listed in this study, there was a gradual trend of DNA methylation values supporting the expression levels of the genes ( Fig. 4a provides DNA methylation levels of the genes as box plot; 4b provides fold change values of the genes, following RNA-Seq analysis, listed in Fig. 4a, and Table 6 provides DNA methylation levels of the top 20 genes).
qRT-PCR confirmation of RNA-seq results. To validate our RNA-Seq data, we selected seven genes (MTRNR2L1, MMP25, CXCL8, G0S2, GZMB, CXCL3 and STAT4) as representative for qRT-PCR analysis. MTRNR2L1, MMP25, CXCL8 and G0S2 were shown to be significantly downregulated and GZMB, CXCL3 and STAT4 significantly up-regulated by RNA-Seq analysis. The entire above mentioned gene expression matched with that of RNA-Seq results as per qRT-PCR data (Fig. 4c).

Discussion
The prevalence of PTSD is high among war veterans and in the general public who experience traumatic events. However, the nature of changes occurring in PTSD, at the molecular level, in the PBMCs, is largely unclear. In the current study, we were therefore interested in exploring the global differences in the gene expression during active PTSD in PBMCs and, most importantly, correlate the difference with any altered epigenetic marks. Consequently, we first obtained global RNA expression pattern in the PBMCs obtained from war veterans diagnosed positive for PTSD during the time of sample collection and identified differentially expressed genes and the related immune system canonical pathways. A novel approach used in the present study is the strategy to combine RNA-Seq data with miRNA array data to simultaneously identify differentially expressed genes and their regulators at the global level. This approach provided broader information on the differentially expressed genes and their regulators by eliminating the limiting factors associated with single gene studies. Thus, we could correlate the differential expression of the target genes with the differential expression pattern of relevant miRNAs. To our knowledge, though preliminarily, this is the first time differentially expressed gene networks specific for the immune system are shown to correlate with global alteration of miRNA expression in the PBMCs of PTSD patients. Furthermore, we correlated the expression of several differentially expressed genes with altered DNA methylations at the corresponding CpG sites of the promoter of the respective genes. Functional enrichment of the differentially expressed genes indicated probable alteration of Th cell differentiation pathway. This pathway was one of the major pathways with three genes (STAT4, TBX21, and HLA-DQA1) present from the significant set of genes. STAT4 and TBX21 have crucial roles in regulating T cell functions. For example, TBX21 is the main transcription factor for the expression of interferon gamma, a pro-inflammatory gene and already reported by our group (Bam et al. 7 ) to be elevated in PTSD. This observation is important because the fate of Th cells decides the outcome of immune cell functions whether to be pro-or anti-inflammatory in nature. It also further supports the report on differential expression of T cell produced pro-inflammatory cytokine(s) in PTSD 6 . Therefore, we conclude that an alteration in the T cell biology is possibly one of the root causes for The top 20 canonical pathways selected for finding the genes common in more than one canonical pathway (overlap). Many of the differentially expressed genes are present in multiple pathways related to immune system biology. Table 2 has the list of genes from our dataset that are present in all the canonical pathways in the list. (b) T helper cell differentiation canonical pathway with genes differentially expressed in PTSD. Red and green colors indicate up-and downregulated genes, respectively, in PTSD. The pathway was generated by analyzing all the differentially expressed genes in IPA. (c) The differentially expressed genes were analyzed on Panther pathways analysis tool. The Panther pathway with highest number of genes (10) from the dataset was "inflammation mediated by chemokine and cytokine signaling pathway" (P00031). the underlying inflammation seen during PTSD. Altogether, our observation corroborates well with previous PTSD reports employing RNA-Seq technique with RNA obtained from peripheral blood leukocytes 8 . The authors reported that several genes involved in the innate immune system network were differentially expressed in PTSD. Similarly, in the present study, the top canonical pathways with differentially expressed genes were from the innate immune system in addition to some disease specific pathways. For example, the top canonical pathway with the largest number of differentially expressed genes from our dataset was 'agranulocyte/granulocyte adhesion and diapedesis' . This pathway describes the stages involved in the movement or migration of leukocytes out of the circulatory system to the site of tissue damage or infection during an inflammatory response [32][33][34] . Chemotactic molecules (chemokines) like those secreted by monocytes, macrophages and other immune cells play an important role during this process 35,36 . Function of chemokines is mainly to bring about migration (homing) of leukocytes in the respective sites during homeostasis and inflammatory processes 37 . Other functions of chemokines are seen during different processes like maturation, activation and differentiation for different types of leukocytes [38][39][40] . We observed that expression of many of the chemokines and their receptors (CCL4, CCL5, CXCL1, CXCL2, CXCL3, CXCL6, CXCL8, CXCR1, and CXCR2) were altered in PTSD patients. Chemokines like CCL4, CCL5, CXCL1, CXCL2, CXCL3, CXCL6 and CXCL8 are considered to be pro-inflammatory in their function 41 . Another example of a canonical pathway with several differentially expressed genes from our dataset was "dendritic cell maturation", which plays a critical role in antigen processing and presentation.
Micro-RNAs are critical regulators of gene expression and their interaction with target mRNAs can lead to destabilization of the mRNA in most (> 90%) of the cases in mammals 23 . Expression of genes of the immune system is also known to be controlled by these small RNA molecules 42 . Interestingly, we found that several up-regulated genes in PTSD were either known or predicted to be a target of the many downregulated miRNAs from our dataset. For example, hsa-miR-125a and hsa-miR-193a-5p, which target Interferon gamma (IFNG) 6 and Interleukin-12B (IL12B) 7 respectively, were downregulated in PTSD as reported from our group previously. Both these genes are pro-inflammatory in their functioning 43,44 . As another example, we observed that TBX-21 and STAT4, the key genes in the Th cell differentiation pathway, were upregulated in PTSD. TBX-21 is one of the main transcription factors involved in the differentiation and functioning of Th1 cells and also plays an important role in the functioning of other cells of the immune system 45,46 . STAT4 is induced in response to signaling via the IL-12 pathway, leading to induction of IFN in Th1 type CD4+ cells 47,48 . There were several downregulated miRNAs from our dataset that were predicted to target TBX-21 and STAT4. These observations suggested that several of the up-regulated genes in the PBMCs of PTSD patients could be resulting from the decreased presence of miRNAs, as a result of their downregulation, that target the genes post transcription. Further analysis of the miRNA dataset indicated that, up-regulated chemokines CXCL2, CXCL3, CCL4 and CCL5 are targets of numerous downregulated miRNAs from our dataset. Altogether, these observations suggest that there is a breakdown in the miRNA-mediated gene regulation in the PBMCs during PTSD and strikingly, it includes several of the up-regulated genes having pro-inflammatory properties. In contrast to the above observations, we found only 7 upregulated miRNAs and very few genes from our dataset were shown to be their targets. Moreover, these genes were not present in the top canonical pathways from our analysis. We observed that a higher percentage (80.36%) of the differentially expressed genes in PTSD patients were downregulated. Correlating with the downregulated gene expression, there was a gradual but clear trend of higher DNA methylation level at CpG sites. Higher DNA methylation in the promoter region of a gene correlates with lower transcription of the gene and vice versa [27][28][29] . For example, our group has shown that IL-12 transcription is increased in PTSD which correlated with lower DNA methylation at its promoter region 7 . Thus, it is possible that the genes in the PBMCs of PTSD patients are differentially expressed, at least in part, because of an altered DNA methylation. As another example, CSRNP1 (aka AXUD1), a tumor suppressor, is one of the significantly downregulated genes in PTSD and it has higher level of DNA methylation in PTSD (82.08%) than control (79.09%). Downregulation of CSRNP1 correlates with progression of cancer 49 . This property can be extrapolated to reason that CSRNP1 downregulation can lead to higher cellular proliferation. Previously, our group has reported that CD4+ T cell, CD8+ T cell and B cell numbers are higher in PTSD patients 6 , and this is an indication that the proliferation of these cells is higher during PTSD. Thus, it is possible that the higher cell numbers of lymphocyte subsets seen in PTSD may result from lower CSRNP1 expression which probably is because of the higher DNA methylation at its promoter region. Furthermore, we found that STAT4 is among the significantly up-regulated genes and has correspondingly lower DNA methylation trend at the CpG island. STAT4, induced by IL-12, leads to activation and proliferation of Th1 type CD4+ cells, which correlates well with our observation of increased T cells and alteration in the Th cell differentiation pathway in PTSD 6 . Thus, we hypothesize that the DNA methylation of immune system genes is altered during PTSD resulting in differential expression of a section of the genes.
In summary, the present work has identified differentially expressed genes and miRNAs and the related canonical immune system pathways in the PBMCs of PTSD inflicted war veterans. Furthermore, we provide evidence that many genes have altered DNA methylation at their CpG islands and the expression of the associated genes inversely correlate in PTSD patients. Taken together, the present and previous reports from our lab, and from other research groups, clearly indicate that miRNAs and DNA methylation play a critical role in the modulation of the immune system, with a special emphasis on chronic inflammation seen in PTSD. Most importantly, our findings, although preliminary, open future directions for studies in a pathway specific manner and targeting specific gene regulators to develop novel management strategies and therapies to control the inflammatory response seen during PTSD in war return veterans and the general population.

RNA-sequencing (RNA-Seq).
For RNA-Seq, five controls and five PTSD patients were analyzed. RNA-Seq libraries were constructed using Illumina TruSeq RNA Sample Preparation kit. Briefly, total RNA was purified from PBMCs using the Qiagen RNA easy kit. The oligo-dT beads were added to 1 μ g of total RNA to isolate mRNA. The purified mRNA was fragmented to 200-400 bases. The RNA fragments were then reverse transcribed into double stranded cDNA fragments. The DNA fragments were repaired to generate blunt ends using T4 DNA polymerase, Klenow polymerase and T4 polynucleotide kinase. After DNA fragments were purified using Qiagen PCR purification kit (Qiagen catalogue #28004), an "A" base was added to the 3′ end of the blunt DNA fragment by Klenow fragment. Sequencing adapters were ligated to the ends of DNA fragments using DNA ligase. The libraries were then amplified by limited PCR (15 cycles) using primers provided by the kit. The PCR products were then separated by 2% agarose gel electrophoresis and fragments with sizes ranging from 250 bp to 400 bp were excised and purified using the QIAquick Gel Extraction Kit (Qiagen catalogue #28704). The concentration and distribution of the library were determined by a NanoDrop spectrophotometer (Thermo Scientific, Wilmington, DE). The library was sequenced by Illumina HiSeq 2000 at Tufts University Genomic core facility. During analysis, we trimmed three nucleotides from the 5′ end. Raw sequencing reads (50 bp single-end) were mapped to human genome build hg19 using Tophat 2 53 . We used the default parameters (TopHat2) present in Galaxy for the mapping. The accepted hits were used for assembling transcripts and estimating their abundance using Cufflinks. The differentially expressed gene, promoter usage and splicing form were determined by Cuffdiff and Cuffcompare 54 . We selected genes keeping a cutoff value of 2 or more fold change plus a p value < 0.005. The heat maps and links were generated using Circos software 55 . The RNA-Seq data is now available in NCBI's GEO database (Accession# GSE83601).
Micro-RNA microarray. For the miRNA microarray, four controls and eight PTSD patients were included in the study. Microarray for the miRNAs was performed by Johns Hopkins Memorial Institute (Deep Sequencing and Microarray Core Facility), Baltimore. Total RNA, including mRNA, miRNA and other small RNA molecules, were isolated from PBMC samples as described above. Next, total RNA samples were used in the analysis of miRNA differential expression by miRNA array hybridization assay using the Affymetrix miRNA-v1 gene chip. Array data normalization and quality control was performed as described previously by Zhou et al. 56 . Linear fold-changes in miRNA up-regulation or down-regulation were calculated to compare the differences of all the miRNAs expressed between PTSD patients and controls. A linear fold-change of at least plus or minus 1.5 was used as a cut off value for the inclusion of a miRNA. Moreover, only the miRNAs which were significant on the basis of p value (less than or equal to 0.05) calculated using student's t test, were included for the analysis. We call these miRNAs differentially expressed. We used two tailed Student's t test to get the p values. The miRNA array data is now available in ArrayExpress (Accession# E-MTAB-4880).  Data analysis tools and functional enrichment of the genes. We employed Ingenuity Pathway Analysis (IPA, http://www.ingenuity.com/, QIAGEN, CA), Panther 57 and DAVID 58,59 for functional analysis of our datasets. IPA has tools to define interactions between miRNAs and target genes. It also has tools to identify the canonical pathways in which a given set of genes are involved. We also took advantage of the IPA Expression pairing tool to identify miRNAs and target genes from our datasets that are known or predicted to interact. Panther is a bioinformatics-based pathway analysis tool which can be used for functional enrichment of a given set of genes. Similarly, DAVID is also a bioinformatics-based tool to perform functional annotation or gene ontology identification and other categorization of a given set of genes. Both Panther and DAVID cannot be used for miRNA-gene interactions however, they both have an exhaustive list of functional annotation data. Functional enrichment of a set of genes can help to identify the biological pathways in which the genes in a dataset are involved and also other genes with which a certain gene interacts, thereby providing a more meaningful understanding of the data. Consequently, we performed functional enrichment analysis in IPA, Panther pathways and DAVID with the RNA-Seq data which included all the genes with a fold-change of 2 or more and p value less than or equal to 0.005 (we call these genes differentially expressed). The three bioinformatics tools are similar in that they search for evidence of enrichment of genes in particular list of genes. However, not all databases have  -miRs-145-5p, -15b-5p, -149-3p, -23b-3p, -30c-2-3p, -30c-1-3p, -342-3p   SCGB3A1  2.54  hsa-miRs-423-5p, -663a, -625-5p, -30e-3p   CXCL3  2.3  hsa-mirs-425-5p, let-7c-3p, -532-3p, -584-5p, -1207-5p, -132-3p, -181-a-5p, -181b-5p, -181c-5p, -181d-5p,  -150-5p, -194-5p   USP9Y  2  hsa-miRs-132-3p, -130b-3p, -130a-3p, -140-5p, -28-3p, -92b-3p, -92a-3p, -181a-5p, -181b-5p, -181c-5p,  -181d-5p, -23a-3p   CHDH  2 hsa-miRs-455-3p, -342-5p, -1231, -140-3p, -28-5p, -29b-2-5p, -324-3p, -505   complete information for all the genes and their interactions. Thus, a more complete analysis is obtained if the same dataset is analyzed using different tools.
Expression pairing of miRNAs and genes. In this part of the analysis, the differentially expressed genes and the miRNAs were analyzed together to see whether there are gene(s) that are targets of the miRNAs from our dataset obtained from similar samples in the present study. This type of analysis can help to simultaneously identify the differentially expressed genes and their regulators (miRNAs, in this case) in the same cells. To do this, we used a tool called Expression pairing, available in the Ingenuity's microRNA target filter analysis section. This tool is useful to determine the miRNA(s) and gene(s) that interact with each other from a given set of miRNAs and genes when provided simultaneously as input datasets. To perform expression pairing, both the miRNA array followed by the RNA-Seq dataset was analyzed in IPA simultaneously. The tool then paired miRNAs with the genes from the dataset based on Ingenuity's knowledge base for miRNAs and their targets. The genes obtained after expression pairing are either known or predicted to be a target of miRNA(s) from the list. Consequently, we got a list of genes from the input dataset and also miRNAs that could be their possible regulator for the observed differential expression. Subsequently, the genes obtained after expression pairing were used for another round of functional enrichment analysis for which, only IPA was used. This functional enrichment was performed to obtain the pathways in which the genes targeted by the miRNAs are involved.
Micro-RNA target gene digging. In our miRNA dataset, most of the miRNAs from PTSD patients were downregulated. Based on our experience, IPA did not cover all of the miRNAs and their targets. So, we performed miRNA-target search to identify additional differentially expressed miRNAs and genes in our dataset, which possibly could be targeted by the miRNAs. We used the publicly available website 60 , www.targetscan.org for performing this analysis. One at a time, the genes from our dataset were analyzed in www.targetscan.org. After obtaining the list of miRNAs that target the genes, we manually searched for the miRNAs that are present in our dataset. As we had a very long list of differentially expressed genes, we listed miRNAs for only the top up-or downregulated genes as per the RNA-Seq analysis.
DNA methylation analysis. For identifying differences in the DNA methylation level in specific genes, we used the publicly available Gene Expression Omnibus (GEO) datasets (GSE21282) 52   dataset, the DNA methylation of CpG sites was obtained from whole blood cells. The observations in the above mentioned dataset are different from ours in that our observations are collected from PBMCs rather than whole blood cells. However, using this dataset and our own earlier Methylated DNA immunoprecipitation (MeDIP) data (not provided), we could previously identify and report differential expression of IL12 as a result of alteration in DNA methylation and other epigenetic mechanisms in PBMCs from PTSD patients 7 . We obtained the average β values for the true CpG sites around TSS of specific genes from the datasets. For our purpose, we included all the β values without using a cutoff value. Moreover, in a case where there was more than one CpG site, we report here the β values of only the CpG sites with higher average β values. We compared the DNA methylation levels for all the differentially expressed genes from our RNA-Seq dataset which had true CpG islands and listed in the Illumina's probe list.
Quantitative Real Time PCR validation of the RNA-Seq data. For validation purpose of the RNA-Seq data, we selected seven genes for analysis. Complementary DNA (cDNA) was prepared from 0.5 μ g of total RNA using miScript RT II kit (Qiagen, Valencia, CA) in a 20 μ l system and used 7 ng of the original amount for qRT-PCR. Quantitative RT-PCR in triplicate wells was performed in an Applied Biosystem ViiA TM 7 Real-Time PCR system (Life Technologies, Carlsbad, CA) real time PCR instrument. As an internal control, 18S rRNA and GAPDH message was quantified along with the genes. The expression level of the genes was expressed as the relative abundance value to the internal control. For the validation of gene expression by employing qRT-PCR, we included 24 controls and 24 PTSD samples.