Differential and unique patterns of synaptic miRNA expression in dorsolateral prefrontal cortex of depressed subjects

Abstract

Altered synaptic plasticity is often associated with major depressive disorder (MDD). Disease-associated changes in synaptic functions are tightly correlated with altered microRNA (miRNA) expression. Here, we examined the role of miRNAs and their functioning at the synapse in MDD by examining miRNA processing machinery at synapse and sequencing miRNAs and analyzing their functions in synaptic and total tissue fractions obtained from dorsolateral prefrontal cortex (dlPFC) of 15 MDD and 15 matched non-psychiatric control subjects. A total of 333 miRNAs were reliably detected in the total tissue fraction. Multiple testing following the Benjamini-Hochberg false discovery rate [FDR] showed that 18 miRNAs were significantly altered (1 downregulated 4 up and 13 downregulated; p < 0.05) in MDD subjects. Out of 351 miRNAs reliably expressed in the synaptic fraction, 24 were uniquely expressed at synapse. In addition, 8 miRNAs (miR-215-5p, miR-192-5p, miR-202-5p, miR-19b-3p, miR-423-5p, miR-219a-2-3p; miR-511-5p, miR-483-5p showed significant (FDR corrected; p < 0.05) differential regulation in the synaptic fraction from dlPFC of MDD subjects. In vitro transfection studies and gene ontology revealed involvement of these altered miRNAs in synaptic plasticity, nervous system development, and neurogenesis. A shift in expression ratios (synaptic vs. total fraction) of miR-19b-3p, miR-376c-3p, miR-455-3p, and miR-337-3p were also noted in the MDD group. Moreover, an inverse relationship between the expression of precursor (pre-miR-19b-1, pre-miR-199a-1 and pre-miR-199a-2) and mature (miR-19b-3p, miR-199a-3p) miRNAs was found. Although not significantly, several miRNA processing enzymes (DROSHA [95%], DICER [17%], TARBP2 [38%]) showed increased expression patterns in MDD subjects. Our findings provide new insights into the understanding of the regulation of miRNAs at the synapse and their possible roles in MDD pathogenesis.

Fig. 1. miRNA expression volcano plot, heatmap, gene ontology prediction, and IPA analysis in total fraction.

Normalized values of 333 detectable miRNAs were plotted in an expression volcano plot and heatmap with hierarchical clustering. a miRNAs with high expression are shown with green color on the map; miRNAs with low expression are shown in red. For the clustering purposes, the Euclidean method was used to measure the distance, and the average linkage algorithm was applied to calculate the average pairwise distance between all pairs of points. Following the average linkage clustering algorithm, the dendrogram was constructed to demonstrate the expression similarities. b Volcano plot of genes differentially expressed between control and MDD groups. The y-axis corresponds to the significance level represented with log₁₀P value, and the x-axis displays the log2 (FC) value. The red dots represent the significantly (p ≤ 0.05) overexpressed genes in MDD (FC ≥ 1.3); The blue dots represent the significantly (p ≤ 0.05) under expressed genes (FC ≤ 1.3) in MDD; the green dots represent the genes whose expression levels did not reach statistical significance (p ≥ 0.05) but expression level was higher (FC ≥ 1.3) in MDD group. c Significantly altered miRNAs in the total fraction of dlPFC from MDD subjects. Out of 333 miRNAs, 4 miRNAs and 18 miRNAs were significantly up- or downregulated in MDD subjects, respectively. d GO analysis for biological process conducted with predicted target genes by significantly altered miRNAs showing significant enrichment of terms in various categories associated with neuronal functions. The lower p value is shown as blue color, and the circle size means the number of gene counts in each GO term. IPA analysis was performed with predicted target genes for significantly up- and downregulated miRNAs separately for canonical pathway (e) and disease and function (f). The results from up- and downregulated miRNA are shown as blue and orange colors respectively.

Fig. 2. miRNA expression heatmap and GO prediction in the synaptic fraction.

a Normalized values of 351 miRNAs were plotted in an expression heatmap with hierarchical clustering. miRNAs with high expression are shown with green color on the map whereas miRNAs with low expression are shown in red. For the clustering purpose, the Euclidean method was used to measure the distance, and the average linkage algorithm was applied to calculate the average pairwise distance between all pairs of points. Following the average linkage clustering algorithm, the dendrogram was constructed to demonstrate the expression similarities. b Significantly altered miRNAs in the synaptic fraction of dlPFC from MDD subjects. Out of 351 miRNAs, 8 miRNAs were significantly altered in MDD subjects (6 upregulated: miR-215-5p, miR-192-5p, miR-202-5p, miR-19b-3p, miR-423-5p, miR-219a-2-3p; 2 downregulated: miR-511-5p, miR-483-5p). c GO analysis for biological process showing significant enrichment of terms in various categories associated with neuronal morphogenesis, growth and differentiation. d GO-based functional network using cellular component as predictor, demonstrating significant enrichment of gene sets central to synaptic morphology, function and regulation. The prediction analysis was done following an FDR corrected p value cutoff 0.05 to determine the gene set enrichment in biological process and cellular component category separately. In each category the enriched terms were used to create networks where nodes are presented with terms and connected with edges. As shown in the graphs, if two nodes are connected, then they share 20% (default) or more genes. Bigger nodes represent larger gene sets. Thicker edges represent more overlapped genes.

Fig. 3. GO prediction based on 24 uniquely expressed miRNAs in synaptosomes.

a Table shows 24 miRNAs that were uniquely expressed in the synaptic fraction. b GO-based functional network using predicted targets of 24 miRNAs recruiting neuro-related terms. Closely connected networks presented with nodes and edges show enrichment of terms related to neuron projection, neuron development, neuron fate commitment and dendritic development. c GO-based functional network using predicted targets of 24 miRNAs recruiting transcriptional regulation related terms. Closely connected networks presented with nodes and edges showing enrichment of terms central to transcriptional regulation. A consensus list of predicted targets was used in standalone Cytoscape program to perform the GO analysis. Display pathways were selected with p values ≤ 0.05. Clustering was done based on the common functionality of genes enriched for specific term. The kappa score was set at 0.4 to define the term-term relationship (i.e., edges between the nodes). Genes involved in more than one function were represented with multiple color combination.

Fig. 4. Validation of target genes using in vitro cell model.

Relative quantification of target gene expression was done in SH-SY5Y neuroblastoma cell line transiently transfected with vehicle (n = 6) or hairpin inhibitors (n = 6). a Schematic diagram of in vitro study. b Mimic miR-19b-3p overexpression oligo (n = 6). Overall group differences in the three groups are as follows: CISD: df = 2; f = 35.2, p < 0.01; CHST7: df = 2; f = 9.9, p = 0.002; CHP1: df = 2; f = 0.8, p = 0.451; CYB56D1: df = 2; f = 1.1, p = 0.37; FUT9: df = 2; f = 0.5, p = 0.61; N6AMT1: df = 2; f = 8.1, p = 0.004; SEL1L3: df = 2; f = 17.8, p = 0.001. c Mimic miR-483-5p overexpression oligo (n = 6). Overall group differences in the three groups are as follows: CCDC9: df = 2; f = 7.3, p = 0.006; CX3CL1: df = 2; f = 10.8, p = 0.001; C5AR1: df = 2; f = 22.7, p < 0.001; FOXO3: df = 2; f = 12.6, p = 0.001; ELK1: df = 2; f = 23, p = 0.001; HBGEF: df = 2; f = 0.5, p = 0.60; IRF1: df = 2; f = 2.6, p = 0.10; NFAM1: df = 2; f = 3.6, p = 0.05; MAP2K3: df = 2; f = 16.9, p < 0.001; TMEM98: df = 2; f = 4.5, p = 0.03. d Mimic miR-511-5p overexpression oligo (n = 6). Overall group differences in the three groups are as follows: CD68: df = 2; f = 0.6, p = 0.54; DISC1: df = 2; f = 1.1, p = 0.355; ELK1: df = 2; f = 23, p < 0.001; IL17RA: df = 2; f = 0.6, p = 0.57; IRF2: df = 2; f = 0.4, p = 0.65; PHLDB1: df = 2; f = 8.1, p = 0.004; TAB2: df = 2; f = 17.8, p = 0.001. The average differences of target gene expression among vehicle, mimic, and hairpin inhibitor were assessed by one-way ANOVA with post hoc Bonferroni correction. GAPDH was used as normalizer for gene expression. Values denote average ± SEM. ‘a’ and ‘b’ denote statistical significance ‘between vehicle and mimic’ and ‘in mimic compared to vehicle and hairpin inhibitor’, respectively. Inhibitor, hairpin inhibitor.

Fig. 5. Pre-/mature-miRNA expressions in synaptosomes.

Scatter plots represent the relative quantification of mature and their respective precursor miRNAs in synaptosomes. The differences between two groups are as follows: miR-19b-3p: df = 8; t = −0.949, p = 0.35; pre-miR-19b-1: df = 19; t = 0.57, p = 0.57; pre-miR-19b-2: df = 21; t = −0.79, p = 0.43; miR-199a-3p: df = 26; t = 0.69, p = 0.497; pre-miR-199a-1: df = 22; t = −2.22, p = 0.036; pre-miR-199a-2: df = 26; t = −1.62, p = 0.11; miR-455-3p: df = 26; t = 0.65, p = 0.51; pre-miR-455: df = 21; t = 2.223, p = 0.037; miR-215-5p: df = 27; t = 0.37, p = 0.71; pre-miR-215: df = 21; t = 2.33, p = 0.03. The average differences were assessed by student’s t test. U6 and geometric means of GAPDH, ACTB, and ribosomal 18S RNA were used as normalizers for mature and precursor miRNAs, respectively. (n = 15/group). Ct, control; miRNA, microRNA; MDD, major depressive disorder.