Identification of an Immune-Related Biomarker Model Based on the CircRNA-Associated Regulatory Network for Esophageal Carcinoma

Esophageal carcinoma (ESCA) is one of the most frequent types of malignant tumor that has a dismal prognosis. This research applied datasets aimed from the GEO and TCGA to create a prognostic signature for forecasting the clinical outcome of ESCA patients on the basis of a circRNA-associated regulatory network. Methods. A regulatory network associated with ESCA was established based on transcriptome data of circRNAs, miRNAs, and mRNAs. Functional annotations were implemented to further explore the mechanism of ESCA. Cox relative regression method was applied to create a risk signature. Besides, the immune microenvironment of the signature was investigated by utilizing the CIBERSORT algorithm. Results. Based on 27 DEcircRNAs, 65 DEmiRNAs, and 780 DEmRNAs, the circRNA-miRNA-mRNA network was finally set up. Functional enrichment unearthed that the regulatory network might participate in phosphorylation negative regulation, MAPK pathway, and PI3K/AKT pathway. This study established a risk scoring signature based on the seven immune-related genes (IRGs) (MARP14, RASGR1, PTK2, HMGB1, DKK1, RARB, and IGF1R), which was validated for its reliability. A stable and accurate nomogram combining immune-related risk scores with clinical features was constructed. Furthermore, we observed that the risk model was also related to the immunocyte infiltration. Conclusion. Our study successfully created a circRNA-associated regulatory network and further developed an immune-related model to forecast the clinical outcome of ESCA patients as well as to assess their immune status.


Introduction
Esophageal cancer (ESCA) is the most frequent type of malignancy, which ranked the sixth for the mortality rate of all tumors [1]. According to GLOBOCAN 2018 data, it was estimated that around 572,000 new cases of ESCA and 508,000 fatalities occurred in that year [2], and the 5-year overall survival (OS) rate was merely 15-25% [3]. Despite improvements in the diagnosis and treatment of ESCA, most patients are not diagnosed until the middle or late stages because of the insidious nature of early symptoms. erefore, there is an urgent need for more precise and effective targets to aid in the diagnosis and personalized therapy of ESCA patients.
CircRNAs, defined as endogenous noncoding RNAs, have a single-stranded, covalently closed-loop structure. eir closed structure is produced by reverse splicing the 5′ and 3′ ends of linear RNAs, making them more resistant to nucleic acid exonucleases, as well as more stable in tissues and plasma compared to linear RNAs [4,5]. In 2011, the ceRNA hypothesis was first proposed, with the implication that circRNAs could negatively regulate the expression of miRNAs and influence their interactions with mRNAs [6].
In recent years, numerous studies showed the carcinogenic roles of circRNAs. CircRNA 100146, for example, could regulate SF3B3 through modulating miR-361-3p, thereby promoting the advancement of lung cancer [7]. CircNRIP1, a tumor promotor, may influence AKT1 expression in gastric cancer by binding to miR-149-5p [8]. Similarly, CircSETD3 was reported to be significantly overexpressed in nasopharyngeal carcinoma, which may dramatically enhance cancer cell invasion and migration [9]. Numerous studies indicate that certain circRNAs have cancer-suppressive effects. CircTADA2A-E6 suppresses breast cancer development by targeting miR-203a-3p [10]. Furthermore, circCUL2 has the potential to limit the progression of gastric cancer via autophagy activation mediated by miR-142-3p/ ROCK2 [11].
Currently, the 8th edition of ESCA TNM staging is the primary tool used clinically to assess the prognosis of patients with ESCA [12]. However, this approach is general and limited by individual differences. erefore, there is a pressing need to exploit new methods that can accurately evaluate the survival of ESCA. Recently, mRNA gene signatures based on immune-related genes have become a research hotspot for predicting the prognosis of patients. For example, Chen et al. developed an immune signature to forecast head and neck cancer patients' prognosis [13]. Moreover, the prognostic index model of 11 immune-related genes was constructed in hepatocellular carcinoma [14].
Although previous research has demonstrated the performance of the immune-related model in many malignancies, few have considered using immune genes to create a prognostic risk model for ESCA.
In this present study, we successfully constructed an esophageal cancer-associated ceRNA network using GEO datasets and the TCGA database. We further developed an immune-related risk signature by selecting seven immunerelated genes from the ceRNA network. Finally, we further detected the expression patterns of signature genes by in vitro analysis. Our study not only provides an optimal prognostic model to predict ESCA patients' OS but also opens up a new insight for ESCA treatment.

Development of the CircRNA-Associated Regulatory
Network.
e basic information of circRNAs was downloaded from circBase (http://circbase.org/). e CSCD online tool (http://gb.whu.edu.cn/CSCD/) was performed to obtain targeted miRNAs of the circRNAs, and these potential miRNAs were determined by DEmiRNAs. Next, the relationship of miRNA and mRNA was generated by the TargetScan, miRDB, and miRTarBase prediction tools [15][16][17]. Afterward, these target mRNAs were intersected with the determined DEmRNAs. Ultimately, a circRNAassociated regulatory network was set up by Cytoscape.

Functional Enrichment Analysis.
Gene ontology (GO) is a bioinformatics method designed to measure gene function. KEGG, a database that integrates genomic, chemical, and biological systems' information, has been widely used to integrate and interpret datasets generated by highthroughput experimental techniques. To uncover the functional mechanism of the regulatory network, GO and KEGG methods were performed, and p < 0.05 was considered significant.

Construction of the CircRNA-Based Prognostic Signature.
We first collected immune-relevant genes (IRGs) from the ImmPort database, which is a public database providing an up-to-date genes' list and functional information for IRGs. en, the overlapping genes between IRGs and the identified mRNAs retrieved from the circRNA regulatory network were utilized to determine differentially expressed IRGs. To develop a favorable risk model, the GSE53625 set was divided into the training cohort and test cohort at 6 : 4 ratio. e prognostic IRGs were first identified by the univariate regression method. To shrink the model of overfitting, LASSO regression was conducted. Finally, multivariate Cox method was conducted to create a signature to predict OS for ESCA cases. e risk score for patients was generated by the following equation: risk score � gene 1 expression × coef 1 + gene 2 expression × coef 2. . . + gene n expression × coef n, where coef n represents the regression coefficient. e test cohort, entire cohort, and TCGA set were utilized to confirm the performance of the model.

Establishment of the Nomogram.
To estimate the independent capacity of the prognostic model, univariate and multivariate Cox methods were carried out. A nomogram to predict OS of ESCA was developed by combining the clinical traits and prognostic model. To validate the predictive power of this signature, we plotted calibration maps and ROC curves.

Immunocyte Infiltration Analysis.
CIBERSORT is a bioinformatics method to identify the infiltration of immune cells. In this case, we input the reference gene profiles and analyzed the infiltrating proportions of each cell type among 22 types of tumor-infiltrating immune cells (TICs).

Gene Set Enrichment Analysis.
We applied GSEA to examine the underlying biological function related to highrisk patients in the signature. e gene sets of Hallmark and KEGG were obtained from the Molecular Signatures Database. Gene sets with a normalized P value less than 0.05 were deemed significant.
2.9. Cell Culture. Two human ESCA cell lines (KYSE-150 and Eca109), as well as one normal esophageal epithelial cell line (HEEC), were provided from the Chinese Academy of Sciences (Shanghai, China). ESCA cells and HEEC were cultured in the RPMI-1640 medium with 10% fetal bovine serum (FBS) and 1% penicillin/streptomycin with 5% CO 2 at 37°C.

RNA Extraction and qRT-PCR.
Total RNA was isolated from cell samples utilizing RNA-easy isolation reagent (Vazyme Biotech), which was then reversed into cDNA by PrimeScript Mix reagent (Takara). SYBR Green Premix (Vazyme Biotech) was prepared for use in the PCR system. e value of individual genes was finally standardized to the GAPDH expression level. e primer sequences for indicated genes are displayed in Supplementary Table S1. 2.11. Statistical Analysis. Kaplan-Meier method was employed to evaluate the clinical outcome differences between two risk groups. In addition, receiver operating characteristic (ROC) analysis was performed to examine the reliability of the model. All statistical data were analyzed using GraphPad 8.0 and R software version 4.0.

Identification of DEcircRNAs, DEmiRNAs, and DEmRNAs.
A total of 28 DEcircRNAs, including 15 highly expressed circRNAs and 13 lowly expressed circRNAs, were identified in the GSE131969 dataset ( Figure 1). We eliminated hsa_circ_0001336 since it was not found in the CSCD database. After performing the CSCD online tool, we identified 1140 miRNAs that might target these 27 circR-NAs. After crosschecking with DEmiRNAs retrieved from the TCGA database, only 65 miRNAs remained. We then predicted potential mRNAs targeted by these 65 DEmiRNAs using TargetScan, miRDB, and miRTarBase. ese candidate mRNAs were intersected with DEmRNAs selected from GSE53625, resulting in 780 genes shared.

Construction of the ceRNA Network. A total of 27
DEcircRNAs, 65 DEmiRNAs, and 780 DEmRNAs were finally determined to establish a circRNA-associated network by Cytoscape 3.8 software (Figure 2). Figure 3(a), the significant function enrichments were Wntrelated pathway, negative regulation of phosphorylation, and negative regulation of protein phosphorylation. KEGG revealed that tumor pathways were greatly enriched, including the MAPK pathway, PI3K-AKT pathway, and Ras pathway (Figure 3(b)).

Construction of the CircRNA-Based Signature in the
Training Set. We randomly divided GSE53625 containing 179 patients into training and internal validation series at a 6 : 4 ratio. e IRGs were overlapped with DEmRNAs from the ceRNA network mentioned above, resulting in 64 differentially expressed IRGs remained. ese genes were further subjected to the univariate Cox method, LASSO Cox, and multivariate Cox analysis consecutively to establish a seven-IRG signature (Table 1).
As expected, the clinical outcome of the high-risk group was dramatically dismal than that of the low-risk group (Figure 4(a)). e prognosis accuracy of our model was assessed by the areas under the curve (AUC) of the ROC (Figure 4(e)). In addition, we explored the distribution of the risk value and survival outcome of each patient (Figure 4(i)).

Validation of the CircRNA-Based
Signature. Our proposed signature is verified to be efficient when tested in the internal cohort. We observed that high-risk ESCA patients showed worse OS than its counterpart (Figure 4(b)). e similar result was also confirmed in the entire set and TCGA cohort (Figure 4(c)). e AUC of ROC was calculated to estimate the reliability of the model in the above cohorts (Figures 4(f )-4(h)). At final, the distribution of the risk score and clinical outcome of each ESCA case are shown in Figures 4(j)-4(l)).

Establishment of the Nomogram Based on the Risk Model.
To evaluate the independent power of the prognostic signature, Cox relative hazard regression methods were conducted (Figures 5(a) and 5(b)). We found that only the risk score was greatly meaningful for forecasting OS in univariate and multivariate Cox methods. In addition, a nomogram was developed by combining the risk model and clinical parameters for forecasting the clinical outcome of ESCA cases ( Figure 5(c)). Calibration curves showed excellent Type Type

Correlation Analysis of Immunocyte Infiltration and Risk
Score. By applying the CIBERSORT algorithm, a bar plot ( Figure 6(a)) was performed to represent the discrepancy in the infiltration level of 22 types of TICs between two risk groups. To further confirm the relationship of the risk signature and immunocyte infiltration in ESCA, the Pearson correlation analysis was performed. e result showed that M0 macrophage, resting memory CD4 + T cell, and activated NK cells had a positive association with the risk score, while CD8 + T cells and plasma cells displayed a negative correlation with the risk score ( Figure 6(b)).

Prognostic
Performance of the CircRNA-Based Signature in the Immunosuppressive Microenvironment. Recently, the cancer-immunity cycle, which depicts the immune surveillance on tumor cells, has become a research highlight in cancer immunotherapy. It involves a series of procedures including cancer cell antigen release [18]. We selected 54 negative regulatory genes involved in the cancer-immunity cycle from the Tracking Tumor Immunophenotype website for further analysis. Next, we analyzed the relationship between negative regulatory genes and our signature. e results suggested that DNMT1, NECTIN3, and SMC3 had a positive correlation with the risk score, while NCR3 presented the opposite trend ( Figure 7).
3.9. GSEA. As shown in Figure 8(a), angiogenesis, IL2/ STAT5 pathway, inflammatory response, oxidative phosphorylation, PI3K/AKT/MTOR pathway, and TGF-beta pathway were dramatically enriched in cases with high-risk scores. We also observed that six KEGG pathways were activated in the high-risk group, including the MTOR pathway and TGF-beta pathway, which are in line with the result of Hallmark (Figure 8(b)).

Detection of Signature Gene Expressions in Cell
Lines by the qRT-PCR Assay. Finally, a circular RNA-associated ceRNA subnetwork was constructed based on our immunerelated risk signature ( Figure 9). Subsequently, we further examined the expression levels of seven signature genes in HEEC, Eca109, and KYSE-150 using qRT-PCR. e result showed that HMGB1, PTK2, and IGF1R had significantly higher expression levels in Eca109 and KYSE-150 and lower expression levels in HEEC. On the contrary, RASGRP1,  MAPK14, RARB, and DKK1 were lowly expressed in Eca109 and KYSE-150 and highly expressed in HEEC ( Figure 10).

Discussion
ESCA ranks eighth and sixth among cancers according to morbidity and mortality, respectively, with high rates of metastasis and recurrence [1]. Consequently, it is urgent to exploit powerful biomarkers for ESCA. Accumulating evidence suggests that circRNAs could interact with miRNA, resulting in a corresponding decrease in miRNA activity [19]. miRNAs act primarily on mRNAs, but the circRNA competition process impedes the regulation of target gene expression by miRNAs, which further affects tumor formation, development, and metastasis. For example, hsa_-circ_0006168 could boost ESCC cell viability and metastasis by binding with miR-100 [20]. As discovered by Meng et al., CIRS-7 could accelerate the development of ESCC by binding with miR-876-5p and upregulating the MAGE-A family gene [21]. In addition, circGSK3β is reported to act as an indicator for the prognosis of ESCA [22]. us, circRNA may become a new molecular target for ESCA treatment.

Journal of Oncology
We identified 27 DEcircRNAs, 65 DEmiRNAs, and 780 DEmRNAs from GEO, TCGA, and CSCD databases, which in turn establish an ESCA-associated ceRNA network. According to the immune-related signature, we also created a circRNA-miRNA-mRNA subnetwork. Wang and his colleagues proved that hsa_circ_0005654 had higher expression in early gastric carcinoma specimens than in normal cases [23]. hsa_circ_0001313, also named circCCDC66, has been shown to have a cancer-promoting role in stomach cancer, colon cancer, kidney carcinoma, and lung cancer [24][25][26][27]. However, other six circRNAs have not been studied. Hence, the role of these circRNAs in tumors, especially in ESCA, needs to be unearthed in future research studies.
To uncover the underlying mechanism of the ceRNA network, we implemented a functional annotation of mRNA from the network. e results showed that the Wnt-related pathway, negative regulation of phosphorylation, and negative regulation of protein phosphorylation were greatly enriched. Enrichment of the KEGG signaling pathway suggested that tumor-related pathways were activated including the MAPK pathway, PI3K-AKT pathway, and Ras pathway. MAPK signaling pathway has the basic function of regulating cell growth, survival, and differentiation. As one of the clearest signaling pathways in tumor biology research, its oncogenic effect when abnormally activated has been demonstrated in many tumors, including esophageal cancer [28]. e mechanisms by which the PI3K-AKT signaling pathway facilitates cancer development may include the following: cell proliferation, metabolic reprogramming, suppressing autophagy, and promoting EMT [29].
In the present study, Cox relative regressions and LASSO regression were performed to establish a signature of seven IRGs, including HMGB1, PTK2, RASGRP1, MAPK14, IGF1R, RARB, and DKK1. HMGB1, a chromatin-associated protein, has been reported to regulate ESCC cell proliferation and radiosensitivity [36]. Di and his colleagues observed that overexpression of HMGB1 could promote cell radioresistance after irradiation. PTK2, also known as FAK, was found to facilitate metastasis of ESCC cells through the miR-4324/FAK pathway [37]. As suggested by Ma et al., high IGF1R expressions mediate cell proliferation and apoptosis, indicating that it could act as a crucial oncogene in esophageal cancer [38]. Moreover, Lyros et al. discovered that the upregulation of DKK1 in esophageal adenocarcinoma enhances tumor growth and progression through AKT-phosphorylation and the Wnt axis [39].
Recently, immune cells have been proved to play a central part in the tumor environment [40,41]. Moreover, infiltration of immune cells has been characterized as a reliable prognostic factor in esophageal cancer. Previous reports suggested that high infiltration levels of M2 macrophages [42,43] and Tregs [44] were strong evidence for tumor development and metastasis. Conversely,    Schumacher et al. revealed that increased enrichment of CD8 + T cells was a favorable prognostic indicator for the clinical outcome [45]. We observed that CD8 + T cells are negatively associated with risk scores, consistent with the hypothesis that a high abundance of CD8 + Tcells may indicate a better clinical outcome. In a word, our study highlights the importance of maintaining certain immune cells in tumor immunotherapy. Moreover, we found that the risk score was also related to the expression of CIC negative regulatory markers. DNMT1 belonged to the DNA methyltransferase family, and its overexpression is identified in human T-cell, B-cell, and myeloid malignancies, indicating that DNMT1 may play a crucial role in tumor maintenance [46,47]. Nectin-3, a nectin family member, took part in regulating the formation of adherent junctions and had been reported to be a novel biomarker in tumorigenesis. It has also been demonstrated that Nectin-3 could contribute to lymphocyte extravasation by interacting with Nectin-2 [48]. e function of SMC3 (structural maintenance of chromosomes 3) was described to induce tumorigenesis, and its expression can be downregulated in CD8 + T cells after immunotherapy [49].
However, this study still had some limitations. Firstly, the construction of the ceRNA network was mainly based on bioinformatics analysis, and its regulatory mechanism required further experiments. In addition, the result of the correlation analysis between the signature and immunocyte infiltration and CIC negative regulatory genes needs to be validated in additional cohorts.

Conclusion
is study determined a range of DEcircRNAs, DEmiRNAs, and DEmRNAs and set up an ESCA-associated network in ESCA. Furthermore, we successfully developed an immunerelated prognostic model which could accurately forecast clinical outcomes, mirror the immune microenvironment, and provide individualized treatment for ESCA patients.

Data Availability
All datasets generated for this study can be found in e Cancer Genome Atlas (https://portal.gdc.cancer.gov/) and the NCBI Gene Expression Omnibus (GSE131969 and GSE53625).

Conflicts of Interest
e authors declare that they have no conflicts of interest.

Authors' Contributions
ZH and SH designed and wrote the manuscript. JX and XC analyzed the data. JT and KZ reviewed and revised the manuscript. All authors approved the final manuscript.

Supplementary Materials
Supplementary Table 1