This study was conducted at the Guangdong Provincial People’s Hospital from August 2020 to January 2021. It involved nine patients with ischemic stroke within 4 h of onset and two healthy individuals as a healthy control group, all of whom were recruited for plasma exoLR sequencing. The other 10 healthy individuals from our previous study, whose data are available from the Gene Expression Omnibus database (Accession No. GSE159657, https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE159657), were used as part of the healthy control group, and the expression profiles were analyzed [14].

The inclusion criteria for the ischemic stroke group were patients aged more than 18 years and with a definite diagnosis of ischemic stroke. The diagnosis was performed according to World Health Organization guidelines and validated with computed tomography or magnetic resonance imaging [15]. Conversely, individuals in the healthy control group were required to display normal renal and liver function, with no history of malignancy, rheumatologic disorders, recent cardiovascular or cerebrovascular events, chronic heart failure, diabetes, acute or chronic infectious disease, myocarditis, pulmonary embolism, aortic dissection, or pericarditis. The exclusion criteria for both groups included patients older than 80 years or diagnosed with acquired immunodeficiency syndrome. The study was conducted according to the guidelines of the Declaration of Helsinki and approved by the ethics committee of the Guangdong Provincial People’s Hospital (approval No. GDREC2019443H). Written informed consent was obtained from all the study participants.

The blood samples were taken within an hour of the patients’ hospital admission and prior to thrombolytic therapy. Venous blood (2 mL) was collected from each patient in tubes containing ethylenediaminetetraacetic acid. For extracting plasma, the blood was centrifuged for 10 min at 3000 g and 4 °C, and the supernatant was carefully transferred to a new tube for repeated centrifugation. Then, the plasma was stored in cryogenic vials at –80 °C until further use. Total RNA in exosomal vesicles was extracted using an exoRNeasy serum/plasma kit (77023, QIAGEN, Hilden, Germany) following the manufacturer’s protocol [16]. The RNA-sequencing libraries were constructed using SMART technology (Clontech Laboratories, Mountain View, CA, USA) [17]. RNA sequencing was performed by Guangzhou Epibiotek Co., Ltd. (Guangzhou, China) on an Illumina Nova-Seq 6000 System. The filter sequences and adapters were trimmed using Cutadapt (version 2.5, Technical University of Dortmund, Ruhr region, Germany). Using the “-rna-strandness RF” parameter, the remaining reads were aligned against the human Ensemble genome GRCh38 [18]. Based on featureCounts (version 1.6.3, The University of Melbourne, Victoria, Australia), reads were mapped to the genome [19]. Messenger RNAs and lncRNAs were annotated using the GENCODE database (version 28, https://www.gencodegenes.org/) [20]. R software (version 4.1.2, R Foundation for Statistical Computing, Vienna, Austria) was used to generate principal component analysis (PCA) plots from fragments per kilobase per million mapped fragments. Unmapped reads were circular RNAs identified using the pipeline (https://code.google.com/p/acfs/) [21]. The Gene Expression Omnibus database (Accession No. GSE186844, https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE186844) contains raw data that supports the results of this study.

The morphological characteristics of the exosomes were assessed using transmission electron microscopy (JEM-1200EX, Japan Electron Optics Laboratory, Tokyo, Japan). The sizes and distribution of the exosomes were measured using a NanoSight NS500 instrument developed by NanoSight Ltd. (Amesbury, UK). Exosomal tumor susceptibility gene 101 (TSG101), CD9, and CD63 protein markers were detected using western blotting [22].

The expression levels of plasma exoLRs were estimated based on high-throughput sequencing through differential expression analysis. The variance percentile rank was filtered at 15% or lower, and at least four counts per million were required to retain an RNA. As a result of filtering the data, the counts were normalized using the M-value normalization method. The COMBAT algorithm was then used to remove the batch effect. The raw counts from the sequence datasets were analyzed using the “DESeq2” R package (version 1.16.1), with design models considering batch effects [23], to obtain differentially expressed exoLRs. A log₂ fold change $\ge$$|$1.0$|$ at a statistical significance threshold of p 0.05 was applied. These steps decreased the influence of batch effects on experimental results. The expression heat maps were generated using the “pheatmap” package (version 1.0.10) in R.

CIBERSORT (version 3.6.2) was used to quantify immune cell abundance in plasma across study groups. The standardized gene expression data were uploaded to a publicly accessible online database (http://cibersort.stanford.edu/). An analysis of 547 marker genes was performed to determine the plasma distribution of 22 human immune cells [24].

We predicted the targeted transcripts of differentially expressed long noncoding RNAs (DElncRNAs) based on antisense, cis-acting, and trans-acting analyses to explore the roles of exosomal lncRNA in ischemic stroke. For antisense DElncRNA analysis, RNAplex (version 2.5.0, https://www.tbi.univie.ac.at/RNA/ViennaRNA/doc/html/man/RNAplex.html) [25] was used to predict the complementary correlations between antisense DElncRNAs and differentially expressed messenger RNAs (DEmRNAs). For cis-acting analysis, lncRNAs located less than 100 kb upstream/downstream of a gene were considered potential cis-regulators. For lncRNA trans-acting analysis, the correlation coefficient between DElncRNAs and DEmRNAs was determined, and a value $>$0.95 was considered to indicate trans-action.

Gene Ontology (GO) enrichment analysis for biological processes regulated by the differentially expressed exoLRs, including the DEmRNAs, target transcripts of DElncRNAs, parental genes of differentially expressed circular RNAs (DEcircRNAs), and DEmRNAs in the competing endogenous RNA (ceRNA) networks, was performed using the ClusterProfiler package (version 3.10.1, https://bioconductor.org/packages/release/bioc/html/clusterProfiler.html) [26] and Metascape (version 3.5, https://metascape.org/gp/index.html) [27]. As previously described, the pathway enrichment analysis of the Kyoto Encyclopedia of Genes and Genomes (KEGG) was used to identify the pathways regulated by differentially expressed exoLRs [28], with a statistical significance threshold of p 0.05.

Based on the ceRNA hypothesis, ceRNA networks of DElncRNA/DEcircRNA-microRNA-DEmRNA were built by predicting microRNA-binding RNA. Thus, micro RNAs regulated by DElncRNAs and DEmRNAs were predicted and selected to construct the ceRNA network. First, the microRNAs targeted by DElncRNAs/DEcircRNA were predicted using the MiRcode database (http://www.mircode.org) [29]. The potential messenger RNAs targeted by the micro RNAs were retrieved from miRDB (http://www.mirdb.org/index.html) [30] and TargetScan (http://www.targetscan.org/vert_71/) [31]. The potential messenger RNAs were intersected with the entire dataset of DEmRNAs. DElncRNAs/DEcircRNAs and targeted DEmRNAs in the opposite expression pattern were removed from the ceRNA network. Cytoscape (version 3.9.0, Institute for Systems Biology, Washington, WA, USA) [32] was applied to visualize the constructed network.

For analyzing the clinical characteristics, continuous variables were expressed as the median (interquartile range). A PCA plot was generated using R software by transforming the data into a log2 scale and plotting it using the plot PCA function. Count variables were analyzed using Fisher’s exact probability test, whereas continuous variables were analyzed using the Kruskal-Wallis rank-sum test. A p value 0.05 indicated a statistically significant difference. All data were analyzed in R software (version 4.1.2, R Foundation for Statistical Computing, Vienna, Austria).

The clinical characteristics of the study participants are summarized in Table 1. In the sequencing group, no differences were observed in age, alcohol consumption, lipid profiles, white blood cell counts, neutrophils (%), and C-reactive protein levels between patients with ischemic stroke and healthy individuals. In the validation group, no differences were observed in sex, lipid profiles, and C-reactive protein levels. Overall, the proportions of participants with hypertension, diabetes, or smoking habits were higher in the ischemic stroke group. Furthermore, an increase in monocyte count was also observed in the ischemic stroke group.

Transmission electron microscopy indicated the membrane-enclosed structures of exosomes (Fig. 1a). The original transmission electron microscope image is shown in Supplementary Fig. 1. The average diameter of the isolated exosomes was 79.11 $\pm$ 18.55 nm (Fig. 1b). Western blotting verified the presence of TSG101, CD9, and CD63 membrane markers on the exosomal membrane (Fig. 1c). The raw data for all western blots are provided in Supplementary Fig. 2. The messenger RNAs (protein coding), pseudogenes, circular RNAs, lncRNAs, and antisense RNAs constituted approximately 68%, 8%, 7%, 6%, and 7% of the total mapped reads, respectively (Fig. 2a). Based on the screening criteria of log₂ fold change $\ge$$|$1.0$|$ and p 0.05, 508 exosomal messenger RNAs consisting of 321 up- and 187 downregulated RNAs (Supplementary Table 1), 40 exosomal lncRNAs consisting of 31 up- and 9 downregulated RNAs (Supplementary Table 2), and 115 exosomal circular RNAs consisting of 67 up- and 48 downregulated RNAs (Supplementary Table 3) were screened as differentially expressed exoLRs. The results of hierarchical cluster analysis revealed that the expression patterns of messenger RNAs, lncRNAs, and circular RNAs in circulating plasma were distinguishable (Fig. 2b–d). The abundance of exosomal circular RNAs and lncRNAs was relatively low, and their expression levels showed individual differences even within the same group.

Fig. 1.

Characterization of exosomes isolated from plasma. (a) Representative transmission electron microscopy image showing the morphology of exosomes isolated from the plasma and vesicle size (scale bar = 100 nm). (b) NanoSight nanoparticle tracking analysis was conducted to determine the size distribution of isolated exosomes. Most isolated extracellular vesicles had a size distribution consistent with exosomes (30–100 nm). (c) Exosome marker proteins (CD63, TSG101, and CD9) detected by western blotting (S1 to S3: ischemic stroke specimens; H1 and H2: control specimens). TSG101, tumor susceptibility gene 101.

Fig. 2.

Overview of the expression profile of exoLRs. (a) Distribution of mapped reads to genes, with annotations showing the proportions of different types of RNA. (b–d) Heat maps showing clustered expression patterns of differentially expressed exosomal messenger RNAs (b), long noncoding RNAs (c), and circular RNA (d) between normal and ischemic stroke specimens. In the heat maps, light yellow indicates low expression, whereas brown indicates higher expression. Abbreviations: circRNA, circular RNA; DE, differentially expressed; lncRNA, long noncoding RNA.

Based on the exosomal sequencing data, we assessed the proportions of 22 circulating immune cells using an established computational resource (CIBERSORT) (Fig. 3a). The proportions of exhausted and gamma delta T cells were statistically higher in patients with ischemic stroke (Fig. 3b). PCA of the immunological profiles showed a nonuniform distribution (Fig. 3c). The expression of exoLRs might be closely associated with the T cell-mediated inflammatory responses.

Fig. 3.

exoLR reflecting relative fractions of different immune cell types. (a) Relative fraction of immunocytes identified using the CIBERSORT algorithm, which estimated the relative subsets of 22 types of immunocytes from 12 normal and 9 ischemic stroke specimens. (b) Comparison of the fractions of exhausted (p = 0.008) and gamma delta T cells (p = 0.016) between ischemic stroke specimens and normal controls. (c) Two-dimensional region-based PCA plot obtained by including all studied samples and all types of immune cells. Each dot indicates one sample. The red plot refers to normal control specimens, whereas the blue plot indicates ischemic stroke specimens. Abbreviations: PCA, principal component analysis.

An analysis and visualization of exosomal messenger RNA functional profiles was performed using ClusterProfiler. A bar plot of the functional profile analysis results (Supplementary Table 4) is shown in Fig. 4a. Among the top 20 biological process enrichment terms, cell adhesion-associated biological processes (cell adhesion, biological adhesion, cell-cell adhesion, and regulation of cell adhesion), metabolic processes (sphingomyelin metabolic process, phospholipid metabolic process, and ammonium ion metabolic process), and T cell regulation-associated processes (T cell co-stimulation and positive regulation of T cell proliferation) were significantly enriched. The top 20 pathway analysis enrichment terms are shown in Fig. 4b, and further details are provided in Supplementary Table 5. Beyond metabolic pathways (ko01100) and cell adhesion molecules (ko04514), immune disease pathways, including ko05320, ko05330, ko05332, ko05310, ko05323, and ko05322, accounted for the majority of the significantly enriched human disease pathways (Fig. 4b).

Fig. 4.

Biological process and pathway enrichment analyses of differentially expressed messenger RNA (DEmRNAs). (a) Bar plot ranking the top 20 biological processes enrichment terms of DEmRNAs. A lower q-value indicates significant enrichment, as indicated by the blue scale bar (high: light, low: dark). (b) Circular plot of the top 20 pathway enrichment terms for DEmRNAs. In the red scale bar, the higher the –log10 (q-value) the greater the enrichment score (high: dark, low: light). GO, Gene Ontology.

Based on the DEmRNAs and DElncRNAs, we obtained 375 trans-acting pairs, including 51 DEmRNAs and 18 DElncRNAs (Supplementary Table 6). However, no antisense or cis-acting pairs of DElncRNAs and DEmRNAs were found. These 51 DEmRNAs were considered the potential targets of DElncRNAs. The results of biological process enrichment analyses (Supplementary Table 7) are shown as a bar plot (Fig. 5a). Among the top 20 biological processes, metabolic processes (response to fructose, response to vitamins, nor-spermidine metabolic process, and cellular carbohydrate metabolic process) were significantly enriched. The results of the pathway enrichment analysis are shown in Supplementary Table 8. Beyond metabolic pathways and human diseases pathways, cellular processes pathways, including phagosome (ko04145), focal adhesion (ko04510), and ferroptosis (ko04216), might be closely associated with the regulation of the DElncRNAs in trans (Fig. 5b).

Fig. 5.

Biological process and pathway enrichment analyses of the target transcripts of differentially expressed long non-coding RNAs (DElncRNAs). (a) Bar plot ranking the top 20 biological processes enrichment terms for DElncRNAs. A lower q-value indicates significant enrichment, as indicated by the blue scale bar (high: light, low: dark). (b) Circular plot of the top 20 pathway enrichment terms for DEmRNAs. In the red scale bar, the higher the –log10 (q-value) the greater the enrichment score (high: dark, low: light).

The enrichment of the parental genes of DEcircRNAs in GO biological processes was assessed (Supplementary Table 9), and catabolism-associated processes, particularly protein catabolic processes, were among the most significantly enriched GO terms (Fig. 6a). The function of DEcircRNAs might be closely associated with protein catabolic processes in ischemic stroke. KEGG pathway enrichment analysis revealed that metabolic pathways were closely associated with the expression of DEcircRNAs (Fig. 6b and Supplementary Table 10).

Fig. 6.

Biological process and pathway enrichment analyses of the parental genes of the differentially expressed circular RNAs (DEcircRNAs). (a) Bar plot ranking the top 20 biological processes enrichment terms for DEcircRNAs. A lower q-value indicates significant enrichment, as indicated by the blue scale bar (high: light, low: dark). (b) Circular plot of the top 20 pathway enrichment terms for DEmRNAs. In the red scale bar, the higher the –log10 (q-value) the greater the enrichment score (high: dark, low: light).

A total of 158 DElncRNA predicted microRNA pairs (Supplementary Table 11) and 441 predicted microRNA-DEmRNA pairs (Supplementary Table 12) were obtained and combined to construct the lncRNA-microRNA-messenger RNA-ceRNA regulatory network (Supplementary Table 13). The lncRNA-related ceRNA networks involving the top 10 up- and downregulated messenger RNAs are shown in Fig. 7a. In the same way, the relationships between 69 DEcircRNA-predicted microRNA pairs (Supplementary Table 14) and 318 predicted microRNA-DEmRNA pairs (Supplementary Table 15) were obtained and combined to construct a circular RNA-microRNA-messenger RNA-ceRNA regulatory network (Supplementary Table 16). The circular RNA-related ceRNA networks involving the top 10 up- and downregulated messenger RNAs are shown in Fig. 8a.

Fig. 7.

LncRNA-microRNA-messenger RNA regulatory network, and biological process and pathway enrichment analyses of the differentially expressed messenger RNAs involved in the networks. (a) ceRNA regulatory network involving the top 10 up- and downregulated messenger RNAs. Red dots represent upregulated messenger RNA, blue dots represent downregulated messenger RNA, yellow dots represent upregulated lncRNA, green dots represent downregulated lncRNA, and beige dots represent microRNA. (b) Bar plot ranking the top 20 biological processes and signaling pathways, based on enrichment scores [–log₁₀ (p value)] from pathway enrichment analysis of the messenger RNAs involved in the networks. ceRNA, competing endogenous RNA; lncRNA, long-noncoding RNA; mRNA, messenger RNA.

Fig. 8.

Circular RNA-microRNA-messenger RNA regulatory network, and biological process and pathway enrichment analyses of the differentially expressed messenger RNAs involved in the networks. (a) ceRNA regulatory network involving the top 10 up- and downregulated messenger RNAs. Red dots represent upregulated messenger RNA, blue dots represent downregulated messenger RNA, yellow dots represent upregulated circular RNA, green dots represent downregulated circular RNA, and beige dots represent microRNA. (b) Bar plot ranking the top 20 biological processes and signaling pathways, based on enrichment scores [–log₁₀ (p value)] from pathway enrichment analysis of the messenger RNAs involved in the networks.

GO biological process and KEGG pathway enrichment analyses of the messenger RNAs involved in the networks were investigated and are presented in Fig. 7b and Fig. 8b, respectively. The messenger RNAs involved in both ceRNA networks were significantly enriched in the functions of cell-substrate adhesion, wound healing and spreading of cells, extracellular matrix-receptor interaction, and positive regulation of T cell proliferation. Both ceRNA networks included the lncRNA-microRNA-messenger RNA and circular RNA-microRNA-messenger RNA networks, sharing 108 common DEmRNAs. We performed GO biological process and KEGG pathway enrichment analyses on the 108 common DEmRNAs to investigate the co-regulatory mechanism. Collagen formation and cell-substrate adhesion were significantly enriched terms. The results are shown in Supplementary Fig. 3.