Blood samples of three pairs of patients with primary MIA and matched healthy individuals were collected at the Department of Neurosurgery, Longyan First Hospital Affiliated to Fujian Medical University (Longyan, China) between April 2018 and June 2018. Patients in the MIA group were selected from patients with onset within 24 h who were confirmed to have one ruptured IA by head computed tomography angiography (CTA) or digital subtraction angiography (DSA). The patients received microsurgical clipping or interventional therapy within the acute phase (24 h) after one of the aneurysms was ruptured. Those patients who had cerebrovascular diseases other than IA, systemic malignant tumor or severe complications were excluded. Individuals in the matched group were selected among individuals who were attending a regular health check-up and underwent head CTA/DSA to exclude IA. To ensure the quality of samples, the present study excluded patients with factors that may have affected the state of peripheral blood, including pregnancy, chemotherapy and fever (≥37.3˚C). Individuals in the MIA and matched groups were matched based on age, sex and past medical history, such as hypertension and smoking. Baseline information and clinical characteristics of the two groups are presented in Table I. Another 20 individuals were recruited as the validation cohort following the same criteria as those for the initial study group used for RNA sequencing (RNA-seq). They were divided into 10 pairs of matched patients with primary MIA as a test group and healthy individuals as a control group. The baseline information and clinical features of the two groups are presented in Table II. All volunteers included provided written informed consent and the present study was approved by the Ethics Committee of Longyan First Hospital Affiliated to Fujian Medical University (Longyan, China).

Total RNA was extracted from peripheral blood mononuclear cell samples using TRlzol reagent (Takara Bio, Inc.) according to the kit's instructions, as described previously (10). The NanoDrop-1000 (Thermo Fisher Scientific, Inc.) was used to quantify total RNA, and optical density at 260 nm (OD260)/OD280 ratios between 1.8 and 2.1 were considered acceptable. Agarose gel electrophoresis and the Nanodrop spectrophotometer were used to check the quantity of RNA prior to sequencing and reverse transcription-quantitative PCR (RT-qPCR).

RNase R (Epicentre; Illumina, Inc.) was used to degrade linear RNA and enrich circRNAs in total RNA. Subsequently, the circRNA library was constructed according to the manufacturer's protocols (NEBNext Ultra Directional RNA Library Prep kit; New England BioLabs, Inc.). The Library Quantification kit (Kapa Biosystems; Roche Diagnostics) was used to identify and quantify sequences using the Illumina HiSeq 4000 system (Illumina, Inc.). By using paired-end sequencing, the nucleotide length of sequencing was 150 bp. Single-stranded DNA was generated using 0.1 M NaOH and its loading concentration was 8 pM, which was measured by qPCR (11,12). Sequencing quality was evaluated by FastQC v0.11.7 software (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/). Adaptors and poor-quality bases were trimmed with Cutadapt v1.14 (https://cutadapt.readthedocs.io/en/stable/). In order to ensure the quality and reliability of data analysis, it was necessary to filter the original data. The filtering comprised the following: Reads with adapter were removed; reads with a ratio of indeterminable base information of >0.002 were removed; when the number of low-quality bases in a single-ended read exceeded 50% of read length, this paired read was removed. Furthermore, genes with low expression levels were filtered in the analysis. The threshold values were as follows: 1.5-fold change (FC), P≤0.05 and mean value of fragments per kilobase of exon model per million reads mapped ≥0.5.

Base calling and image processing were performed using Solexa pipeline v1.8 (Off-Line Base Caller software, v1.8). Quantile normalization and subsequent data processing were performed using the Ballgown package v2.8.4 (http://www.bioconductor.org/packages/release/bioc/html/ballgown.html) in R software v3.4.1 (www.r-project.org). Differential expression patterns of circRNAs in the samples were visualized using hierarchical clustering and a scatter plot. Significantly differentially expressed circRNAs (FC >2 and P<0.05) were identified using volcano plot filtering.

To validate the changes in circRNA expression identified by RNA-sequencing, 10 circRNAs, including 5 upregulated and 5 downregulated circRNAs, were selected. They met the following requirements: i) Length of circRNA between 200 and 3,000 bp; ii) |log2 FC|>2; iii) P<0.05; and iv) exon-related circRNA. The expression levels of circRNAs were evaluated by fluorescence real-time PCR (C1000; Applied Biosystems; Thermo Fisher Scientific, Inc.) in 10 pairs of matched patients with primary MIA as a test group and healthy individuals as a control group using SYBR Green as the probe, as previously described (13). The housekeeping gene β-actin was used as an internal control. The results were calculated using the 2^-ΔΔCq method (14).

GO enrichment analysis (https://david-d.ncifcrf.gov/) was conducted to explore the potential biological functions of the target genes of the circRNAs identified. KEGG pathway analysis (http://www.genome.jp/kegg/) revealed the signaling pathways the target genes of identified circRNAs were involved in at the molecular level. Genes with a corrected P<0.05 were considered to be enriched.

miRNA binding sites of circRNAs were predicted using TargetScan v7.2(15), circBank v2014(16) and miRanda v3.3a. Information on miRNA-mRNA regulation was obtained using miRTarBase (17) and TargetScan. The potential effects of the differentially expressed circRNAs associated with ‘leukocyte transendothelial migration’ were further predicted by constructing a competing endogenous RNA (ceRNA) network of circRNA-miRNA-mRNA interactions. Cytoscape v3.6.1 was used to construct the graph of the circRNA-miRNA-mRNA network.

Normally distributed data were analyzed using a two-tailed t-test and a Mann-Whitney U test was used for skewed data. All statistical analyses were performed using SPSS v19.0 (IBM Corp.). P<0.05 (two-tailed) was considered to indicate a statistically significant difference.

The normalized expression profiles are presented in Fig. 1. Hierarchical clustering revealed the difference in circRNA expression profiles between the two groups (Fig. 1A). A total of 3,328 differentially expressed circRNAs between the MIA group and the matched group were identified by circRNA sequencing analysis. The distribution of circRNAs in three pairs of matched samples is presented in a scatter plot (Fig. 1B). Compared to the matched group, 60 circRNAs were revealed to be significantly differentially expressed in the MIA group (|log2 FC|≥2; P<0.05). Among them, 20 were upregulated and 40 were downregulated (Fig. 1C). Furthermore, 85.79% of the circRNAs were previously identified in circbase (18) and 14.21% were novel (Fig. 1D). Most circRNAs had a predicted spliced length <1,000 nt (Fig. 1E). These circRNAs were widely distributed among all chromosomes except the Y chromosome, since only samples from female patients were collected (Fig. 1F). The top 10 most upregulated and downregulated circRNAs are listed in Table III.

A total of 10 differentially expressed circRNAs, including 5 upregulated [Homo sapiens circRNA (hsa_circ)_0008911, hsa_circ_0074837, hsa_circ_0131628, hsa_circ_0083229 and hsa_circ_0135895] and 5 downregulated (hsa_circ_0018168, hsa_circ_0002100, hsa_circ_0000690, hsa_circ_0000711 and hsa_circ_0000682) circRNAs, were selected to confirm the sequencing data using RT-qPCR. As presented in Fig. 2A, the changes observed in all selected circRNAs were similar between circRNA sequencing and the RT-qPCR assays, which confirmed the accuracy of the sequencing. As indicated in Fig. 2B, hsa_circ_0135895, hsa_circ_0000690 and hsa_circ_0000682 were significantly differentially expressed between samples from 10 patients with MIA and 10 healthy individuals (P<0.05).

GO enrichment and KEGG pathway analysis were performed with genes of a total of 3,328 differentially expressed circRNAs (Fig. 3). The top 10 most enriched GO terms by the upregulated and downregulated genes of identified circRNAs in the categories biological process (BP), cellular component and molecular function (MF) are presented in Fig. 3A and B. The top 10 GO terms for the downregulated genes of identified circRNAs in the category BP were mainly involved in metabolic processes, while those in the category MF mainly involved different types of binding functions. The KEGG pathway analysis revealed that the upregulated genes of identified circRNAs mapped to 8 pathways, including ‘thiamine metabolism’, ‘pyrimidine metabolism’ and ‘DNA replication’ (Fig. 3C and E). Furthermore, the downregulated genes of identified circRNAs were mapped to 144 pathways, including ‘leukocyte transendothelial migration’ in Fig. 4A, ‘natural killer cell mediated cytotoxicity’ in Fig. 4B, ‘T cell receptor signaling pathway’ and ‘NF-κB signaling pathway’ (Fig. 3D and F). The 60 significantly differentially expressed circRNAs were also analyzed through GO enrichment and KEGG pathway analysis. The GO enrichment analysis was unsuccessful due to the low number of circRNAs. Only five pathways were obtained in the KEGG pathway analysis of the 60 host genes of the circRNAs identified, namely ‘leukocyte transendothelial migration’, ‘natural killer cell mediated cytotoxicity’, ‘VEGF signaling pathway’, ‘axon guidance’ and ‘Wnt signaling pathway’.

TargetScan, circBank and miRanda were utilized to predict miRNAs targeted by the differentially expressed circRNAs. The top 5 miRNAs predicted to be targeted by each of the three selected significantly differentially expressed circRNAs, which were all contained in the three databases, are presented in Table IV. The miRNA-mRNA regulatory relationships were further identified using miRTarBase and TargetScan. The miRNAs targeted by the circRNAs associated with ‘leukocyte transendothelial migration’ and their regulated mRNAs were further selected for ceRNA network construction. In this network, an upregulated circRNA (hsa_circ_0135895) and two downregulated circRNAs (hsa_circ_0000682 and hsa_circ_0000690) regulate protein tyrosine kinase 2 (PTK2), protein kinase Cβ (PRKCB) and integrin subunit αL (ITGAL), respectively, by sponging a number of different miRNAs (Fig. 5). hsa-miR-4778-3p may combine with both hsa_circ_0135895 and hsa_circ_0000690, while hsa-miR-543 may combine with both hsa_circ_0135895 and hsa_circ_0000682. Several genes may combine with >20 predicted miRNAs, including ITGAL, IL2 inducible T-cell kinase, phosphatidylinositol-4,5-bisphosphate 3-kinase catalytic subunit α, phospholipase C γ1, PRKCB, PTK2, protein tyrosine phosphatase non-receptor type 11, RAP1A member of RAS oncogene family, ras homolog family member H, rho associated coiled-coil containing protein kinase 1 and signal-induced proliferation-associated 1. Among the three circRNAs, hsa_circ_0135895 sponged the most mRNAs.