Prognostic plasma exosomal microRNA biomarkers in patients with substance use disorders presenting comorbid with anxiety and depression

Psychiatric disorders such as anxiety and depression precipitated by substance use occurred during both use and withdrawal. Exosomes play significant roles in biological functions and regulate numerous physiological and pathological processes in various diseases, in particular substance use disorders (SUDs) and other psychiatric disorders. To better understand the role of exosomal miRNAs in the pathology of symptoms of anxiety and depression in patients with SUDs, we first isolated circulating exosomes from heroin-dependent patients (HDPs) and methamphetamine-dependent patients (MDPs) and identified exosomal miRNAs that were differentially expressed between patients and healthy controls (HCs). Furthermore, the correlations between exosomal DE-miRNAs and symptoms of anxiety and depression which were measured using Hamilton-Anxiety (HAM-A)/Hamilton-Depression (HAM-D) Rating Scales in the participants. Notably, the expression level of exosomal hsa-miR-16-5p, hsa-miR-129-5p, hsa-miR-363-3p, and hsa-miR-92a-3p showed significantly negative correlations with HAM-A scores in both HDPs and MDPs. But all of the 4 DE-miRNAs lost significant correlations with HAM-D scores in HDPs. Functional annotation analyses showed that the target genes of the DE-miRNAs were mainly enriched for “synapse”, “cell adhesion”, “focal adhesion” and “MHC class II protein complex”. Our study suggests that a set of circulating exosomal miRNAs were associated with anxiety and depression in SUD patients and may have clinical utility as diagnostic and prognostic biomarkers.


Results
Demographics and hematological parameters of the study participants. The demographics of the study participants are displayed in Table 1, the age, gender and years of education of the MDPs and HDPs were matched to those of the HCs. A total of 20 hematological parameters were examined (Supplementary Table S1). Compared with the HCs, four of these parameters in both HDPs and MDPs, related to white blood Table 1. Demographics characterisristics of study participants. M mean, SD standard deviation, MDPs methamphetamine-dependent patients, HDPs heroin-dependent patients, HCs healthy controls, BMI Body Mass Index.
Behavioral characteristics of the study participants. Data for all three behavioral scales aforementioned were collected for all the study participants ( Table 2). The degree of drug craving was assessed by VAS and the scores among the three age-, gender-, and education-matched groups are shown in Table 2. Patients with SUDs were associated with higher VAS scores (MDPs: p = 2.19e−03; HDPs: p = 5.5e−03) compared with HCs. No significant differences were found between MDPs and HPDs (p > 5e−03). As shown in Table 2, the intergroup analysis of the HAM-A scale revealed that both MDPs and HDPs presented significantly higher scores than the HCs (both p < 1e−02). A comparison of the median scores of the HAM-D scale among the groups indicated that significant differences existed between the HCs and HDPs (p = 8.89e−03), but not between the HCs and HDPs (p > 5e−02). No differences in HAM-A or HAM-D scale scores were observed between the MDPs and HDPs (p > 5e−02).
Characteristics of exosomes obtained from the plasma of the study participants. We performed Nanoparticle Tracking Analysis (NTA), Transmission Electron Microscopy (TEM), and immunoblot analysis of protein markers specific for plasma-derived EVs to confirm that the nanoparticles we obtained were indeed exosomes. TEM showed that exosomes from HCs, HDPs, and MDPs all exhibited a cup-shaped morphology with diameters ranging from 50-150 nm (Fig. 1B), the median diameters were 99.4 ± 26.35, 100.8 ± 23.5, 99.5 ± 25.14, respectively. And the representative exosomal diameter was shown in Fig. 1A. There are no differences between health control and SUDs in both exosome size distribution ( Figure S2A) and numbers (Fig-Table 2. Behavioral characteristics of study participants. M mean, SD standard deviation, VAS Visual-Analogue Craving Scale, MDPs methamphetamine-dependent patients, HDPs heroin-dependent patients, HCs healthy controls. HCs (M ± SD, n = 10) HDPs (M ± SD, n = 10) Adjusted p value (HC vs MDPs) MDPs (M ± SD, n = 10) www.nature.com/scientificreports/ ure S2B). Exosomes were also identified by western blotting for three commonly used exosomal markers (CD63, TSG101, and ALIX), and calnexin was used as a negative control (Fig. 1C). All these data indicated that the nanoparticles obtained were indeed exosomes.
Overview of small RNA sequencing data and small RNA biotype mapping. A total of 630,546,931 raw reads were obtained after sequencing (Supplementary Table S1), which were generated from 30 small RNA libraries (10 replicates per group) derived from HCs, HDPs, and MDPs. After removing contaminant reads from the original data, a total of 240,754,247 clean reads were obtained from the 30 sequencing samples (Supplementary Table S2). Reads obtained from sequencing were used for alignment and mapping to the human genome after adapter clipping and quality filtering. A total of 1299 distinct miRNAs were identified from all three groups from the combined raw read count data (Supplementary Table S2). The clean reads were then mapped to noncoding RNA databases, including the Rfam and miRBase databases, to annotate the classifications of the small RNAs. As shown in Supplementary Figure S3A, miRNAs were the most abundant exosomal RNA species in all the groups (64.27 ± 5.88%, 70.91 ± 5.28%, and 73.55 ± 4.82% in HCs, HDPs, and MDPs, respectively).
Length distribution analysis of the miRNAs showed that 20-24 nucleotides (nt) was the most frequent read length, with the peak at 22 nt, indicating that mature miRNAs were well enriched during the sequencing library preparation process (Supplementary Figure S3B). After identification of the conserved miRNAs, additional filtering was performed using miRDeep2 to identify potentially novel miRNAs. A total of 133 novel miRNAs were detected across the 30 libraries (Supplementary Table S2), and all displayed typical miRNA features at the genomic level. A representative readout of the predicted novel miRNAs is shown in Supplementary Figure S3C.

DE-miRNAs and RT-qPCR validation in HCs and
SUDs. Differential expression analysis was performed using linear contrast in the DESeq2 RNA-seq, and controlled for 5% false-discovery rate (FDR) using the Benjamini-Hochberg method for each pairwise comparison. MiRNAs with an adjusted p-value < 5e−02, as determined by DESeq2, were assigned as being differentially expressed. We found 34 differentially expressed miRNAs between HCs and MDPs, 7 of which were significantly upregulated and 27 significantly downregulated ( Fig. 2A, Supplementary Table S3a). Nineteen miRNAs were identified as DE-miRNAs between HCs and HDPs ( Fig. 2A), 7 of which were markedly upregulated and 12 significantly downregulated (p < 5e−02) (Supplementary Table S3b). We also compared miRNAs expression between MDPs and HDPs, and identified 16 miRNAs with p value < 5e−02, but lost significance after FDR-adjusted. Fifteen miRNAs were identified as being differentially expressed in both the HCs vs MDPs and HCs vs MDPs comparisons.
Hierarchical clustering analysis was performed for these 15 DE-miRNAs ( Fig. 2B) while volcano plots were used to assess the overall distribution of these DE-miRNAs in HCs vs HDPs (Fig. 2C) and HCs vs MDPs (Fig. 2D) comparisons.
Next, we recruited validation sets and used 30 samples for each group to isolate exosomes and validated the selected miRNAs identified by RNA-seq. We selected five of the DE-miRNAs for verification by qPCR assay, using 5′ nuclease probes, and normalized with U6 ( Fig. 3), details of primer sequence were shown in Supplementary  Table S4. Four of these showed a consistent trend with the RNA sequencing results for the MDP samples, and all were downregulated in MDPs when compared with HCs: hsa-miR-143-3p (Fig. 3A), hsa-miR-200a-3p (Fig. 3B), hsa-miR-363-3p (Fig. 3C), and hsa-miR-125b-5p ( Fig. 3D) (all p < 5e−02); the change in hsa-miR-141-3p expression lost statistical significance in the qPCR validation (p > 5e−02). Three of the five hsa-miRNAs showed a consistent trend with the RNA sequencing results for the HDP samples: hsa-miR-143-3p, hsa-miR-363-3p, and hsa-miR-125b-5p (all p < 9.0e−04); the changes in hsa-miR-200a-3p and hsa-miR-141-3p expression level lost statistical significance in the qPCR validation (p > 5.0e−02).  Hierarchical clustering analysis of miRNAs with altered expression among the three groups (p < 5e−02, fold change > 2). Red strip, high relative expression; blue strip, low relative expression; white strip, no change in gene expression. Color intensity reflects the degree of expression increase or decrease. (C) Volcano plot showing DE-miRNAs from HCs vs MDPs with various fold changes and p-values. Vertical line, fold change = 2 (log2 transformed); horizontal line, p = 5e−02 (−log10 transformed). Red dots, p < 5e−02, fold change > 2; green dots, change fold > 2, p < 5e−02; blue dots, p > 5e−02, fold change > 2; gray dots, insignificantly changed miRNAs. (D) Volcano plot showing DE-miRNAs from HCs vs HDPs with various fold changes and p-values. Vertical line, fold change = 2 (log2 transformed); horizontal line, p = 5e−02 (−log10 transformed). Red dots, p < 5e−02, fold change > 2; green dots, change fold > 2, p < 5e−02; blue dots, p > 5e−02, fold change > 2 ; gray dots, insignificantly changed miRNAs. MDPs Exosomes obtained from methamphetamine-dependent patients, HDPs Exosomes obtained from heroin-dependent patients, HCs Exosomes obtained from healthy controls, SUDs substance use disorders. Potential regulatory roles of the identified DE-miRNAs. We then performed GO and KEGG enrichment analyses to predict overlapping target genes and analyze the potential pathway involved in SUD-related anxiety and depression. For the comparison between HCs and HDPs (Fig. 5A), GO analysis indicated that the predicted targets were mainly enriched in nervous system (synapse and sympathetic nervous system development), cell mobility and proliferation (cell junction; cell adhesion; facal adhesion), and immune system (MHC class II protein complex; TAP complex) (Fig. 5B). The top 20 enriched GO terms of biological process, cellular component, and molecular function are displayed in Fig. 5C, which were predicted to play crucial roles in single-organism process, biological regulation, metabolic process. KEGG pathway analysis indicated that axon guidance, antigen processing and presentation, focal adhesion, and Notch signaling pathway were enriched in both MDPs and HDPs, implying that nervous system, cell mobility, and immune system werepresented among the enriched KEGG terms (Fig. 5D).  Table S5a) between HCs and HDPs were connected by their respective mRNA targets, illustrating the network influenced by significantly differently regulated miRNAs (Fig. 6A). The respective mRNA targets of the top 10 DE-miRNAs   Table S5b) between HCs and MDPs are depicted in Fig. 6B. The targeted 5299 genes were carried out overlapping with the aforementioned database on DE-miRNAs between HCs and patients with SUDs.

Discussion
In the present study, we compared the expression profiles of plasma exosomal miRNAs derived from SUD patients and HCs. The expression level of 34 and 19 miRNAs were significantly altered in MDPs and HDPs, respectively, when compared with those of HCs. To further explore the functions of DE-miRNAs identified, GO and KEGG pathway analysis were performed for the predictive target genes. The results showed that target genes of the DE-miRNAs were mainly enriched in nervous system, cell mobility and proliferation, and immune system both in MDPs vs HCs and HDPs vs HCs. There are 16 DE-miRNAs were identified between MDPs and HDPs, but all these 16 miRNAs lost significant difference when adjusted by FDR. Among the DE-miRNAs identified, part of them were reported to be involved in anxiety and depression. MiR-144-5p and miR-16-3p are involved in the response to mood stabilizer treatment 41 and stress responses 42 . www.nature.com/scientificreports/ MiR-451a was altered in both serum and CSF from MDD patients 43 . And some other DE-miRNAs have been reported to be associated with neurological diseases, such as miR-181a and miR-484 were reported as mitomiRs involved in modulating mitochondrial function 44 , which played a vital role in neurodegenerative and psychological disease 45,46 .
We then compared the symptoms of anxiety and depression between SUD patients and HCs using HAM-A and HAM-D measures, and found that the HDPs had significantly higher scores in both the HAM-A and HAM-D measures, but MDPs lost statistical significance in the HAM-D. Notably, we also found statistically significant correlations between some miRNAs (hsa-miR-92a-3p, hsa-miR-363-3p, hsa-miR-16-5p, and hsa-miR-129-5p) and the total HAM-A scores, even after adjusting for years in education, in an ordinal logistic regression analysis. Hsa-miR-92a-3p, hsa-miR-363-3p, hsa-miR-16-5p, and hsa-miR-129-5p showed a significant correlation with total HAM-D scores in MDPs, but not in HDPs. It's worth noting that the expression of all the 4 DE-miRNAs aforementioned in MDPs displayed significantly negative correlations to HAM-D scores, while the HAM-D scores in the MDPs showed no difference when comparing to HCs. These findings implied that these four identified DE-miRNAs may have independent function in different type of SUDs, and have potential as biomarkers to assess the severity of anxiety and depression in SUD patients who are not diagnosed by scales.
Researches have shown that many of the targeted genes, such as those associated with neurotransmission/ neuromodulation, including GABAergic, glutamatergic, dopaminergic, and serotonergic signaling, or with HPA axis functioning, are known to be aberrantly regulated in anxiety and depression 47 . Our study focused on exosomal miRNA profiles of SUD patients. We found that those 4 DE-miRNAs aforementioned were closely related to anxiety and depression. Gheysarzadeh et al. reported that the serum level of miR-16 was significantly downregulated in patients with anxiety and depression 42,48 , which was in agreement with our results. Shao et al. found the CSF miR-16, which has been suggested to be a biomarker for anxiety 49 , has a role in the pathogenesis of depression via affecting raphe SERT expression 50 , and play a role in the development of depressive-like behaviors in the hippocampus by targeting and downregulating BNDF expression 51 via PI3K/Akt/mTOR pathway 52 . Furthermore, miR-16-5p was predicted to target to BDNF, NPY4R, GLUD1, and FKBP5, which are reported to be involved in the pathophysiology of anxiety and depression [53][54][55][56][57] . MiR-129-5p was predicted to negatively regulate neuropeptide FF receptor 2 (NPFFR2), which has been shown to activate the HPA axis and trigger anxiety-and depressive-like behaviors 58 . Nonetheless, to date, no systematic study has been carried out on circulating exosomal miRNA expression profiles of patients with SUDs and differential exosomal miRNAs in SUD patients may provided a plausible regulatory mechanism how peripheral molecular affecting central nerve system and behaviors.
In the present study, the circulating exosomal miRNA profiles of both MDPs and HDPs were investigated for the first time, leading to the identification of DE-miRNAs in MDPs and HDPs presenting with comorbid anxiety and depression. We explored the possible mechanism underlying the comorbidity of anxiety and depression associated with SUDs. Notably, GO and KEGG analyses consistently showed that these target genes were enriched for categories related to psychological disorder-associated processes, cell mobility, and the immune system. Molecular enrichment results have revealed that serine/threonine-protein kinases are involved in depressivelike behavior through the phosphorylation of eIF4E (eukaryotic initiation factor 4E) and consequent regulation of 5-HT neurotransmission 59 . Several cell adhesion molecules (CAMs), such as Np65 and NLGN-4, have been identified at neuromuscular junctions and synapses, where they recruit scaffolding proteins, neurotransmitter receptors, and synaptic vesicles 60,61 . They were indicated to be involved in regulating anxiety-and depressive-like behaviors, and proposed to have potential as transdiagnostic biomarkers. A different study showed that anxiety disorders are mediated through the modulation of miRNAs associated with the regulation of genes involved in axonal guidance 62 , a term that was also enriched in our biological analysis. DNA-Binding/Differentiation proteins, which was enriched in molecular analysis, have been shown to be connected with depression 63 .
Combined, our findings suggest that the identified DE-miRNAs may contribute to the biological functions of SUDs-derived exosomes via the regulation of the relevant pathways. Nevertheless, several limitations should be considered when interpreting our results. First, the power of this study was limited because of the relatively small sample size. To validate the sensitivity and specificity of the miRNAs as biomarkers for SUD-associated anxiety and depression, further research should be carried out with a larger sample size and in different centres. Second, the cross-sectional design of the study does not allow causal interpretations. Further investigation of DE-miRNAs and target pathways in longitudinal clinical studies or animal studies will help to validate the effects of these biomarkers.

Conclusion
We identified exosomal miRNAs differentially expressed in patients with SUDs. The dysregulated miRNAs might be involved in the underlying pathophysiology of SUDs through biological pathways. Moreover, our results highlighted the importance of exosomal miRNAs as potential biomarkers for SUD patients suffering from depression and anxiety. Four DE-miRNAs were associated with anxiety in both MDPs and HDPs, but only associated with depression in MDPs, suggesting that these DE-miRNAs may be potential biomarkers for the diagnosis and treatment of anxiety and depression in SUDs.

Materials and methods
Ethics statement. All protocols and recruitment procedures described in this study were approved by the www.nature.com/scientificreports/ Edition of the Diagnostic and Statistical Manual of Mental Disorders (DSM-5) criteria. Healthy individuals were recruited from community sites. All the participants were male smokers with HIV(−), HBV(−) and HCV(−), then they were divided into a test set and a validation set. In the test set, 10 HCs and 20 gender-and age-matched SUD patients (10 HDPs and 10 MDPs) were recruited, poly drugs abusers were excluded. Thirty participants in each group were validated in the validation set.
Sample collection and preparation. The demographics and substance use characteristics of the participants included in the study are presented in Table 1. Blood samples were collected from fasted SUD patients and HCs at 08:00 ~ 10:00 AM and using vacutainer blood collection tubes containing EDTA anticoagulant. The anticoagulant-treated blood samples were mixed by inverting the tube several times. Blood samples (10 ml were then centrifuged at 3000 rpm for 10 min at 20 °C. The upper layer containing plasma was transferred to a 2-ml eppendorf tube and centrifuged at 3000 × g for 15 min at 4 °C. Plasma from each sample was collected into a 4-ml eppendorf tube, immediately placed on dry ice, and then stored at − 80 °C until use. Scale administration. The Visual-Analogue Craving Scale (VAS) is used to measure the degree of drug craving using a 10 point scale, with 0 indicating "none at all" and 10 indicating "very much". The scale records the subjective craving experienced by the study subjects at a specific time.
The HAM-D scale is a 24-item questionnaire administered by an interviewer that measures the severity of depressive symptoms. A total score above 20 is considered indicative of major depression 65 .
Isolation of plasma exosomes. The exosomes was isolated by SEC (size exclution chromatography) methods as described previously with minor modifications 66  After incubation with a specific secondary anti-mouse or -rabbit horseradish peroxidase-conjugated antibody (1:5000 dilution; GeneTex, USA) at room temperature for 1 h, the protein bands were detected using a chemiluminescence detection kit (Millipore Co., MA, USA) and scanned using the iBright FL1500 chemiluminescence imaging system (Thermo Fisher Scientific, USA).
Transmission electron microscopy. For transmission electron microscopy (TEM) analysis, 10 μl of exosomal suspension was loaded on a copper mesh, incubated for 10 min at room temperature, and washed with distilled water. After adsorption, each sample was negatively stained with 10 μl of a 2% (w/v) uranyl acetate solution at room temperature for 1 min, and excess liquid was blotted with filter paper. Then, the samples were observed and photographed using a Hitachi H-7650 transmission electron microscope (Hitachi, Tokyo, Japan) operating at 80 kV. Nanoparticle tracking analysis. Vesicle suspensions with concentrations between 1 × 10 7 -10 9 /mL were examined using the ZetaView PMX 110 (Particle Metrix, Meerbusch, Germany) equipped with a 405 nm laser to determine the size and quantity of particles isolated 67 . A video of 60-s duration was taken with a frame rate of 30 frames/sec, and particle movement was analyzed using NTA software (ZetaView 8.02.28, https:// www. parti cle-metrix. de).

RNA preparation and library construction.
Total RNA in the exosomes was extracted using the Qiagen miRNeasy Mini kit (Qiagen, CA, USA) following the manufacturer's protocol. The concentration and quality of the RNA were determined by Agilent 2100 Bioanalyzer. Small RNA libraries were prepared using the QIAseq miRNA Library Kit (Qiagen, Frederick, MD) according to the manufacturer's protocols. For each library, 1-500 ng of total RNA from each sample was used in all experimental procedures and index codes were added to attribute sequences to each sample. At last, library quality was assessed on the Agilent Bioanalyzer 2100 and qPCR. The clustering of the index-coded samples was performed on acBot Cluster Generation System using TruSeq PE Cluster Kitv3-cBot-HS (Illumina, San Diego, CA, USA) according to the manufacturer's instructions. After cluster generation, the library preparations were sequenced on an Illumina Hiseq platform and paired-end reads were generated. www.nature.com/scientificreports/ Quantification and differential expression analysis of miRNA. Use Bowtie tools soft, The Clean Reads respectively with Silva database, GtRNAdb database, Rfam database and Repbase database sequence alignment, filter ribosomal RNA (rRNA), transfer RNA (tRNA), small nuclear RNA (snRNA), small nucleolar RNA (snoRNA) and other ncRNA and repeats. The remaining reads were used to detect known miRNA and new miRNA predicted by comparing with known miRNAs from miRbase and Human Genome (GRCh38), respectively. Read count for each miRNA was obtained from the mapping results, and Transcripts Per Million (TPM) was calculated.
Prediction and functional analysis of miRNA-targeted genes. The mRNA targets of the DE-miRNAs were predicted using two databases (miRanda and RNAhybrid). DE-miRNAs were analyzed using a Venn diagram to identify overlapping genes. Only the miRNA-target gene interactions that showed overlap in two databases (miRanda and RNAhybrid) aformentioned were considered candidate targets. Gene-gene and a miRNA-target gene interaction networks were visualized using Cytoscape 3.0. The shared DE-miRNA-target interactions were also visualized using Cytoscape 3.0 (https:// cytos cape. org). GO and KEGG (www. kegg. jp/ kegg/ kegg1. html) pathway enrichment analyses were performed to assess the potential functions of the identified target genes using topGO 3. and GraphPad Prism 8.0 (GraphPad Software, San Diego, CA, USA, https:// www. graph pad. com). Quantitative data were expressed as means ± standard deviation (SD). Exosome concentrations were analyzed using the one-way ANOVA among three groups. An adjusted p value < 5e−02 was considered to indicate a statistically significant difference.