Comprehensive circRNA Expression Prole and Transcriptomics-Based Molecular Pathway Analysis of the Sodium-Glucose Co-Transporter-2 Inhibitor Dapagliozin in Diabetic Tubular Epithelial Cells

Diabetic kidney disease (DKD) is a serious diabetes complication. Sodium-glucose cotransporter 2 inhibitors (SGLT2i) are novel anti-diabetes drugs and have clinical renal protection. However, the molecular mechanisms involved remain unclear. Here, human proximal tubular epithelial cells (PTECs) were treated with normal glucose (NG), high glucose (HG), and three types of anti-diabetes agents including SGLT2i (dapagliflozin), metformin, and dipeptidyl peptidase-4 inhibitor (DPP4i, vildagliptin) to perform microarray analysis. A total of 2,710 differentially expressed (DE) circRNAs were identified in PTECs. Network pharmacology, Gene Ontology (GO) and the Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses showed that the effects of dapagliflozin on PTECs primarily involved lipid metabolism, Rap1, and MAPK signaling pathways. The effects of metformin were mainly on the AMPK and FOXO signaling pathways, whereas vildagliptin were on insulin secretion and the HIF-1 signaling pathway. Furthermore, circRNA-miRNA-mRNA networks revealed that dapagliflozin might regulate the primary bile acid biosynthesis, arginine and proline metabolism, and vesicular transport through the hsa_circ_1586-miR-4739/hsa-miR-7851-3p/miR-1273g-3p networks, and regulate the MAPK, ErbB, and insulin resistance signaling pathways through the has_circ_012448-hsa-miR-378g/hsa-miR-29b-5p networks. Overall, our study elucidates circRNA expression profile in PTECs treated with dapagliflozin for the first time, providing novel clues for exploring the molecular mechanisms of dapagliflozin on DKD. diabetes 1 . Considering the high incidence rate (40%) of new-onset CKD in patients with diabetes and the limited efficacy of renin angiotensin system inhibitors 2,3 , new therapeutic strategies are urgently needed to improve the renal prognosis of patients with diabetes. Tubular atrophy, interstitial fibrosis, and rarity of tubular capillaries are closely related to renal function decline in patients with DKD. The sodium-glucose cotransporter 2 inhibitors (SGLT2i) are new oral anti-diabetes agents, that can specifically bind to SGLT2 in the proximal tubular, reduce glucose reabsorption, promote urinary glucose excretion, and effectively reduce blood glucose 4 . Multicenter clinical trials have demonstrated that SGLT2i considerably hinders DKD progression (EMPA-REG 5 , CANVAS 6 , and Declare-TIMI 7 ). Previous studies have indicated that SGLT2i can protect the kidneys by regulating tubuloglomerular feedback. However, renal outcome trial (ROT) studies, including CREDENCE 8 and DAPA-CKD 9 , have confirmed that canagliflozin and dapagliflozin effectively reduce the risk of end stage renal disease (ESRD) and cardiovascular death in patients with or without DKD, suggesting the existence of other mechanisms. However, the exact molecular backgrounds remain unclear. addition, a new report indicated that metformin promoted autophagy in the kidneys of individuals with chronic kidney disease, which may contribute to the actions of SGLT2i increasing autophagy and muting inflammation in diabetic kidneys 21 . Here, we performed functional analysis of circRNAs and mRNAs co-regulated by three drugs. GO and KEGG analyses revealed several important pathways related to cell communication, cell death, transcription factor, MAPK, p53, ErbB, and HIF-1 signaling pathways, which are reported to participate in the pathophysiological process of DKD. In addition, we identified the hub genes co-regulated by three drugs. Notably, t has been demonstrated that the hub genes VAV1 22 and CDKN1A 23 are involved in cell proliferation, differentiation and cell cycle regulation. VEGFA plays a key role in the progression of renal fibrosis 24 . DDIT3 is involved in ischemic injury and oxidative stress responses 25 . These findings suggested that the combined use of the three drugs may improve renal tubular injuries by regulating PTEC proliferation, apoptosis, and fibrosis, providing clues for further exploration of the mechanisms of multi-drug combination in vivo . indicated that SGLT2i is involved in DKD tubular cell apoptosis, tubulointerstitial and fibrosis. Our previous study also has reported that dapagliflozin improves tubulointerstitial fibrosis in DKD. Another study indicated that molecular pathways targeted by dapagliflozin and associated with DKD were related to energy metabolism, mitochondrial and endothelial functions 26 . Based on the CNC network constructed by the most significant DEcircRNAs, we revealed that the target genes of circRNAs upregulated by HG were associated with calcium ion transport, cell apoptosis, intercellular junctions, and mitophagy pathways. In contrast, the target genes of circRNAs downregulated by HG were associated with gene expression, oxoacid metabolic process, amino acid metabolism, and lysosome pathways. However, these target genes of DEcircRNAs induced by dapagliflozin were most enriched in gene expression, glucose metabolic process, extracellular exosome, epigenetics, and amino acid metabolism pathways. Strikingly, extracellular exosome signaling provides us with a new hypothesis, that is, whether the exosomes secreted by renal tubules induced by dapagliflozin affect the glomerulus, thereby improving the glomerular function and restoring the renal function. Collectively, these studies indicated that dapagliflozin protects DKD through multiple signaling pathways are feasible. In particularly, some signaling pathways have not been previously reported in the field of DKD and may be the starting point for future studies. We selected hsa_circRNA_001586 and hsa_circRNA_012448 to construct ceRNA network and revealed a new interaction between circRNA, miRNA and mRNA related to dapagliflozin. The hub miRNAs found in this study include hsa-miR-4739, hsa-miR-7851-3p, hsa-miR-1273g-3p, hsa-miR-378g and hsa-miR-29b-2-5p. Previous studies have reported hsa-miR-4739 is elevated in the urinary exosomes of patients with


Introduction
Diabetic kidney disease (DKD) is one of the serious microvascular complications of diabetes 1 . Considering the high incidence rate (40%) of new-onset CKD in patients with diabetes and the limited efficacy of renin angiotensin system inhibitors 2,3 , new therapeutic strategies are urgently needed to improve the renal prognosis of patients with diabetes. Tubular atrophy, interstitial fibrosis, and rarity of tubular capillaries are closely related to renal function decline in patients with DKD. The sodium-glucose cotransporter 2 inhibitors (SGLT2i) are new oral anti-diabetes agents, that can specifically bind to SGLT2 in the proximal tubular, reduce glucose reabsorption, promote urinary glucose excretion, and effectively reduce blood glucose 4 . Multicenter clinical trials have demonstrated that SGLT2i considerably hinders DKD progression (EMPA-REG 5 , CANVAS 6 , and Declare-TIMI 7 ). Previous studies have indicated that SGLT2i can protect the kidneys by regulating tubuloglomerular feedback. However, renal outcome trial (ROT) studies, including CREDENCE 8 and DAPA-CKD 9 , have confirmed that canagliflozin and dapagliflozin effectively reduce the risk of end stage renal disease (ESRD) and cardiovascular death in patients with or without DKD, suggesting the existence of other mechanisms. However, the exact molecular backgrounds remain unclear.
Circular RNAs (circRNAs) feature a stable structure formed by special loop splicing 10 . These circRNAs are highly conserved and abundantly expressed in mammalian tissues 11 . Due to the lack of terminal 5′-and 3′-ends, circRNAs exhibit higher stability and resistance against RNA exonucleases 12 . Recent studies have revealed that circRNAs regulate gene expression via multiple mechanisms, such as acting as microRNA (miRNA) sponges 13 , forming RNA-protein complexes 14 , and even being transcribed into proteins 15 . However, limited studies on circRNAs have focused on DKD pathogenesis. Particularly, the relationship of SGLT2i and circRNAs in DKD has not been studied.
In this study, circRNA and mRNA microarrays were performed in PTECs (HK-2 cell) treated with NG, HG, and SGLT2i (dapagliflozin) to systematically identify differentially expressed circRNAs and mRNAs. We identified, constructed, and functionally analyzed the circRNA-associated mRNA and competing endogenous RNA (ceRNA) networks in HK-2 cells. In addition, metformin and dipeptidyl peptidase-4 inhibitors (DPP4i), as classical anti-diabetes drugs, have been reported to have renal-protective effects, but the mechanisms remain unclear. To further clarify the difference and superiority of SGLT2i in renal protection, we also analyzed the circRNA and mRNA microarray in HK-2 cells treated with metformin and vildagliptin. Our study provides new insights and a theoretical basis for further exploration in the molecular mechanism of action of dapagliflozin in treating DKD.

Network pharmacology on the effect of dapagliflozin, metformin and vildagliptin on DKD
Using network pharmacology analysis, we predicted a total of 47, 37, and 18 targets for dapagliflozin, metformin, and vildagliptin on STITCH (http://stitch.embl.de/), respectively. GO term and KEGG pathway enrichment analyses showed that fatty acid degradation and metabolism, the PPAR signaling pathway, and ion transport participate in the biological processes of dapagliflozin on DKD. FOXO and AMPK signaling pathways participate in the biological processes of metformin, and insulin secretion is involved in the biological processes of vildagliptin ( Fig. 1a-f). Based on the DKD targets predicted by the DisGeNET (https://www.disgenet.org/) database and Venny (http://bioinfogp.cnb.csic.es/tools/venny/index.html) analysis, we obtained hub genes of different drugs on DKD using the cytoHubba app in Cytoscape ( Supplementary Fig. 1). KEGG analysis indicated that the primary mechanisms by which dapagliflozin improve DKD might be related to fatty acid degradation and metabolism and the PPAR signaling pathways, while FOXO, AMPK, neuroactive ligand-receptor interaction, cAMP, and insulin secretion signaling pathways play key roles in the mechanisms of renal protection by metformin and vildagliptin ( Fig. 1g-i).

CircRNA and mRNA expression profile in HK-2 cells treated with different drugs (dapagliflozin, metformin, vildagliptin)
To further clarify the precise mechanisms of action of different drugs on DKD, we conducted RNA microarray assays to understand the role of circRNA and mRNA profiles in HK-2 cells treated with NG, HG, dapagliflozin, metformin and vildagliptin. A total of 2,710 significantly differentially expressed circRNAs were identified (fold change ≥ 1.5, P < 0.05). Compared with the NG group, 119 circRNAs were upregulated and 434 circRNAs were downregulated in the HG group. A total of 1,406 circRNAs were upregulated, and 1,151 were downregulated (185 upregulated and 49 downregulated circRNAs in the dapagliflozin group, 901 upregulated and 736 downregulated circRNAs in the metformin group, and 320 upregulated and 366 downregulated circRNAs in the vildagliptin group) in the drug groups compared with the HG group ( Fig. 2a-d). Moreover, a heatmap of the significant DEcircRNAs was constructed ( Fig. 2e-h). To validate the microarray analysis results, six circRNAs were randomly chosen for validation using qRT-PCR. Five circRNAs were consistent with microarray data, including three upregulated circRNAs in HG group compared with the NG group, and two upregulated circRNAs in the metformin group compared with the HG group ( Fig. 2i and j). The qRT-PCR results supported the reliability of the microarray data.
To provide a comprehensive understanding of circRNAs, the parental linear transcripts for circRNAs were annotated and evaluated by GO and KEGG pathway analyses. The results showed that GO terms were related to gene expression and RNA binding in the HG group (Fig. 3a), gene expression and RNA binding in the dapagliflozin group (Fig. 3b), and protein modificationprocesses and gene expression in the metformin and vildagliptin groups, respectively ( Fig. 3c and d). KEGG pathway analyses showed that the cell cycle was the most involved in the dysregulated circRNAs in the HG group (Fig. 3e). However, in the mechanism action of dapagliflozin in protecting HK-2 cells, the role of circRNAs was closely related to the Rap1 signaling pathway (Fig. 3f), whereas metformin is associated with the MAPK and PI3K-Akt signaling pathways (Fig. 3g). Among the protective mechanisms of vildagliptin, the HIF-1 signaling pathway changed most significantly (Fig. 3h). Collectively, these biological processes and molecular functions may provide clues to further explore the molecular mechanisms of action of different drugs on DKD.

GO and KEGG Analyses of mRNAs co-regulated by dapagliflozin, metformin and vildagliptin
To further compare the mechanisms of different drugs, we used Venny analysis to investigate the circRNAs and mRNAs co-regulated by the three types of drugs. The results showed that dapagliflozin, metformin, and vildagliptin jointly regulated 54 circRNAs and 43 mRNAs, with 20 upregulated circRNAs and 38 upregulated mRNAs, and 34 downregulated circRNAs and five downregulated mRNAs ( Fig. 4a-d). Furthermore, GO and KEGG analyses were performed on these mRNAs. Cell communication, cell death, transcription factor, MAPK, p53, ErbB, and HIF-1 signaling pathways were closely related to the mechanisms of the three drugs in improving HK-2 cell injury induced by HG ( Fig.  4e and f). In addition, the hub genes of these key pathways were VAV1, CDKN1A and VEGF, which are related to the cell cycle, proliferation, differentiation, and glycolysis ( Fig.  4g). To some extent, these drugs may alleviate HK-2 cell injury by regulating similar pathways.

Exploration of relationships between the key circRNAs and their target genes
Next, we focused on dapagliflozin-related circRNAs. We chose seven significantly DEcircRNAs, including top four upregulated circRNAs and top three downregulated circRNAs in the HG group compared with NG group (Table 2). These circRNAs were validated by qRT-PCR and the results showed a similar trend to the microarray data ( Fig.  5a and c). Notably, there was one circRNAs that was induced by HG but repressed by dapagliflozin and two circRNAs that were repressed by HG but induced by dapagliflozin, indicating that dapagliflozin might alleviate HK-2 cell injury by regulating these circRNAs (Fig. 5e). To identify the functions of circRNAs, we constructed coding and noncoding coexpression (CNC) networks according to the degree of correlation of circRNA target genes ( Fig. 5b, d, f). GO and KEGG analyses showed that the target genes of upregulated circRNAs were mainly enriched in calcium ion transport, cell apoptosis, intercellular junctions, and mitophagy pathways in the HG group, whereas the target genes of downregulated circRNAs were mainly enriched in gene expression, oxoacid metabolic process, amino acid metabolism, and lysosome pathways. However, the target genes of dapagliflozin-related circRNAs primarily focused on gene expression, glucose metabolic process, extracellular exosome, epigenetics, amino acid metabolism, and lysosome pathways (Fig. 6), providing a theoretical basis for further exploring the mechanistic action of dapagliflozin.

Establishment of the PPI network and identification of hub genes
Using the STRING database (http://string-db.org/), we constructed protein-protein interaction (PPI) network of mRNAs involved in the ceRNA network ( Fig. 8a and e). Using the cytoHubba app in Cytoscape, we defined the hub genes of the hsa_circRNA_001586-miRNA network, including CCND1, CCND2, CDK6, CDK4, CDKN1A, CDK2, CCNA1, CCNE1, GSK3B, and SNAP25.which are involved in the cell cycle (Fig. 8b). The hub genes of hsa_circRNA_012448-miRNA were VEGFA, IGF1, KIT, AREG, PRF1, TNFRSF9, TNFRSF18, GSK3B, CDKN1A, and ERBB3, which are related to insulin signaling and PI3K-Akt signaling pathways (Fig. 8f). GO and KEGG analyses showed that the hsa_circRNA_001586-miRNA-mRNA axis was primarily involved in gene expression, cell adhesion and junctions, potassium channel regulator activity, amino acid metabolism, Oglycan biosynthesis, and vesicle transport ( Fig. 8c and d). However, the hsa_circRNA_012448-miRNA-mRNA axis was mainly involved in nitrogen metabolism, vesicle secretion, MAPK, ErbB, and insulin secretion signaling pathways ( Fig. 8g and h). Functional annotation showed that dapagliflozin participated in the pathophysiology of DKD through different circRNA-miRNA networks.

Discussion
Currently, the treatment of DKD lacks effective drugs. As a new anti-diabetes drug, SGLT2i (dapagliflozin) improves renal outcomes partly beyond its anti-hyperglycemic effect. However, the underlying mechanisms remain unclear. The recent rise of high throughput sequencing and microarray analysis have shed new light on studies of multiple diseases including DKD. Especially, some noncoding RNAs may be potential diagnosis biomarkers, and provide new perspectives to study the potential molecular mechanisms and regulatory targets of DKD. In this study, we used network pharmacology and RNA microarrays to detect the differential expression profiles of circRNA/mRNA in HK-2 cells treated with NG, HG, and different drugs, and explored the biological functions of circRNA/mRNA. We then predicted the different mechanism actions of dapagliflozin, metformin, and vildagliptin using GO and KEGG analyses. Furthermore, through the CNC networks and ceRNA networks, we identified the potential regulatory network of the circRNAs-miRNAs and signaling pathways of dapagliflozin in renal protection.
Compared with the HG group, we identified 234 DEcircRNAs in the dapagliflozin group, 1,637 in the metformin group, and 686 in the vildagliptin group, indicating that circRNAs play an important regulatory role in HK-2 cells treated with different drugs. Moreover, our data showed that linear transcripts from the same parent could express multiple circRNA subtypes in HK-2 cells, which is consistent with a previous study where some "hot spot" genes produced a variety of different circRNA subtypes 16 . Strikingly, although the predicted the targets functions of dapagliflozin and vildagliptin by network pharmacology were the lipid metabolism process and insulin secretion signaling pathways, GO and KEGG pathway analyses of the parental genes of DEcircRNAs showed that the pathways for dapagliflozin were involved with the Rap1, transcription regulation, and protein processing in endoplasmic reticulum pathways, while the pathways for vildagliptin were mainly the HIF-1, focal adhesion, and VEGF signaling pathways. Reportedly, defective fatty acids oxidation in PTECs is a key mechanism of kidney fibrosis development 17 , and high oxidative stress induced by fatty acids could trigger endoplasmic reticulum (ER) stress 18 . These studies further confirmed that the importance of lipid metabolism in dapagliflozin protecting PTECs. However, the main pathways for metformin were the FOXO and AMPK signaling pathways, consistent with the predicted results 19 . These results indicated that the role of dapagliflozin and vildagliptin may have tissue and cell specificity, and extend the molecular mechanisms in the pathophysiological process of HK-2 cells.
Generally, the treatment of diabetes and its complications involves a multi-drug combination. Some studies have explored the superiority of combination therapies (SGLT2i along with metformin/DDP-4i) in treatment-naïve patients with diabetes with combination therapies as an initial approach to observe targeting kidney outcomes. A previous study reported that combined SGLT2i (dapagliflozin) and DPP4i (saxagliptin) reduced the activation of the NLRP3 inflammasome and delayed the DKD progression in mice with type 2 diabetes 20 . In addition, a new report indicated that metformin promoted autophagy in the kidneys of individuals with chronic kidney disease, which may contribute to the actions of SGLT2i increasing autophagy and muting inflammation in diabetic kidneys 21 . Here, we performed functional analysis of circRNAs and mRNAs coregulated by three drugs. GO and KEGG analyses revealed several important pathways related to cell communication, cell death, transcription factor, MAPK, p53, ErbB, and HIF-1 signaling pathways, which are reported to participate in the pathophysiological process of DKD. In addition, we identified the hub genes co-regulated by three drugs. Notably, t has been demonstrated that the hub genes VAV1 22 and CDKN1A 23 are involved in cell proliferation, differentiation and cell cycle regulation. VEGFA plays a key role in the progression of renal fibrosis 24 . DDIT3 is involved in ischemic injury and oxidative stress responses 25 . These findings suggested that the combined use of the three drugs may improve renal tubular injuries by regulating PTEC proliferation, apoptosis, and fibrosis, providing clues for further exploration of the mechanisms of multi-drug combination in vivo.
Most circRNAs are epigenetic regulatory factors whose expression levels and functions are significantly related to the expression of protein-coding genes. Currently, studies in vivo and in vitro have indicated that SGLT2i is involved in DKD tubular cell apoptosis, tubulointerstitial inflammation and fibrosis. Our previous study also has reported that dapagliflozin improves tubulointerstitial fibrosis in DKD. Another study indicated that molecular pathways targeted by dapagliflozin and associated with DKD were related to energy metabolism, mitochondrial and endothelial functions 26 . Based on the CNC network constructed by the most significant DEcircRNAs, we revealed that the target genes of circRNAs upregulated by HG were associated with calcium ion transport, cell apoptosis, intercellular junctions, and mitophagy pathways. In contrast, the target genes of circRNAs downregulated by HG were associated with gene expression, oxoacid metabolic process, amino acid metabolism, and lysosome pathways. However, these target genes of DEcircRNAs induced by dapagliflozin were most enriched in gene expression, glucose metabolic process, extracellular exosome, epigenetics, and amino acid metabolism pathways. Strikingly, extracellular exosome signaling provides us with a new hypothesis, that is, whether the exosomes secreted by renal tubules induced by dapagliflozin affect the glomerulus, thereby improving the glomerular function and restoring the renal function. Collectively, these studies indicated that dapagliflozin protects DKD through multiple signaling pathways are feasible. In particularly, some signaling pathways have not been previously reported in the field of DKD and may be the starting point for future studies.
We selected hsa_circRNA_001586 and hsa_circRNA_012448 to construct ceRNA network and revealed a new interaction between circRNA, miRNA and mRNA related to dapagliflozin. The hub miRNAs found in this study include hsa-miR-4739, hsa-miR-7851-3p, hsa-miR-1273g-3p, hsa-miR-378g and hsa-miR-29b-2-5p. Previous studies have reported hsa-miR-4739 is elevated in the urinary exosomes of patients with DKD 27 . Hsa-miR-4739 in plasma can be used as a potential biomarker for severe limb ischemia in diabetic patients 28 . In addition, hsa-miR-378p not only participates in fat metabolism in obese mice, but also alleviates renal tubulointerstitial fibrosis 29 . Thus, these hsa_circRNA_001586, hsa_circRNA_012448 and miRNAs may be prospective targets for the prevention or treatment of dapagliflozin in diabetic patients. The KEGG pathway analysis showed that hsa_circRNA_001586-miRNAs network induced by dapagliflozin were functionally involved in the cell adhesion and junctions, potassium channel regulator activity, amino acid metabolism, O-glycan biosynthesis, and vesicle transport. While the hsa_circRNA_012448-miRNAs networks were most enriched nitrogen metabolism, vesicle secretion, MAPK, ErbB, and insulin secretion signaling pathways. At present, newer strategies based on dapagliflozin are still in the experimental stages and limited. Thus, these ceRNA networks may be new biomarkers and therapeutic targets for DKD treatment. These mechanisms need to be further verified in vivo and in vitro experiments as well as large-scale clinical trials in the future.
Previous studies based on RNA microarray analysis have mostly analyzed noncoding RNA profiling in DKD, but not conducted between anti-diabetes drugs, circRNAs, and signaling pathways. In this study, we made a perfect combination of circRNAs-genessignaling pathways and identified several important circRNAs-miRNAs networks in DKD such as hsa_circ_1586-miR-4739/hsa-miR-7851-3p/miR-1273g-3p networks, and has_circ_012448-hsa-miR-378g/hsa-miR-29b-5p networks. Furthermore, a few novel signaling pathways were found to be involved in renal protection of dapagliflozin. However, we only identified candidate aberrantly circRNAs and predicted some potential miRNAs-genes regulatory network in HK-2 cells based on bioinformatics analysis. Further molecular biological experiments are necessary to confirm the function of the identified cirRNAs as well as the interaction between dapagliflozin and circRNAs in DKD. In future studies, we also need to analyze more clinical samples to validate these results.

CircRNA and mRNA microarrays
Three pairs of HK-2 cells from each group were used for circRNA and mRNA microarray analysis by KangChen Biotech Inc (Shanghai, China). Sample labeling and microarray hybridization were conducted using Arraystar Human circRNA Array v2.0 (Arraystar). Briefly, total RNA was digested with RNase R (Epicenter, Madison, WI, USA). Then, linear RNAs were removed and circular RNA were concentrated. Subsequently, the concentrated circRNAs were amplified and transcribed into fluorescent RNA using the Arraystar Super RNA Labeling Kit using the random priming method. Thereafter, the labeled cRNAs were hybridized, washed, and scanned on Arraystar Human circRNA Array v2.0 (Arraystar), using an Agilent Scanner G2505C (Agilent Technologies, Santa Clara, CA, USA). Differentially expressed circRNAs were identified based on a P-value < 0.05 with a fold change > 1.5. For mRNAs, we used Arraystar Human mRNA arrays (Arraystar) following the manufacturer's instructions and quantile standardization was performed on the original data using GeneSpring GX v12.1 software (Agilent Technologies). Differentially expressed mRNAs were identified based on a P-value < 0.05 with a fold change > 2.0. All raw data were deposited in the Gene Expression Omnibus under the accession code GSE179226 and GSE158534.

RNA isolation and Validation by qRT-PCR
Total RNA was isolated from HK-2 cells using TRIzol reagent (Takara Bio Inc, Japan). The concentrations were assessed using a NanoDrop ND-1000 spectrophotometer (Thermo Scientific, Wilmington, DE, USA). Total RNA was reverse transcribed to complementary DNA using a reverse transcription kit (Sangon Biotech, Shanghai, China). CircRNA expression were assessed by qRT-PCR using SGExcel FastSYBR Mixture (with ROX) (Sangon Biotech) in an ABI 7500 system (Applied Biosystems, CA, USA). Ten predicted miRNAs were reverse transcribed into cDNA using miRNA First Strand cDNA Synthesis (Tailing Reaction) (Sangon Biotech). miRNA expression was assessed by qRT-PCR using a MicroRNA qPCR Kit (Sangon Biotech). β-actin and U6 were served as the internal normalization control using the 2 -ΔΔCt cycle threshold. The primer sequences are listed in Table 1.

Construction of coding and noncoding co-expression analysis
The co-expression analysis was based on calculating the Pearson correlation coefficient between coding and noncoding genes according to their expression levels. The absolute value of PCC ≥ 0.95, P-value < 0.05 was recommended and retained for further analysis.

Construction of PPI network and identification of hub genes
The protein-protein interaction (PPI) network of the DEmRNAs was established using the STRING database 33 (http://string-db.org/), and then visualized using Cytoscape 3.7.1 software. Subsequently, the cytoHubba app in Cytoscape was used to determine the hub genes.

GO and KEGG pathway analyses of mRNAs in the network
To understand functional enrichment, GOs and KEGG pathway analyses of DEmRNAs were performed using the Database for Annotation, Visualization, and Integration Discovery (DAVID) 34 . The GO terms and pathways were considered statistically significant at P < 0.05.

Statistical Analysis
Comparisons between two groups were analyzed with Student's t-test and multiple comparisons were calculated using one-way analysis of variance (ANOVA) followed by Bonferroni's multiple comparisons test using GraphPad Prism software (version 9.0). Data were presented as the mean ± S.D. The statistical significance was defined as P-value < 0.05.

Acknowledgements
This work was supported by the National Natural Science Foundation of China (grant numbers: 81770812 and 81974110).

Author Contributions
Guijun Qin and Yi Song concepted and designed the study. Yi Song and Feng Guo administrated the project and collected data. Yi Song, Feng Guo, and Yifan Liu validated and analyzed these data. Yanyan Zhao, Fengjuan Huang and Lin Zhao analyzed bioinformatics. Guijun Qin, Yi Song write this manuscript. Guijun Qin, Feng Guo, XunjieFan and Yi Song modified and discussed this manuscript.