The Transcriptional and Immunological Roles of Six2 in Clear Cell Renal Cell Carcinoma
1The Ministry of Education Key Laboratory of Laboratory Medical Diagnostics, College of Laboratory Medicine, Chongqing Medical University, Chongqing, 400016, China
2Department of Laboratory Medicine, Women and Children’s Hospital of Chongqing Medical University, Chongqing, 400016, China
*Corresponding Authors: Tao Song. Email: CQSFYsongtao@163.com; Qianyin Li. Email: firstname.lastname@example.org
Received: 28 March 2022; Accepted: 25 May 2022
Abstract: Background: Six2, a transcription factor, exerts an oncogenic role in clear cell renal cell carcinoma (ccRCC). Increased Six2 expression could enhance cancer metastasis. However, the regulatory mechanism of Six2 in promoting metastasis remains unclear. The purpose of this study is to analyze the regulatory pattern of Six2 and the potential role of Six2 in the tumor immune microenvironment. Materials and Methods: Firstly, transcriptional data in TCGA-KIRC cohorts was used to analyze the relationship between Six2 expression and clinical information. Secondly, we detect the association between Six2 and the tumor immune microenvironment in ccRCC. Then, we analyzed Six2-related differentially expressed genes (DEGs) and constructed a prognostic model using the Lasso-Cox algorithm by integrating Six2 ChIP data and co-expressed genes. Next, we analyzed the clinical significance and immunotherapy sensitivity of this model. Results: Six2 was overexpressed in RCC cells compared with normal kidney cells and upregulated Six2 was positively linked with clinical stage, grade, T stage, M stage, and poor survival. And Six2 was correlated with the remodeling of the tumor microenvironment. Potential downstream effectors and biological functions regulated by Six2 were identified using in silico analysis. Meanwhile, a risk model based on 8 Six2 target genes was established to classify ccRCC patients into high-and-low risk groups. This risk model showed a reliable ability to forecast the overall survival of ccRCC patients and predict the susceptibility to immunotherapy. Conclusions: Our findings provide a promising prognostic indicator for ccRCC patients and help better understand the transcriptional and immunological role of Six2 in ccRCC.
Keywords: Six2; ccRCC; TCGA; tumor immune environment
Renal cell carcinoma (RCC) is one of the 10 most common cancers in both sexes, accounting for about 4% of human malignancies. An expected 76,080 new cases of kidney & renal pelvis cancer and 13780 deaths from these diseases were predicted in the United States in 2021 . RCC is classified into different pathological groups according to genomic and clinical features, with clear cell RCC (ccRCC) being the most common subtype . Although the therapy landscape is improving by introducing targeted treatment and immunotherapy, approximately one-quarter of patients receiving nephrectomy finally show disease relapse and metastasis. However, few reliable genetic biomarkers indicate treatment outcomes in both localized and metastatic tumors [3,4]. Therefore, RCC-specific diagnostic biomarkers for early detection and tumor progression monitoring are highly needed in the clinical management of RCC patients .
Transcription factor Six2 plays a pivotal role in maintaining the mesenchymal progenitor population as undifferentiated (MET) [6,7]. Recent studies found that Six2 promotes metastasis by decreasing E-cadherin to facilitate epithelial-mesenchymal transition (EMT) and upregulating SOX2 and NANOG expression to enable metastatic tumor colonization [8,9]. Apart from the metastatic ability of carcinoma cells, it is also well-established that cancer metastasis frequently accompanies the deposition and remodeling of the extracellular matrix, an important component of the tumor microenvironment . However, the potential role of Six2 in the immune microenvironment remains unclear. Although several downstream targets of Six2 were identified, there is still an urgent need to better understand the additional target genes regulated by Six2.
In this study, we discovered that Six2 expression showed a significant increase in ccRCC cells, and abnormal expression of Six2 was correlated with worse clinical outcomes and immune environment. Then, Gene Ontology (GO), Kyoto Encyclopedia of Genes (KEGG), and Gene set enrichment analysis (GSEA) of Six2-related DEGs were performed. Next, we integrated Six2-associated DEGs, ChIP-seq data from the GTRD database, and co-expression genes of Six2 to generate a prognostic model. After validation of this model, we assessed its susceptibility to immunotherapy. Our data provided a comprehensive analysis of the potential regulatory network and immunological role of Six2. The workflow diagram is depicted in Fig. S1.
2 Materials and Methods
2.1 Patient Sample Collection and Processing
Transcriptome data (FPKM) and clinical and pathological information for KIRC were downloaded from The Cancer Genome Atlas (TCGA) public database, which includes 72 normal samples and 539 tumor samples. All RNA-seq data were normalized before analysis.
The relationship between the Six2 expression level and clinic-pathological information, including TNM stage and Fuhrman grade, was analyzed and visualized by R. The survival data and survival curve were analyzed by R software (version 3.6.3) using survival package. Cancer samples were separated into two groups according to Six2 median mRNA level: Low (lower than median expression level) and high (higher than median expression level). R packages (survival, survminer) were utilized to establish KM survival curves.
2.2 Six2-Related Differentially Expressed Genes (DEGs)
RNA expression data was divided into two degrees by Six2 median expression levels: High (263) and low (263). Then DEGs between high and low groups were identified by package limma in Rstudio. The adjust p < 0.05 and |FC| > 1.2 were seen as DEGs. Packages ggpubr, ggthemes mapped volcano plot and heatmap of Six2-related DEGs.
2.3 Immune Microenvironment Analysis
The stromal score, immune score, and estimate score were analyzed by the Estimate R package. CIBERSORT (https://cibersort.stanford.edu/) was used to quantify the relative proportion of human immune cells. The correlation analysis detected the potential immune cells affected by Six2.
2.4 Enrichment Analysis
Subsequently, gene ontology (GO), including biological process (BP), cellular component (CC), and molecular function (MF), KEGG pathway, Gene set enrichment analysis (GSEA) for DEGs were performed in R using R packages “clusterprofiler.” R package ggplot2 was used for graphical plotting.
2.5 CHIP Data Collection and Screening Core Genes
The target genes of Six2 were projected based on The Gene Transcription Regulation Database (GTRD; http://gtrd.biouml.org/) . Cytoscape was used to visualize the network of Six2 target genes. We obtained the highly correlated genes by analyzing the Pearson correlation coefficient in TCGA KIRC samples. Then the shared genes among DEGs (|FC| > 1.2), co-expressed genes (|correlation coefficient| > 0.2, p < 0.05), and Six2 potential regulating genes were selected to further model construction.
2.6 Six2-Related Risk Model Construction and Validation
The 48 core genes were used for establishing a risk score model. The lasso-Cox regression algorithm was employed to narrow down the number of candidate genes. R “glmnet” package was utilized to decide the Lasso-Cox coefficients of constructed nine-gene model. TCGA datasets were divided into two groups equally for training sessions and testing sessions to examine the robustness of this model. The risk score of each ccRCC patient was calculated based on our equation. Then, the KIRC samples were assigned to high or low-risk groups according to the median risk score. KM survival and ROC curves were done in R software. R package “rms” was utilized to generate a scoring nomogram to assess the risk status of patients. The Cancer Immunome Atlas (TCIA, https://tcia.at/) provides the results of immunophenoscore (IPS) of KIRC patients.
2.7 Cell Culture and Antibodies
HEK293T, HK2, and 786-O cells were cultured in Dulbecco’s Modified Eagle Medium (DMEM), DMEM/F12 -Dulbecco’s Modified Eagle Medium: F-12, RMPI 1640 medium (Gibco, Carlsbad, CA, USA), respectively. 10% fetal bovine serum (ExCell Bio, China) and 1% penicillin/streptomycin (PS) (Invitrogen, Grand Island, NY, USA) were added to these media. The cells mentioned above were grown at 37°C in an incubator with 5% CO2.
2.8 Western Blot Analysis
Cells were lysed in 1% SDS lysis buffer and transferred into 1.5 ml eppendorf tubes for 10-min boiling at 95°C. Lysates were centrifuged at 12,000 rpm, at RT for 10 min, and the supernatants were collected. Then BCA Kit was used to determine the concentration of each sample. The proteins (10 μg each sample) were separated on 10% SDS-PAGE. Gels were electrically transferred to the PVDF membrane and blocked with 5% skim milk diluted in TBST. All blots were incubated with prim y antibodies anti-Six2 (Proteintech) and anti-Tubulin (TRANSGEN) at 4°C overnight. After washing with TBST, the corresponding HRP-conjugated secondary antibodies were applied to detect proteins. The signals were visualized by ECL Reagents (Merck, Billerica, MA, USA).
2.9 The Human Protein Atlas(HPA)
The Human Protein Atlas (https://HumanProteinAtlasproteinatlas.org/) provides a large amount of omics data by integrating proteomics, transcriptomics, and pathology data. This database contains pathology information from 20 different common cancer types. In this study, the protein expression of CKAP4, CUBN, IQGAP2, SLC12A8, and SLC47A1 were downloaded from the Human Protein Atlas.
2.10 Statistical Analysis
Chi-squared test was applied to compare the differences in clinical information. Differences in overall survival (OS), disease-specific survival (DSS), and progression-free interval (PFI) were determined by Kaplan–Meier and the log-rank test. R “timeROC” package plotted the receiver operating characteristic (ROC) curves. The effectiveness of the risk score for predicting OS was assessed by the section under the ROC curve. Two-sided p-values with p < 0.05 were viewed as statistically significant. Data processing, analysis, and plot were performed in R software and GraphPad Prism 8.
3.1 Correlation Between Six2 Expression and Clinicopathological Characteristics in ccRCC
To determine the Six2 expression pattern in ccRCC, we detected Six2 protein expression levels in 293T, HK2, and 786-O cells. Higher Six2 expression was observed in 293T and 786-O cells compared with that of HK2 cells (Fig. 1A). Then the impact of Six2 on prognosis in TCGA cohorts was evaluated. It suggested that elevated Six2 expression remarkably predicted poor overall survival, disease-specific survival, and progression-free interval (Figs. 1B–1D).
Next, we assessed the relations between Six2 expression and clinicopathological characteristics in TCGA datasets. It was found that the progression of the clinical stage, grade, T stage, and M stage of the tumor was in line with the Six2 expression increase (Figs. 1E–1H), which may indicate Six2 can be a potential biomarker of ccRCC.
3.2 Six2-Related Immune Cell Infiltration Analysis
To comprehensively understand the role of Six2, we also analyzed the role of Six2 in the tumor microenvironment. The results showed that higher Six2 expression was correlated with higher stromal scores, immune scores, and estimate scores among ccRCC patients (Figs. 2A–2C). To be more specific, we then analyzed the immune cells’ differences in different Six2 expression groups (Fig. S2). We observed higher activated CD4 memory T cells, regulatory T cells, and M0 macrophages in the Six2 high expression group and lower M1 macrophages, resting Dendritic cells, and resting Mast cells in the Six2 low expression group (Fig. 2D).
3.3 Identification of DEGs Based on the Six2 Expression Level
To find genes closely related to Six2, the tumor samples in TCGA-KIRC datasets were analyzed with limma package, and 796 DEGs (498 up-regulated genes and 298 down-regulated genes) were identified from Six2 high and low expression groups. The heatmap and volcano map displaying DEGs wereshown in Figs. 3A, 3B. Then, we performed GO and KEGG pathway analysis to clarify the function of Six2 in ccRCC. GO analysis suggested that these DEGs mainly participated in extracellular structure organization, extracellular matrix, and leukocyte migration. For molecular function (MF), Six2 was significant in functions associated with extracellular matrix structural constituent, glycosaminoglycan binding, and heparin binding. The following terms related to cellular component (CC) were mainly enriched in collagen-containing extracellular matrix, apical part of cell, and endoplasmic reticulum lumen (Fig. 3C). Furthermore, the following KEGG enrichment results included PI3K-Akt signaling pathway, cytokine-cytokine receptor interaction, tight junction, and ECM-receptor interaction (Fig. 3D). GSEA analysis for Six2 showed that Six2 was involved in angiogenesis, apical junction, epithelial-mesenchymal transition, hypoxia, ECM regulators, extracellular matrix organization, ECM-receptor interaction, and cytokine-cytokine receptor interaction (Fig. 3E).
3.4 Construction and Validation of an Integrated-Gene Signature Model
ChIP-seq is a technique for understanding transcriptional regulation. To fully understand the regulatory patterns of Six2, we obtained Six2 ChIP-seq data from GTRD (Fig. 4A) and used R studio to analyze co-expressed genes of Six2 (Fig. S3A). Then, 48 genes were selected for Lasso-Cox regression analysis to produce the optimal model (Fig. 4B). An eight-gene risk signature was generated, and the risk score of every sample was figured out based on mRNA levels of eight Six2-related genes: Risk score = 0.0871348 * CLMP + 0.0420739 * CKAP4 + 0.0962748 * COL22A1 + 0.1474064 * SLC12A8–0.3146996 * MOB3B-0.0329607 * IQGAP2-0.0853117 * CUBN−0.0955748 * SLC47A1. And the survival curves, sitecount (TF binding sites in promoter region), and immunochemistry results of these genes were depicted in Figs. S3B, S3C, and Fig. S4A. KIRC samples were classified into low-and high-risk groups formulated on the median value of the risk score. Figs. 4C–4E displayed the expression profile of prognostic genes, risk patterns, and survival status between the high-and low-risk groups. It was also presented that patients in the high-risk group had a worse clinical outcome than the low-risk groups (Figs. 4G, 4H). To further assess the prognostic value of the risk subgroup in KIRC patients, we performed univariate and multivariate analyses and found out the risk score could be an independent risk factor in ccRCC (Fig. 5A). Furthermore, the clinicopathologic information was also analyzed between high-and low-risk groups. It was shown that the low-risk group has better OS than that of the high-risk group in terms of age, gender, stage, and tumor stage (Fig. 5B).
3.5 Predicting Risk and Response of Immunotherapy for Individual Patients
To quantify the risk evaluation of each KIRC patient, a nomogram was constructed to predict 1-, 3- and 5-year survival probability with eight parameters (Fig. 6A). After scoring each patient, we examined several representative immune checkpoints, and found that the low-risk group showed higher PD-L1 expression, which may provide hints for potentially applicable drugs for ccRCC patients (Fig. S5). Then, we applied TCIA to predict the sensitivity of patients to immunotherapy (Fig. 6B). We discovered that the low-risk group showed better IPS than the high-risk group, which means that the low-risk group might be more susceptible to immune checkpoint inhibitors.
Initially found to be highly expressed in nephron progenitor cells, Six2 has been found to be highly expressed in nearly every cancer type examined, such as colorectal cancer and lung cancer [12,13]. Silenced Six2 in adult kidneys was found reactivated in renal carcinoma, and its misexpression has been linked to the etiology of renal cancer . It was demonstrated that Six2 promotes distant metastasis through downregulating E-cadherin levels and upregulating Sox2 expression by directly binding with the Srr2 enhancer of Sox2. Undoubtedly, a high expression level of Six2 could promote EMT and enhance the metastatic capacity of cancer cells. Beyond cell property change of malignant cells, the tumor immune microenvironment also contributed to cancer metastasis . However, the downstream effectors and immunological role of Six2 in kidney cancer have not been comprehensively examined.
In the present study, we analyzed the potential downstream targets of Six2 using transcriptomic data in TCGA and ChIP data in GTRD . Firstly, an elevated expression of Six2 was detected in RCC cell lines and was positively correlated with malignant behavior. Then, the role of Six2 in the tumor microenvironment was explored. Thirdly, we divided TCGA patients into two groups by Six2 expression and identified DEGs and their related function by GO and KEGG analysis. In addition, we further constructed a Six2-related risk model based on the overlapping gene signature among DEGs, ChIP data, and co-expressed genes. After evaluating the model, we established a personalized risk assessment nomogram and predicted response to immunotherapy of KIRC patients.
RCC affects 400,000 more people worldwide every year. Although it can be effectively treated with nephrectomy in the early stage, more than a third of patients finally develop into metastatic RCC, which has high mortality rates and shows different biological features from the localized disease . Our data revealed that patients with higher Six2 levels showed advanced pathological stage, and metastasis status, indicating that Six2 may influence cancer metastasis. It was suggested that cancer-cell invasion is the first step of metastasis cascades, whereas the EMT program has been reported to play a pivotal role in tumor invasion [17,18]. Interestingly, previous studies hadshown the promoting role of Six2 in the EMT process. Wan et al. found that Six2 was upregulated and boosted EMT under TGF-β treatment . Next, we found that Six2 expression was correlated with the proportion of M1 macrophage and regulatory T cells. Regulatory T cells are one of the primary immunosuppressive cells in the tumor microenvironment, whereas M1 macrophages produce high level of inflammatory cytokines and play an anti-tumor role in cancer [20,21].
Our data also noticed that most of the DEGs shown in the heatmap showed a pro-metastatic role in various cancers, including CTHRC1, CXCL5, LOX, etc. [22–24]. Further studies could analyze the relationship or correlating role between these genes and Six2 in EMT development. Then, GO, KEGG, and GSEA methods were used to analyze Six2-related DEGs. PI3K-Akt signaling pathway, focal adhesion, leukocyte migration, ECM-receptor interaction, and cytokine-cytokine receptor interaction were enriched in our analysis. It has been shown that EMT companies with loss and reconstruction of cell-cell junctions, cell-matrix interactions, and cytoskeleton. PI3K-Akt signaling pathway has also been verified to participate in the EMT process . Six2 may also have a latent function in ECM remodeling and further regulating tumor immune microenvironment since hypoxia and ECM-regulators were enriched. All of these implicated that Six2 could be an important regulator in EMT and immune microenvironment, which are two important factors in triggering tumor metastasis .
It was noticed that dysregulation of gene expression programs controlled by transcription factors can cause a wide range of diseases, including cancer. Mutation and amplification in transcription factors have been known to facilitate carcinogenesis . However, the direct inhibition of transcription factors is untractable due to the difficulties of targeting protein-DNA interactions . Thus, identifying novel downstream effectors or pathways may provide ideal therapeutic targets controlled by Six2. In the present study, we combined the transcriptomic data, ChIP data, and co-expression analysis data in the public database to analyze potential downstream effectors of Six2. Moreover, a novel Six2-related gene risk model was constructed to predict the OS of KIRC patients in training and testing sessions. As cancer immunotherapy successfully treats multiple neoplasms by activating the immune response, we also tested the potential sensitivity of patients to immunotherapy and found out that the risk score model may be an available predictor of immunotherapy [29–32].
There are still several limitations in our study. Firstly, using TCGA transcriptomic data cannot well simulate the gene expression profile acquired from Six2 overexpression or knockdown treatment. Secondly, the ChIP-seq data of Six2 is mainly based on human fetal kidney tissues rather than kidney cancer tissues. Although there are some common features between embryogenesis and tumorigenesis, the data cannot fully reflect the overall regulatory change of Six2 in cancer . Besides, detailed experiments should be done to get a comprehensive understanding of the immunological function of Six2 in ccRCC.
Our study revealed that elevated Six2 expression was closely related to unfavorable survival and cancer progression. Six2 also played an essential part in the tumor immune microenvironment. Novel potential downstream targets of Six2 were selected by integrating transcriptomic data and ChIP-seq data of Six2, and a Six2-related prognostic model was established. Our findings provide insight into the significant transcriptional and immunological role of Six2 in ccRCC.
Acknowledgement: The authors would like to appreciate data availability in the TCGA database.
Authorship: Dayu Tian and Yang Shi have contributed equally to this work. Tao Song and Qianyin Li conceptualized and designed this study. Lei Li, Xiangmin Qiu wrote the first draft of the manuscript.
Ethics Approval and Informed Consent Statement: TCGA is a public database. The information of patients in the database has obtained ethical approval. Users have free access to the data for research and publishing relevant articles. Our study analyzed open-source data, so there are no ethical issues.
Availability of Data and Materials: The datasets in this study are available in the Cancer Genome Atlas (https://portal.gdc.cancer.gov/) and GTRD (http://gtrd.biouml.org/). Further inquiries can be directed to the corresponding authors.
Funding Statement: This study was funded by the National Natural Science Foundation of China (Grant No. 81873932), the Natural Science Foundation of Chongqing (Grant No. cstc2021jcyj-msmX0019), the Science and Technology Research Project of Chongqing Municipal Education Commission of China (Grant No. KJQN202000438), and Chongqing Science and Technology Bureau Technology Innovation and Application Development Special Key Projects (Grant No. cstc2019jscx-dxwtBX0032).
Conflicts of Interest: There is no conflict of interest. The authors have no relevant financial or non-financial interests to disclose.
|This work is licensed under a Creative Commons Attribution 4.0 International License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.|