ScRNA-seq Data Reveal Gene Upregulation and Downregulation in Oxygen-Induced Retinopathy

Introduction

Retinopathy of prematurity (ROP) is a vascular disease in premature infants characterized by pathological retinal neovessels, retinal vascular dilation, and hemorrhage, which could lead to retinal detachment and blindness [1]. Its uniqueness lies in the concurrent physiological and pathological angiogenesis in the developing retina, providing a unique model for studying angiogenesis regulation [2–4]. Pathological angiogenesis is a core feature of ROP, and understanding its molecular mechanisms is vital for developing treatment strategies for various diseases [2–4]. The most widely used oxygen-induced retinopathy (OIR) model has enabled the study of the pathophysiological processes occurring in ROP, uncovering both the mechanisms behind ROP and potential therapeutic interventions [5].

Current ROP screening methods require substantial resources in terms of equipment and ophthalmologist expertise, while biomarkers with high sensitivity and specificity offer a more convenient approach for early detection [6]. The interaction of different cell types and the regulation of various genes contribute to the development of pathological angiogenesis of ROP [4,6]. Therefore, exploring the gene regulations in ROP can aid in more precisely targeted therapies.

Single-cell RNA sequencing (scRNA-seq) has revolutionized the study of gene expression at a cellular resolution [7]. The advantage of scRNA-seq lies in its ability to discover rare cell types and track cell differentiation pathways, thereby enhancing the understanding of complex biological systems and disease mechanisms [7,8]. Some studies have used scRNA-seq to characterize the signal pathways or gene regulations associated with retinopathy and its relative therapies. This includes investigating how these pathways affect retinal neovascularization and promote retinal neurodegeneration [9,10].

However, research on gene regulation in OIR using scRNA-seq remains limited. This study investigated gene regulation in OIR rat models using data obtained from the Gene Expression Omnibus (GEO) platform. The key genes identified may serve as valuable risk biomarkers to enhance ROP screening.

Material and Methods

DATA DOWNLOAD:

The data set GSE150703 and annotation file GPL19057 were downloaded from the GEO data platform. GSE150703 contains 12 samples of retinal tissues obtained from mice under normoxic (NORM) conditions and OIR conditions on postnatal days 14 and 17. In this bioinformatics analysis, the data of the animals retrieved from the database were divided into 4 groups: P14-NORM (n=2), the NORM control group on postnatal day 14; P14-OIR (n=5), the OIR group on postnatal day 14; P17-NORM (n=2), the NORM control group on postnatal day 17; and P17-OIR (n=3), the OIR group on postnatal day 17.

This research is a bioinformatics analysis and ethical approval certification is not required. The relevant data for the research was obtained from the database and analyzed.

QUALITY CONTROL:

The R package DropletUtils [11] was used to determine the expression levels in each cell, and barcodes were filtered. The cells were further filtered according to the number of unique molecular identifiers in each cell. Then, the scater [12] package was used to quantify the expression of genes in cells, and the cells were filtered based on the proportion of mitochondrial gene reads (<10%), as shown in Figure 1A–1C. We retained those cells whose proportion of mitochondrial gene-encoded unique molecular identifiers (percent.mt) was less than 20% and whose detected number of genes (nFeature_RNA) was greater than 200 and less than 6000, as well as those with a total unique molecular identifier count (nCount_RNA) greater than 500 for further analysis.

DATA PREPROCESSING:

The NormalizeData function in the Seurat [13] package was used to normalize the expression matrix of each sample after filtering.

PRINCIPAL COMPONENT ANALYSIS:

The 2000 genes with the greatest differential expression between cells were identified by using the FindVariableFeatures function in the Seurat package. Then, the ScaleData function in the Seurat package was used to linearly scale the expression data. Finally, the RunPCA function in the Seurat package was used for linear dimensionality reduction analysis, as shown in Figure 1D–1G.

CELL CLUSTERING AND ANNOTATION:

We selected the principal component with the highest standard deviation and used the FindNeighbors and FindClusters functions in the Seurat package for cell cluster analysis. Then, the RunUMAP function in the Seurat package was used for nonlinear dimensionality reduction analysis, as shown in Figure 1H and 1I.

IDENTIFICATION OF MARKER GENES:

The FindMarkers function in the Seurat package was used to identify the differentially expressed genes (DEGs) between the cells in each cluster and all other cells (log2 fold change [FC] ≥2; P value ≤0.05) to identify marker genes (the top 500 ranked by logFC). The cells were labelled according to the existing marker [14] genes, and a cluster diagram was generated, as shown in Figure 1J and 1K. We used the FindMarkers function in the Seurat package for differential expression analysis, and the minimum expression ratio in the cell population was set at 0.25.

GENE ENRICHMENT ANALYSES:

Based on the GO [15] database and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway database [16], functional enrichment analyses were performed with the candidate genes. GO enrichment analysis was performed to discover the biological processes enriched with these DEGs and demonstrated the important functions of the differential gene expression among the groups. GO enrichment analysis consists of 3 categories: biological processes, cellular components, and molecular functions.

PROTEIN–PROTEIN INTERACTION NETWORK PREDICTION:

The protein–protein interactions (PPIs) of the proteins encoded by the candidate genes were analyzed using the STRING [17] online tool. In this study, a required confidence level (combined score) of >0.7 was selected as the PPI threshold. By connectivity degree analysis of the network statistics, the hub proteins, important nodes involved in the protein interactions in the PPI network, were identified.

Results

IDENTIFICATION OF CELL TYPES:

The cell type in each cluster was determined based on the marker gene for each cell and the marker gene for each cell cluster. We identified 10 kinds of cells: bipolar cells, neurons, Muller glia, amacrine cells, cones, perivascular cells, microglia, retinal ganglion cells, vascular endothelial cells, and horizontal cells, as shown in Figure 2A and 2B. Then, the marker genes of all these cell types were reanalyzed to identify new marker genes, as shown in Figure 2C.

ANALYSIS OF CELL PROPORTIONS:

The proportions of the types of cells in each sample with respect to the corresponding total number of cells were determined, as shown in Figure 2D. Then, the t test was used to determine the significance of the differences in the corresponding comparisons, as shown in Figure 2E. Figure 2D and 2E showed that there were significant differences among 3 kinds of cells: bipolar cells (P<0.001), Muller glia (P=0.027), and retinal ganglion cells (P=0.043). According to the results of the above analysis, we selected the most significantly different cells, bipolar cells, as the core cells for subsequent analyses. Bipolar cells act as relay stations for retinal signal transmission and sensitive nodes for energy metabolism for OIR. Bipolar cells have mitochondrial dysfunction, imbalance in lactate metabolism, and changes in synaptic plasticity. These characteristics make them both the early targets for hypoxia-induced damage and the key links in amplifying retinal dysfunction and structural degeneration.

DIFFERENTIAL EXPRESSION ANALYSIS:

Differential expression analysis of bipolar cells with the FindMarkers function in the Seurat package revealed 355 upregulated genes and 5 downregulated genes, as shown in Tables 1 and 2. These genes are also illustrated in a scatter plot in Figure 3A.

GO ENRICHMENT ANALYSIS OF DEGS:

GO enrichment analysis showed that the top 9 biological process terms were RNA splicing, RNA processing, messenger (m)RNA processing, gene expression, RNA metabolic process, mRNA metabolic process, organelle organization, regulation of RNA metabolic process, and regulation of organelle organization, as shown in Figure 3B–3D. The top 9 cellular component terms were nucleus, nucleoplasm, nuclear body, nuclear lumen, nuclear speck, neuron projection, nuclear periphery, spliceosomal complex, and intracellular membrane-bounded organelle, as shown in Figure 3E–3G. The top 9 molecular function terms were RNA binding, DNA binding, mRNA binding, kinase binding, mRNA 3′ untranslated region binding, histone deacetylase binding, single-stranded DNA binding, lysine-acetylated histone binding, and nucleoside-triphosphatase activity, as shown in Figure 3H–3J.

In GO analysis, the upregulated genes with high expression included Ywhab, heterogeneous nuclear ribonucleoprotein U (Hnrnpu), casein kinase 1 alpha 1 (Csnk1a1), Y-box-binding protein 1 (Ybx1), vesicle-associated membrane protein 2 (Vamp2), purb, and Chchd10, and the downregulated genes included the long noncoding RNA (lncRNA) X-inactive specific transcript (lncRNA Xist) and mitochondrial DNA (mtDNA)-encoded cytochrome c oxidase subunit I (mt-Co1).

KEGG PATHWAY ENRICHMENT ANALYSIS OF DEGS:

To characterize the coordinative relations between genes and the roles of the genes in biological functions, the DEGs were subjected to KEGG pathway enrichment analysis, in which biochemical metabolic and signal transduction pathways were found to be enriched. As shown in Figure 4A, the top 15 pathways are listed in the bubble plot: spliceosome, RNA transport, phototransduction, mineral absorption, endocrine and other factor-regulated calcium reabsorption, protein processing in endoplasmic reticulum, lysine degradation, thyroid hormone signalling pathway, thyroid hormone synthesis, salivary secretion, Alzheimer disease, HIF-1 signaling pathway, glutamatergic synapse, proximal tubule bicarbonate reclamation, and central carbon metabolism in cancer.

In KEGG analysis, the upregulated genes with high expression included Hnrnpu, phosphofructokinase 1 (Pfkp), Rho, Tpr, Mt1, Vamp2, and Csnk1a1; and the downregulated gene was nicotinamide adenine dinucleotide ubiquinone oxidoreductase subunit 5 (Ndufs5), as shown in Figure 4B–4D. The analysis revealed significantly enriched KEGG pathways associated with RNA transport, phototransduction, thyroid hormone signaling, salivary secretion, mineral absorption, and the spliceosome, as presented in Figures 5A–E and 6.

PPI NETWORK ANALYSIS:

Based on the above DEGs, PPI network analysis was performed using the STRING online tool. Then, the hub genes were selected based on their connectivity degree, as shown in Table 3 and Figure 6. We identified the genes with a connectivity degree of 25 or greater as hub genes, which included serine/arginine splicing factor 1 (Srsf1), Srsf11, Sf3b1, Hnrnpu, hnrnph1, and heterogeneous ribonucleoprotein-C (Hnrnpc).

Discussion

LIMITATIONS:

This study identified hub genes related to OIR but has limitations. The study relied on bioinformatics without experimental validation, and the small sample size may limit the generalization of the results. The study also overlooked epigenetic and protein modification interactions and did not fully consider genetic, environmental, and lifestyle factors. Future research should validate findings experimentally and integrate multi-omics data for a more comprehensive understanding of the pathogenesis of ROP.

Conclusions

This study identifies key hub genes in OIR through enrichment and PPI network analyses. These genes, which have been implicated in cancer and immune diseases, are closely associated with the pathological processes of OIR in the present study. These genes may be potential biomarkers for ROP; however, further in vivo experiments or clinical trials are required.