Identification of the circRNA-miRNA-mRNA Prognostic Regulatory Network in Lung Adenocarcinoma

Background: Numerous studies have identified that circular RNAs (circRNAs) can serve as competing endogenous RNAs (ceRNAs) to regulate tumor progression. However, there are still a large number of circRNAs to be deciphered. Objective: The purpose of this study was to reveal novel circRNAs and their potential role in lung adenocarcinoma (LUAD). Methods: To unveil LUAD-related circRNAs, microRNA (miRNAs), and messenger RNA (mRNA) and elucidate their possible molecular mechanisms, we employed a strategy combining extensive data mining and bioinformatics methods. According to the results of bioinformatics workflow analysis, a novel circRNA-miRNA-mRNA network was constructed. Results: Ten circRNAs with different expressions were acquired from four Gene Expression Omnibus (GEO) microarray datasets. Seven Prognostic-related differential miRNAs of LUAD were gained from The Cancer Genome Atlas (TCGA). Simultaneously, the miRNA reaction components corresponding to the ten circRNAs were predicted. Two circRNA–miRNA interactions including two circRNAs (hsa_circ_0008234 and hsa_circ_0002360) and two miRNAs (hsa-miR-490-3p and hsa-miR-1293) were identified above. Then, target genes of the two miRNAs and differently expressed genes (DEGs) from TCGA on LUAD were collected. Three hub-genes (ADCY9, NMUR1, SYT1) were determined according to prognosis in patients with LUAD ulteriorly. Conclusions: hsa_circ_0008234/hsa-miR-490-3p/SYT1 and hsa_circ_0002360/hsa-miR-1293/ (ADCY9, NMUR1) networks were established, and identified molecules may be involved in pathogenesis and prognosis in patients with LUAD.


Introduction
Lung cancer is the first and most frequent cause of cancer-related death worldwide, partly due to its poor response to existing conventional target treatments or immune therapies [1]. Thus, the 5-year survival rate of patients and median overall survival (OS) remain poor depending on the different TNM stage of the tumor. Lung cancer is typically classified into about 85% non-small cell lung cancer (NSCLC) and about 15% small cell lung cancer (SCLC). In these subtypes, NSCLC were mainly refers to lung adenocarcinoma (LUAD) and squamous cell carcinoma (LUSC) [2]. Surgical excision is still the most effective treatment strategy for patients in the early stage. Combined with radiotherapy and chemotherapy, the postoperative survival rate of lung cancer patients can be significantly improved [3]. However, clinically successful treatment therapies are not very promising due to the lack of early diagnosis, postoperative tumor relapse, metastasis, and the development of drug resistance. Chemotherapy and targeted biological treatment were still the best for advanced lung cancer patients, especially LUAD patients. Therefore, exploring the genetic and epigenetic regulation of LUAD and finding new tumor molecular markers and therapeutic drugs play an important role in ameliorating the therapeutic effect and diminishing recurrence rate and mortality rate.
Circular RNAs (circRNAs) are a kind of endogenous non-coding RNA, which originates from the exon or intron region of a gene and forms a covalent closed-loop structure The microarray datasets providing circRNA expression profiles for human samples in lung cancer were obtained from the GEO database. We selected GSE101684, GSE 101586, GSE112214, and GSE158695. The GSE101684 dataset contains four LUAD tissues derived from the inferior lobe of the right lung and four para-cancerous tissues. The GSE 101586 has five LUAD tissues and five normal lung tissues. The GSE 112214 and GSE 158695 include three NSCLC samples and three normal lung tissues. The details of four circRNA microarray data are shown in Table 1. MiRNA and mRNA expression profile data, and clinical information were obtained from TCGA, including 521 LUAD tissue samples and 46 adjacent normal tissues. LUAD tissue with missing clinical information were excluded, and 500 LUAD samples were shortlisted for further analysis.

Differential circRNA Expression and Venn Analysis
We used GEO2R online analysis tool on the NCBI website and Venn diagram software to obtain the mutual DECs in the four GEO datasets above. The criteria for the significant differentially expressed were |log2(FC)| > 1 and adjusted p < 0.05. In this study, we selected circRNAs expression profiles that were co-highly expressed in three or four GEO databases by Venn analysis. Meanwhile, the commonly significant lower expression circRNAs in the three GEO datasets were obtained and the simultaneously lower expressed were screened out through the Venn analysis.

Definition of the miRNAs-Related Prognostic Model
To verify the association between the survival time and target miRNAs from the TCGA database, we carried out the Univariate Cox, LASSO, and multivariate Cox regression analyses via the glmnet (version 2.0.18) and survival (version 2.44.1.1) R packages. First, univariate Cox regression was performed to determine the correlation between miRNAs expression levels and overall survival. A p value less than 0.05 is considered a criterion, which indicates that the difference was statistically significant. Then, the LASSO regression analysis was performed on these miRNAs that met this criterion. In this way, we can further screen for candidate prognostics miRNAs, the minimum λ value was used as the inclusion criteria and can represent this model's most appropriate number of variables. The miRNAs screened by the above methods were further analyzed by multivariate Cox regression to evaluate the independent contribution of each miRNA to prognosis and hazard ratios (HR), and 95% confidence intervals (CI) for the key miRNAs were calculated. According to the formula: Σ (expmiRNAn × βmiRNAn), we used the Cox regression coefficient (β) and expression levels of miRNAs to calculate prognosis-related risk scores. The LUAD samples were divided into high-risk and low-risk groups based on risk score.

Evaluation of the Risk Model and Identification of Potential Biomarkers
Time-dependent receiver operator characteristic (ROC) curves were established, and the area under the curve (AUC) was calculated to assess the performance of the prognostic model above. DEMs survival times were plotted by K m survival analyses. Multivariate Cox regression were performed to ensure which miRNAs differentially expressed could serve as potential prognostic molecular markers. The above indicators with statistically significant differences were p values < 0.05.

Prediction of miRNA Binding Sites and Target mRNAs
The potential target miRNA of the circRNAs were predicted using the circbank (http: //www.circbank.cn/searchMiRNA.html accessed on 1 May 2021). At the same time, we selected the DEMs through the above prognosis model analysis. Next, we conducted a Venn diagram analysis; overlapping predicted miRNAs were chosen for further investigation and research. Furthermore, we used miRwalk (http://mirwalk.umm.uni-heidelberg.de/ accessed on 1 May 2021) to predicate target mRNAs of these selected miRNAs, which can bind to the target mRNA 3 -UTR. The differentially expressed genes (DEGs) in LUAD were gained from the TCGA database. Similarly, overlapping predicated mRNAs were selected for the following bioinformatics analysis.

GO and KEGG Pathway Analyses
We used DAVID (https://david.ncifcrf.gov/ accessed on 1 May 2021) bioinformatics resources to carry out GO, and KEGG pathway analysis of the overlapping predicted differential expressed mRNAs. The criteria for the statistically significant differential expressed is a p value < 0.05.

Establishment of PPI Network and Identification of Hub-Genes
We constructed a PPI network for targeted DEGs via using the STRING (https://cn. string-db.org/ accessed on 1 May 2021) online website. Then, PPI network was embellished through Cytoscape software. The Cytoscape plug-in cytoHubba was used to screen hub genes with the node degree. The Cytoscape plug-in MCODE to screen out the essential modules in the PPI network with the selected criteria as follows: degree cut-off = 2, node score cut-off = 0.2, Max depth = 100, and k-score = 2.

Differentially Expressed and OS Analysis of the Hub Genes
The differential expression data of these hub genes in LUAD tissues and normal tissues were gained using the TCGA. The effect of these hub genes on OS rate in lung cancer was identified using the Gene Expression Profiling Interactive Analysis (GEPIA, http://gepia2.cancer-pku.cn/#index accessed on 3 May 2021), a web-based tool derived from the TCGA.

Correlation Analysis of Hub Gene and Identification of the Association between Hub Gene and Immune Infiltration in LUAD
Based on the GEPIA, we determined the correlation of the hub genes mentioned above with each other. Meanwhile, the Tumor Immune Estimation Resource (TIMER, http://timer.cistrome.org/ accessed on 3 May 2021) website was used to analyze the correlation between the expression of hub genes and the infiltration of several immune cells in LUAD.

Differentially Expressed circRNAs Analysis in the NSCLC
With a p value < 0.05, |log2(FC)| > 2 as the standard of statistical difference, we selected a total of 813 DECs from four microarray datasets (GSE101684, GSE101586, GSE112214, GSE158695). The principal information for the four datasets is presented in Table 1. A total of 411 DECs were identified in the GSE101684 dataset, 175 circRNAs were upregulated and 236 circRNAs were downregulated ( Figure 1A); 68 DECs with 21 were upregulated and 47 were downregulated and identified from the GSE101586 dataset ( Figure 1B); 149 DECs were determined in the GSE112214 dataset, 133 of which were upregulated and 16 were downregulated ( Figure 1C); 101 upregulated circRNAs and 84 downregulated circRNAs were found derived from the GSE158695 dataset ( Figure 1D). Then, we used the Venn diagram analysis tool to identify the overlapped circRNAs from the four microarray datasets. The upregulated circRNAs shared in three or four data sets were selected for further analysis (hsa_circ_103415, hsa_circ_101066, hsa_circ_104513, hsa_circ_103134, hsat_circ_100395, hsa_circ_101213, hsa_circ_102046, hsa_circ_102442) ( Figure 1E). Similarly, circRNAs downregulated in all three data sets were selected for further analysis (has_circ_002178, has_circ_103123) ( Figure 1F). The details concerning the selected circR-NAs were listed in Table 2.

Prognostic-Related Differential miRNAs from TCGA
Based on the miRNA expression profile obtained from the TCGA database and the cut-off criterion of |log2(FC)| > 2, p < 0.05, we screened out 298 significantly DEMs in the LUAD patient sample (Figure 2A). Among these selected miRNAs, including 241

Prognostic-Related Differential miRNAs from TCGA
Based on the miRNA expression profile obtained from the TCGA database and the cutoff criterion of |log2(FC)| > 2, p < 0.05, we screened out 298 significantly DEMs in the LUAD patient sample (Figure 2A). Among these selected miRNAs, including 241 upregulated and 57 downregulated miRNAs, these miRNAs were selected for predicting prognosis analysis. Firstly, we used the univariate Cox regression to evaluate the correlation between DEMs and patients' OS. 41 miRNAs have an obvious correlation with OS of patients. Then, LASSO regression analysis was performed on these miRNAs. Ultimately, 25 miRNAs were determined for further analysis when the λ value was at a minimum ( Figure 2B,C). At the same time, we performed a multivariate Cox regression analysis to acquire HR values and 95% CIs for selected miRNA; the results showed that seven prognostic-related differential miRNAs of LUAD were obtained (hsa-miR-1293, hsa-miR-142-3p, hsa-miR-490-3p, hsa-miR-543, hsa-miR-548au-5p, hsa-miR-548v, hsa-miR-5571) ( Figure 2D). All specimens were divided into two groups: high-risk group and low-risk group ( Figure 2E). We drew a time-dependent ROC curve according to the above median risk score values. The AUCs for 3-and 5-year survival were 0.754 and 0.784, respectively, demonstrating that the risk score model has a stable performance ( Figure 2F).

Construction of circRNA-miRNA-mRNA Network
Firstly, a total of 5388 DEGs were obtained from the TCGA database, including 3852 upregulated and 1721 downregulated genes ( Figure 4E). Additionally, 714 target mRNAs of hsa-miR-490-3p and 884 target mRNAs of hsa-miR-1293 were obtained from the miR-Walk. We identified 96 upregulated and 89 downregulated target genes that play an essential role in LUAD by cross-predicting the target gene and DEGs, respectively ( Figure  4F,G). Subsequently, the circRNA-miRNA interplay, and miRNA-mRNA interplay were integrated to establish a circRNA-miRNA-mRNA network ( Figure 5A), which primarily hsa-miR-490-3p and 884 target mRNAs of hsa-miR-1293 were obtained from the miRWalk. We identified 96 upregulated and 89 downregulated target genes that play an essential role in LUAD by cross-predicting the target gene and DEGs, respectively ( Figure 4F,G). Subsequently, the circRNA-miRNA interplay, and miRNA-mRNA interplay were integrated to establish a circRNA-miRNA-mRNA network ( Figure 5A), which primarily shed light on the correlation between the DECs (hsa_circ_0008234 and hsa_circ_0002360), miRNAs (hsa-miR-490-3p and hsa-miR-1293) and the 185 mRNAs. The PPI network was constructed using STRING and Cytoscape to gain insight into the biological interactions of 185 genes and to define the circRNA-miRNA-mRNA regulatory network ( Figure 5A). Then, we used the plug-in of CytoHubba in Cytoscape to screen out the top 10 node degrees to represent the central genes of the PPI network, including NMUR, ADCY9, CXCL16, SSTR1, S1PR1, OPRK1, ELAVL3, SYT1, CNTN2 and KCNC1( Figure 5B).

Functional Enrichment Analysis
GO and KEGG analyses were performed by DAVID online database to evaluate the function of the target mRNAs of these two DEMs in the subnetwork. The top six enriched terms under the three common categories: biological process (BP), molecular function (MF), and cellular component (CC), were shown ( Figure 6A). Results from the first six BPenriched terms revealed that the target mRNAs were primarily associated with "calcium icon-regulated exocytosis of neurotransmitters," "positive regulation of gene expression," "brain development," "neuron differentiation," "cell adhesion" and "positive regulation of potassium ion transmembrane transport." The top six KEGG pathways showed the target mRNAs most likely to be involved in the "cGMP-PKG singling pathway", which provides more evidence in support of this assumption that DECs may play vital roles in cancer initiation and progression, and metastasis.

Functional Enrichment Analysis
GO and KEGG analyses were performed by DAVID online database to evaluate the function of the target mRNAs of these two DEMs in the subnetwork. The top six enriched terms under the three common categories: biological process (BP), molecular function (MF), and cellular component (CC), were shown ( Figure 6A). Results from the first six BP-enriched terms revealed that the target mRNAs were primarily associated with "calcium icon-regulated exocytosis of neurotransmitters", "positive regulation of gene expression", "brain development", "neuron differentiation", "cell adhesion" and "positive regulation of potassium ion transmembrane transport". The top six KEGG pathways showed the target mRNAs most likely to be involved in the "cGMP-PKG singling pathway", which provides more evidence in support of this assumption that DECs may play vital roles in cancer initiation and progression, and metastasis.

Identification of Hub Genes
The relative expression level of the top ten hub genes and the effect on OS were obtained from the TCGA. Only three hub genes (ADCY9, NMUR1, SYT1) are related to the overall survival of patients. Specifically, LUAD patients with overexpressed SYT1 ( Figure  7A,D) have a low OS rate, which indicates it can be a candidate for poor prognostics. Furthermore, the downregulated hub genes ADCY9 (Figure 7B,E) and NMUR1 ( Figure 7C,F) in the LUAD negatively affected the OS rate, which may be a favorable prognostics factor. Moreover, using the GEPIA, we examined the correlation of three differentially expressed hub genes and found that the expressions of ADCY9, NMUR1 and SYT1 were strongly associated with each other in LUAD ( Figure 7G-I).

Identification of Hub Genes
The relative expression level of the top ten hub genes and the effect on OS were obtained from the TCGA. Only three hub genes (ADCY9, NMUR1, SYT1) are related to the overall survival of patients. Specifically, LUAD patients with overexpressed SYT1 (Figure 7A,D) have a low OS rate, which indicates it can be a candidate for poor prognostics. Furthermore, the downregulated hub genes ADCY9 (Figure 7B,E) and NMUR1 ( Figure 7C,F) in the LUAD negatively affected the OS rate, which may be a favorable prognostics factor. Moreover, using the GEPIA, we examined the correlation of three differentially expressed hub genes and found that the expressions of ADCY9, NMUR1 and SYT1 were strongly associated with each other in LUAD ( Figure 7G-I).

The Hub Genes Expression Was Correlated yo Immune Cell Infiltration in LUAD
The analysis of the TIMER website indicated that the expression of three hub genes (ADCY9, NMUR1, SYT1) was significantly associated with the infiltration of several immune cells in LUAD, including B cells, CD4+ T cells, CD8+ T cells, neutrophils, macrophages, and DCs. The result expounded that the expression of ADCY9 ( Figure 8A) and NMUR1 ( Figure 8B) were positively correlated with B cells, CD4+ T cells, CD8+ T cells, neutrophils, macrophages, and DCs. Thus, we speculated that the high expression of ADCY9 and NMUR1 might facilitate the anti-tumor immunity in the LUAD microenvironment. This result also proved that patients with high expression of ADCY and NMUR1 genes had a high survival rate. On the contrary, we found no significant correlation between the SYT1 ( Figure 8C) gene expression and immune cells infiltration. Therefore, SYTI may be a proto-oncogene that may inhibit the anti-tumor immune response and is negatively correlated with patient survival time.

The Hub Genes Expression Was Correlated yo Immune Cell Infiltration in LUAD
The analysis of the TIMER website indicated that the expression of three hub genes (ADCY9, NMUR1, SYT1) was significantly associated with the infiltration of several immune cells in LUAD, including B cells, CD4+ T cells, CD8+ T cells, neutrophils, macrophages, and DCs. The result expounded that the expression of ADCY9 ( Figure 8A) and NMUR1 ( Figure 8B) were positively correlated with B cells, CD4+ T cells, CD8+ T cells, neutrophils, macrophages, and DCs. Thus, we speculated that the high expression of ADCY9 and NMUR1 might facilitate the anti-tumor immunity in the LUAD microenvironment. This result also proved that patients with high expression of ADCY and NMUR1 genes had a high survival rate. On the contrary, we found no significant correlation between the SYT1 ( Figure 8C) gene expression and immune cells infiltration. Therefore, SYTI may be a proto-oncogene that may inhibit the anti-tumor immune response and is negatively correlated with patient survival time.

Discussion
Human genes can be divided into protein-coding genes and non-coding genes. Although more than 90% of human genes can be actively transcribed, the known proteincoding genes only account for about 1.1% of the human genome, while the majority of the gene transcription variants may be classified as ncRNAs [12,13]. In the past few years, many ncRNAs have been identified in human tissues, body fluids and cells to benefit from the development of high-throughput sequencing and genomic analysis platforms. Numerous studies have suggested that these ncRNAs have various biological functions and possibly participate in multiple metabolic pathways [14]. Compared with known ncRNAs, circRNA is a novel research hotspot and has gained popularity among researchers. In addition, due to the biological properties of circRNAs being covalently closed loops, circRNAs can be more stable existing in tissue, cells, or plasma than other ncRNAs [15]. It has been reported that circRNAs play various essential roles in biological functions, including microRNA sponges, RBP-binding molecules, transcriptional regulation factors, or protein translation templates [16,17]. Recent studies have indicated that a large number of circRNAs are differentially expressed in specific tumor cells and tissues, which may be involved in occurrence and development of tumors. In these studies, one of the biological functions of circRNA that has been studied extensively is that it can act as a sponge of miRNA to regulate gene expression to promote or inhibit the occurrence and progression of tumors [18]. For example, the upregulated circRNA circSEPT9 in TNBC could release the inhibition of leukocyte inhibitor factors by downregulating mir-637 and accelerate the carcinogenesis, development and metastasis of triple-negative breast cancer [19]. Furthermore, some circRNAs may serve as a tumor suppressor to inhibit tumor growth, invasion, and metastasis, so the downregulation of circRNAs may exert opposite effects. For example, circ-HuR was determined to be a tumor suppressor circRNA, can suppress HuR expression and gastric cancer progression, and can be a candidate for a potential therapeutic target for gastric cancer [20]. Circular RNA hsa_circ_0000326, a novel circRNA found in

Discussion
Human genes can be divided into protein-coding genes and non-coding genes. Although more than 90% of human genes can be actively transcribed, the known proteincoding genes only account for about 1.1% of the human genome, while the majority of the gene transcription variants may be classified as ncRNAs [12,13]. In the past few years, many ncRNAs have been identified in human tissues, body fluids and cells to benefit from the development of high-throughput sequencing and genomic analysis platforms. Numerous studies have suggested that these ncRNAs have various biological functions and possibly participate in multiple metabolic pathways [14]. Compared with known ncRNAs, circRNA is a novel research hotspot and has gained popularity among researchers. In addition, due to the biological properties of circRNAs being covalently closed loops, circRNAs can be more stable existing in tissue, cells, or plasma than other ncRNAs [15]. It has been reported that circRNAs play various essential roles in biological functions, including microRNA sponges, RBP-binding molecules, transcriptional regulation factors, or protein translation templates [16,17]. Recent studies have indicated that a large number of circRNAs are differentially expressed in specific tumor cells and tissues, which may be involved in occurrence and development of tumors. In these studies, one of the biological functions of circRNA that has been studied extensively is that it can act as a sponge of miRNA to regulate gene expression to promote or inhibit the occurrence and progression of tumors [18]. For example, the upregulated circRNA circSEPT9 in TNBC could release the inhibition of leukocyte inhibitor factors by downregulating mir-637 and accelerate the carcinogenesis, development and metastasis of triple-negative breast cancer [19]. Furthermore, some circRNAs may serve as a tumor suppressor to inhibit tumor growth, invasion, and metastasis, so the downregulation of circRNAs may exert opposite effects. For example, circ-HuR was determined to be a tumor suppressor circRNA, can suppress HuR expression and gastric cancer progression, and can be a candidate for a potential therapeutic target for gastric cancer [20]. Circular RNA hsa_circ_0000326, a novel circRNA found in lung cancer patient microarray, can act as a miR-338-3p sponge and alter the function of miR-338-3p to facilitate LUAD proliferation and metastasis [21].
In this study, we investigated the circRNA expression profile in LUAD tissues and para-cancerous tissues from four GEO databases (GSE101684, GSE112214, GSE 101586, and GSE158695). We focused on differentially expressed circRNAs in lung cancer tissues, which were remarkably upregulated or downregulated in cancer tissues, and may be significantly related to tumor occurrence, development, invasion, and metastasis. Here, we selected eight upregulated circRNAs (hsa_circRNA_103415, hsa_circRNA_101066, hsa_circRNA_104513, hsa_circRNA_103134, hsa_circRNA_100395, hsa_circRNA_101213, hsa_circRNA_102046, and hsa_circRNA_102442) and two downregulated circRNAs (hsa_circRNA_002178 and hsa_circRNA_103123) for further analysis. As highly conserved covalently closed RNAs, these selected circRNAs may contain abundant miRNA binding sites, indicating that they can serve as a sponge to adsorb the corresponding miRNA and thus function as ceRNAs to regulate correlated protein gene expression. Herein, we found the possible miRNAs binding sites of these ten dysregulated expressed circRNAs by circBank. At the same time, according to the data derived from the TCGA database, survival analysis and prognosis assessment of DEMs in LUAD were performed to construct a prognosis model and identify a potential miRNA biomarker. Finally, we determined hsa-mir-490-3p and hsa-mir-1293 in overlapped miRNAs for further analysis. Our prognosis model proved that the upregulated mir-490-3p could be associated with a good prognosis in LUAD. Previous studies have revealed that mir-490-3p can serve as tumor suppressor miRNA, and it is related to the occurrence and progression of various tumors. Zhiyong et al. demonstrated that miR-490-3p inhibited the Malignant Progression of LUAD cells by downregulating the Wnt/β-catenin signaling pathway [22]. Kang et al. reported that mir-490-3p was downregulated and mir-490-3p can promote cell proliferation, metastasis, and invasion via targeting HMGA2 [23]. In contrast, mir-1293 is considered a marker of poor prognosis and adverse to patient survival time. Although, Takagawa et al. observed that mir-1293 can suppress tumor progression by inhibiting DNA repair pathways and suggested that miR-1293 is a candidate for developing miRNA-based cancer therapeutics [24]. The molecular mechanism of miR-1293 in LUAD has not been reported yet. Next, the target genes of mir-490-3p and mir-1293 were identified through the miRwalk and TCGA database. Functionally, GO and KEGG pathway analysis was conducted for the target mRNAs to further investigate the underlying biological pathway of these target genes. A total of 185 target mRNAs were enriched in BP, MF, and CC terms, including "positive regulation of gene expression", "cell junction", "cell adhesion", and "RNA polymerase II core promoter proximal region sequence-specific binding", The results of KEGG pathway analysis indicated that "cGMP-PKG signaling pathway" was significantly enriched. Subsequently, we constructed a PPI network, screening ten hub genes from the PPI network. We investigated the expression of these genes in lung cancer tissues from the TCGA database and analyzed the association between both hub genes and patient survival time. The results showed that the downregulation of ADCY9 and NMUR1 in LUAD can inhibit the growth and aggressiveness of cancer, while the upregulation of SYT1 had the contrary effects. Additionally, the expressions of ADCY9 and NMUR1 in LUAD showed a positive correlation, while the expressions of SYT1 and ADCY9, SYT1 and NMUR1 showed a negative correlation. As for the functions of the selected hub genes, mutations in SYT1, the master switch responsible for allowing the human brain to release neurotransmitters, could lead to a rare neurodevelopmental disorder [25]. Recent studies have shown that SYT1 is also involved in cancer regulation [26,27]. ADCY9 is a widely distributed adenylate cyclase that catalyzes formation of cyclic AMP from ATP [28]. In addition, studies have demonstrated that ADCY9 takes to participate in the regulation of cellular function in certain cancers and drug responses [29,30]. NMUR1 is widely distributed in various organs of the human body and participates in the regulation of activation of phospholipase C activity, chloride transport, second messenger-mediated signaling and other physiological functions, and also plays an important role in the occurrence and development of various cancers [31]. However, there is little research on these hub genes in LUAD. Further experiments are needed to prove that hub genes participate in tumorigenesis. The microenvironment of a tumor consists of heterogeneous populations, including cancer cells themselves, infiltrating immune cells, and stromal cells; tumor-infiltrating immune cells play a vital role in tumorigenesis, progression, and metastasis, which not only inhibit tumor progression by attacking and killing cancer cells but also promote tumor progression by changing the tumor microenvironment [32,33]. With the development and application of new technologies such as single-cell RNA sequencing and mass cytometry, emerging evidence has shown that immune cells within the tumor microenvironment may possess various functions besides conventional tumor-antagonizing functions [34,35]. Hence, we assess the association between three hub genes and immune cells, respectively. The results demonstrated the positive correlation between the expression of ADCY9 and NMUR1 and immune cells. Therefore, the dysregulated expression and functional impairment of ADCY and NMUR1 could affect the anti-tumor immune reaction of the body's immune system. However, the SYT1 gene showed the opposite result to the two genes mentioned above. Based on the above analysis results, we speculated that these hub genes might become new targets for immunotherapy or potential tumor immunomarkers.

Conclusions
In this study, a regulatory network of circRNA-miRNA-mRNA was successfully constructed. Our work demonstrated that some novel circRNAs, miRNAs, and hub genes may have clinical application value as prognostic markers or novel biomarkers in LUAD. Although the current research has certain limitations, the combination of molecular signaling pathway mechanism research and bioinformatics analysis to identify potential cancer diagnostic biomarkers and therapeutic targets still has extensive advantages. Thus, the above evidence indicated that regulatory network analysis of circRNAs could be a powerful tool to explore the molecular mechanism of circRANAs, miRNAs, and hub genes in LUAD development. Additionally, we found that the hub gene may affect the immune response of the body, and it can become a candidate gene for immunotherapy targets or potential tumor immune markers. Nevertheless, further validation of the biological function and mechanism of these circRNAs, miRNAs, and hub genes need to be carried out to assess whether they can serve as novel biomarkers or therapeutic targets in LUAD patients to improve the diagnosis and treatment of LUAD.