Integrated bioinformatics analysis of retinal ischemia/reperfusion injury in rats with potential key genes

Abstract

The tissue damage caused by transient ischemic injury is an essential component of the pathogenesis of retinal ischemia, which mainly hinges on the degree and duration of interruption of the blood supply and the subsequent damage caused by tissue reperfusion. Some research indicated that the retinal injury induced by ischemia-reperfusion (I/R) was related to reperfusion time.

In this study, we screened the differentially expressed circRNAs, lncRNAs, and mRNAs between the control and model group and at different reperfusion time (24h, 72h, and 7d) with the aid of whole transcriptome sequencing technology, and the trend changes in time-varying mRNA, lncRNA, circRNA were obtained by chronological analysis. Then, candidate circRNAs, lncRNAs, and mRNAs were obtained as the intersection of differentially expression genes and trend change genes. Importance scores of the genes selected the key genes whose expression changed with the increase of reperfusion time. Also, the characteristic differentially expressed genes specific to the reperfusion time were analyzed, key genes specific to reperfusion time were selected to show the change in biological process with the increase of reperfusion time.

As a result, 316 candidate mRNAs, 137 candidate lncRNAs, and 31 candidate circRNAs were obtained by the intersection of differentially expressed mRNAs, lncRNAs, and circRNAs with trend mRNAs, trend lncRNAs and trend circRNAs, 5 key genes (Cd74, RT1-Da, RT1-CE5, RT1-Bb, RT1-DOa) were selected by importance scores of the genes. The result of GSEA showed that key genes were found to play vital roles in antigen processing and presentation, regulation of the actin cytoskeleton, and the ribosome. A network included 4 key genes (Cd74, RT1-Da, RT1-Bb, RT1-DOa), 34 miRNAs and 48 lncRNAs, and 81 regulatory relationship axes, and a network included 4 key genes (Cd74, RT1-Da, RT1-Bb, RT1-DOa), 9 miRNAs and 3 circRNAs (circRNA_10572, circRNA_03219, circRNA_11359) and 12 regulatory relationship axes were constructed, the subcellular location, transcription factors, signaling network, targeted drugs and relationship to eye diseases of key genes were predicted. 1370 characteristic differentially expressed mRNAs (spec_24h mRNA), 558 characteristic differentially expressed mRNAs (spec_72h mRNA), and 92 characteristic differentially expressed mRNAs (spec_7d mRNA) were found, and their key genes and regulation networks were analyzed.

In summary, we screened the differentially expressed circRNAs, lncRNAs, and mRNAs between the control and model groups and at different reperfusion time (24h, 72h, and 7d). 5 key genes, Cd74, RT1-Da, RT1-CE5, RT1-Bb, RT1-DOa, were selected. Key genes specific to reperfusion time were selected to show the change in biological process with the increased reperfusion time. These results provided theoretical support and a reference basis for the clinical treatment.

Introduction

As one of the highest oxygen-consuming tissues in the body, the retina has a delicate and complex structure and a vigorous metabolism [1]. The retina is a transparent membrane that is divided into ten layers from the outside, namely, the retinal pigment epithelium, the cone and rod layer of the retinal neurosensory layer, the outer limb, the outer nuclear layer, the outer plexiform layer, the inner nuclear layer, the inner plexiform layer, the ganglion cell layer composed of ganglion cell nuclei, the nerve fiber layer composed of nerve cell axons, and the inner limb [2, 3]. Retinal tissue has a dual blood supply system for adequate energy supply and oxygen consumption: the outer five layers of the retina are supplied by the choroidal blood vessels (referred to as the choroidal circulation) [4,5,6], whereas the inner five layers are supplied by the central retinal artery (referred to as the retinal circulation) [7, 8]. The choroidal circulation is much denser than the retinal circulation, and the sparse retinal circulation is more favorable for light passage but makes the retina more susceptible to vascular disease [9]. Any pathological damage or retinal ischemia and hypoxia caused by retinal vascular obstruction can lead to infarction of the retinal tissue cells, which results in the loss of the function of receiving and transmitting light stimuli from the external environment. After the blood supply is restored, the damage to the structure and function of the retina is further aggravated, resulting in retinal ischemia /reperfusion (RIR) injury [10, 11]. A series of scientific studies have confirmed that the pathogenesis of RIR injury also includes a series of deleterious events, including complement system activation and leukocyte recruitment, endoplasmic reticulum stress, calcium overload, decreased oxidative phosphorylation, increased free radicals, vascular endothelial cell dysfunction, apoptosis signaling, necrosis and autophagy [12,13,14]. It is the main cause and pathogenesis of retinal thinning and atrophy, retinal ganglion cell (RGC) death, and visual function impairment [15].

Non-coding RNA (ncRNA) regulation is one of the essential epigenetic regulatory mechanisms. ncRNAs are functional RNA molecules that cannot be translated into proteins and have regulatory roles, mainly including microRNA (miRNA), long non-coding RNA (lncRNA), and circular RNA (circRNA) [16, 17]. In recent years, competing endogenous RNA (ceRNA) networks consisting of ncRNA and mRNA have gradually attracted the attention of researchers. The nature of the ceRNA network is that lncRNA or circRNA competes with corresponding mRNA and binds to the same miRNA, thereby participating in the regulation of gene expression at the post-transcriptional level [18, 19]. The primary ceRNA networks are circRNA-/lncRNA-miRNA-mRNA networks. Findings suggest that ceRNA networks may have critical regulatory roles in the development of RIR processes in multiple organs, including RIR [20,21,22]. For example, lncRNA Ttc3-209 is significantly up-regulated after RIR injury, which up-regulates the translation of Wnt8a mRNA through sponging miR-484, thereby promoting RGC cell apoptosis [20]. However, screening critical genes and their specific mechanism and functions involved in the development of RIR injury still needs more investigation.

The primary purpose of this study was to screen and identify the key genes in RIR injury by integrating mRNA, lncRNA and circRNA sequencing data, to provide theoretical support and reference basis for clinical treatment. In this study, we screened the differentially expressed circRNAs, lncRNAs and mRNAs between the control group and the model group, and within different reperfusion times (24h, 72h, 7d) using whole-transcriptome sequencing, and obtained the expression trends of the time-varying mRNAs, lncRNAs, and circRNAs by time-sequencing analysis. Then, candidate circRNAs, lncRNAs, and mRNAs were obtained by intersection of differentially expressed genes and tendency change genes. key genes whose expression changed with the prolongation of reperfusion time were selected by the importance scores of the genes. The biological pathways and potential regulatory mechanisms of the key genes were analyzed by bioinformatics, and the drugs associated with them as well as the specific molecular regulatory mechanisms during different periods of reperfusion were investigated. In addition, characteristic differentially expressed genes specific to the reperfusion time were analyzed, and key genes specific to the reperfusion time were screened to show the changes in biological processes with the prolongation of reperfusion time. The workflow of the study was listed in Fig. 1.

Fig. 1

The Workflow of this study

Material and methods

Construction of high intraocular pressure (IOP) model

High IOP model was constructed for whole transcriptome sequencing. 12 Sprague Dawley (SD) rats (Half male and half female, 220±20 g) were purchased from Beijing Silaikedake Laboratory Animal Technology Co., LTD. Rats were fed and managed in strict compliance with the Vision and Ophthalmology Research Association Statement and approved by the Ethics Committee of the First Affiliated Hospital of Kunming Medical University. SD rats were randomly divided into a control (rats was maintained at normal IPO levels), model_24h group (subjected to high IOP process and followed with reperfusion for 24h), model_72h group (subjected to high IOP process and followed with reperfusion for 72h), and model_7d group (subjected to high IOP process and followed with reperfusion for 7 days), with three rats in each group. As mentioned before [23, 24], all groups of rats were injected intraperitoneally with 100 mg/kg ketamine and 5 mg/kg xylazine to anesthetize rats for experimental injury, and 0.5% Alcaine ophthalmic solution was used to anesthetize the corneas of rats topically, and 1% tropicamide was used to dilate the pupils. Subsequently, rats in the model_24 h, model_72 h, and model_7 d groups were subjected to corneal cannulation to elevate IOP to 110 mmHg for 60 min to induce the ischemic process. Corneal cannulation in the control group of rats was maintained at normal IOP levels. Further, the eyes of rats in the model_24 h, model_72 h, and model_7 d groups performed reperfusion for 24 h, 72 h, and 7 d, respectively. Finally, the retinas from each group were collected for whole transcriptome sequencing.

RNA extraction and quality assessment

Total RNA was extracted from the retinal tissues of each group of rats using the mirVana miRNA Isolation Kit (Ambion, 1561) with reference to the manufacturer’s protocol. Subsequently, total RNA integrity was qualified and assessed by NanoDrop 2000 spectrophotometer (Thermo Fisher Scientific, MA, USA) and Agilent 2100 Bioanalyzer (Agilent Technologies, CA, USA). Samples with RNA integrity number ≥7 were subjected to subsequent analysis.

Library construction and high-throughput sequencing

RNA library construction is the process of converting RNA into double-stranded DNA recognizable by a second-generation sequencer through processes such as reverse transcription and splice junctions. Construction of RNA libraries from qualified RNA samples was complicated using TruSeq Stranded TotalRNA with Ribo-Zero Gold (Illumina, CA, USA). For lncRNA sequencing, rRNA was removed by Ribo-Zero Gold, and then libraries of lncRNA were prepared. For circRNA sequencing, rRNA was removed, and RNA was treated with RNase R (Epicentre, Madison, WI, USA) to remove linear transcripts and was then fragmented to approximately 200 bp. The purified RNA fragments were subjected to first- and second-strand cDNA synthesis following adaptor ligation according to the NEB Next Ultra™ RNA Library Prep Kit manual for Illumina (NEB, USA). Subsequently, these libraries were sequenced using an Illumina Hiseq 2500 sequencing platform (Illumina), and 150 bp paired-end reads were generated. The raw reads in FastQ format were further quality filtered by Trimmomatic [25] software for adapter removal and filtering. Clean reads were quality inspected by fastqc [26] and then aligned to the human reference genome using hisat2 [27]. The sequencing reads of each sample was aligned with the mRNA sequence by transcript4sequences, aligned with the known lncRNA sequences, and lncRNA prediction sequences by bowtie2. Data was imported to eXpress to make a gene quantitative analysis. The FPKM value and counts value were obtained. For circRNAs, we used BWA software [28] to align the sequencing reads of each sample with the reference genome. CIRI2 was used to scan the circRNA-paired chiastic clipping signal [29, 30]. The expression values were normalized by the reads per million (RPM) algorithm [31], and the number of junction reads counts and fold changes were normalized by DEseq2 in the R script. False discovery rate correction (adj. P<0.05) was used in the analysis.

Differentially expressed gene screening

The objective of differential expression analysis is to find genes that have significant differences in expression levels between samples. These samples can represent different biological states, such as drug-treated versus control, diseased versus healthy individuals, different tissues or different developmental stages. Principal component analysis (PCA) was a widely used data dimensionality reduction algorithm to obtain the overall degree of similarity or difference in gene expression values between samples, which was accomplished with scatterplot3d package (v0.3-43). The differentially expressed circRNAs, lncRNAs, and mRNAs between model_24h group and control group, model_72h group and control group, model_7d group, and control group, and model group (the combination of samples in model_24h group, model_72h group, and model_7d group) and control group were screened using the DEseq2 R package. The threshold value was set as |log2 (fold change, FC) | > 0.5 and P adj. < 0.05 for lncRNAs and mRNAs, |log2 (fold change, FC) | > 0.5 and P < 0.05 for circRNAs. The ggplot2 R package was utilized to draw volcano plots of the DEGs. The top 10 up-regulated genes and the top 10 down-regulated genes are labeled in the volcano plot sorted according to the ploidy of difference log2FC Heat maps were constructed using the heatmap package in the R script.

Identification of Chronological Expression Analysis

In order to identify the circRNAs, lncRNAs, and mRNAs that changed in the same trend with the increase of reperfusion time, the gene expression in each group was clustered using the R package Mfuzz based on fuzzy c-means clustering (FCM), with a desired number of clusters of 10. Mfuzz is a software tool that can be used to analyze gene changes over time. It can also help to understand how genes interact with each other in biological processes and how they affect each other's functions. The gene expression in control (with the reperfusion time as 0h), 24 h reperfusion, 72 h reperfusion, and 7 d reperfusion were set as point1, point2, point3 and point4. All the genes in the sequencing results with the same expression trend in different periods were analyzed by clustering, and then the trend genes were merged to obtain the trend change genes.

Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment

The GO database provides specialized terminology to define the properties of gene products, which encompasses biological processes (BP), cellular components (CC), and molecular functions (MF) [32]. BP represents a process of molecular activity events, including a functional collection of cells, tissues, organs, and species. CC indicates the cell or the external environment it is in. MF describes the active component of the gene product at the molecular level. KEGG is a database that systematically analyzes gene function and links genomic and functional information [33]. In this study, GO and KEGG pathway enrichment analysis of DEmRNAs was performed by the clusterProfiler package [34]. P <0.05 is considered statistically significant.

Protein-protein interaction network (PPI)

The PPI of candidate mRNAs and time-specific mRNAs were obtained by String (https://string-db.org/). Cystoscope was utilized to visualize the interconnections between differential genes using the degree of node importance (Degree) as a metric, which indicates the number of neighbors of the node, and highly connected nodes indicated a more critical role in the network. The CytoHubba was used to calculate the importance of the genes in the modules using the topology algorithm MCC.

Gene Set Enrichment Analysis (GSEA) enrichment analysis for candidate genes

GSEA was performed to detect the biological function, chromosomal location, or regulation of genes [35]. GSEA is a gene function analysis method that is used to reveal differences in gene expression under two different conditions. GSEA focuses on the co-regulatory patterns of an entire gene set, rather than on the expression levels of individual genes. Through statistical methods, GSEA is able to identify the enhancement or inhibition of a biological pathway or function, thereby revealing the underlying biological mechanism. Here, the correlation between key genes and other genes in the transcriptome sequencing was determined by Pearson correlation analysis to rank genes. The CP: KEGG dataset under C2 in the rat (Rattus norvegicus) species was extracted as a background gene set by the R package msigdbr and analyzed for GSEA enrichment. Finally, GSEA enrichment analysis was performed using the GSVA R package. The enrichplot R package was used to visualize the Top 5 pathways.

Regulatory networks and prediction

The competitive endogenous RNA (ceRNA) regulatory network consists of mRNAs, coding pseudogenes, long-chain noncoding RNAs, and miRNAs. ceRNA networks form a highly interconnected regulatory system by mediating the competitive binding of miRNAs between different RNA molecules. The significance of studying ceRNA networks lies not only in expanding the knowledge of gene regulatory mechanisms, but also in their potential applications in various biological and medical fields. To construct ceRNA networks, the miRNAs targeted to key genes were predicted by miRWalk (Ver. 2.0, http://mirwalk.umm.uni-heidelberg.de/). And the binding relationship of candidate lncRNAs and candidate circRNAs with miRNA was predicted by miRanda (http://www.microrna.org/ microrna/home.do). To find other genes associated with the key genes, GeneMANIA (http://www.genemania.org/) was used to predict the genes associated with the biomarker function and the functions involved. Subcellular localization pertains to the specific site where a protein or its expression product is found within a cell, such as the nucleus, cytoplasm, or cell membrane. Increasing evidence indicates that RNAs residing in various subcellular organelles exhibit distinct functionalities during biological processes, highlighting the significance of subcellular localization in unraveling the intricate biological functions of RNAs. For subcellular location, the base sequences of five key genes (Cd74, RT1-Da, RT1-CE5, RT1-Bb, RT1-DOa) were queried through the NCBI website (https://www.ncbi.nlm.nih.gov/). The base sequences of the genes were then entered into the mRNALocater (http://bio-bigdata.cn/mRNALocater) database to obtain predicted scores for the five subcellular localizations, with the highest score being used for the final specific, accurate localization. X2K (eXpression2Kinases, https://amp.pharm.mssm.edu/X2K/) is used to identify upstream regulators that may be responsible for the observed genome-wide gene expression patterns by integrating information based on transcriptome data and existing signal transduction databases. To explore the potential regulatory mechanisms of key genes, proteins that interact with transcription factors and upstream kinases that regulate key genes were analyzed by eXpression2Kinases. SIGNOR (http://signor.uniroma2.it/) is a manually annotated database of causal relationships between human proteins, biologically relevant chemicals, stimuli, and phenotypes. In order to analyze the signaling relationships of critical genes further, key genes were annotated through the SIGNOR database to construct an essential gene signaling information network. The DGIdb (https://www.dgidb.org) database was utilized to search for drugs targeting key genes for potential drug discovery of biomarkers and to construct a drug-gene interaction network. Since the database is for human-targeted genes, we used the R package "homologene" to convert rat key genes into homologous human genes and then searched the DGIdb website. To analyze the role of key genes in other eye diseases, the Comparative Toxicogenomics Database was used (CTD, https://ctdbase.org/) to analyze the relationship between ocular diseases and key genes. Bar charts were drawn to show the top five diseases with the highest predicted scores for each key gene respectively.

Real-tine quantitative reverse transcription (RT-qPCR) assay

The mRNA expression of key genes was detected by RT-qPCR assay. The extracted total RNA was used as a template, and reverse transcription was carried out using an RNA first-strand cDNA synthesis kit. The relative expression levels of mRNA were determined by the qPCR method. qPCR reaction conditions were as follows: 95℃ 3 min, 95℃ 15s, 60℃ 35s, 72℃ 30 s, a total of 40 cycles. GAPDH was used as the internal reference of mRNA, and the relative expression level was calculated by the 2^-ΔΔCt method. The primers were listed in Table 1.

Table 1 Primer sequences for qPCR assay

Statistical analysis

GraphPad Prism 8.0 was used to analyze the data and graph. Data were presented as mean ± standard deviation ($\overline{x }\pm SD$) and Turkey’s test was used for comparison between the two groups. P<0.05 represents the difference is significantly significant.

Results

Differential expressed genes screening

To clarify the gene expression changes in RIR injury at the transcriptome level, we constructed a high IOP model (24 h, 72 h, and 7 d of reperfusion, respectively) and performed whole transcriptome sequencing. The principal component analysis (PCA) results revealed that the samples from the high IOP model with different reperfusion time exhibited significant differences in mRNA, lncRNA, and circRNA expression. (Supplementary Figure 1). As a result, 3020 differentially expressed mRNA (DEmRNA1, 919 upregulated and 2102 downregulated mRNAs), 862 differentially expressed lncRNAs (DElncRNA1, 280 upregulated and 582 downregulated lncRNAs), and 44 circRNAs (DEcircRNA1, 20 upregulated and 22 downregulated circRNAs) were screened between model_24 h and control group (criteria: |log2FC|>1 and adj.P<0.05 for mRNA and lncRNAs, and |log2FC|>1 and P<0.05 for circRNA). The gene expression distribution of differentially expressed genes was distributed as heat maps (Fig. 2A-F, Supplementary Tables 1, 2 and 3). Also, 2334 differentially expressed mRNAs (DEmRNA2, 563 upregulated and 1771 downregulated mRNAs), 369 differentially expressed lncRNAs (DElncRNA2, 110 upregulated and 259 downregulated lncRNAs), and 30 circRNAs (DEcircRNA2, 18 upregulated and 12 downregulated circRNAs) were screened between model_72 h and control group. The gene expression distribution of differentially expressed genes was distributed as heat maps (Fig. 2G-L, Supplementary Tables 4, 5 and 6). 1195 differentially expressed mRNAs (DEmRNA3, 46 upregulated and 1149 downregulated mRNAs), 170 differentially expressed lncRNAs (DElncRNA3, 33 upregulated and 137 downregulated lncRNAs), and 36 circRNAs (DEcircRNA3, 13 upregulated and 23 downregulated circRNAs) were screened between model_7d and control group. The gene expression distribution of differentially expressed genes was distributed as heat maps (Fig. 2M-R, Supplementary Tables 7, 8 and 9). To identify differential mRNAs, lncRNAs, and circRNAs between retinal tissue from the model group versus those from the control group, the model group was set as a collection of models in different reperfusion time. 2098 differentially expressed mRNAs (DEmRNA4, 1906 upregulated and 192 downregulated mRNAs), 389 differentially expressed lncRNAs (DElncRNA4, 349 upregulated and 40 downregulated lncRNAs), and 76 circRNAs (DEcircRNA4, 63 upregulated and 13 downregulated circRNAs) were screened between model and control group. The gene expression distribution of differentially expressed genes was distributed as heat maps (Fig. 2S-T, Supplementary Tables 10, 11 and 12).

Fig. 2

Differentially expressed mRNAs, lncRNAs, circRNA in high IOP model. A The volcano plot revealed DEmRNAs between model_24 h and control group. B The heatmap revealed DEmRNAs between model_24 h and control group. C The volcano plot revealed DElncRNAs between model_24 h and control group. D The heatmap revealed DElncRNAs between model_24 h and control group. E The volcano plot revealed DEcircRNAs between model_24 h and control group. F The heatmap revealed DEcircRNAs between model_24 h and control group. G The volcano plot revealed DEmRNAs between model_72 h and control group. H The heatmap revealed DEmRNAs between model_72 h and control group. I The volcano plot revealed DElncRNAs between model_72 h and control group. J The heatmap revealed DElncRNAs between model_72 h and control group. K The volcano plot revealed DEcircRNAs between model_72 h and control group. L The heatmap revealed DEcircRNAs between model_72 h and control group. M The volcano plot revealed DEmRNAs between model_ 7d and control group. N The heatmap revealed DEmRNAs between model_ 7d and control group. O The volcano plot revealed DElncRNAs between model_ 7d and control group. P The heatmap revealed DElncRNAs between model_ 7d and control group. Q The volcano plot revealed DEcircRNAs between model_ 7d and control group. R The heatmap revealed DEcircRNAs between model_ 7d and control group

Chronological analysis of differentially expressed genes

On the other side, chronological changes in gene expression with increasing reperfusion time were detected by the Mfuzz algorithm. The R package Mfuzz was used to perform cluster analysis of the model and control groups to identify the RNAs with trend changes in samples from different reperfusion time. All the RNAs obtained from the sequencing results with the same expression trend in different reperfusion periods (control group, model_24h group, model_72h group, and model_7d group) were clustered. Then, the trend RNAs were concatenated to obtain the trend change RNAs. As a result, mRNA expression in cluster 1 and 7 were increased with the increase of reperfusion time, and therefore, the genes in the clusters were merged, and 3733 trend mRNAs were obtained (Fig. 3A, Supplementary Tables 13 and 14). The lncRNA expression in cluster 5 and 9 were increased, while lncRNA expression in cluster 3 were decreased with the increase of reperfusion time, and 4584 trend lncRNAs were obtained (Fig. 3B, Supplementary Tables 15 and 16). The circRNA expression in cluster 7 and 9 were increased, while circRNA expression in cluster 8 and 10 were decreased with the increase of reperfusion time, and 4635 trend circRNAs were obtained (Fig. 3C, Supplementary Tables 17 and 18). Then, 316 candidate mRNAs, 137 candidate lncRNAs, and 31 candidate circRNAs were obtained by the intersection of differentially expressed mRNAs, lncRNAs, and circRNAs with trend mRNAs, trend lncRNAs and trend circRNAs (Fig. 4A-C, Supplementary Tables 20 and 21). To evaluate the biological functions that candidate genes enriched, GO annotation and KEGG enrichment was preformed. Candidate mRNAs were enriched in 686 GO terms (629 in BP terms, 21 in CC terms, and 36 in MF terms), such as positive regulation of cell activation, adaptive immune response, positive regulation of leukocyte activation, external side of plasma membrane, extracellular matrix, external encapsulating structure, immune receptor activity, cytokine binding, peptide binding, involving genes such as A2m, Acp5, C1qc, Acta2, Bgn, Cd180, Cd1d1, Cd74, Clu, Also, candidate mRNAs were enriched in 57 KEGG terms, such as Cell adhesion molecules, Cytokine-cytokine receptor interaction, Epstein-Barr virus infection, Human T-cell leukemia virus 1 infection, Coronavirus disease-COVID-19, including genes such as RT1-Da, Icoslg, RT1-Db1, RT1-Ba, RT1-Bb, Icos (Fig. 4D-F, Table 2, Supplementary Tables 22 and 23).

Fig. 3

chronological analysis of differentially expressed genes. A mRNA expression results at different times. B lncRNA expression results at different times. C circRNA expression results at different times. Note: Red and purple represent genes with high MEM.SHIP values, and yellow and green color represent genes with low MEM.SHIP values. The horizontal axis indicates the rat retinal tissue samples at different reperfusion time, and the vertical axis indicates the changes in gene expression

Fig. 4

Candidate gene selection and their function enrichment. A Candidate mRNAs were obtained by the intersection of differentially expressed mRNAs, with trend mRNAs. B Candidate lncRNAs were obtained by the intersection of differentially expressed lncRNAs, with trend lncRNAs. C Candidate circRNAs were obtained by the intersection of differentially expressed circRNAs, with trend circRNAs. D Scatterplot showed GO enrichment results of candidate mRNAs. Horizontal coordinate is the enrichment factor, vertical coordinate is the name of the enriched pathway, dot size indicates the number of differential genes enriched into the pathway, color indicates the range of p.adjust. E Chord diagrams of KEGG terms of candidate mRNAs. The color of the left gene ribbon represents the logFC of the gene, and different ribbons on the right represent different pathways. F scatterplot showed KEGG enrichment results of candidate mRNAs. Horizontal coordinate is the enrichment factor, vertical coordinate is the name of the enriched pathway, dot size indicates the number of differential genes enriched into the pathway, color indicates the range of p.adjust

Table 2 Top 10 GO enrichment and Top 5 KEGG enrichment of Candidate mRNAs

Key gene identification and the prediction of regulation

In order to understand the interactions between proteins encoded by candidate genes, 316 candidate mRNAs were uploaded to the STRING database for PPI network construction (interaction score >0.7), the protein-protein interaction network that contained 297 nodes and 281 edges was obtained (average node degree: 1.89). Then, PPI was imported to Cytoscape (Ver. 3.7.2), and the importance scores of the genes were obtained by the topology algorithm MCC to identify key genes. 5 key genes (Cd74, RT1-Da, RT1-CE5, RT1-Bb, RT1-DOa) were identified as a result (Fig. 5, Supplementary Tables 24 and 25).

Fig. 5

The protein- protein interaction of candidate mRNAs. Left: key genes obtained by MCC algorithm; Right: The protein- protein interaction of candidate mRNAs

Then, GSEA was used to evaluate the function of key genes. Cd74 was found to be positively related to ribosome, antigen processing, and presentation, regulation of actin cytoskeleton (Fig. 6a, Supplementary Table 34). RT1-Bb was found to be positively related to ribosome, antigen processing, and presentation (Fig. 6b, Supplementary Table 26). RT1-CE5 was found to be positively related to ribosome (Fig. 6C, Supplementary Table 27). RT1-Da was found to be positively related to ribosome, antigen processing, and presentation, regulation of actin cytoskeleton (Fig. 6D, Supplementary Table 28). RT1-DOa was found to be negatively related to ribosome (Fig. 6E, Supplementary Table 29).

Fig. 6

GSEA result for key genes. A Cd74 gene enrichment signaling pathway. B RT1-Bb gene enrichment signaling pathway. C RT1-CE5 gene enrichment signaling pathway. D RT1-Da gene enrichment signaling pathway. E RT1-Doa gene enrichment signaling pathway

To predicted the targeted miRNAs and ncRNA with key genes, ceRNA networks was constructed. A network included 4 key genes (Cd74, RT1-Da, RT1-Bb, RT1-DOa), 34 miRNAs and 48 lncRNAs and 81 regulatory relationship axes (32 related with Cd74, 12 related with RT1-Da, 23 related with RT1-Bb, 14 related with RT1-DOa) was constructed (Fig. 7A, Supplementary Table 31). Also, a network included 4 key genes (Cd74, RT1-Da, RT1-Bb, RT1-DOa), 9 miRNAs and 3 circRNAs (circRNA_10572, circRNA_03219, circRNA_11359) and 12 regulatory relationship axes (4 related with Cd74, 1 related with RT1-Da, 5 related with RT1-Bb, 2 related with RT1-DOa) was constructed (Fig. 7B, Supplementary Table 32). RT1-CE5 was not included in the networks since no miRNA targeted was predicted in the miRWalk database. GeneMANIA was utilized to predict the genes associated with the key genes and the functions involved to find the genes associated with the key genes. 20 related genes with 1144 total links were predicted (Physical Interactions accounted for 77.64%, Co-expression 8.01%, Predicted 5.37%, Co-localization 3.63%, Genetic Interactions 2.87%, Pathway 1.88%. Shared protein domains accounted for 0.6%). These related genes were enriched in 98 functions, including antigen processing and presentation of peptide antigen, antigen processing and presentation, MHC protein complex, antigen binding, antigen processing and presentation of exogenous peptide antigen, antigen processing and presentation of exogenous antigen, antigen processing and presentation of peptide antigen via MHC class II (Fig. 7C, Supplementary Tables 33 and 34).

Fig. 7

Regulation prediction of key genes. A Key mRNA-miRNA -lncRNA Networks. Red dots are mRNAs, yellow triangles are miRNAs, and blue squares are lncRNAs. B Key mRNA-miRNA-circRNA network. Red dots are mRNAs, yellow triangles are miRNAs, and green squares are circRNAs. C GeneMANIA Network. The middle is the key gene, the outer circle is the related genes with similar functions to the key gene, different colors connecting lines indicate different networks, and different colors in the pie chart indicate different functions of the genes. D Subcellular localization of key genes. Horizontal coordinates are genes, vertical coordinates are scoring at different sites, with higher scores indicating a higher likelihood of being at that site. E Regulatory network of transcription factors, kinases and related proteins. Transcription factors in red, kinases in green, and proteins in gray. F SIGNOR-based signaling information network for key genes. G Gene-Drug Network. Genes are shown in red and drug names in blue. H Gene-Disease Prediction. Horizontal coordinate is the prediction score, vertical coordinate is the predicted disease name, and different color bars indicate different key genes

mRNALocater was used to predict the subcellular location of key genes. Based on the 5 highest scores for subcellular localization, RT1-CE5 localized to the nucleus, Cd74, RT1-Da, RT1-Bb, and RT1-DOa localized to the cytoplasm (Fig. 7D, Supplementary Table 35). In order to explore the potential regulatory mechanisms of key genes, transcription factors, proteins that interact with transcription factors, and upstream kinases that regulate key genes were analyzed by eXpression2Kinases. The results showed that a total of 10 transcription factors (ZNF778, AEBP2, GCM2, ADNP, ESR2, ZNF577, NR3C2, ZNF516, MESP2, HES1), 28 proteins interacting with transcription factors, and 10 upstream kinases (BUB1B, CDK2 CSNK2A1, AKT1, HIPK2. PINK1, AURKB, ABL1, CDK7, PDPK1) were predicted (Fig. 7E, Supplementary Table 36). In order to analyze the signaling relationship of key genes further, the key genes were annotated by the SIGNOR database to construct the key gene signaling network. The results showed that the proteins CIITA and RFX5 were transcriptional activators of RT1-Da, the complex RFX was also transcriptional activator of RT1-Da, the protein NFX1 was transcriptional repressor of RT1-Da, the proteins Marchf1 and March 8 were repressor of RT1-Da, and the complex EBVgH:gL: gp42 was activator of RT1-Da. gp42 activates RT1-Da protein. At the same time, RFX also activates RT1-Bb and RT1-DOa proteins (Fig. 7F, Supplementary Table 37).

Also, DGIdb predicted the drugs that targeted key genes, and the network was constructed. A total of 27 drugs were predicted to be targeted with 4 key genes (Cd74, RT1-Da, RT1-CE5, RT1-Bb), Cd74 was predicted as the target of 1 drug (MILATUZUMAB), RT1-Da was predicted as the target of 6 drugs, RT1-CE5 was predicted as the target of 15 drugs, RT1-Bb was predicted as the target of 12 drugs. Also, 2 drugs (CLAVULANIC ACID, AMOXICILLIN) were predicted to be common targeted to RT1-Da, RT1-CE5, and RT1-Bb. 2 drugs (CARBAMAZEPINE, TICLOPIDINE) were predicted to be common targeted to RT1-CE5, RT1-Bb. 1 drug (FLOXACILLIN) was predicted to be common targeted to RT1-Da RT1-Bb (Fig. 7G, Supplementary Table 38).

To analyze the role of key genes in eye diseases, the relationship between eye diseases and key genes was analyzed using the Comparative Toxicogenomics Database, and bar graphs were plotted to show the top five diseases with the highest predicted scores for each key gene, respectively. The results showed that the Cd74 gene was mainly associated with diseases such as Retinal Diseases, Eye Abnormalities, Cataract, and Vision Disorders. RT1-Bb Gene was mainly associated with Cataract, Retinal Degeneration, Eye Abnormalities, and Dry Eye Syndromes. RT1-CE5 genes were mainly associated with Eye Abnormalities, Pathologic Nystagmus, Vision Disorders, Retinoschisis, and Horner Syndrome. The RT1-Da gene was mainly associated with Vision Disorders, Optic Nerve Diseases, Cataract, Retinal Diseases, and Eye Abnormalities. The RT1-Da gene was mainly associated with Vision Disorders, Optic Nerve Diseases, Cataract, Retinal Diseases, and Eye Abnormalities. The RT1-DOa gene was mainly associated with Retinal Diseases, Lens Subluxation, Night Blindness, Uveitis, and Anterior and Corneal Injuries. In conclusion, key genes were mainly associated with Retinal diseases, Eye abnormalities, Cataract, and Vision disorders (Fig. 7H, Supplementary Table 39). In order to observe the changes in the expression of key genes over time, vertical scatter plots were plotted to show the expression levels of the key genes in the samples from transcriptome data from self-sequencing (mRNA) among different groups.

The expression of the key genes was detected in transcriptome data and verified by qPCR to observe the changes in the expression among different reperfusion time. The results showed that the expression of Cd74, RT1-Da, and RT1-CE5 was all increased after reperfusion for 24h, 72h, and 7d compared to the control group. RT1-Bb was significantly higher expressed in model_72h group and model_7d compared to the control group. RT1-DOa was significantly higher expressed in model_7d compared to the control group (Fig. 8 A-J).

Fig. 8

Expression of key genes. A-E Expression of key genes across groups in transcriptome data. F-J Expression of key genes across groups were detected in the retinal samples by qPCR

Identification of time-specific genes and their regulatory networks

In order to obtain the characteristic differentially expressed genes specific to the reperfusion time, Venn plots including DEmRNA1 (model_24h group vs. control), DEmRNA2 (model_72h group vs. control), DEmRNA3 (model_7d group vs. control) was plotted, 1370 characteristic differentially expressed mRNAs (spec_24h mRNA) were found in model_24h group, 558 characteristic differentially expressed mRNAs (spec_72h mRNA) were found in model_72h group, 92 characteristic differentially expressed mRNAs (spec_7d mRNA) were found in model_7d group (Fig. 9A, Supplementary Table 40). 612 characteristic differentially expressed lncRNAs (spec_24h lncRNA) were found in model_24h group, 127 characteristic differentially expressed lncRNAs (spec_72h lncRNA) were found in model_72h group, 38 characteristic differentially expressed lncRNAs (spec_7d lncRNA) were found in model_7d group (Fig. 9B, Supplementary Table 41). Also, 36 characteristic differentially expressed circRNAs (spec_24h circRNA) were found in the model_24h group, 22 characteristic differentially expressed circRNAs (spec_72h circRNA) were found in the model_72h group, and 31 characteristic differentially expressed circRNAs (spec_7d circRNA) were found in the model_7d group (Fig. 9C, Supplementary Table 42).

Fig. 9

Time-specific differentially expressed gene identification. A Specific DEmRNA identification. B Specific DElncRNA identification. C Specific DEcircRNA identification

The GO and KEGG enrichment evaluated the function of characteristic differentially expressed mRNAs. For spec_24h mRNA, genes were enriched in 385 GO terms (241 in BP terms, 84 in CC terms, and 33 in MF terms), such as ribosome, ribosomal subunit, camera-type eye development, retina development in camera-type eye, cellular response to chemical stress, structural constituent of ribosome, signaling receptor activator activity, including genes such as Rps23, Rpl3, Rpl5, Pax4, Lpcat1. Also, spec_24h mRNAs were enriched in 12 KEGG terms such as Coronavirus disease-COVID-19, Ribosome, MAPK signaling pathway, Cytokine-cytokine receptor interaction, Focal adhesion, involving genes such as Rpl7, Rpl3, Rpl5, Rpl4, Rps15a, (Fig. 10A-C, Table 3, Supplementary Tables 43 and 44). Then, a PPI was constructed to investigate the interaction of the model. For spec_72h mRNA, genes were enriched in 264 GO terms (144 in BP terms, 88 in CC terms, and 62 in MF terms), such as regulation of membrane potential, cell junction assembly, regeneration, synaptic membrane, intrinsic component of synaptic membrane, ion channel activity, channel activity, passive transmembrane transporter activity, involving genes such as Gabrb2, Cacna1c, Sparcl1, Spock3, Kctd16, Dcc. Also, 25 KEGG terms were enriched to spec_72h mRNA, such as Neuroactive ligand-receptor interaction, Retrograde endocannabinoid signaling, cAMP signaling pathway, Calcium signaling pathway, Serotonergic synapse, including genes such as Kcnq3, Gad2, Gad1, Pde3a, Creb5, Ppp2r2c (Fig. 10D-F, Table 4, Supplementary Tables 45 and 46). For spec_7d mRNA, genes were enriched in 89 GO terms (73 in BP terms, 12 in CC terms, and 4 in MF terms), such as adaptive immune response, lymphocyte differentiation, mononuclear cell differentiation, external side of the plasma membrane, extracellular matrix, external encapsulating structure, extracellular matrix structural constituent, cytokine receptor activity, MHC class II protein complex binding, including genes such as C4a, RT1-DOa, Cd3g, Pou2af1, Ikzf3, Tectb, Dcn, Zp2, RT1-DMa. 25 KEGG terms were found enriched to spec_7d mRNA, such as Th1 and Th2 cell differentiation, Th17 cell differentiation, Epstein-Barr virus infection, Intestinal immune network for IgA production, Hematopoietic cell lineage, including genes such as Cd2, Cd79a, Tap1, Cd28, Tnfrsf13b (Fig. 10G-I, Table 5, Supplementary Tables 47 and 48).

Fig. 10

Function enrichment of time-specific differentially expressed gene. A Scatterplot showed GO enrichment results of 24h_specific mRNAs. Horizontal coordinate is the enrichment factor, vertical coordinate is the name of the enriched pathway, dot size indicates the number of differential genes enriched into the pathway, color indicates the range of p.adjust. B Chord diagrams of KEGG terms of 24h_specific mRNAs. The color of the left gene ribbon represents the logFC of the gene, and different ribbons on the right represent different pathways. C scatterplot showed KEGG enrichment results of 24h_specific mRNAs. Horizontal coordinate is the enrichment factor, vertical coordinate is the name of the enriched pathway, dot size indicates the number of differential genes enriched into the pathway, color indicates the range of p.adjust. D Scatterplot showed GO enrichment results of 72h_specific mRNAs. Horizontal coordinate is the enrichment factor, vertical coordinate is the name of the enriched pathway, dot size indicates the number of differential genes enriched into the pathway, color indicates the range of p.adjust. E Chord diagrams of KEGG terms of 72h_specific mRNAs. The color of the left gene ribbon represents the logFC of the gene, and different ribbons on the right represent different pathways. F scatterplot showed KEGG enrichment results of 72h_specific mRNAs. Horizontal coordinate is the enrichment factor, vertical coordinate is the name of the enriched pathway, dot size indicates the number of differential genes enriched into the pathway, color indicates the range of p.adjust. G Scatterplot showed GO enrichment results of 7d_specific mRNAs. Horizontal coordinate is the enrichment factor, vertical coordinate is the name of the enriched pathway, dot size indicates the number of differential genes enriched into the pathway, color indicates the range of p.adjust. H Chord diagrams of KEGG terms of 7d_specific mRNAs. The color of the left gene ribbon represents the logFC of the gene, and different ribbons on the right represent different pathways. I scatterplot showed KEGG enrichment results of 7d_specific mRNAs. Horizontal coordinate is the enrichment factor, vertical coordinate is the name of the enriched pathway, dot size indicates the number of differential genes enriched into the pathway, color indicates the range of p.adjust

Table 3 Top 10 GO enrichment and Top 5 KEGG enrichment of characteristic differentially expressed mRNA in 24h reperfusion (spec_24h mRNA)

Table 4 Top 10 GO enrichment and Top 5 KEGG enrichment of characteristic differentially expressed mRNA in 72h reperfusion (spec_72h mRNAs)

Table 5 Top 10 GO enrichment and Top 5 KEGG enrichment of characteristic differentially expressed mRNA in 7d reperfusion (spec_7d mRNAs)

Then, PPI was constructed using STRING and visualized by Cytoscape to investigate the interactions in time-specific genes in 24h as reperfusion time. A network with 1150 nodes and 8020 edges was constructed (Fig. 11A, Supplementary Table 49). The top 10 key genes (Rps14, Rack1, Rpl10a, Rpl24, Rpl19, Rps23, Rps16. Uba52, Rpl27, Rpl18) in spec_24h mRNAs were selected by MCC algorithm (Fig. 11B, Supplementary Table 50). The targeted miRNAs, lncRNAs, and circRNAs of key genes were predicted, and the interaction networks were constructed. As a result, a network included 2 key genes (Rack1, Rps16), 2 miRNAs and 8 lncRNAs, and 8 regulatory relationship axes (7 related with Rack1, 1 related with Rps16) were constructed (Fi. 11C, Supplementary Table 51). Also, a network that included 2 key genes (Rack1, Rps16), 2 miRNAs, 6 circRNAs, and 6 regulatory relationship axes (5 related with Rack1, 1 related with Rps16) was constructed (Fig. 11D, Supplementary Table 52). The results of GSEA showed that all the key genes were to be positively related to ribosome, eukaryotic translation elongation, nonsense-mediated decay, selenoamino acid metabolism, cytoplasmic ribosomal proteins (Supplementary Figure 2A, Supplementary Tables 53, 54, 55, 56, 57, 58, 59, 60, 61 and 62).

Fig. 11

Analysis of the regulatory mechanism between the model_24h group and the control group of the model. A The protein- protein interaction of 24h_specific mRNAs. B key genes obtained by MCC algorithm. C Key mRNA-miRNA -lncRNA Networks. Red dots are mRNAs, yellow triangles are miRNAs, and green squares are lncRNAs. D Key mRNA-miRNA-circRNA network. Red dots are mRNAs, yellow triangles are miRNAs, and blue squares are circRNAs

A PPI of time-specific genes in 72h as reperfusion time was also constructed, and the network included 396 nodes and 1077 edges (Fig. 12A, Supplementary Table 63). Top 10 key genes (Cdca3, Ttk, Nuf2, Cenpf, Cenpe, Nusap1, Kif20b, Cenpa, Kif20a, Plk1) in spec_72h mRNA were selected by MCC algorithm (Fig. 12B, Supplementary Table 64). The regulation network was constructed to predict the targets of key genes. A network included 7 key genes (Cdca3, Ttk, Cenpf, Nusap1, Kif20b, Cenpa, Plk1), 21 miRNAs and 385 lncRNAs and 616 regulatory relationship axes (4 related with Cdca3, 3 related with Ttk, 43 related with Cenpf, 35 related with Nusap1, 53 related with Kif20b, 473 related with Cenpa, 5 related with Plk1) was constructed (Fig. 12C, Supplementary Table 65). As for circRNAs, a network included 5 key genes (Cdca3, Cenpf, Nusap1, Kif20b, Cenpa), 18 miRNAs and 345 circRNAs and 387 regulatory relationship axes (6 related with Cdca3, 33 related with Cenpf, 26 related with Nusap1, 34 related with Kif20b, 288 related with Cenpa,) was constructed (Fig. 12D, Supplementary Table 66). The key genes were positively related to ribosome, cell adhesion molecules, complement and coagulation cascades, cytokine-cytokine receptor interaction, and focal adhesion (Supplementary Figure 2B, Supplementary Tables 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77 and 78).

Fig. 12

Analysis of the regulatory mechanism between the model_72h group and the control group of the model. A The protein- protein interaction of 72h_specific mRNAs. B key genes obtained by MCC algorithm. C Key mRNA-miRNA -lncRNA Networks. Red dots are mRNAs, yellow triangles are miRNAs, and green squares are lncRNAs. D Key mRNA-miRNA-circRNA network. Red dots are mRNAs, yellow triangles are miRNAs, and blue squares are circRNAs

As for time-specific genes in 7 days as reperfusion time, a network that included 45 nodes and 83 edges was constructed (Fig. 13A, Supplementary Table 77). The top 10 key genes (Cd3e, Zap70, Cd2, Cd28, Cd3g, Cd247, Cd79a, Cxcr5, Cxcr3, Pou2af1) in spec_7d mRNA were selected by MCC algorithm (Fig. 13B, Supplementary Table 78). A network included 6 key genes (Cd3e, Zap70, Cd28, Cd247, Cxcr3, Pou2af1), 37 miRNAs and 379 lncRNAs and 752 regulatory relationship axes (17 related with Cd3e, 174 related with Zap70, 29 related with Cd28, 272 related with Cd247, 87 related with Cxcr3, 173 related with Pou2af1) was constructed (Fig. 13C, Supplementary Table 79). As for circRNAs, a network included 6 key genes (Cd3e, Zap70, Cd28, Cd247, Cxcr3, Pou2af1), 37 miRNAs and 306 circRNAs and 426 regulatory relationship axes (19 related with Cd3e, 84 related with Zap70, 12 related with Cd28, 150 related with Cd247, 47 related with Cxcr3, 114 related with Pou2af1) was constructed (Fig. 12D, Supplementary Table 80). The key genes were mainly positively related to cytokine-cytokine receptor interaction, ribosome, JAK/STAT signaling pathway, leishmania infection, and hemopoietic cell lineage (Supplementary Figure 2C, Supplementary Tables 81, 82, 83, 84, 85, 86, 87, 88, 89 and 90).

Fig. 13

Analysis of the regulatory mechanism between the model_7d group and the control group of the model. A The protein- protein interaction of 7d_specific mRNAs. B key genes obtained by MCC algorithm. C Key mRNA-miRNA -lncRNA Networks. Red dots are mRNAs, yellow triangles are miRNAs, and green squares are lncRNAs. D Key mRNA-miRNA-circRNA network. Red dots are mRNAs, yellow triangles are miRNAs, and blue squares are circRNAs

Discussion

Retinal ischemia/reperfusion (RIR) injury is a remarkably complex pathophysiological process that is widely seen in a variety of ocular diseases, such as retinal vascular occlusions, glaucoma, diabetic retinopathy, and retinopathy of prematurity, which can lead to blindness [15, 36]. The nature of RIR is that the blockage and the subsequent restoration of blood flow to the tissues induces a series of oxidative stress and inflammatory effects, which ultimately leads to damage to retinal nerve cells, especially retinal ganglion cells (RGCs) [11, 15, 37]. Recently, the pathogenesis of RIR has been poorly understood, and no clinically approved drugs can effectively rescue ischemic retinal nerve cells. Currently, treatments for RIR damage include the counteraction of oxygen radicals and oxidative stress, inhibition of calcium overload, inhibition of apoptosis, inhibition of inflammatory responses and reduction of retinal edema, and counteraction of neurotoxicity of nitric oxide and excitatory amino acids [11, 38,39,40]. Despite the many methods of treating RIR damage, treatment outcomes remain suboptimal. The tissue damage caused by transient ischemic injury is an essential component of the pathogenesis of retinal ischemia, which mainly hinges on the degree and duration of interruption of the blood supply and the subsequent damage caused by tissue reperfusion [41]. Some research indicated that the retinal injury induced by ischemia/reperfusion (I/R) was related to reperfusion time. Wang et al. found that retinal edema was seen in the early stage and followed by retina atrophied gradually in 72 h and 144 h as reperfusion time [42]. Zhang et al. found loss of cells in the retinal ganglion cell layer was apparent 2 days after I/R injury, and the number of degenerated capillaries increased greatly by 7 to 8 days after the injury [43]. These available studies show that the pathological changes in the retina after different reperfusion times vary considerably. This may indicate that the treatment strategy for retinal ischemia should be adapted according to the reperfusion time.

In this manuscript, we screened the differentially expressed circRNAs, lncRNAs, and mRNAs between the control and model groups and at different reperfusion time (24h, 72h, and 7d) with the aid of whole transcriptome sequencing technology, and the trend changes in time-varying mRNA, lncRNA, circRNA were obtained by chronological analysis. Then, candidate circRNAs, lncRNAs, and mRNAs were obtained as the intersection of differentially expression genes and trend change genes. The candidate mRNAs were mainly enriched in immune-related terms such as lymphocyte-mediated immunity, complement and coagulation cascades, Hematopoietic cell lineage, antigen processing, and expression. The blood-retinal barrier consists of the inner tight junction between retinal capillary endothelial cells and the outer tight junction between the RPE [44, 45]. On the one hand, due to the structure of the blood-retinal barrier, the macromolecular antibodies in the retinal vessels and choroidal vessels cannot play their functions. On the other hand, there are no lymphoid tissues in the retina, so the antigens do not cause the clonal proliferation of specific T cells or B cells. The retina, therefore, has long been recognized as a privileged site for immunity [11]. However, the immunologic response to various stress cues has been found to play a pivotal role in the retina [13]. When ischemia and reperfusion occur, the permeability of the blood-retinal barrier is changed, and the activation of microgila [46, 47] and the increase of Treg cells [48, 49] are observed. Toll signaling activation was also found and induced inflammasome formation [50]. Our results of candidate gene enrichment also showed that immune response was crucial in the process of RIR.

Importance scores of the genes selected the key genes whose expression changed with the increase of reperfusion time. As a result, 5 key genes, Cd74, RT1-Da, RT1-CE5, RT1-Bb, and RT1-DOa, were selected. Cd74 is a receptor for the cytokine macrophage migration inhibitory factor (MIF) [51]. CD74 regulates T-cell and B-cell development, dendritic cell (DC) motility, macrophage inflammation, and thymic selection [52]. Cd74, RT1-Ba, RT1-Bb, RT1-Da, and RT1-Db1 were referred to as major histocompatibility complex (MHC) class II members [53]. The MHC-II molecule is a central molecule in the protein presentation pathway. It binds to processed short peptides and presents them to T lymphocytes, activating them to become effector T cells [54]. Abcouwer et al. [55] Minocycline was particularly effective in decreasing the appearance of MHCII+ inflammatory leukocytes and reduced leukocyte adhesion and invasion, as well as vascular permeability in RIR. The result of GSEA showed that key genes were found to play vital roles in antigen processing and presentation, regulation of the actin cytoskeleton, and the ribosome. The research from Honjo [56] showed that the cytokines and growth factor in the aqueous humor activate Rho after the ischemia happened, and the Rho/ROCK signal transduction participates in RIR injury via rearrangement of the actin cytoskeleton that was attributed to improved outflow. With the increase in reperfusion time, key gene expression also increased. Also, key genes were mainly associated with Retinal diseases, Eye abnormalities, Cataract, and Vision disorders. Clavulanic acid and Amoxicillin were predicted to be commonly targeted to RT1-Da, RT1-CE5, and RT1-Bb, which further proved that the key genes might participate in the process of RIR injury by involving the adaptive immune response.

Also, the characteristic differentially expressed genes specific to the reperfusion time were analyzed. Key genes specific to reperfusion time were selected to show the change in biological process with the increased reperfusion time. Intriguingly, the specific genes in 24h as reperfusion time were mainly enriched in the pathways that related to retinal function and cellular response to external stress, such as camera-type eye development, retina development in camera-type eye, visual perception, detection of abiotic stimulus, detection of light stimulus. However, the specific genes in 72h as reperfusion time were enriched in the pathways related to membrane potential and neurodevelopment, such as regulation of postsynaptic membrane potential, neurotransmitter receptor activity, and ion channel activity. The specific genes in 7 days as reperfusion time were enriched in the pathways related to immune responses, such as adaptive immune response, lymphocyte differentiation, and positive T cell selection. The retina is a highly specialized neural tissue that continues the central nervous system. The existing research on the novel drug for RIR injury was mainly targeted to the retinal ganglion cell damage and inflammation. Guan et al. [57] Puerarin can ameliorate RIR injury by suppressing apoptosis and TLR4/NLRP3 inflammasome activation in RGCs. Lee et al. [58] proved Nicotinamide mononucleotide significantly suppressed retinal functional damage, as well as inflammation. Our results indicated that the biological process in different reperfusion time seems to be different. This may help reveal the mechanisms of the onset and progression of RIR injury and offer a novel aspect of its treatment.

In summary, we screened the differentially expressed circRNAs, lncRNAs, and mRNAs between the control and model groups and at different reperfusion time (24h, 72h, and 7d). 5 key genes, Cd74, RT1-Da, RT1-CE5, RT1-Bb, RT1-DOa, were selected. Key genes specific to reperfusion time were selected to show the change in biological process with the increased reperfusion time. These results provided theoretical support and a reference basis for the clinical treatment.