Identification of hsa_circRNA_100632 as a novel molecular biomarker for fulminant type 1 diabetes

Objective Circular RNAs (circRNAs) are associated with diabetes, but their role in fulminant type 1 diabetes (FT1D) is unclear. Thus, we characterized the role of circRNAs in FT1D. Research design and methods CircRNA expression profiles were detected in peripheral blood mononuclear cells (PBMCs) of five FT1D patients and five controls using a circRNA microarray. An independent cohort comprised of 40 FT1D cases, 75 type 1 diabetes (T1D) cases, and 115 controls was used to verify the circRNAs using quantitative real-time polymerase chain reaction (qRT-PCR). Spearman’s correlation analysis and receiver operating characteristic (ROC) curve analysis were performed to determine the clinical diagnostic capability of circRNAs. Bioinformatics was used to identify potential biological functions and circRNA–miRNA–mRNA interactions. Results There were 13 upregulated and 13 downregulated circRNAs in PBMCs of patients with FT1D. Five circRNAs were further verified in a second cohort. Hsa_circRNA_100632 was significantly upregulated in the FT1D and T1D groups. Hsa_circRNA_100632 was differentiated between patients with FT1D and controls [area under the curve (AUC) 0.846; 95% CI 0.776–0.916; P<0.0001] as well as between patients with FT1D and patients with T1D (AUC 0.726; 95% CI 0.633–0.820; P<0.0001). Bioinformatics analysis showed that hsa_circRNA_100632 may be involved in 47 circRNA–miRNA–mRNA signaling pathways associated with diabetes. Conclusions CircRNAs were aberrantly expressed in PBMCs of patients with FT1D, and hsa_circRNA_100632 may be a diagnostic marker of FT1D.


Introduction
Fulminant type 1 diabetes (FT1D) is a novel subtype of type 1 diabetes (T1D) first identified by Imagawa et al. in 2000 (1, 2). FT1D is characterized by aggressive disease onset, dramatic loss of Cpeptide concentration, nearly normal levels of hemoglobin A1c (HbA1c) at onset, and diabetic ketoacidosis (DKA) (3,4). Failure to accurately diagnose FT1D can result in death (5,6). Diagnosis of FT1D requires measurement of plasma glucose (PG) at onset, measurement of HbA1c, and onset of clinical symptoms of DKA within 1 week after onset. Then, fasting C-peptide (FCP) and 2 h postprandial C-peptide (PCP) must be measured after complete remission of ketoacidosis, which typically occurs within 1 or 2 weeks, to confirm the diagnosis (7,8). These criteria make the diagnosis of FT1D difficult and complex, preventing precise intervention and worsening prognosis. Therefore, it is of great clinical significance to explore markers for identification of FT1D to allow for simple or early diagnosis and effective intervention.
Circular RNAs (circRNAs), a new class of non-coding RNAs, are characterized by a closed-loop structure, resistance to degradation, and excellent stability (9,10). Several studies have shown that circRNAs are involved in many physiological and pathological disease processes, and circRNAs can serve as diagnostic and prognostic molecular biomarkers (10,11). Our previous study showed that several circRNAs are differentially expressed in the peripheral blood of T1D patients and may contribute to T1D disease progression through interactions with miRNA and circRNA-derived peptides (12). Other studies have shown that circRNAs promote the onset and progression of T1D (13,14). However, the expression profiles of circRNAs and the value of circRNAs as biomarkers of FT1D have not been characterized.
In the present study, we investigated the circRNA expression profiles in peripheral blood mononuclear cells (PBMCs) from patients with FT1D using microarray data, analyzed their biological functions using bioinformatics, and validated their expression levels using quantitative real-time polymerase chain reaction (qRT-PCR). As a result, a circRNA that can be used for diagnosis and differential diagnosis of FT1D was identified using a receiver operating characteristic (ROC) curve. The present study demonstrated that circRNAs have potential to aid in the diagnosis and treatment of patients with FT1D.

Participants
A total of 240 individuals were placed into two cohorts, and Table 1 shows the clinical and demographic data of these individuals. We selected five patients with new onset FT1D (median diabetes duration one month) and five age-and sexmatched control individuals for analysis in the first cohort. In the second cohort, 115 patients with FT1D (N = 40) or T1D (N = 75) and 115 controls were recruited to validate differentially expressed circRNAs identified in the first cohort.
Patients with diabetes met the diagnostic criteria of the American Diabetes Association (15). Individuals with T1D were recruited based on the following criteria: (1) acute onset of ketosis or ketoacidosis requiring insulin replacement therapy; (2) at least six months of continuous insulin dependence after diagnosis; and (3) positive glutamic acid decarboxylase antibodies (GADAs). Diagnosis of FT1D was performed according to the criteria of the Committee of the Japan Diabetes Society in 2012 (16). The inclusion criteria for FT1D were as follows (1): occurrence of diabetic ketosis or ketoacidosis within a week after onset of hyperglycemic symptoms; (2) initial PG level > 16.0 mmol/L and HbA1c < 8.7% (72 mmol/mol) at onset; and (3) serum FCP < 100 pmol/L (0.3 ng/ml) and PCP < 170 pmol/L (0.5 ng/ml) at disease onset. Controls had fasting plasma glucose (FPG) < 5.6 mmol/L and HbA1c < 6.1%. The exclusion criteria for control subjects were infectious diseases, pregnancy, malignant disease, or a family history of diabetes.
All data were collected from patients from the Second Xiangya Hospital of Central South University (Changsha, Hunan, China). The study was approved by the Ethical Committee of the Second Xiangya Hospital of Central South University (LYF 2021100), and all subjects provided written informed consent prior to inclusion in the study.

Data collection and clinical measurements
Sex, age, body height, and weight were recorded for all subjects.
Fasting venous blood samples were tested for triglycerides (TG), total cholesterol (TC), high density lipoprotein-cholesterol (HDL-C), low density lipoprotein-cholesterol (LDL-C), FPG, and HbA1c using standard procedures in the hospital clinical laboratory. GADAs were detected as previously described (17).

Peripheral blood mononuclear cell collection
Fasting peripheral blood was collected from patients and controls into EDTA blood tubes, and PBMCs were isolated using density gradient centrifugation. Ficoll-Paque PLUS (GE Healthcare, USA) was added to the blood samples with an equal amount of phosphate buffered saline (PBS; Invitrogen, USA). Finally, PBMCs were stored in TRIzol reagent (Invitrogen, USA) at -80°C.

Total RNA extraction and quality control
Total RNA was extracted from PBMCs of each participant using TRIzol reagent (Invitrogen, USA) according to the manufacturer's instructions. The RNA purity and concentration were determined using a NanoDrop 2000 instrument (Thermo Scientific, USA). The range of mean A260/A280 was 1.8 to 2.0.

CircRNA microarray expression profiling
Sample preparation and microarray hybridization were performed according to the manufacturer's instructions (Arraystar Human circRNA Array V2.0, 8 ×15K). Total RNA was digested with Rnase R (Epicentre, USA) to remove of linear RNA and enrich circRNA. The enriched circRNAs were then amplified and transcribed into fluorescent cRNA according to the Arraystar Super RNA Labeling protocol (Arraystar, USA). Finally, the hybridized arrays were washed, fixed, and scanned using an Agilent Scanner G2505C. The details of the method have been previously described (12). Agilent Feature Extraction software (version 11.0.1.1) was used to analyze the acquired array images. Quantile normalization and subsequent data processing were performed using the limma package in R software. Significantly differentially expressed circRNAs between the FT1D group and control group were identified using volcano plot filtering. CircRNAs with absolute fold differences ≥ 1.5 and P < 0.05 were considered significantly differentially expressed. Hierarchical clustering was performed to show the distinguishable circRNA expression pattern among samples.

GO and KEGG pathway enrichment analyses
OmicShare tools was used to perform Gene Ontology (GO) and Kyoto Encyclopedia Genes and Genomes (KEGG) analyses. GO analysis was conducted based on the parent genes of circRNAs to investigate the properties of aberrantly expressed circRNAs. Pathways were analyzed based on differentially expressed circRNAs using KEGG analysis.

QRT-PCR assay
One microgram of total RNA was converted into cDNA using a GoScript ™ Reverse Transcription System (Promega, USA) according to the manufacturer's instructions. Target gene expression levels were determined using a SYBR Green kit (Promega, USA) in 10-ml samples using an ABI PRISM Step One Sequence Detection System (Applied Biosystems, USA). The thermocycler conditions were as follows: 95°C for 10 min followed by 45 cycles of 95°C for 15 s and 60°C for 1 min. Internal standard was b-actin, and all reactions were performed in triplicate. The designed and optimized primers for circRNA transcripts are shown in Supplemental Table S1.

CircRNA-miRNA-mRNA network analysis
CircRNA can sponge miRNA, and mRNA can be targeted by miRNA, which may play critical roles in disease onset and progression. To elucidate the relationship between circRNAs and miRNAs, TargetScan and miRanda prediction software were used to predict the interactions of circRNAs-miRNAs-mRNAs. The circRNA-miRNA-mRNA network was constructed using the Cytoscape bioinformatics software.

Statistical analysis
Data are presented as the mean ± SD (standard deviation) or as the median (25th-75th percentile). The 2 -DDCT method was used to calculate the relative expression levels of selected circRNAs analyzed using qRT-PCR. Normality of the data was checked using the Shapiro-Wilk test. One-way analysis of variance (ANOVA) was performed to compare groups if the data were normally distributed, and a rank test (Mann-Whitney U test or Kruskal-Wallis H test) was used if the data were not normally distributed. The Chi-square test was used to analyze categorical data among groups. Generalized linear regression model analysis was used to identify the associations between circRNA expression levels and FT1D after adjusting the confounders. Spearman's correlation coefficient was used to evaluate the relationship between circRNAs and clinical characteristics of FT1D. A ROC curve was used to assess the diagnostic and differential diagnostic performance of circRNA for FT1D. Statistical analyses were performed using SPSS 21.0 (IBM SPSS, Inc., Chicago, IL, USA) and GraphPad Prism 5.01 (GraphPad Software, San Diego, CA, USA). P < 0.05 was considered statistically significant.

Identification of aberrantly expressed circRNAs in PBMCs of patients with FT1D
We detected 13,563 circRNAs in PBMCs from patients with FT1D. The classification of circRNAs according to expression level was performed by hierarchical clustering ( Figure 1A). In the volcano plots ( Figure 1B), separately expressed circRNAs were grouped by absolute fold changes ≥ 1.5 and P < 0.05, resulting in 26 differentially expressed circRNAs in patients with respect to controls, of which 13 were upregulated and 13 were downregulated ( Table 2).

Bioinformatics analysis of differentially expressed circRNAs
To predict the potential biological functions of the 26 differentially expressed circRNAs, we used their parent genes for GO enrichment and KEGG pathway analyses. The GO terms for parent genes of differentially expressed circRNAs were classified and summarized according to biological process (BP), cellular component (CC), and molecular function (MF). The terms with the most genes in BP, CC, and MF were cellular process, binding, and cell, respectively (Figure 2A). For BP, the term with the highest enrichment score was negative regulation of phosphate metabolic process ( Figure 2B), and the term with the highest enrichment score for CC was lamellar body membrane and alveolar lamellar body membrane ( Figure 2C). For MF, the term with the highest enrichment score was D-ribulokinase activity ( Figure 2D).
The enriched pathways of the parent genes of the 26 differentially expressed circRNAs were analyzed using KEGG pathway annotation. The non-homologous end-joining and 2oxocarboxylic acid metabolism terms had the highest correlations with the enrichment scores ( Figure 3A). Furthermore, KEGG pathway annotation showed that most genes were involved in cell growth and death as well as the immune system ( Figure 3B).

Validation of differentially expressed circRNAs
To further validate the microarray data, five differentially expressed circRNAs, including three upregulated circRNAs (hsa_circRNA_ 100246, hsa_circRNA_100245, and hsa_circRNA_100632) and two downregulated circRNAs (hsa_circRNA_005528 and hsa_ circRNA_406299), were selected for qRT-PCR verification. The following four inclusion criteria were used: (1) circRNA with a high differential expression fold; (2) circRNA length between 200 bp and 1000 bp; (3) circRNA average raw intensity > 100; and (4) exonic-circRNAs. As shown in Figure 4, patients with FT1D and T1D had significantly increased expression of hsa_circRNA_100632 (P<0.0001) and significantly decreased expression of hsa_circRNA_005528 (P<0.01) compared to controls. Moreover, patients with FT1D showed higher expression of hsa_circRNA_100632 than patients with T1D (P<0.0001) and hsa_circRNA_100632 has a significant A B difference in FT1D and controls (P<0.0001) regardless of before or after adjusting for age (Table S2).

Association of significantly differentially expressed circRNAs with clinical parameters
To analyze whether significantly differentially expressed circRNAs are associated with clinical parameters, Spearman's correlation analysis was performed in the FT1D group, the T1D group, the control group, and for all subjects. In patients with T1D, the expression levels of hsa_circRNA_100632 were positively correlated with FCP (P<0.05), and the expression levels of hsa_circRNA_005528 were positively associated with sex (P<0.01). In all subjects, hsa_circRNA_100632 was negatively correlated with age (P<0.05), BMI (P<0.01), TG (P<0.05), TC (P<0.01), and LDL-C (P<0.01), but hsa_circRNA_100632 was positively associated with FPG (P<0.01) and HbA1c (P<0.01). In addition, hsa_circRNA_005528 was negatively correlated with HbA1c in all subjects (P<0.01) ( Table 3).

Diagnostic value of hsa_circRNA_100632 as a biomarker of FT1D
To further evaluate whether has_circRNA_100632 has molecular diagnostic value for FT1D, we generated a ROC curve to evaluate its sensitivity and specificity. As shown in Figure 5, the area under the curve (AUC) for hsa_circRNA_100632 was 0.846 (95% CI 0.776-0.916) with a sensitivity of 75.0% and specificity of 79.1% between patients with FT1D and controls (P<0.0001). Furthermore, the AUC of hsa_circRNA_100632 was 0.726 (95% CI 0.633-0.820) with a sensitivity of 82.5% and specificity of 58.8% between patients with FT1D and patients with T1D (P<0.0001).

Discussion
In the present study, we profiled circRNA expression in PBMCs of patients with FT1D. We verified that hsa_circRNA_100632 was differentially upregulated and that hsa_circRNA_005528 was differentially downregulated among patients with FT1D, patients with T1D, and controls. We found that hsa_circRNA_100632 levels were associated with islet b-cell destruction in patients with T1D. Previous studies have suggested that circRNAs may play a role in T1D or type 2 diabetes (T2D) through regulation of islet b-cell dysfunction (18,19). CircGlis3, a b-cell-derived exosomal circRNA, QRT-PCR analysis of five differentially expressed circRNAs among patients with fulminant type 1 diabetes, patients with type 1 diabetes, and controls. (fulminant type 1 diabetes, N = 40; type 1 diabetes, N = 75; controls N = 115). ** P < 0.01 and ****P < 0.0001. has been shown to play a role in lipotoxicity-induced b-cell disorder and development of diabetes through suppression of insulin secretion and cell proliferation (20). In the present study, we identified hsa_circRNA_100632 as an important circRNA in FT1D. Has_circRNA_100632 is derived from sterile alpha motif domain containing 8 (SAMD8), which is a ceramide phosphoethanolamine synthase in the endoplasmic reticulum. This circRNA may act as a ceramide sensor, and its N-terminal SAM structural domain controls endoplasmic reticulum ceramide levels and inhibits ceramide-induced mitochondrial apoptosis (21). High glucose can lead to mitochondrial dysfunction through oxidative stress (22). We hypothesized that hsa_circRNA_100632 may promote progression of T1D through induction of b-cell failure. Future studies that include functional experiments using immune cells and animal models are needed.
The present study indicated that hsa_circRNA_100632 in PBMCs may be a diagnostic marker of FT1D. Increasing evidence has indicated that circRNAs may be diagnostic indicators of T1D and T2D (23). Hsa_circ_0054633, a diagnostic biomarker of prediabetes and T2D, shows diagnostic capability with AUC values of 0.751 and 0.793, respectively (24). In addition, hsa_circ_0063425 and hsa_circ_0056891 may be valuable markers for early detection Receiver operating characteristic (ROC) plot of hsa_circRNA_100632 in fulminant type 1 diabetes. The blue line indicates the ROC curve for hsa_circRNA_100632 distinguishing fulminant type 1 diabetes from controls, and the black line represents the ROC curve for hsa_circRNA_100632 distinguishing fulminant type 1 diabetes from type 1 diabetes. Hsa_circRNA_100632-miRNA-mRNA network. Blue ellipses, pink hexagons, and red octagons represent mRNAs related to diabetes, diabetes-associated miRNAs, and hsa_circRNA_100632, respectively. The orange solid lines indicate the relationship between miRNAs and mRNAs, and the dotted lines represent the relationship between circRNAs and miRNAs.
of T2D (25). A recent study has suggested that the expression level of hsa_circ_0060450 is upregulated in patients with T1D and may represent a novel therapeutic target for T1D (14). Furthermore, our previous study found that hsa_circ_0072697 is involved in 50 signaling pathways related to diabetes, indicating that it may be a diagnostic indicator of T1D (12). However, the role of circRNA in FT1D has not been characterized. The present study is the first to indicate that circRNAs may play a key role in FT1D, suggesting that circRNAs may be valuable diagnostic and differential diagnostic markers for FT1D. ROC analysis showed that hsa_circRNA_100632 had a sensitivity of 75.0% and specificity of 79.1% for distinguishing between patients with FT1D and control subjects (AUC = 0.846). Moreover, hsa_circRNA_100632 distinguished between FT1D and T1D with 82.5% sensitivity and 58.8% specificity (AUC = 0.726). Several diagnostic markers have been previously explored in smaller cohorts, and models have been built from data derived from multiple clinical studies (5,(26)(27)(28). For example, the glycated albumin (GA)/HbA1c ratio has been verified as a sensitive marker for glucose excursion of FT1D in a cohort comprised of 56 outpatients (27). The fulminant index shows an ability to identify FT1D based on DKA (5). In addition, the serum 1,5anhydroglucitol/GA index has been shown to be a suitable indicator for early differential diagnosis between FT1D and T1D when HbA1c < 8.7% with an optimal cut-off point of 0.3 (28). In contrast, the present study focused on the expression levels of hsa_circRNA_100632, which may simplify and reduce the cost associated with diagnosis of FT1D. In addition, KEGG analysis showed that differentially expressed circRNAs may be associated with cell growth and death as well as the immune system. FT1D has a rapid onset, severe symptoms, and a high mortality rate (29,30). Initially, FT1D was considered related to nonautoimmune pathogenesis due to the absence of islet autoantibodies, such as GADAs (1). Over the past 20 years, FT1D has gained increasing attention. Studies have shown that immunity plays a critical role in the initiation and progression of FT1D (29,31). Cellular and humoral immunity may play an important role in FT1D (29). Viral infection and drug induction may lead to b-cell loss in FT1D by triggering an immune response (32,33). The present results indicated that circRNAs were abnormally expressed in peripheral immune cells in FT1D, and bioinformatics analysis indicated that circRNAs may be involved in the immune response. Furthermore, the present results showed that hsa_circRNA_100632 may be involved in 47 diabetes-associated circRNA-miRNA-mRNA interaction signaling pathways. In the hsa_circRNA_100632 network, hsa_circRNA_100632 may sponge hsa-miR-27a-3p and modulate the nuclear factor of activated T cells 5 (NFAT5) target gene. NFAT5 activates various immune cells, especially T lymphocytes (34). These findings suggest that hsa_circRNA_100632/hsa-miR-27a-3p/NFAT5 may play a key role in FT1D by regulating the immune response. However, it is not known whether the upregulated expression levels of hsa_circRNA_100632 are due to increased expression of specific leukocytes in PBMCs of FT1D patients or whether the increasing hsa_circRNA_100632 can affect specific leukocyte expression. Further studies are needed to verify this finding.
The major strength of the present work was that it included a large validation cohort considering the rarity of FT1D. Furthermore, the protocols for recruitment and examination of subjects were highly standardized. However, the present study had several limitations. The sample size of the discovery cohort was small and the multiple-testing correction was not performed in the discovery analysis, indicating that the relevance of hsa_circRNA_100632 to the pathogenesis of FT1D requires further immune cell and animal studies. In addition, the present study focused on FT1D in a Chinese population, and further studies are needed to determine if these results regarding hsa_circRNA_100632 can be generalized to other ethnic groups.
In summary, we characterized the circRNA transcriptome in PBMCs derived from patients with FT1D for the first time. We then predicted the biological functions of differentially expressed circRNAs using GO and KEGG pathway analyses. The results indicated that hsa_circRNA_100632 may be a molecular biomarker for FT1D. Longitudinal studies are needed to validate hsa_circRNA_100632 as a molecular biomarker of FT1D progression.

Data availability statement
The original contributions presented in the study are publicly available. This data can be found here: https://ngdc.cncb.ac.cn/ omix.(OMIX repository, accession number OMIX002351).

Ethics statement
The studies involving human participants were reviewed and approved by the ethical committee of the Second Xiangya Hospital of Central South University. The patients/participants provided their written informed consent to participate in this study. Written informed consent was obtained from the individual(s) for the publication of any potentially identifiable images or data included in this article.

Author contributions
WY contributed to the experiments and data analysis and wrote the first draft of the manuscript. SL and ZZho proposed the project, designed the study, provided critical suggestions, and revised the manuscript. JQ, ZXia, and ZZha collected the study data and samples. JQ, ZXie, and XL reviewed the manuscript and contributed to the discussion. All authors contributed to the article and approved the submitted version.