Involvement of CircRNA Expression Profile in Diabetic Retinopathy and Its Potential Diagnostic Value.

Background: Circular RNAs (circRNAs), a class of non-coding and undegradable RNAs, play many pathological functions by acting as miRNA sponges, interacting with RNA-binding proteins, and others. The recent literature indicates that circRNAs possess the advanced superiority for the early screening of diabetic retinopathy (DR).
Methods: CircRNA sources of peripheral blood mononuclear cells (PBMCs) from healthy controls (n = 4), diabetes mellitus patients (DM) (n = 4), and DR patients (n = 4) were extracted for circular RNA microarray analysis. Enriched biological modules and signaling pathways were analyzed by Gene Ontology Enrichment and Kyoto Encyclopedia of Genes and Genomes analysis, respectively. Real-time quantitative reverse transcription PCR (RT-qPCR) was performed to validate differentiated levels of several circRNAs (fold change ≥2, p < .05) in different groups of healthy control subjects (n = 20), DM patients (n = 60), and DR patients (n = 42). Based on our clinical data from DR, the diagnostic performance of candidate circRNAs was measured by operating characteristic curves (ROCs). Subsequently, their circRNA-miRNA networks were constructed by bioinformatics analysis.
Results: Circular RNA microarray analysis was performed, and 2,452 and 289 circRNAs were screened with differential expression in DR patients compared to healthy controls and DM patients, respectively. Enrichment analyses showed that circRNAs in DR patients were enriched in extracellular matrix (ECM)-receptor interaction and focal adhesion pathways. The top 5 differential circRNAs in circRNA microarray analysis were subsequently quantified and verified by RT-qPCR. Consistently, a significant 2.2-fold reduction of hsa_circ_0095008 and 1.7-fold increase in hsa_circ_0001883 were identified in DR patients compared to DM patients. Meanwhile, the area under curves of hsa_circ_0095008 and hsa_circ_0001883 were 0.6710 (95% CI, 0.5646-0.7775) (p = 0.003399) and 0.6071 (95% CI, 0.4953-0.7189) (p = 0.06644), respectively, indicating a good diagnostic value.
Conclusion: Our study provided a new sight for the pathological mechanism of DR and revealed the potential value of hsa_circ_0095008 and hsa_circ_0001883 as diagnostic biomarkers for the early diagnosis of DR patients.

Introduction

The incidence of diabetes mellitus (DM) has increased dramatically worldwide in recent years, ranking as the ninth leading cause of death. Of these, 90% of the patients had type 2 diabetes mellitus (T2DM) (Zheng et al., 2018). Unhealthy dietary habits, lifestyle, and genetic factors were involved in the development of T2DM (Kautzky-Willer et al., 2016). It is reported that the vast majority of DM have at least one complication, including diabetic nephropathy, cardiovascular disease, and diabetic retinopathy (DR) (Naqshbandi et al., 2008; Gu et al., 2020). As a common complication of DM, DR is the main cause of impaired vision in diabetic patients. Retinal microvascular leakage and obstruction were attributed to fundus lesions, macular edema, and others (Gonzalez-Casanova et al., 2021). Currently, vascular endothelial growth factor A (VEGFA) inhibitors are the only drugs in clinical therapy that can effectively treat DR (Antonetti et al., 2021). Furthermore, vascular endothelial growth factor A (VEGFA) inhibitors are not effective in all patients with DR (Funatsu et al., 2009).

Circular RNAs (circRNAs) have become a new hot spot in the field of non-coding RNA research apart from microRNA (miRNA) and long non-coding RNA (Cai et al., 2019). CircRNAs are widely distributed in eukaryotic cells and participate in the pathogenesis and development of multiple types of diseases, including cancers, neurological diseases, and others (Altesha et al., 2019; Mengxue Xu et al., 2020; Lin et al., 2020; Wu et al., 2020; Fang et al., 2021). CircRNAs regulate various cellular activities via affecting RNA polymerase prolongation, acting as miRNA sponges to regulate target gene expression, and interacting with RNA-binding proteins to regulate the translation process (Zhang et al., 2013; Ebbesen et al., 2017; Hsiao et al., 2017). Increasing evidence shows that circRNAs are closely related to a variety of human diseases, such as tumors and DR (Wang et al., 2017; Zhang et al., 2017). For example, circHIPK3, circRNA cZNF609, and hsa_circ_0005015 have been proven to play a vital role in the progression of DR by regulating the growth, proliferation, migration, and tube formation of retinal vascular endothelial cells (Liu et al., 2017; Shan et al., 2017). Furthermore, previous studies identified that circRNA such as circular RNA-ZNF532 and circ-PSEN1 regulate DR progression by serving as miRNA sponges (Jiang et al., 2020; Zhu et al., 2021).

CircRNAs are universally expressed and conserved in human and vertebrate neural retina (Sun et al., 2019; Meng-Lan Li et al., 2021). In mouse and rat retinas, circRNA population increases during development, with significant developmental stage specificity (Han et al., 2017; George et al., 2019; Mellough et al., 2019; Kaining Chen et al., 2021; Gang Chen et al., 2021). CircRNAs are aberrantly expressed in retina-related diseases (Wang et al., 2018; Cao et al., 2019; Chen et al., 2020a; Sun et al., 2020; Sun et al., 2021). Furthermore, expression abnormalities of circRNAs appear earlier than the disease onset in a retinal degeneration model (Chen et al., 2020b).

Based on the prevailing biological functions of circRNAs, it suggests that circRNAs may be an ideal molecular marker for DR diagnosis and therapeutic targets. To investigate the circRNAs associated with DR occurrence, we analyzed circRNA expression profiles in healthy controls, DM, and DR patients, followed by the verification of differential expressions of circRNAs by RT-qPCR.

Materials and Methods

Peripheral Blood Mononuclear Cell Collection

Peripheral blood mononuclear cells (PBMCs) were obtained from 60 cases of DM patients, 42 cases of DR patients and 20 cases of healthy individuals by using anticoagulation tubes. Ficoll-Paque PLUS (GE Healthcare, United States) was added in blood samples with an equal amount of PBS. The mixture was centrifuged at 1,500 rpm for 40 min, and then in the middle, PBMCs were washed with PBS. Finally, PBMCs were kept at –80°C in TRIzol (Sigma-Aldrich, United States) for subsequent use. All participants agreed and signed the informed consent. The study was approved by the Research Ethics Board (REB) of the Affiliated People’s Hospital of Ningbo University (approval number 2019–048).

Total RNA Extraction and Quantitative Real-Time PCR

Total RNA was extracted from PBMCs using TRIzol (Sigma-Aldrich, United States) and then reverse transcribed into cDNA RNA using PrimeScriptTM RT reagent (Takara, United States). CircRNAs were quantified by qRT-PCR using Power SYBR Green Master Mix (Applied Biosystems, United States) on 7500 Fast Real-Time PCR System (Applied Biosystems, United States), according to manufacturer’s instructions. Primer sequences of all circRNAs are listed in Table1.

TABLE 1

Primers of validated circRNAs in RT-qPCR.

circRNA ID	Forward (5–3′)	Reverse (5–3′)
hsa_circ_0095008	ATG​CGA​CCA​TCC​ACC​TCA​AAG	ACA​TCA​CAC​ACA​ATC​ACG​GCA
hsa_circ_0005062	TCA​TCA​GCA​CCC​TGT​CGT​CT	CTG​CTT​TTC​CTG​TGA​TTT​TAC​CCA
hsa_circ_0001883	AGA​GAG​TAC​CAG​ACC​CGA​CA	GCA​AGT​GAG​CGA​AAT​GCT​CTT
hsa_circ_0040707	GCT​CTT​TGC​AGG​GTC​GAC​AA	AGT​GGT​TTT​TGG​GGC​CGT​TG
hsa_circ_0002031	GTG​ATC​GTT​GGC​GGA​CAT​TT	ATG​CTG​CTG​TCA​TGT​GCT​TCT
GAPDH	ATG​GAA​ATC​CCA​TCA​CCA​TCT​T	CGCCCCACTTGATTTTGG

Circular RNA Microarray Analysis

For circRNA microarray analysis, the total RNA sample was isolated using TRIzol (Sigma-Aldrich, United States) and then purified by an miRNA Isolation Kit (MACHEREY- NAGEL, Cat#740955, Germany), according to the manufacturer’s instructions. The extracted RNA was subsequently amplified and labeled with an Ambion WT Expression Kit (Cat#740955, Ambion, United States). Labeled samples were dissolved in a hybridization solution to load onto a Capital Bio Technology Human CircRNA Array v2 microarray (Agilent, United States) overnight. The circRNA microarray results were analyzed by Agilent GeneSpring software. The circRNAs with fold change ≥2 and p ≤ 0.05 were considered as upregulated or downregulated in circRNA microarray analysis.

CircRNA–miRNA Network Prediction

miRanda-3.3 software was used to predict circRNA-targeting miRNAs based on the degree of sequence complementarity between miRNAs and circRNAs. These circRNA–miRNA pairs were combined at entropy values below 20 and then constructed into networks using the open source bioinformatics software Cytoscape (v3.19.0, Institute of Systems Biology, United States).

Statistical Analysis

SPSS version 13.0 (SPSS Inc, IL, United States) was utilized for statistical analysis of all data in this study. The significant difference of RT-qPCR results was analyzed by using Student’s t test. For the comparison of each group, p < 0.05 was considered statistical significant.

Results

CircRNA Expression Profiles in PBMCs of Different Groups

PBMC samples were obtained and extracted from healthy controls and DM patients (with or without DR). Human CircRNA microarray v2 (Capital Bio Technology) was performed to detect the profile of circRNA expression, which revealed significant differences in DR patients compared to healthy controls (Figure 1A) or DM patients (Figure 1B) by hierarchical clustering. Volcano plot filtering was used to represent significant changes in differential circRNAs (FC ≥ 2 and p ≤ 0.05) between two groups. The results confirmed that a total of 104 circRNAs were significantly upregulated and 185 circRNAs downregulated in DR patients compared to DM patients (Figure 1C). In addition, 1,106 circRNAs were significantly upregulated and 1,346 circRNAs downregulated in DR patients, compared to healthy controls (Figure 1D).

FIGURE 1

Profile of circRNA level in DR samples. (A) Identified circRNAs between DR and control groups by heat map analyzes. (B) Identified circRNAs between DR and DM groups by heat map analyzes. (C) Differentially expressed circRNAs between DR and control groups by Volcano map analysis; (D) The differently expressed circRNAs between DR and DM groups by Volcano map analysis. The red dots represent significant upregulation, and green dots represent significant downregulation.

Functional Analysis of Differential circRNAs by Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) Analysis

We analyzed the parent genes of differential circRNAs by using KEGG and GO analyses to predict their biological functions. The top 10 pathways enriched in DR patients were mainly related to cellular components, molecular functions, and biological processes. Those enriched signaling pathways were associated with cell periphery, plasma membrane, extracellular matrix component, proteinaceous extracellular matrix, and extracellular matrix (Figure 2A). In terms of molecular functions, guanyl nucleotide exchange factor activity and Ras guanyl nucleotide exchange factor activity were the most enriched aspects (Figure 2B). The most enriched biological functions included multicellular organismal process, single-multicellular organism process, anatomical structure morphogenesis, and the movement of cell or subcellular components, all of which were associated with cell growth and proliferation (Figure 2C).

FIGURE 2

Differentially expressed circRNAs in DR samples by GO analysis. (A) DR vs control group, (B) DR vs DM group.

In addition, KEGG pathway analysis showed the top 30 enriched pathways in DR patients compared to DM patients (Figure 3A) and to healthy controls (Figure 3B). Specifically, those differential circRNAs were mainly focused on ECM–receptor interaction and focal adhesion.

FIGURE 3

Differentially expressed circRNAs in DR patients by KEGG analysis. (A) DR vs control group, (B) DR vs. DM group.

Validation of circRNA Level by RT-qPCR

To further identify biomarkers for DR patients, the top five upregulated or downregulated circRNAs (p ≤ 0.05 and raw processed signal ≥100) between DR and either DM or control groups are listed in Table 2. RT-qPCR technology was performed in an independent cohort (healthy controls, n = 20; DM patients without DR, n = 60; DR patients, n = 42) to detect these five candidate circRNAs (Figure 4). Detailed clinical characteristics of the study population are shown in Table 3. Trends of hsa_circ_0095008 and hsa_circ_0001883 were consistent with the results of circRNA arrays. hsa_circ_0095008 expression level was decreased by 2.18-fold (p = 0.013) and 2.47-fold (p = 0.001) in DR patients compared to DM patients and healthy controls, respectively (Figure 4A). hsa_circ_0001883 expression was increased by 1.73-fold (p = 0.025) and 3.04-fold (p = 0.015) in DR patients compared to DM and healthy controls, respectively (Figure 4B). A striking increase in hsa_circ_0040707 (p = 0.022) was identified in DR compared to healthy controls (Figure 4C). Moreover, the expression of hsa_circ_0005062 was significantly elevated (fold change = -1.953, p = 0.037) in DR patients, compared to DM patients (Figure 4D). However, there was no significant difference in hsa_circ_0002031 expression among these groups (Figure 4E).

TABLE 2

List of differentially expressed circRNAs in patients with DR.

circRNA ID	DR vs controls			DR vs DM
	Fold change	p-value	Regulation	Fold change	p-value	Regulation
hsa_circ_0095008	3.302807	0.000866	Down	2.057623	0.040505	Down
hsa_circ_0005062	2.764961	0.009428	Down	4.017943	0.000528	Down
hsa_circ_0001883	3.217526	0.004683	Up	2.126805	0.041181	Up
hsa_circ_0040707	3.140882	0.019335	Up	3.287805	0.016502	Up
hsa_circ_0002031	2.379722	0.005157	Down	3.232438	0.042273	Down

FIGURE 4

Validation of five candidate circRNAs in DR samples (n = 42) by RT-qPCR. (A) hsa_circ_0095008; (B) hsa_circ_0001883; (C) hsa_circ_0040707; (D) hsa_circ_0005062; (E) hsa_circ_0002031), compared with DM samples (n = 60) and controls (n = 20).*p < 0.05, **p < 0.01, ***p < 0.001.

TABLE 3

Detailed clinical parameters of study population.

	Controls	DM	DR
No.	20	60	42
Age (years)	41.35 ± 3.116	63.58 ± 1.344	65.511.271
Sex, female	9	35	22
BMI		24.35 ± 0.3988	23.72 ± 0.5683
HbA1c (%)		7.652 ± 0.1992	7.443 ± 0.2271
SBP (mmHg)		148.0 ± 2.570	144.2 ± 2.309
DBP (mmHg)		76.12 ± 1.965	74.29 ± 2.173
Glucose (mmol/L)		7.246 ± 0.2674	7.586 ± 0.3935
TCHO (mmol/L)		4.759 ± 0.1309	5.156 ± 0.1604
TG (mmol/L)		1.746 ± 0.1433	1.576 ± 0.1369
LDL (mmol/L)		1.255 ± 0.03701	1.275 ± 0.04632
HDL (mmol/L)		2.804 ± 0.1109	3.126 ± 0.1376
Tbil (umol/L)		12.90 ± 0.6937	12.94 ± 0.9465

BMI, body mass index; HbA1c, glycated hemoglobin A1c; SBP, systolic blood pressure; DBP, diastolic blood pressure; TCHO, total cholesterol; TG, triglyceride; LDL, low-density lipoprotein; LDL, low-density lipoprotein; TBil, total bilirubin.

Diagnostic Value of Differentially Expressed circRNAs in PBMCs of DR Patients

Receiver operating characteristic curve (ROC) was performed to analyze the sensitivity and specificity of circRNAs. The area under curve (AUC), cutoff value, sensitivity, and specificity of ROC analysis are listed in Table 4. These results showed that the AUCs of hsa_circ_0095008 were 0.671 (95% CI, 0.565–0.778) (p = 0.003) between DR and DM patients, and 0.804 (95% CI, 0.681–0.926) (p = 0.0001) between DR patients and healthy controls, respectively (Figure 5A). The AUCs of hsa_circ_0001883 were 0.607 (95% CI, 0.495–0.719) (p = 0.066) between DR and DM patients, and 0.725 (95% CI, 0.595–0.855) (p = 0.004) between DR patients and healthy controls (Figure 5B).

TABLE 4

AUC characteristics of hsa_circ_0095008 and hsa_circ_0001883.

Group	CircRNA ID	AUC	Sensitivity	Specificity	PPV	NPV
		(95% CI)	% (95%CI)	% (95%CI)	% (95%CI)	% (95%CI)
DR vs DM	hsa_circ_0095008	0.671	52.38	81.67	66.7	71
		0.565–0.778	36.42–68.00	69.56–90.48	50.6–82.8	60.3–81.7
	hsa_circ_0001883	0.607	50.00	70.00	52.5	66.1
		0.495–0.719	34.2–65.8	56.8–81.2	37–68	54.3–77.9
DR vs controls	hsa_circ_0095008	0.804	71.43	85	90.9	58.6
		0.681–0.926	55.42–84.28	62.11–96.79	81.1–100.7	40.7–76.5
	hsa_circ_0001883	0.725	61.90	75	81.25	46.7
		0.595–0.855	45.6–76.4	50.9–91.3	67.7–94.8	28.8–64.5

AUC, area under curve; PPV, positive predictive value; NPV, negative predictive value.

FIGURE 5

ROC analysis of hsa_circ_0095008 and hsa_circ_0001883 in DR samples. Diagnostic value of (A) hsa circ 0095008 and (B) hsa circ 0001883 for DR. Blue lines represent ROC curves of circRNAs for distinguishing DR patients from healthy controls ,while red lines represent ROC curves of circRNAs for distinguishing DR patients from DM patients without DR.

We further considered whether these circRNAs could serve as diagnostic biomarkers to distinguish DR from DM. The optimal cutoff value of hsa_circ_0095008 was 0.294, with a sensitivity of 52.38%, a specificity of 81.67%, a positive predictive value (PPV) of 66.7%, and a negative predictive value (NPV) of 71%. The cutoff value of hsa_circ_0001883 was 0.066, with a sensitivity of 50%, a specificity of 70%, a PPV of 52.5%, and an NPV of 66.1% for DR patients.

To diagnose DR from healthy controls, the sensitivity and specificity of hsa_circ_0095008 were 71.43 and 85%, respectively (cutoff value = 0.447, PPV = 90.9%, NPV = 58.6%); the sensitivity and specificity of hsa_circ_0001883 were 61.9 and 75%, respectively (cutoff value = 1.744, PPV = 81.25%, NPV = 46.7%). Taken together, hsa_circ_0095008 and hsa_circ_0001883 have diagnostic potential for DR.

Construction of circRNA–miRNA Network

It is well established that circRNAs have potential to regulate mRNA levels as miRNA sponges (Lyu and Huang, 2017). CircRNAs contain many miRNA binding sites and act as miRNA sponges to regulate the expression of miRNA targets (Ebbesen et al., 2017). For example, circRNA ciRS-7 contains 70 conserved binding sites for miR7. Thus, miR-7 levels can be increased by inhibiting the expression of circRNA ciRS-7 (Hansen et al., 2013). To explore whether hsa_circ_0095008 and hsa_circ_0001883 were miRNA sponges, we predicted their binding miRNAs using miRanda-3.3 software combined with entropy values below 20. Bioinformatics software Cytoscape was utilized to construct circRNA–miRNA networks. The prediction results demonstrated that they had more than 100 binding miRNAs, respectively, suggesting that hsa_circ_0095008 and hsa_circ_0001883 might regulate miRNA levels by acting as miRNA sponges in PBMCs (Figure 6). Among them, only 2 miRNAs had more than one binding site for hsa_circ_0095008, whereas there were 100 miRNAs with more than one binding site for hsa_circ_0001883. The current study also demonstrated that beyond acting as miRNA sponges, circRNAs could also serve as protein sponges or encode proteins to perform regulatory functions (Yang et al., 2018; Zhang et al., 2018; Huang et al., 2020; Wu et al., 2021). Thus, the additional roles of hsa_circ_0095008 and hsa_circ_0001883 in DR progression are still urgently needed to be further investigated.

FIGURE 6

Target miRNA prediction of hsa_circ_0095008 and hsa_circ_0001883 in DR samples. The CircRNA/miRNA regulatory networks of (A) hsa_circ_0095008 and (B) hsa_circ_0001883 were predicted by miRanda- 3.3 software.The red and blue dots represent circRNAs and miRNAs, respectively.

Discussion

As one of the most common diabetic complications, DR is more difficult to diagnose than others, considering that its early symptoms are not obvious and the clinical presentation is not very specific (Stitt et al., 2016). Generally, DR patients could be diagnosed only when irreversible eye damage occurred, including blurred vision and eventual blindness (Wong et al., 2016; Sheen et al., 2020). The increasing incidence of retinal neovascularization is the core factor contributing to DR (Li et al., 2018). However, the current therapeutic option for DR is still limited due to severe side effects (Sun and Jampol, 2019; Crabtree and Chang, 2021; Everett and Paulus, 2021). Therefore, it is crucial to improve the early diagnosis of DR.

CircRNAs, a class of circular non-coding RNAs, play vital biological functions in many physiological processes. CircRNAs were extensively explored as diagnostic and predictive biomarkers for various diseases, including cancer, neurological disease, inflammatory bowel disease, skin diseases, and Alzheimer’s disease (Ye et al., 2019; Zhang et al., 2019; Wang et al., 2020; Wu et al., 2020; Shengnan Li et al., 2021). To explore more sensitive biomarkers for DR screening, many samples extracted from different optical tissues, including retinal, vitreous humor, and serum, were studied (Gu et al., 2017; Jiang et al., 2020; He et al., 2020). He et al. studied the profile of circRNA expression in the vitreous humor of patients with proliferative diabetic retinopathy (PDR) (He et al., 2020). Gu et al. testified about many dysregulated circRNAs in serum from DR patients. hsa_circ_063981, hsa_circ_404457, hsa_circ_100750, hsa_circ_406918, hsa_circRNA_104387, hsa_circ_103410, and hsa_circ_100192 were significantly upregulated in DR patients compared to DM patients and healthy controls (Gu et al., 2017). Li et al. identified the profile of exosomal circRNAs in serum of PDR patients to explore their possible pro-angiogenic role (Xinsheng Li et al., 2021).

As it is difficult to obtain RNA from vitreous fluidand and serum exosome, and exon circRNA in serum are unstable, PBMCs are the ideal material for clinical diagnosis of DR (Wen et al., 2020). Hence, our study first identified the profile of differential circRNA in the PBMCs of healthy controls, DM patients, and DR patients. According to our outcomes represented in Figure 3, the differential circRNAs in PBMCs of DR patients had an important role in cell migration based on enrichment in ECM–receptor interaction and focal adhesion pathways. It was also confirmed by RT-qPCR in Figures 4A,B that hsa_circ_0095008 significantly decreased and hsa_circ_0001883 significantly upregulated in the DR compared to DM and control groups, which is consistent with the results of microarray analysis. It was summarized that these two abnormal circRNAs were closely associated with the occurrence of DR.

Interestingly, orosomucoid-1 (ORM1), the host gene of hsa_circ_0001883, is an acute phase response (APR) protein that regulates angiogenesis by stimulating VEGF (Luo et al., 2015). Increasing evidence suggests that retinal neuronal degeneration started in the early stages of DR (Srinivasan et al., 2017; Rolev et al., 2021), while hsa_circ_0095008 is spliced from neural cell adhesion molecule 2 (NCAM2) that modulates neuronal morphogenesis and differentiation (Parcerisas et al., 2021), which also indicates that circRNAs play a potential role in the early stages of DR. Moreover, circRNAs mainly exert their regulatory functions as miRNA sponges. They can affect mRNA stability, regulate host gene transcription, bind to RNA-binding proteins, and even directly translate proteins (Zhou and Kuang, 2021). It is reported that circSMARCA5 inhibits host gene transcription by interacting with the host gene (Xiaolong Xu et al., 2020).

However, our study has several limitations. First, more clinical samples are warranted to inspect the potential role of abnormal circRNAs as early diagnostic biomarkers of DR. Second, both the region and Asian species of samples also limit their representation and further diagnostic application. Last, the mechanism of these functional circRNAs in DR progression is still unclear.

In conclusion, our study provided a new sight for the pathological mechanism of DR and revealed the potential value of hsa_circ_0095008 and hsa_circ_0001883 as non-invasive biomarkers for the early diagnosis of DR.