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 sex-matched 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.

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).

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 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.

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.

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.

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-μl 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 β-actin, and all reactions were performed in triplicate. The designed and optimized primers for circRNA transcripts are shown in Supplemental Table S1 .

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.

Data are presented as the mean ± SD (standard deviation) or as the median (25th–75th percentile). The 2^-ΔΔCT 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.

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 ).

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 2-oxocarboxylic 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 ).

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 difference in FT1D and controls (P<0.0001) regardless of before or after adjusting for age ( Table S2 ).

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 ).

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).

CircRNAs can function as competing endogenous RNA that can sponge corresponding miRNAs and affect mRNA expression in various diseases, including diabetes. Hsa_circRNA_100632 was predicted by circBank, a comprehensive database of human circRNA, to sponge 70 miRNAs. Five miRNAs related to diabetes were selected (hsa-miR-1-3p, hsa-miR-27a-3p, hsa-miR-221-3p, hsa-miR-222-3p, and hsa-miR-503-5p) for miRNA–mRNA prediction. According to TargetScan and miRanda software, these 5 miRNAs may bind to 333 mRNAs. We next selected 41 mRNAs related to diabetes to form a specific circRNA–miRNA–mRNA network related to diabetes. The results showed that hsa_circRNA_100632 may be involved in 47 circRNA–miRNA–mRNA signaling pathways associated with diabetes ( Figure 6 ). For example, hsa_circRNA_100632 sponges hsa-miR-221-3p and hsa-miR-222-3p, resulting in modulation of the thrombospondin 1 (THBS1), solute carrier family 40 member 1 (SLC40A1), and aldehyde dehydrogenase 1 family, member A1 (ALDH1A1) target genes.