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. Supplementary Information The online version contains supplementary material available at 10.1186/s12864-024-10288-0. Keywords: Retinal ischemia/reperfusion (RIR), ncRNAs, Whole transcriptome sequencing Introduction As one of the highest oxygen-consuming tissues in the body, the retina has a delicate and complex structure and a vigorous metabolism [[35]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 [[36]2, [37]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) [[38]4–[39]6], whereas the inner five layers are supplied by the central retinal artery (referred to as the retinal circulation) [[40]7, [41]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 [[42]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 [[43]10, [44]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 [[45]12–[46]14]. It is the main cause and pathogenesis of retinal thinning and atrophy, retinal ganglion cell (RGC) death, and visual function impairment [[47]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) [[48]16, [49]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 [[50]18, [51]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 [[52]20–[53]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 [[54]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. [55]1. Fig. 1. [56]Fig. 1 [57]Open in a new tab 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 [[58]23, [59]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 [[60]25] software for adapter removal and filtering. Clean reads were quality inspected by fastqc [[61]26] and then aligned to the human reference genome using hisat2 [[62]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 [[63]28] to align the sequencing reads of each sample with the reference genome. CIRI2 was used to scan the circRNA-paired chiastic clipping signal [[64]29, [65]30]. The expression values were normalized by the reads per million (RPM) algorithm [[66]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) [[67]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 [[68]33]. In this study, GO and KEGG pathway enrichment analysis of DEmRNAs was performed by the clusterProfiler package [[69]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 ([70]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 [[71]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, [72]http://mirwalk.umm.uni-heidelberg.de/). And the binding relationship of candidate lncRNAs and candidate circRNAs with miRNA was predicted by miRanda ([73]http://www.microrna.org/ microrna/home.do). To find other genes associated with the key genes, GeneMANIA ([74]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 ([75]https://www.ncbi.nlm.nih.gov/). The base sequences of the genes were then entered into the mRNALocater ([76]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, [77]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 ([78]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 ([79]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, [80]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 [81]1. Table 1. Primer sequences for qPCR assay Gene Sequence Product length GAPDH Forward 5’-CTCCTCGAAGTACCCTGTGC-3’ 353 bp Reserve 5’-CATGGTGCAGCGATGCTTTA-3’ Cd74 Forward 5’-GTGACACTGGGCTACTCGTC-3’ 92 bp Reserve 5’-GAGCAGTCAAGCCCTCCATT-3’ RT1-Da Forward 5’-GTCAGCCCGGGGTACAGC-3’ 943 bp Reserve 5’-AGCCAGCGATAGTCTGTGGT-3’ RT1-Bb Forward 5’-AGCCAGCGATAGTCTGTGGT-3’ 943 bp Reserve 5’-TGTACCCCGGGCTGACGAT-3’ RT1-CE5 Forward 5’-TACCATGCTGGAGTTGGTGG-3’ 359 bp Reserve 5’-AAATCTTCCACACAGATCCCC-3’ RT1-DOa Forward 5’-TCCCCGTCAAATTCGTGTGT-3’ 96 bp Reserve 5’-CATCCCTCACCCGATACAGC-3’ [82]Open in a new tab s Statistical analysis GraphPad Prism 8.0 was used to analyze the data and graph. Data were presented as mean ± standard deviation ( [MATH: x¯±SD :MATH] ) 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 [83]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. [84]2A-F, Supplementary Tables [85]1, [86]2 and [87]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. [88]2G-L, Supplementary Tables [89]4, [90]5 and [91]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. [92]2M-R, Supplementary Tables [93]7, [94]8 and [95]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. [96]2S-T, Supplementary Tables [97]10, [98]11 and [99]12). Fig. 2. [100]Fig. 2 [101]Open in a new tab 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. [102]3A, Supplementary Tables [103]13 and [104]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. [105]3B, Supplementary Tables [106]15 and [107]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. [108]3C, Supplementary Tables [109]17 and [110]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. [111]4A-C, Supplementary Tables [112]20 and [113]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. [114]4D-F, Table [115]2, Supplementary Tables [116]22 and [117]23). Fig. 3. [118]Fig. 3 [119]Open in a new tab 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. [120]Fig. 4 [121]Open in a new tab 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 ID Description Counts P Value qvalue Category GO:0050867 positive regulation of cell activation 45 2.39264E-23 2.14E-47 GO_BP GO:0002250 adaptive immune response 45 3.8263E-23 1.62E-50 GO_BP GO:0002696 positive regulation of leukocyte activation 43 9.36341E-23 2.56E-44 GO_BP GO:0002443 leukocyte mediated immunity 42 7.99351E-20 5.00E-39 GO_BP GO:0002683 negative regulation of immune system process 41 1.45064E-17 5.98E-48 GO_BP GO:1903131 mononuclear cell differentiation 41 2.26455E-17 1.32E-35 GO_BP GO:0002697 regulation of immune effector process 40 7.77577E-18 4.62E-41 GO_BP GO:0045785 positive regulation of cell adhesion 40 2.91266E-16 3.47E-45 GO_BP GO:0032103 positive regulation of response to external stimulus 39 3.31556E-16 8.34E-33 GO_BP GO:0002253 activation of immune response 38 1.06302E-21 7.46E-31 GO_BP GO:0009897 external side of plasma membrane 49 1.76287E-27 4.45358E-25 GO_CC GO:0031012 extracellular matrix 25 2.97486E-09 1.38221E-07 GO_CC GO:0030312 external encapsulating structure 25 3.28275E-09 1.38221E-07 GO_CC GO:0000323 lytic vacuole 20 5.42278E-05 0.000856229 GO_CC GO:0005764 lysosome 20 5.42278E-05 0.000856229 GO_CC GO:0043235 receptor complex 18 5.16901E-05 0.000856229 GO_CC GO:0062023 collagen-containing extracellular matrix 17 2.28492E-07 6.41381E-06 GO_CC GO:0001772 immunological synapse 12 1.7685E-12 2.2339E-10 GO_CC GO:0005581 collagen trimer 12 8.66691E-11 5.47384E-09 GO_CC GO:0005770 late endosome 12 0.00057603 0.006929687 GO_CC GO:0140375 immune receptor activity 18 2.17255E-12 8.96465E-10 GO_MF GO:0019955 cytokine binding 17 1.18871E-10 1.635E-08 GO_MF GO:0042277 peptide binding 16 5.03693E-05 0.000989713 GO_MF GO:0033218 amide binding 16 0.000544452 0.008640692 GO_MF GO:0030246 carbohydrate binding 15 1.41925E-05 0.000344487 GO_MF GO:0004896 cytokine receptor activity 14 2.47717E-10 2.55539E-08 GO_MF GO:0003823 antigen binding 13 2.461E-11 5.07742E-09 GO_MF GO:0005539 glycosaminoglycan binding 12 0.000172029 0.003086288 GO_MF GO:0008528 G protein-coupled peptide receptor activity 11 1.90402E-05 0.000436477 GO_MF GO:0001653 peptide receptor activity 11 3.32433E-05 0.000685862 GO_MF rno04514 Cell adhesion molecules 20 1.45E-11 1.20E-09 KEGG rno04060 Cytokine-cytokine receptor interaction 20 1.90E-08 2.53E-07 KEGG rno05169 Epstein-Barr virus infection 16 9.66E-07 8.92E-06 KEGG rno05166 Human T-cell leukemia virus 1 infection 16 3.45E-06 2.74E-05 KEGG rno05171 Coronavirus disease - COVID-19 16 7.21E-05 0.00042816 KEGG [122]Open in a new tab 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. [123]5, Supplementary Tables [124]24 and [125]25). Fig. 5. [126]Fig. 5 [127]Open in a new tab 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. [128]6a, Supplementary Table [129]34). RT1-Bb was found to be positively related to ribosome, antigen processing, and presentation (Fig. [130]6b, Supplementary Table [131]26). RT1-CE5 was found to be positively related to ribosome (Fig. [132]6C, Supplementary Table [133]27). RT1-Da was found to be positively related to ribosome, antigen processing, and presentation, regulation of actin cytoskeleton (Fig. [134]6D, Supplementary Table [135]28). RT1-DOa was found to be negatively related to ribosome (Fig. [136]6E, Supplementary Table [137]29). Fig. 6. [138]Fig. 6 [139]Open in a new tab 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. [140]7A, Supplementary Table [141]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. [142]7B, Supplementary Table [143]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. [144]7C, Supplementary Tables [145]33 and [146]34). Fig. 7. [147]Fig. 7 [148]Open in a new tab 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. [149]7D, Supplementary Table [150]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. [151]7E, Supplementary Table [152]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. [153]7F, Supplementary Table [154]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. [155]7G, Supplementary Table [156]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. [157]7H, Supplementary Table [158]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. [159]8 A-J). Fig. 8. [160]Fig. 8 [161]Open in a new tab 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. [162]9A, Supplementary Table [163]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. [164]9B, Supplementary Table [165]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. [166]9C, Supplementary Table [167]42). Fig. 9. [168]Fig. 9 [169]Open in a new tab 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. [170]10A-C, Table [171]3, Supplementary Tables [172]43 and [173]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. [174]10D-F, Table [175]4, Supplementary Tables [176]45 and [177]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. [178]10G-I, Table [179]5, Supplementary Tables [180]47 and [181]48). Fig. 10. [182]Fig. 10 [183]Open in a new tab 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) ID Description Counts P Value qvalue Category GO:0001654 eye development 58 3.17E-07 0.000141998 GO_BP GO:0150063 visual system development 58 3.93E-07 0.000161268 GO_BP GO:0048880 sensory system development 58 5.58E-07 0.000211499 GO_BP GO:0043010 camera-type eye development 56 4.01E-08 3.13E-05 GO_BP GO:0007015 actin filament organization 55 1.62E-06 0.00049719 GO_BP GO:0010959 regulation of metal ion transport 55 1.09E-05 0.001672816 GO_BP GO:0042391 regulation of membrane potential 55 2.06E-05 0.002473344 GO_BP GO:0062197 cellular response to chemical stress 54 2.52E-07 0.00012416 GO_BP GO:0001654 eye development 58 3.17E-07 0.000141998 GO_BP GO:0150063 visual system development 58 3.93E-07 0.000161268 GO_BP GO:0005840 ribosome 62 5.19E-18 8.69E-16 GO_CC GO:0098984 neuron to neuron synapse 58 3.73E-07 1.56E-05 GO_CC GO:0044391 ribosomal subunit 56 2.41E-19 6.05E-17 GO_CC GO:0014069 postsynaptic density 52 3.73E-06 0.000103985 GO_CC GO:0032279 asymmetric synapse 52 5.95E-06 0.000142259 GO_CC GO:0031252 cell leading edge 52 8.78E-06 0.000191603 GO_CC GO:0099572 postsynaptic specialization 52 3.34E-05 0.000492938 GO_CC GO:0031012 extracellular matrix 51 2.18E-06 6.89E-05 GO_CC GO:0030312 external encapsulating structure 51 2.52E-06 7.44E-05 GO_CC GO:0031253 cell projection membrane 50 4.36E-08 3.13E-06 GO_CC GO:0003735 structural constituent of ribosome 52 5.40E-17 4.64E-14 GO_MF GO:0030546 signaling receptor activator activity 52 3.10E-05 0.004616756 GO_MF GO:0048018 receptor ligand activity 51 3.31E-05 0.004616756 GO_MF GO:0015267 channel activity 49 0.000285193 0.015291589 GO_MF GO:0022803 passive transmembrane transporter activity 49 0.000285193 0.015291589 GO_MF GO:0003779 actin binding 45 0.000144 0.010574374 GO_MF GO:0046873 metal ion transmembrane transporter activity 45 0.000210388 0.012892223 GO_MF GO:0005216 ion channel activity 44 0.00059622 0.021312244 GO_MF GO:0022836 gated channel activity 35 0.001195484 0.033083858 GO_MF GO:0003735 structural constituent of ribosome 52 5.40E-17 4.64E-14 GO_MF rno05171 Coronavirus disease - COVID-19 60 8.55E-15 1.15E-12 KEGG rno03010 Ribosome 51 5.24E-16 1.41E-13 KEGG rno04010 MAPK signaling pathway 35 0.000346598 0.012770921 KEGG rno04060 Cytokine-cytokine receptor interaction 33 0.000217877 0.011515421 KEGG rno04510 Focal adhesion 29 2.51E-05 0.001682982 KEGG [184]Open in a new tab Table 4. Top 10 GO enrichment and Top 5 KEGG enrichment of characteristic differentially expressed mRNA in 72h reperfusion (spec_72h mRNAs) ID Description Counts P Value qvalue Category GO:0042391 regulation of membrane potential 36 1.10E-08 3.62E-05 GO_BP GO:0007409 axonogenesis 30 8.09E-06 0.00231256 GO_BP GO:0034329 cell junction assembly 29 4.63E-06 0.00167271 GO_BP GO:0031099 regeneration 28 1.88E-08 3.62E-05 GO_BP GO:0006816 calcium ion transport 27 5.63E-05 0.005831814 GO_BP GO:0006875 cellular metal ion homeostasis 25 0.000117623 0.008538448 GO_BP GO:0070588 calcium ion transmembrane transport 24 1.48E-05 0.003109349 GO_BP GO:0007611 learning or memory 23 2.05E-05 0.003380869 GO_BP GO:0021537 telencephalon development 23 4.92E-05 0.005433553 GO_BP GO:0072507 divalent inorganic cation homeostasis 23 8.34E-05 0.006976169 GO_BP GO:0097060 synaptic membrane 44 3.56E-14 3.56E-14 GO_CC GO:0099240 intrinsic component of synaptic membrane 34 9.58E-16 9.58E-16 GO_CC GO:0031012 extracellular matrix 34 1.05E-09 1.05E-09 GO_CC GO:0030312 external encapsulating structure 34 1.19E-09 1.19E-09 GO_CC GO:0099699 integral component of synaptic membrane 33 5.27E-16 5.27E-16 GO_CC GO:0045211 postsynaptic membrane 33 6.12E-12 6.12E-12 GO_CC GO:0043235 receptor complex 32 1.61E-08 1.61E-08 GO_CC GO:1902495 transmembrane transporter complex 30 2.10E-08 2.10E-08 GO_CC GO:1990351 transporter complex 30 8.39E-08 8.39E-08 GO_CC GO:0099572 postsynaptic specialization 30 2.25E-06 2.25E-06 GO_CC GO:0005216 ion channel activity 37 6.51E-11 1.83E-08 GO_MF GO:0015267 channel activity 37 1.27E-09 7.95E-08 GO_MF GO:0022803 passive transmembrane transporter activity 37 1.27E-09 7.95E-08 GO_MF GO:0022836 gated channel activity 33 1.54E-11 8.65E-09 GO_MF GO:0005261 cation channel activity 31 4.54E-10 4.20E-08 GO_MF GO:0046873 metal ion transmembrane transporter activity 27 8.98E-06 0.000229803 GO_MF GO:0005543 phospholipid binding 24 0.0007170 0.010205534 GO_MF GO:0015276 ligand-gated ion channel activity 20 1.60E-10 3.01E-08 GO_MF GO:0022834 ligand-gated channel activity 20 2.40E-10 3.29E-08 GO_MF GO:0005244 voltage-gated ion channel activity 20 1.10E-07 3.65E-06 GO_MF rno04080 Neuroactive ligand-receptor interaction 32 6.58E-10 1.51E-07 KEGG rno04723 Retrograde endocannabinoid signaling 17 9.88E-08 4.54E-06 KEGG rno04024 cAMP signaling pathway 17 1.87E-05 0.00039076 KEGG rno04020 Calcium signaling pathway 17 9.93E-05 0.00189851 KEGG rno04726 Serotonergic synapse 16 3.06E-08 2.34E-06 KEGG [185]Open in a new tab Table 5. Top 10 GO enrichment and Top 5 KEGG enrichment of characteristic differentially expressed mRNA in 7d reperfusion (spec_7d mRNAs) ID Description Counts P Value qvalue Category GO:0002250 adaptive immune response 14 2.07E-09 0.000156495 GO_BP GO:0030098 lymphocyte differentiation 12 2.72E-07 2.95E-05 GO_BP GO:1903131 mononuclear cell differentiation 12 1.05E-06 0.000161552 GO_BP GO:0042113 B cell activation 11 5.82E-08 0.000176511 GO_BP GO:0002449 lymphocyte mediated immunity 10 1.22E-06 0.000490106 GO_BP GO:0002460 adaptive immune response based on somatic recombination of immune receptors built from immunoglobulin superfamily domains 10 1.48E-06 0.001039951 GO_BP GO:0070661 leukocyte proliferation 10 5.36E-06 0.001713198 GO_BP GO:0002443 leukocyte mediated immunity 10 1.54E-05 0.001039951 GO_BP GO:0002683 negative regulation of immune system process 10 3.74E-05 0.001092719 GO_BP GO:0046651 lymphocyte proliferation 9 1.60E-05 0.000156495 GO_BP GO:0009897 external side of plasma membrane 11 6.57E-06 0.000209643 GO_CC GO:0031012 extracellular matrix 9 8.32E-05 0.000918831 GO_CC GO:0030312 external encapsulating structure 9 8.63E-05 0.000918831 GO_CC GO:0043235 receptor complex 9 8.63E-05 0.000918831 GO_CC GO:0098802 plasma membrane signaling receptor complex 7 1.66E-05 0.000397478 GO_CC GO:0062023 collagen-containing extracellular matrix 6 0.000843948 0.008084135 GO_CC GO:0098636 protein complex involved in cell adhesion 5 1.78E-06 8.53E-05 GO_CC GO:0042101 T cell receptor complex 4 1.71E-07 1.64E-05 GO_CC GO:0042613 MHC class II protein complex 3 2.31E-05 0.000443124 GO_CC GO:0042611 MHC protein complex 3 6.49E-05 0.000918831 GO_CC GO:0005201 extracellular matrix structural constituent 5 1.64E-05 0.002732036 GO_MF GO:0004896 cytokine receptor activity 4 0.000827139 0.03874492 GO_MF GO:0023026 MHC class II protein complex binding 3 2.92E-05 0.002732036 GO_MF GO:0023023 MHC protein complex binding 3 0.00012028 0.00751221 GO_MF rno04658 Th1 and Th2 cell differentiation 8 3.10E-09 1.99E-07 KEGG rno04659 Th17 cell differentiation 8 1.08E-08 3.46E-07 KEGG rno05169 Epstein-Barr virus infection 7 4.34E-05 0.000397816 KEGG rno04672 Intestinal immune network for IgA production 6 2.34E-08 5.02E-07 KEGG rno04640 Hematopoietic cell lineage 6 2.06E-06 2.65E-05 KEGG [186]Open in a new tab 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. [187]11A, Supplementary Table [188]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. [189]11B, Supplementary Table [190]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. [191]11C, Supplementary Table [192]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. [193]11D, Supplementary Table [194]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 [195]2A, Supplementary Tables [196]53, [197]54, [198]55, [199]56, [200]57, [201]58, [202]59, [203]60, [204]61 and [205]62). Fig. 11. [206]Fig. 11 [207]Open in a new tab 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. [208]12A, Supplementary Table [209]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. [210]12B, Supplementary Table [211]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. [212]12C, Supplementary Table [213]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. [214]12D, Supplementary Table [215]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 [216]2B, Supplementary Tables [217]67, [218]68, [219]69, [220]70, [221]71, [222]72, [223]73, [224]74, [225]75, [226]76, [227]77 and [228]78). Fig. 12. [229]Fig. 12 [230]Open in a new tab 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. [231]13A, Supplementary Table [232]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. [233]13B, Supplementary Table [234]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. [235]13C, Supplementary Table [236]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. [237]12D, Supplementary Table [238]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 [239]2C, Supplementary Tables [240]81, [241]82, [242]83, [243]84, [244]85, [245]86, [246]87, [247]88, [248]89 and [249]90). Fig. 13. [250]Fig. 13 [251]Open in a new tab 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 [[252]15, [253]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) [[254]11, [255]15, [256]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 [[257]11, [258]38–[259]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 [[260]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 [[261]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 [[262]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 [[263]44, [264]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 [[265]11]. However, the immunologic response to various stress cues has been found to play a pivotal role in the retina [[266]13]. When ischemia and reperfusion occur, the permeability of the blood-retinal barrier is changed, and the activation of microgila [[267]46, [268]47] and the increase of Treg cells [[269]48, [270]49] are observed. Toll signaling activation was also found and induced inflammasome formation [[271]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) [[272]51]. CD74 regulates T-cell and B-cell development, dendritic cell (DC) motility, macrophage inflammation, and thymic selection [[273]52]. Cd74, RT1-Ba, RT1-Bb, RT1-Da, and RT1-Db1 were referred to as major histocompatibility complex (MHC) class II members [[274]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 [[275]54]. Abcouwer et al. [[276]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 [[277]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. [[278]57] Puerarin can ameliorate RIR injury by suppressing apoptosis and TLR4/NLRP3 inflammasome activation in RGCs. Lee et al. [[279]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. Supplementary Information [280]Supplementary Material 1.^ (2.4MB, xlsx) [281]Supplementary Material 2. ^ (1.1MB, xlsx) [282]Supplementary Material 3. ^ (888.4KB, xlsx) [283]Supplementary Material 4. ^ (2.3MB, xlsx) [284]Supplementary Material 5. ^ (954.8KB, xlsx) [285]Supplementary Material 6. ^ (771.2KB, xlsx) [286]Supplementary Material 7. ^ (2.2MB, xlsx) [287]Supplementary Material 8. ^ (980.2KB, xlsx) [288]Supplementary Material 9. ^ (794.8KB, xlsx) [289]Supplementary Material 10. ^ (2.4MB, xlsx) [290]Supplementary Material 11. ^ (1.1MB, xlsx) [291]Supplementary Material 12. ^ (1.2MB, xlsx) [292]Supplementary Material 13. ^ (863KB, xlsx) [293]Supplementary Material 14. ^ (433.2KB, xlsx) [294]Supplementary Material 15. ^ (792.8KB, xlsx) [295]Supplementary Material 16. ^ (315.6KB, xlsx) [296]Supplementary Material 17.^ (512.1KB, xlsx) [297]Supplementary Material 18. ^ (256.7KB, xlsx) [298]Supplementary Material 19. ^ (21.5KB, xlsx) [299]Supplementary Material 20. ^ (17.9KB, xlsx) [300]Supplementary Material 21. ^ (16.3KB, xlsx) [301]Supplementary Material 22. ^ (8.2KB, csv) [302]Supplementary Material 23. ^ (80.4KB, xlsx) [303]Supplementary Material 24. ^ (26.4KB, tsv) [304]Supplementary Material 25. ^ (126B, csv) [305]Supplementary Material 26. ^ (1.4KB, xls) [306]Supplementary Material 27. ^ (922B, xls) [307]Supplementary Material 28. ^ (636B, xls) [308]Supplementary Material 29. ^ (1.4KB, xls) [309]Supplementary Material 30. ^ (631B, xls) [310]Supplementary Material 31. ^ (39.1KB, xlsx) [311]Supplementary Material 32. ^ (23.1KB, xlsx) [312]Supplementary Material 33. ^ (100.7KB, xlsx) [313]Supplementary Material 34. ^ (206.9KB, pdf) [314]Supplementary Material 35. ^ (440B, csv) [315]Supplementary Material 36. ^ (20.7KB, xlsx) [316]Supplementary Material 37. ^ (20.4KB, xlsx) [317]Supplementary Material 38. ^ (19KB, xlsx) [318]Supplementary Material 39. ^ (30.4KB, xlsx) [319]Supplementary Material 40. ^ (240.4KB, xlsx) [320]Supplementary Material 41. ^ (100.4KB, xlsx) [321]Supplementary Material 42. ^ (27.2KB, xlsx) [322]Supplementary Material 43. ^ (3.2KB, tsv) [323]Supplementary Material 44. ^ (90.5KB, tsv) [324]Supplementary Material 45. ^ (4KB, tsv) [325]Supplementary Material 46. ^ (47.4KB, tsv) [326]Supplementary Material 47. ^ (3KB, tsv) [327]Supplementary Material 48. ^ (12.4KB, tsv) [328]Supplementary Material 49. ^ (522.8KB, xlsx) [329]Supplementary Material 50. ^ (16KB, xlsx) [330]Supplementary Material 51. ^ (17.9KB, xlsx) [331]Supplementary Material 52. ^ (18.1KB, xlsx) [332]Supplementary Material 53. ^ (22.1KB, xls) [333]Supplementary Material 54. ^ (19.2KB, xls) [334]Supplementary Material 55. ^ (24KB, xls) [335]Supplementary Material 56. ^ (21.8KB, xls) [336]Supplementary Material 57. ^ (16.2KB, xls) [337]Supplementary Material 58. ^ (26.8KB, xls) [338]Supplementary Material 59. ^ (21.3KB, xls) [339]Supplementary Material 60. ^ (16.1KB, xls) [340]Supplementary Material 61. ^ (18.1KB, xls) [341]Supplementary Material 62. ^ (16.9KB, xls) [342]Supplementary Material 63. ^ (87.1KB, xlsx) [343]Supplementary Material 64. ^ (16KB, xlsx) [344]Supplementary Material 65. ^ (76.4KB, xlsx) [345]Supplementary Material 66. ^ (60.4KB, xlsx) [346]Supplementary Material 67. ^ (18.2KB, xls) [347]Supplementary Material 68. ^ (18.2KB, xls) [348]Supplementary Material 69. ^ (18.6KB, xls) [349]Supplementary Material 70. ^ (18.2KB, xls) [350]Supplementary Material 71. ^ (17KB, xls) [351]Supplementary Material 72. ^ (18.2KB, xls) [352]Supplementary Material 73. ^ (18.4KB, xls) [353]Supplementary Material 74. ^ (18.2KB, xls) [354]Supplementary Material 75. ^ (18.2KB, xls) [355]Supplementary Material 76. ^ (18KB, xls) [356]Supplementary Material 77. ^ (21.7KB, xlsx) [357]Supplementary Material 78. ^ (16KB, xlsx) [358]Supplementary Material 79. ^ (86.7KB, xlsx) [359]Supplementary Material 80. ^ (64.5KB, xlsx) [360]Supplementary Material 81. ^ (20.3KB, xls) [361]Supplementary Material 82. ^ (20KB, xls) [362]Supplementary Material 83. ^ (18.9KB, xls) [363]Supplementary Material 84. ^ (19KB, xls) [364]Supplementary Material 85. ^ (20KB, xls) [365]Supplementary Material 86. ^ (19.4KB, xls) [366]Supplementary Material 87. ^ (18.2KB, xls) [367]Supplementary Material 88. ^ (19.8KB, xls) [368]Supplementary Material 89. ^ (19KB, xls) [369]Supplementary Material 90. ^ (12.9KB, xls) [370]Supplementary Material 91. ^ (2.8MB, tif) [371]Supplementary Material 92. ^ (6.5MB, tif) [372]Supplementary Material 93. ^ (17.3MB, zip) Acknowledgments