Identifying potential therapeutic targets for ischemic stroke through immune infiltration analysis and construction of a programmed cell death‑related ceRNA network.

The present study aimed to uncover the underlying mechanisms and potential intervention targets of ischemic stroke (IS). An immune cell infiltration analysis using CIBERSORT was performed on two stroke-related datasets from the Gene Expression Omnibus database to generate a competitive endogenous RNA (ceRNA) network. Gene Ontology (GO) functional annotation and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses were used to predict potential biological functions of the network. Differentially expressed genes in the ceRNA network and programmed cell death-related genes were intersected to obtain common genes for a stroke ceRNA network related to programmed cell death. These genes were further verified through a middle cerebral artery occlusion rat model using reverse transcription-quantitative PCR. This built a ceRNA regulatory network based on bioinformatic analysis. In addition to the biological functions extracted from the GO and KEGG enrichment analyses, it was discovered that long non-coding (lnc)RNA-mediated ceRNA regulatory pathways were associated with programmed cell death. This included five for apoptosis (lncRNA deleted in lymphocytic leukemia 2 like (DLEU2L)/micro (miR)-4500/sulfite oxidase, lncRNA DLEU2L/miR-4500/transforming growth factor β receptor III, lncRNA DLEU2L/miR-4500/BTB and CNC homology 1 (BACH1), lncRNA DLEU2L/miR-4500/zinc finger and BTB domain containing 5 and lncRNA LINC00266-1/miR-363-3p/zinc finger protein 354B and two for ferroptosis (lncRNA DLEU2L/miR-4500/ribo-nucleotide reductase regulatory subunit M2 and lncRNA DLEU2L/miR-4500/BACH1). Based on the aforementioned results, the present study provided potential approaches for bridging programmed cell death, immune infiltration and ceRNA regulatory networks in IS. The present study may provide novel insights into the clinical diagnosis and treatment of IS, and may improve the knowledge of the regulation of pathophysiological processes for IS. Copyright: © Cui et al.

Introduction

Ischemic stroke (IS) is a common disease among middle-aged and elderly individuals, typically caused by cerebral ischemic injury. IS has high morbidity and mortality worldwide, causing severe disability and negatively affecting quality of life (1). The majority of patients generally receive conservative treatment or interventional therapies, such as thrombolytics (2). New treatment strategies are urgently needed to reduce the mortality and morbidity of patients with IS, and therefore improve patient prognosis and quality of life. Thus, the present study may provide novel directions in which to explore the physiological and pathological mechanisms for IS.

Long non-coding RNAs (lncRNAs), microRNAs (miRNAs/miRs) and circular RNAs (circRNAs) are types of non-coding RNAs (ncRNAs) that play important roles at the biomolecular level, such as transcriptional regulation, post-transcriptional regulation and epigenetics (3). Previous studies have revealed that stroke is influenced by genetic factors and that miRNAs and lncRNAs are specifically linked to IS (4,5). For example, some miRNAs, such as miR-155, miR-125a/b-5p and miR-22, can influence patient recovery during IS therapy by regulating blood pressure (6). High expression of the hyperglycemia-related FAS gene and miRNA has-let-7b-5p are associated with a poor outcome in IS (7). Some abnormal lncRNAs have been shown to have a role in IS, such as metastasis associated lung adenocarcinoma transcript 1(8), maternally expressed 3(9) and H19(10). However, evidence also indicates that lncRNAs may not directly cause these effects, but instead establish vigorous modulatory crosstalk networks through interactions with other ncRNAs by competitively combining with certain ncRNAs. This is known as the ceRNA network theory (11). Studies have stressed the key importance of the ceRNA network; however, there are also reports that a lncRNA-miRNA-mRNA net is an important mechanism for regulating post-transcriptional gene translation (12). Therefore, the present study aimed to build a complete ceRNA network that could improve understanding of cerebral ischemic injury disease treatments.

A previous study demonstrated a relationship between stroke and programmed cell death, which indicated that a strong/dysregulated inflammatory response may aggravate programmed cell death after a stroke t (13). A recent study explored a lncRNA-miRNA-mRNA regulatory network in IS (14). After IS, the disrupted blood-brain barrier infiltrates and initiates the innate phase of the immune response which rapidly activates the immune cells. The release of brain antigens results in the rapid infiltration of dendritic cell precursors into the brain after the innate stage; these cells persist in the brain and are said to determine the outcome of programmed cell death by influencing long-term T cell responses (15). Since the immune system is closely related to programmed cell death, the present study considered whether there is a connection between the ceRNA network and programmed cell death in IS. However, to the best of our knowledge, there have been only a few bioinformatics-related studies on the mechanism of programmed cell death and ceRNA networks (16,17), and therefore it is important to explore the relationship between programmed cell death and the ceRNA network in IS.

The present study focused on constructing a global lncRNA-miRNA-mRNA network based on NCBI GEO datasets for lncRNA, miRNA and mRNA expression profiles. Driving endonuclease genes (DEGs) between patients with IS and healthy controls were extracted and analyzed. The lncRNA-miRNA and miRNA-mRNA pairs were identified, then the lncRNA-miRNA-mRNA network was constructed. The differentially expressed mRNAs in the ceRNA network were analyzed using Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) functional enrichment analyses. Next, the genes related to programmed cell death from specific source datasets were intersected with the DEGs in the ceRNA network. These regulatory pathways were then verified using reverse transcription-quantitative PCR (RT-qPCR). Immune infiltration analysis was also performed to characterize the immune environment resulting from IS (Fig. 1).

Figure 1

Flowchart of construction and analysis of ceRNA networks. GEO, Gene Expression Omnibus; FC, fold change; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; RT-qPCR, quantitative reverse transcription-polymerase chain reaction.

Materials and methods

Dataset collection

Two expression profile datasets were downloaded from the NCBI GEO database (https://www.ncbi.nlm.nih.gov/geo/). The lncRNA and mRNA dataset GSE58294 was based on the GPL570 platform [(HG-U133_Plus_2) Affymetrix Human Genome U133 Plus 2.0 Array] (18), and included 69 cardioembolic stroke samples and 23 healthy control samples. The miRNA array dataset GSE55937 was obtained using the GPL16384 platform [(miRNA-3) Affymetrix Multispecies miRNA-3 Array] (19), and contained data for 24 patients with IS and 24 healthy controls.

Differential expression analysis

The limma package (version 3.52.2) (20) in R software (version 4.1.1, http://www.R-project.org/) (21) was used for differential expression analysis of lncRNAs, miRNAs and mRNAs between patients with IS and healthy controls. The selection standard was set to P<0.05 and |log2 fold change (FC)|>1 as threshold values. The ggplot2 (https://github.com/tidyverse/ggplot2, version 3.3.6) and pheatmap packages (https://www.rdocumentation.org/packages/pheatmap/versions/1.0.12) in R software were used for drawing heat maps and volcano maps, respectively (22).

Immune infiltration analysis

CIBERSORT (https://cibersort.stanford.edu/, version 1.03) was used for immune cell infiltration analysis. The difference in immune cell infiltration between IS and control tissues was analyzed using the corrplot package (https://github.com/taiyun/corrplot, version 0.92) in R software. The barplot and vioplot were produced by pheatmap and vioplot (https://github.com/TomKellyGenetics/vioplot, version 0.3.7) packages in R software, respectively, for visualization.

Establishment of the lncRNA-miRNA-mRNA network

miRNA-lncRNA pairs were identified using the StarBase database, and miRNA-mRNAs interaction pairs were also identified using miRDB, TargetScan (https://www.targetscan.org/vert_80/, version 8.0) and miRTarBase databases (23-28). Pairs were selected according to whether both the lncRNA and mRNA were targeted or co-expressed with a common miRNA in the lncRNA-miRNA-mRNA network. These lncRNA-miRNA-mRNA groups were recognized as co-expression competing triplets, and the corresponding ceRNA regulatory network was constructed. Cytoscape 3.9.0 was employed to visualize the ceRNA regulatory networks (29).

Functional enrichment analysis

GO functional (http://geneontology.org/, version 1.17) annotation and KEGG (https://www.kegg.jp/, version 100.0) analyses were performed, and the clusterprofiler package in R software was applied to predict potential biological functions of the mRNAs underlying the ceRNA network (30). Gene sets with P<0.05 were considered greatly enriched, and bar plots and bubble charts, which were generated using R software, were used to visualize the results.

Programmed cell death-related gene identification in the ceRNA network

The differentially expressed mRNAs in the ceRNA network were compared with the programmed cell death-related genes from GeneCards (https://www.genecards.org/, version 5.6), and with the genes related to apoptosis and ferroptosis. A Venn diagram was generated using R software to visually show the genes common to each set.

Middle cerebral artery occlusion (MCAO) model

Adult male Sprague-Dawley rats (2 months old) were purchased from Qingdao Peng Yue Animal Husbandry Co., Ltd and were raised at the Experimental Animal Center of Qingdao University. Experiments were implemented on the rats (n=12; weight, 220±10 g), which were maintained in a controlled environment (22±2˚C room temperature, 50-60% humidity and dark/light cycle 12 h/12 h) with water and chow ad libitum. Animals were randomized into the following two groups: i) IS rats (n=6), which received MCAO surgery; and ii) sham rats (n=6), which underwent surgery without IS. A MCAO model was established according to Zea Longa's method to evaluate neural function after 24 h (31). Rats were fixed in the supine position and anesthetized using intraperitoneal injection of 10% chloral hydrate solution (350 mg/kg) (32). None of the experiment animals exhibited signs of peritonitis, pain or discomfort. The internal and external carotid arteries of the common carotid were carefully isolated, and the proximal common carotid and distal external carotid arteries were ligated. A nylon bolt was slowly inserted into the internal carotid artery and secured with a fixation wire. After 90 min of blood flow occlusion, the bolts were opened, and reperfusion was performed for 3, 5 and 24 h. Eventually, the wounds were sutured layer by layer, during which the ambient temperature was maintained at 37˚C±0.5˚C, and the mice were monitored for rectal temperature, respiratory rate and heart rate. This animal experiment was reviewed and approved by the Ethics Committee of Qingdao Affiliated Hospital (approval no. 20210125SD3620211228015).

RT-qPCR Assay

Following ischemia-reperfusion, the infarct core area and corresponding cortex were obtained from rats. For RT-qPCR testing of programmed cell death-related ceRNA network mRNA expression, total RNA from each sample was extracted using PureLink™ RNA Mini kit (Invitrogen, Thermo Fisher Scientific, Inc.). RNA purity was tested using a NanoDrop ONEc (Thermo Fisher Scientific, Inc.). The total RNA used Evo M-MLV RT Mix kit with gDNA Clean for qPCR (Hunan Aikerui Biological Engineering Co., Ltd.) for reverse transcription reactions; cDNA synthesis was conducted for 15 min at 37˚C and the reaction was terminated by heating at 85˚C for 5 sec. Subsequently, the pre-amplified cDNA samples were mixed with SYBR Green Premix Pro Taq HS qPCR Kit (Hunan Aikerui Biological Engineering Co., Ltd.), for which the following thermocycling conditions were required: One cycle at 95˚C for 30 sec, 40 cycles at 95˚C for 5 sec and 60˚C for 30 sec. The reaction was conducted on a Quant Studio™ 3 Real-Time PCR Instrument (Thermo Fisher Scientific, Inc.). Experiments were repeated in triplicate. GAPDH was used as an internal reference and the relative expression level was calculated using the 2-ΔΔCq method (Table I) (33).

Table I

Primer sequences of RNAs for reverse transcription-quantitative PCR.

Primer	Species	Primer sequence (5'-3')
GAPDH	Rat	F' GCATCTTCTTGTGCAGTGCC
		R' ACCAGCTTCCCATTCTCAGC
BACH1	Rat	F' GACTGTGAGGTGAAACTGCCTTT
		R' GCAGCGATCCTGTTCTTGCT
SUOX	Rat	F' GACTCACTTAAGCCGCAGCAT
		R' TCCGTGGTGACTCCTGAGAAG
ZNF354B	Rat	F' CATGGATGCTCCTCTCTTGGA
		R' GTCAGGTCTGCAGGATTTTCCT
TGFBR3	Rat	F' AGCGCTTCAGCTTCCTGTTC
		R' TGGTGGCATCGAGAGAAGTG
ZBTB5	Rat	F' ACTTTCCCAGGTCTGGTTTGG
		R' CAGAGATCTGTGAGCTCCTGACA
RRM2	Rat	F' GCATTGAGGAAGAGCCGTTACT
		R' AGTGCTGAATATCCTTGGAAAGGT

RT-qPCR, reverse transcription-quantitative PCR; BACH1, BTB and CNC homology 1; SUOX, sulfite oxidase; ZNF354B, zinc finger protein 354B; TGFBR3, transforming growth factor beta receptor III; ZBTB5, zinc finger and BTB domain containing 5; RRM2, ribonucleotide reductase regulatory subunit M2.

Statistical analysis

GraphPad Prism 9.0 (GraphPad Software, Inc.) software was used for graphing and performing statistical analyses. All data were repeated three times and are presented as the mean ± standard deviation. Both datasets were analyzed for differentially expressed genes using one-way ANOVA with Tukey's post hoc test. P<0.05 was considered to indicate a statistically significant difference.

Results

Volcano plots and heat maps of lncRNAs, miRNAs and mRNAs with differential expression were obtained based on the comparative analysis of lncRNA, miRNA and mRNA expression profiles between patients with IS and healthy controls (Fig. 2). Plot criteria of P<0.05 and |log2 FC|>1 were used as thresholds for significant differential expression. Altogether, 15 lncRNAs (11 upregulated and four downregulated), 43 miRNAs (27 upregulated and 16 downregulated) and 463 mRNAs (347 upregulated and 116 downregulated) were identified as genes differentially expressed between patients with IS and healthy controls (Table II, Table III and Table IV; top 30 mRNAs are listed in Table IV).

Figure 2

Differentially expressed analysis of lncRNAs, miRNAs, and mRNAs. Heatmap of differentially expressed (A) lncRNAs, (B) miRNAs and (C) mRNAs. Volcano plot of (D) 23 differentially expressed lncRNAs, (E) 46 differentially expressed miRNAs and (F) 300 differentially expressed mRNAs. Red represents upregulated genes and green represents downregulated genes. DE, differentially expressed; FC, fold change; lncRNA, long non-coding RNA; miRNA, microRNA.

Table II

Identification of differentially expressed lncRNA in ischemic stroke.

A, Upregulated differentially expressed lncRNA
lncRNA	logFC	P-value
ASAP1-IT2	1.114	<0.001
THUMPD3-AS1	1.097	<0.001
LINC00266-1	1.668	<0.001
DLEU2L	1.428	<0.001
IQCH-AS1	1.066	<0.001
ERVK13-1	2.985	0.002
ATP1A1-AS1	3.766	0.003
LINC00189	3.016	0.005
HYMAI	1.721	0.004
C20orf197	2.158	0.006
SNRK-AS1	1.465	0.011
B, Downregulated differentially expressed lncRNA
IncRNAs	logFC	P-value
LINC00342	-1.706	<0.001
HOXB-AS1	-1.622	<0.001
LINC00565	-1.598	0.004
LINC00920	-4.013	0.003

IS, ischemic stroke; FC, fold change; lncRNA, long non-coding RNA; ASAP1-IT2, ASAP1 intronic transcript 2; THUMPD3-AS1, THUMPD3 antisense RNA 1; LINC00266-1, long intergenic non-protein coding RNA 266-1; DLEU2L, deleted in lymphocytic leukemia 2 like; IQCH-AS1, IQCH antisense RNA 1; ERVK13-1, endogenous retrovirus group K13 member 1; ATP1A1-AS1, ATP1A1 antisense RNA 1; LINC00189, long intergenic non-protein coding RNA 189; HYMAI, hydatidiform mole associated and imprinted (non-protein coding); C20orf197, chromosome 20 open reading frame 197; SNRK-AS1, SNRK antisense RNA 1; LINC00342, long intergenic non-protein coding RNA 342; HOXB-AS1, HOXB cluster antisense RNA 1; LINC00565, long intergenic non-protein coding RNA 565; LINC00920, long intergenic non-protein coding RNA 920.

Table III

Identification of differentially expressed miRNA in ischemic stroke.

A, Upregulated differentially expressed miRNA
miRNA	logFC	P-value
hsa-miR-140-3p	8.843	0.002
hsa-miR-145	10.319	0.002
hsa-miR-3130-5p	4.309	0.004
hsa-miR-99b	9.857	0.004
hsa-miR-330-3p	4.053	0.005
hsa-miR-425	4.690	0.007
hsa-miR-409-3p	9.304	0.009
hsa-miR-550a	5.588	0.009
hsa-miR-363	1.761	0.010
hsa-miR-339-3p	8.747	0.013
hsa-miR-93	9.339	0.018
hsa-miR-625	5.741	0.018
hsa-miR-641	1.500	0.019
hsa-miR-574-3p	6.021	0.019
hsa-miR-103a	5.800	0.020
hsa-miR-671-3p	1.998	0.021
hsa-miR-342-3p	9.402	0.023
hsa-miR-1285	5.294	0.025
hsa-miR-433	1.264	0.027
hsa-miR-491-5p	4.041	0.030
hsa-miR-20b	2.064	0.030
hsa-miR-3200-5p	5.262	0.031
hsa-miR487b	2.319	0.032
hsa-miR345	5.854	0.037
hsa-miR-3152-3p	1.011	0.037
hsa-miR-200c	3.685	0.037
hsa-miR-423-3p	7.780	0.041
B, Downregulated differentially expressed miRNA
miRNA	logFC	P-value
hsa-miR-3135b	-8.480	0.007
hsa-miR-320d	-6.770	0.011
hsa-miR-148a	-3.303	0.011
hsa-miR-320e	-5.395	0.014
hsa-miR-1231	-6.055	0.017
hsa-miR-122	-5.834	0.018
hsa-miR-4429	-5.949	0.021
hsa-let-7i	-11.751	0.036
hsa-miR-3185	-5.591	0.037
hsa-miR-4500	-1.691	0.039
hsa-miR-15a	-8.475	0.041
hsa-miR-4486	-3.999	0.041
hsa-miR-4498	-2.322	0.042
hsa-miR-19a	-3.672	0.044
hsa-let-7c	-8.036	0.048
hsa-miR-4750	-4.902	0.049

IS, ischemic stroke; FC, fold change; miRNA/miR, microRNA.

Table IV

Identification of differentially expressed mRNA (Top 30) in ischemic stroke.

A, Upregulated differentially expressed mRNA
mRNA	logFC	P-value
MAF	8.870	<0.001
BACH1	8.691	<0.001
HTRA1	7.681	<0.001
ZNF354B	7.114	0.001
PTX3	6.322	<0.001
ZSWIM3	6.281	<0.001
NCAPG	6.193	0.004
ZNF174	5.434	0.010
SPR	4.914	0.016
PLB1	4.653	0.008
REX02	4.626	0.002
TGFA	4.298	0.003
GASP2	4.270	0.002
PPFIA1	4.193	<0.001
SLC39A6	4.100	<0.001
FUT7	3.832	<0.001
AMDHD1	3.822	<0.001
TFRC	3.811	0.002
ITPKC	3.775	<0.001
C8orf58	3.578	<0.001
TSPAN7	3.538	<0.001
C9orf84	3.488	0.002
OSBPL1A	3.399	0.008
RBM41	3.328	0.009
COPS8	3.239	0.007
CCDC88A	3.217	<0.001
PRR3	3.197	<0.001
SCAMP1	3.161	0.004
KIF11	3.123	<0.001
TRNT1	3.104	0.003
B, Downregulated differentially expressed mRNA
mRNA	logFC	P-value
ABL1	-1.004	0.016
MYH3	-1.005	0.001
MFSD12	-1.007	<0.001
ZC3H8	-1.008	<0.001
TUT1	-1.016	<0.001
BACH2	-1.018	0.001
HBQ1	-1.026	<0.001
SNX29	-1.029	0.006
NIPSNAP1	-1.033	<0.001
IL11RA	-1.049	<0.001
NOG	-1.051	<0.001
ELP5	-1.053	<0.001
SURF2	-1.057	<0.001
HPS4	-1.058	<0.001
KANK3	-1.064	0.005
PPIL1	-1.065	<0.001
ZBTB5	-1.071	<0.001
MAP2K2	-1.077	<0.001
KAT2A	-1.081	0.014
CCDC85C	-1.081	<0.001
KCNN4	-1.089	<0.001
GGA2	-1.089	0.003
BANK1	-1.102	<0.001
ABCA3	-1.104	0.021
SPSB2	-1.104	0.003
CBFA2T2	-1.110	<0.001
PASK	-1.113	0.020
WDR46	-1.114	<0.001
BCL2L12	-1.115	0.007
GTF3C1	-1.116	<0.001

IS, ischemic stroke; FC, fold change.

Analysis of infiltrated immune cells in ischemic stroke tissues from the GEO expression array dataset

The proportion of infiltrating immune cells in IS tissue was investigated using GEO expression array data. A total of 92 samples were analyzed using CIBERSORT (P<0.05). As presented in Fig. 3, there were 19 subpopulations of immune cells in 69 stroke samples and 23 healthy control samples. Different colors represent different types of immune cells. Compared with normal tissue, the violin plot of the immune cells demonstrated that ‘T cells CD4 memory-activated’, ‘T cells follicular helper’, ‘macrophages M0’ and ‘neutrophils’ infiltrated statistically more, while ‘B cells naive’, ‘T cells CD8’, ‘T cells CD4 naive’ and ‘mast cells resting’ infiltrated statistically less in IS tissue compared with control tissue (Fig. 4).

Figure 3

Proportions of 19 immune cell subsets in IS tissue. X-axis: each sample; Y-axis: the percentage of each immune cell. GSM1406033-GSM1406055 correspond to healthy individuals.

Figure 4

Violin plot showing differences in immune infiltration between normal control tissue and IS tissue. The normal control group is shown in blue, and the IS group is shown in red. IS, ischemic stroke.

In addition, a correlation analysis of immune cells infiltrating IS revealed multiple pairs of immune cells that were either positively or negatively correlated (Fig. 5). The score indicated the degree of relevance: The bigger size of the numbers statistics data represented the more positive or negative the correlation. This result suggested that CD8 T cells and naive B cells had a positive correlation (value=0.42). Conversely, CD8 T cells and neutrophils indicated a significant correlation (value=-0.61).

Figure 5

Co-expression patterns among immune cell branches. Red: positive correlation, blue: negative correlation.

LncRNA-miRNA-mRNA ceRNA network construction

Based on interacting elements, five miRNA-lncRNA pairs and 31 miRNA-mRNA pairs were identified in the ceRNA network. After predicting lncRNA-miRNA and miRNA-mRNA pairs, the interactions were visualized using Cytoscape. This established the lncRNA-miRNA-mRNA ceRNA network, with lncRNA DLEU2L linked to one differentially expressed microRNA (DEmiR) and eight DEGs, and lncRNA LINC00266-1 linked to one DEmiR and 23 DEGs (Fig. 6).

Figure 6

Construction of the ceRNA network. Green represents mRNA, pink represents lncRNA and blue represents miRNA. miR/miRNA, microRNA; ceRNA, competing endogenous RNA.

GO and KEGG functional enrichment analysis

The clusterprofiler package in the R software was adopted for functional enrichment analysis to improve the exploration of the potential biological functions genes expressed in the ceRNA network. The results indicated that differentially expressed mRNAs were partly enriched in the following categories: ‘In utero embryonic development’ (biological process), ‘lysosomal membrane’ (cellular component), ‘lytic vacuole membrane’ (cellular component), ‘protein serine/threonine kinase activity’ (molecular function) and mRNA interactions of the ceRNA network (Fig. 7A and B). Furthermore, KEGG pathways revealed that differentially expressed mRNAs were associated with ‘salmonella infection’, ‘human cytomegalovirus infection’ and the ‘Ras-proximate-1 signaling pathway’ in the ceRNA network (Fig. 7C and D).

Figure 7

Functional enrichment analysis of differentially expressed mRNAs in the ceRNA network. (A) GO biological functional analyses using a barplot. (B) GO biological function analyses using a bubble chart. (C) KEGG pathway analyses using a barplot. (D) KEGG pathway analyses using a bubble chart. BP, biological process; CC, cellular component; MF, molecular function; ceRNA, competing endogenous RNA; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes.

Programmed cell death-related genes in the ceRNA network

In the ceRNA network, 31 differentially expressed mRNAs were intersected with 3,101 programmed cell death-related genes from GeneCards, resulting in six common genes (Fig. 8A; Table V. Next, 31 differentially expressed mRNAs in the ceRNA network interacted with 2,720 apoptosis-related genes in GeneCards and five common genes were obtained (Fig. 8B; Table V). A similar analysis for ferroptosis identified two common genes in the ceRNA network (Fig. 8C; Table V). A common gene was obtained through interacting with 31 differentially expressed mRNAs in the ceRNA network, 2,720 apoptosis-related genes and 253 ferroptosis-related genes (Fig. 8D; Table V).

Figure 8

Identification of the candidate programmed cell death-related genes in the ceRNA network. (A) Venn diagram of mRNAs common between differentially expressed mRNAs in our ceRNA network and programmed cell death-related mRNAs. (B) Venn diagram for identifying mRNAs common between differentially expressed mRNAs in our ceRNA network and apoptosis-related mRNAs. (C) Venn diagram for identifying mRNAs common between differentially expressed mRNAs in our ceRNA network and ferroptosis-related mRNAs. (D) Interactions between differentially expressed mRNAs in our ceRNA network and apoptosis and ferroptosis-related mRNAs. ceRNA, competing endogenous RNA.

Table V

Identified lncRNA-mediated ceRNAs and the intersection between programmed cell death related genes and mRNA in the ceRNA networks.

lncRNA	miRNA	mRNA	Programmed cell death	Apoptosis	Ferroptosis
DLEU2L	hsa-miR-4500	SUOX	X	X
DLEU2L	hsa-miR-4500	SEMA4C
DLEU2L	hsa-miR-4500	RRM2	X		X
DLEU2L	hsa-miR-4500	YOD1
DLEU2L	hsa-miR-4500	TGFBR3	X	X
DLEU2L	hsa-miR-4500	BACH1	X	X	X
DLEU2L	hsa-miR-4500	STRN
DLEU2L	hsa-miR-4500	ZBTB5	X	X
LINC00266-1	hsa-miR-363-3p	GLYR1
LINC00266-1	hsa-miR-363-3p	DUSP10
LINC00266-1	hsa-miR-363-3p	ACOX1
LINC00266-1	hsa-miR-363-3p	ZNF354B	X	X
LINC00266-1	hsa-miR-363-3p	ZNF230
LINC00266-1	hsa-miR-363-3p	PRRG4
LINC00266-1	hsa-miR-363-3p	ANKIB1
LINC00266-1	hsa-miR-363-3p	CXXC5
LINC00266-1	hsa-miR-363-3p	FAM160B1
LINC00266-1	hsa-miR-363-3p	CREB3L2
LINC00266-1	hsa-miR-363-3p	HIPK1
LINC00266-1	hsa-miR-363-3p	CAND1
LINC00266-1	hsa-miR-363-3p	MIER1
LINC00266-1	hsa-miR-363-3p	RGL1
LINC00266-1	hsa-miR-363-3p	SERTAD3
LINC00266-1	hsa-miR-363-3p	CPEB2
LINC00266-1	hsa-miR-363-3p	RBM28
LINC00266-1	hsa-miR-363-3p	MPP7
LINC00266-1	hsa-miR-363-3p	GPD2
LINC00266-1	hsa-miR-363-3p	CPEB4
LINC00266-1	hsa-miR-363-3p	ZFYVE21
LINC00266-1	hsa-miR-363-3p	TBL1XR1
LINC00266-1	hsa-miR-363-3p	HIPK3
C20orf197	hsa-miR-363-3p
SNRK-AS1	hsa-miR-363-3p
ERVK13-1	hsa-miR-363-3p

lncRNA, long non-coding RNA; ceRNA, competing endogenous RNA.

Potential biomarker expression genes in the programmed cell death regulatory pathways by RT-qPCR

A total of six genes were selected for RT-qPCR testing to validate the programmed cell death-related genes. This included the upregulated and downregulated mRNAs in the present study's ceRNA network as follows: Five common apoptosis-related genes (ZBTB5, BACH1, SUOX, ZNF354B and TGFBR3); and two ferroptosis-related mRNAs (BACH1 and RRM2). Comparing IS from rats' MCAO models with healthy normal controls, RT-qPCR results showed that the ferroptosis-related genes BACH1 and RRM2 were confirmed to be differentially upregulated in IS, significantly after 24 h of ischemia-reperfusion (Fig. 9). These results also showed that the apoptosis-related genes BACH1, SUOX, ZNF354B, and TGFBR3 were also markedly upregulated compared with the sham group, while ZBTB5 was significantly downregulated in IS.

Figure 9

Gene validation in programmed cell death regulatory pathways. Expression of the following mRNAs in the ceRNA network; (A) BACH1, (B) TGFBR3, (C) RRM2, (D) SUOX, (E) ZNF354B and (F) ZBTB5. *P<0.05, ***P<0.001 and ****P<0.0001. ceRNA, competing endogenous RNA; I/R, ischemia/reperfusion; BACH1, BTB and CNC homology 1; TGFBR3, transforming growth factor beta receptor III; RRM2, ribonucleotide reductase regulatory subunit M2; SUOX, sulfite oxidase; ZNF354B, zinc finger protein 354B; ZBTB5, zinc finger and BTB domain containing 5.

Discussion

In the present study, a lncRNA-miRNA-mRNA ceRNA regulatory network was constructed based on GEO datasets for IS. GO and KEGG enrichment analysis revealed that mRNAs with differential expression in the ceRNA network were particularly enriched in the ‘in utero embryonic development’ (biological process) and ‘lysosomal membrane’ (cellular component) functions. After the ceRNA network was compared with programmed cell death-related genes, six mRNAs (ZBTB5, BACH1, SUOX, ZNF354B, TGFBR3 and RRM2) were identified. Finally, RT-qPCR testing revealed that several ceRNA regulatory pathways are closely related to programmed cell death in IS pathophysiology. LncRNAs are a special type of transcript family without protein coding functions (34). LncRNAs likely play a ‘bridge’ role in regulating genes levels in the ceRNA network (35). LncRNAs act as natural sponges; they compete with each other to adsorb specific miRNAs, thereby reducing the binding of miRNAs to corresponding target genes and altering the expression of miRNA target genes (36). Research has indicated that certain lncRNAs target miRNAs to regulate signaling pathways that alleviate apoptosis and inflammation. For example, lncRNA GAS5 may regulate IS as a competing endogenous RNA to miR-137 to, in turn, regulate the Notch1 signaling pathway (37), and lncRNA SNHG1 reduces apoptosis and inflammation during an IS by targeting miR-376a and modulating the cystathionine-β-synthase (CBS)/H2S pathway (38). However, it is unclear at the time of IS occurrence whether abnormal lncRNAs act as ceRNAs of specific miRNAs to regulate the expression of target mRNAs, and whether this has specific effects on the arterial wall.

In the present study's comprehensive bioinformatic analysis, lncRNA DLEU2L and lINC00266-1 were two key nodes of the ceRNA network associated with miRNA and mRNA. LncRNA DLEU2L is usually regarded as a prognostic biomarker for cancer cells (39). A recent study showed that DLEU2L can target miR-210-3p to construct a ceRNA aimed at suppressing the proliferation, migration and infiltration of pancreatic cancer cells while promoting apoptosis (39). Additionally, the expression of lncRNA LINC00266-1 and LINC00266-1 encoded miR-548c-3p in osteosarcoma cells suppress the proliferative and metastatic abilities of the cells and promote apoptosis, significantly reducing the growth of osteosarcoma cells in vivo (40). However, the mechanism of lncRNAs in IS is still uncertain. Consequently, there is reason to consider that the present study's ceRNA regulatory network may arise during or post IS.

In addition to the lncRNAs, miRNAs may also play a role in the ceRNA regulatory network of IS. To further explore the biological functions of differentially expressed mRNAs in IS, the present study performed GO and KEGG pathway analyses on the 31 differentially expressed mRNAs in the ceRNA network. These mRNAs were significantly enriched in the ‘in utero embryonic development’ (biological process), especially in the regulation of immune cells. However, how the extracellular matrix organization regulates the regulation of immune cells remains to be further investigated.

This is consistent with the view that cerebral and peripheral immune responses induced by cerebral ischemic injury are involved in the regulation of neural repair after stroke (15). Therefore, it is reasonable to hypothesize that the present study's newly established ceRNA network might play a role in immune infiltration in patients with IS. The combined effects of inflammatory alternative pathways and ceRNA networks may be potential targets for IS pathogenesis.

Programmed cell death is also a biological function extracted from the analysis of the GO and KEGG pathways and is considered essential in the regulation of cerebral injury in IS (41). Programmed cell death, including apoptosis and ferroptosis, represent the primary means by which organisms coordinate the elimination of damaged cells that either present a risk of tumor transformation or are hijacked by microorganisms for pathogen replication (42).

Interactions between significantly differentially expressed mRNAs were identified in the ceRNA network with programmed cell death-related genes. In this network, six mRNAs associated with these two subtypes of programmed cell death were obtained. BACH1, SUOX, ZNF354B, TGFBR3 and ZBTB5 are related to apoptosis and, among these, BACH1 was the most representative gene. BACH1, known as a transcription factor, can negatively regulate various antioxidant genes, such as heme oxygenase-1, that are involved in oxidative stress (43). The circRNA HIPK3, which is enriched in the heart, is associated with oxidative stress-induced ischemia/reperfusion damage (44).

However, the impact of dynamic regulatory crosstalk networks, including the hub nodes on the paths of the IS regulatory network, has not been fully studied. The newly identified ceRNA network not only regulated apoptosis but also programmed cell death of all subtypes, potentially affecting the pathophysiological process of IS. Nonetheless, the detailed mechanism of each step needs further verification.

Furthermore, the present study revealed that BACH1 and RRM2 were related to ferroptosis, which has gained recent attention (44). It is an iron-dependent regulation of cell death characterized by the accumulation of lipid peroxides to lethal levels, resulting in oxidative damage to the cell membrane. A recent study tentatively confirmed the link between ferroptosis and cardiovascular disease (45). Previous studies have shown that ferroptosis is associated with IS, but few studies have thoroughly explored the regulatory mechanisms of ferroptosis in IS (46-48).

In the present study, DLEU2L/has-miR-4500/RRM2 and DLEU2L/has-miR-4500/BACH1 may have been involved in ferroptosis. Unfortunately, the ferroptosis-related hub nodes detected in the present study have previously been explored only in non-IS diseases. For example, RRM2 exerts an anti-ferroptotic effect on liver cancer cells by maintaining GSH synthesis (49). This may be used as a biomarker in serum to evaluate the degree of ferroptosis suppression and improve the diagnostic efficiency of liver cancer (49). BACH1 can repress genes that combat labile iron-induced oxidative stress, and ferroptosis is stimulated at the transcriptional level by BACH1 upon disruption of the balance between the transcriptional induction of protective genes and accumulation of iron-mediated damage (50).

Evidence of crosstalk between ceRNA network-mediated BACH1 regulation and ferroptosis has been reported only in cancer, and not in IS (51,52). Nonetheless, these promising results provide potential research points and information on IS for linking ferroptosis and ceRNA regulatory networks, supporting further exploration of functional variations and specific interactions between ferroptosis and ceRNA regulatory networks.

The present study was based on genetic information downloaded from the GEO database, and it is recommended to additionally perform sequence analysis (such as DNA-seq and RNA-seq) on samples obtained from humans or animals. It is also recommended to use expanded human cohorts or well-designed animal studies to validate these observations and explore integrative mechanisms. Furthermore, the potential binding affinities between the identified lncRNAs, miRNAs and mRNAs require further experimental studies, and the protein expression levels of BACH1, ZBTB5, SUOX, ZNF354B and TGFBR3 could be verified at future studies. Finally, if any of the newly proposed ceRNA regulatory networks for programmed cell death are further validated, the clinical diagnosis and treatment of IS could be improved.

In the present study, a ceRNA regulatory network was built based on bioinformatic analysis. In addition to the biological functions extracted from GO and KEGG enrichment analyses, six newly discovered lncRNA-mediated ceRNA regulatory pathways were identified and validated. These pathways influenced programmed cell death, including apoptosis (lncRNA DLEU2L/miR-4500/SOUX, lncRNA DLEU2L/miR-4500/TGFBR3, lncRNA DLEU2L/miR-4500/BACH1, lncRNA DLEU2L/miR-4500/ZBTB5 and lncRNA LINC00266-1/miR-363-3p/ZNF354B) and ferroptosis (lncRNA DLEU2L/miR-4500/RRM2 and lncRNA DLEU2L/miR-4500/BACH1). The present study's results provided novel insights and potential research options for regarding link between programmed cell death, immune infiltration, and ceRNA regulatory networks in IS.