The prognostic value and immunological role of the small mother against decapentaplegic proteins in kidney renal clear cell carcinoma.

BACKGROUND: The small mother against decapentaplegic proteins (SMADs) are a family of transforming growth factor (TGF)-β signal transduction molecules, playing a vital role in the initiation and development of tumors. This research aimed to determine SMADs' prognostic values and their involvement in immune infiltration.
METHODS: Expression patterns and prognostic values of SMADs were evaluated by pan-cancer analysis in multiple cancer cohorts based on The Cancer Genome Atlas data. cBioPortal database was used for genetic mutation analyses. UALCAN and LinkedOmics databases were applied for the analysis of the methylation level and its correlation with gene expression, respectively. The correlation of gene expression was analyzed by Gene Expression Profiling Interactive Analysis platform. Additionally, we utilized the Tumor Immune Estimation Resource database to explore the correlation between SMAD expressions and the number of tumor-infiltrating immune cells. Functional prediction was performed by Gene Set Enrichment Analysis (GSEA) method.
RESULTS: We found that the expressions of SMAD1, 2, 3, 4, 6 were significantly decreased whereas the expression of SMAD9 was significantly increased in kidney renal clear cell carcinoma (KIRC) tissues than in normal control tissues. And aberrant DNA methylation in the promoter regions may cause the dysregulation of these differentially expressed SMADs. Also, we found that the expressions of SMAD1, 2, 3, 4, 6 decreased significantly with the progression of KIRC tumors, and their high expression level was significantly associated with favorable prognoses of KIRC patients. Genetic mutations analysis using the cBioPortal database found that there were missense mutations in SMAD2 and 4, and truncation mutations in SMAD2 and 3. Further, SMAD1, 2, 3, 4, 6 expressions showed correlations with diverse immune infiltrating cells and immune markers. In particular, SMAD1, 2, 4 expressions were strongly correlated with monocyte, tumor-associated macrophage, M1/M2 macrophage, revealing their potential to regulate the polarity of macrophages. Finally, function prediction by GSEA indicated that SMAD1, 2, 3, 4, 6 were closely involved in immune-related signaling pathways.
CONCLUSIONS: Our findings indicate that SMAD1, 2, 3, 4 and 6 were potent biomarkers for predicting the prognosis and immune cell infiltration of KIRC patients. 2021 Translational Cancer Research. All rights reserved.

Introduction

Renal cell carcinoma (RCC), present in the epithelium of renal tubules, accounts for more than 85% of all renal malignancies and 2% of the global cancer burden. Kidney renal clear cell carcinoma (KIRC), which comprises about 70% to 80% of RCC, is the most common histological subtype and a major cause of cancer-related deaths (1,2). Tumor recurrence occurs in approximately 30–40% of KIRC patients after curative treatment (2,3). Despite the fast development of novel approaches to combating KIRC (4), less than 10% of patients with metastatic KIRC can survive for 5 years after diagnosis (1). Therefore, it is urgently needly to explore potential new prognostic biomarkers or therapeutic targets, to improve the diagnosis, prognosis and individualized treatment of KIRC patients.

Transforming growth factor (TGF)-β family signaling regulates cell proliferation and carcinogenesis via activating the canonical small mother against decapentaplegic (SMAD) dependent as well as the non-canonical SMAD independent pathways (5). The SMAD family comprises transcription factors that mediate the diverse effects of the TGF-β superfamily (6). Eight members of the SMAD family exist in vertebrates, which can be divided into three different subgroups according to their distinct functions in TGF-β signal transduction: receptor-regulated SMADs (SMAD1, 2, 3, 5, 9), common mediator SMAD4, and inhibitory SMADs (SMAD6, 7) (7,8). Studies have suggested that abnormalities of SMAD-related signaling pathways are involved in tumorigenesis in some malignancies, such as cancers of the colon, lung, pancreas, etc. (9-11). Abnormal expression of SMADs has been well documented, and related studies have demonstrated that some SMADs may also be used as prognostic biomarkers in some cancers (12,13). However, the expression of SMADs family genes in patients with KIRC and their prognostic significance is still elusive. Thus, we hypothesized that some genes of SMAD could also have an important prognostic role in KIRC.

Immune cells infiltrated into tumors form a complex ecosystem in the surrounding microenvironment to control tumor development and are closely correlated with the clinical consequences of various cancers including KIRC (14-16). Tumor immune infiltrating cells are generally classified as adaptive immune cells (T and B lymphocytes) and innate immune cells [e.g., macrophages, neutrophils and dendritic cells (DCs)]. These immune cells with different functions communicate with each other directly or indirectly, and jointly affect the growth of tumor cells (17). In recent years, vigorous development of immunotherapy interventions has greatly improved the clinical outcomes in patients with KIRC (18). Therefore, immune infiltrating cells are considered promising targets for KIRC therapy.

However, few studies have deciphered the expression pattern, prognostic value, and immune infiltrating relevance of SMAD family genes in KIRC. The aims of this research were to determine SMADs’ prognostic values and their involvement in immune infiltration. In this study, we performed an integrated and systematic analysis including the expression pattern, prognostic value, and immune infiltrating relevance of SMAD family genes in KIRC by using The Cancer Genome Atlas (TCGA)-based gene expression data and online public platforms, and further explored their potential functions.

We present the following article in accordance with the REMARK reporting checklist (available at https://dx.doi.org/10.21037/tcr-21-178).

Methods

TCGA and cBioPortal database analysis

RNA-Seq data of KIRC patients were obtained from the TCGA database (https://portal.gdc.cancer.gov/), which was extracted for gene expression analysis and Gene Set Enrichment Analysis (GSEA) method. The data of SMADs genetic mutation in KIRC (TCGA, firehose legacy) were extracted from cBioPortal (https://www.cbioportal.org). Overall survival (OS) and disease-free survival (DFS) were set as the clinical endpoint to evaluate the role of SMADs in KIRC. OS means the estimation of time from start of treatment to death due to any cause. DFS indicates the percentage of people who are alive for a given period of time after diagnosis. The study was conducted in accordance with the Declaration of Helsinki (as revised in 2013).

Methylation analysis in UALCAN

UALCAN (http://ualcan.path.uab.edu) is a publicly available and powerful web portal for in-depth analysis of cancer OMICS data (TCGA and MET500) (19). In our study, the methylation level of SMADs was analyzed in KIRC by the “TCGA” dataset.

LinkedOmics

LinkedOmics (http://www.linkedomics.org/login.php), an easy-to-use and publicly available portal, provides an interactive platform for clinicians and biologists to obtain, compare and analyze tumor multi-omics data across 32 TCGA tumor types (20). We used the Linked Finder module to analyze the correlation of gene expression and methylation level of SMADs in KIRC. Pearson correlation coefficient was used in this analysis.

Gene correlation analysis in Gene Expression Profiling Interactive Analysis (GEPIA)

GEPIA (http://gepia.cancer-pku.cn/) is an online database based on data from 8,587 normal tissues and 9,736 tumors from TCGA and the GTEx Project (21). In our study, GEPIA was applied to analyze SMAD gene expression and prognostic value in KIRC. And correlation analysis of gene expression was performed with the “correlation analysis” module in GEPIA. Correlation analysis was performed with Spearman correlation coefficient.

Protein-protein interaction (PPI) network analysis

To understand the interaction relationship and molecular mechanism among differentially expressed SMADs, we constructed PPI networks using the STRING database (https://string-db.org/) and the GeneMANIA database (http://genemania.org).

Immune infiltration analysis in Tumor Immune Estimation Resource (TIMER)

TIMER database (https://cistrome.shinyapps.io/timer/) is used for analyzing immune infiltration based on RNA sequencing data from a variety of tumors (22). We explored the connection between the expressions of SMADs and six tumor-infiltrating immune cells (B cells, CD4+ T and CD8+ T cells, neutrophils, macrophages, and DCs). We also explored the expression correlation between differentially expressed SMAD genes and marker genes of immune cells, including tumor-associated macrophage (TAM), M1/M2 macrophage, DCs, neutrophils, Th1, Th2, Treg and natural killer cells through relevant modules. We specifically assessed the correlation between differentially expressed SMADs and macrophage polarity through the “correlation” module, and macrophage related markers are as follows: CD86/CSF1R for monocyte, IL10/VSIG4/PCGR2A for TAM, NOS2/PTGS2 for M1 macrophages and CD163/NS4A4A/MRC1 for M2 macrophages. Correlation analysis was performed with Spearman correlation coefficient.

GSEA

GSEA focuses on groups of genes that share biological function, common regulation or chromosomal location (23). In this study, GSEA was applied to explore the potential regulatory mechanisms of SMADs expression of KIRC (TCGA).

Statistical analysis

GraphPad Prism 6.0 or SPSS Version 20 software was used to perform the statistical analyses in our study. The Student’s t-test was used to determine the statistical difference in the two groups. The correlation of gene expression (GEPIA) and the connection between gene expression and immune cell infiltration or immune cell marker genes (TIMER) was evaluated by statistical significance as well as Spearman’s correlation.

Results

The expression pattern and methylation level of SMAD family genes in KIRC patients

Firstly, we used 530 KIRC tissues and 72 normal control tissues, as well as 72 paired KIRC tissues and corresponding adjacent normal tissues in TCGA to evaluate the expression pattern of SMADs. Consistently, we found that the expressions of SMAD1, 2, 3, 4 and 6 were significantly reduced in KIRC tissues than in normal control tissues, while expression of SMAD9 was significantly increased in cancer tissues (). And these differentially expressed SMADs were further verified in the GEPIA database (). DNA methylation, especially the excessive methylation of the promoters of tumor suppressor genes, can lead to transcriptional silencing, which has a crucial impact on the development of KIRC (24). In order to explore the methylation levels of these differentially expressed SMADs in KIRC, we used UALCAN to evaluate the methylation status of their promoter regions. As shown in , the promoters of SMAD1, 2, 3, 4, and 6 showed more methylation in KIRC tissues. Further, the correlation between expression patterns and the methylation status of SMADs was explored in KIRC. As presented in , the expressions of all differentially expressed SMAD genes were significantly negatively correlated with their methylation levels. Altogether, these data manifested that the changes in the methylation level of the SMAD promoter regions may be an important cause for regulating the expression pattern of SMADs in KIRC.

Figure 1

The expression pattern and methylation level of SMAD family genes in KIRC patients. (A,B) The expression patterns of SMAD1, 2, 3, 4, 5, 6, 7 and 9 in 530 KIRC patient tissues and 72 normal tissues as well as 72 paired KIRC tissues and corresponding adjacent normal tissues (TCGA). The expression patterns of SMAD1 (C), SMAD2 (D), SMAD3 (E), SMAD4 (F), SMAD6 (G), SMAD9 (H). Promoter methylation status of SMAD1 (I), SMAD2 (J), SMAD3 (K), SMAD4 (L), SMAD6 (M), SMAD9 (N) in 324 KIRC patient tissues and 160 normal tissues (UALCAN). The correlation between the expressions of SMAD1 (O), SMAD2 (P), SMAD3 (Q), SMAD4 (R), SMAD6 (S), SMAD9 (T) with their promoter methylation status (LinkedOmics). P<0.05 was considered significant (*, P<0.05; ***, P<0.001; **** P<0.001). SMAD, small mother against decapentaplegic protein; KIRC, kidney renal clear cell carcinoma.

The expression pattern of these SMADs was also assessed in the main tumor stages of KIRC patients. As shown in , we found that the expressions of SMAD1, 2, 3, 4 and 6 decreased significantly with the progression of KIRC tumors. This suggested that expressions of these SMADs had significant impacts on the inhibition of tumorigenesis and progression of KIRC.

Figure 2

The expression pattern of SMAD genes in the main pathological stages of KIRC patients. The expression patterns of SMAD1 (A), SMAD2 (B), SMAD3 (C), SMAD4 (D), SMAD6 (E), SMAD9 (F) in the main pathological stages of KIRC patients. P<0.05 was considered significant. SMAD, small mother against decapentaplegic protein; KIRC, kidney renal clear cell carcinoma.

The prognostic value of differential expression of SMAD genes in KIRC patients

Next, we investigated the prognostic value of SMADs in KIRC patients. By drawing survival significance maps according to the Cox proportional hazard ratio (HR) using the GEPIA database, we evaluated whether differential expressions of these SMADs have impact on OS and DFS of KIRC patients. We found that low expressions of SMAD1, 2, 3, 4 and 6 were significantly associated with poor OS of KIRC patients (), and low expressions of SMAD2, 3, 4 and 6 were significantly associated with poor DFS of KIRC patients (). The Kaplan-Meier plots further showed that KIRC patients with decreased SMAD1, 2, 3, 4, and 6 expressions had shorter OS (), and KIRC patients with decreased SMAD2, 3, 4, 6 expressions had shorter DFS (). Since SMAD1, 2, 3, 4 and 6 had better prognostic value, the following analysis mainly focused on these 5 SMAD genes.

Figure 3

Prognostic value of differential expression SMAD genes in KIRC patients. (A,B) Survival significant maps of differentially expressed SMADs showed the overall survival (OS) and disease-free survival (DFS) analysis results based on the Cox proportional hazard ratio (HR) through GEPIA. The OS curve of SMAD1 (C), SMAD2 (D), SMAD3 (E), SMAD4 (F), SMAD6 (G), and the DFS curve of SMAD2 (H), SMAD3 (I), SMAD4 (J), SMAD6 (K) in KIRC. The group cutoff choice for OS was the median. A log-rank test was used to estimate the difference in OS and DFS, and a value of P<0.05 was considered statistically significant. SMAD, small mother against decapentaplegic protein; KIRC, kidney renal clear cell carcinoma.

Genetic mutations, expression correlation, and interaction network analyses of SMAD genes in KIRC patients

Then, an overall exploration of the molecular properties of these 5 SMAD genes was performed. The cBioPortal online tool was applied to assess the genetic mutations of these SMADs. In total, 538 KIRC patients (TCGA, firehose legacy) were selected. exhibited that the missense mutations were found in SMAD2 and 4, and the truncating mutation was found in SMAD2 and 3. Next, the potential co-expression of SMADs was analyzed. We found that there was a moderate to high correlation among the expressions of SMAD1, 2, 3 and 4 (). Moreover, a PPI network analysis of these SMADs was conducted by the STRING database to assess the potential interactions among them. As expected, 5 nodes and 10 edges were obtained in the PPI network, indicating the close connection among these SMADs (). In addition, GeneMANIA’s PPI network results further showed that the functions of these SMADs are mainly related to transmembrane receptor protein serine/threonine kinase signaling pathway, BMP signaling pathway, TGF-β signaling pathway, SMAD protein complex assembly and SMAD binding, etc. ().

Figure 4

Genetic mutations, expression correlation, and interaction network analyses of SMAD family genes in KIRC patients. (A) Gene mutation information of SMAD genes in KIRC (cBioPortal). (B) Correlation heat map of SMAD genes in KIRC (GEPIA). (C) Protein-protein interaction network of SMAD genes (STRING). (D) The interaction network and function prediction of SMAD genes (GeneMANIA). SMAD, small mother against decapentaplegic protein; KIRC, kidney renal clear cell carcinoma.

SMAD1/2/3/4/6 expressions were correlated with immune cell infiltration in KIRC patients

Infiltrating immune cells in the tumor microenvironment is critical to modulate KIRC progression and has shown promising prognostic value (15). We further utilized the TIMER database to explore the relevance between SMAD1, 2, 3, 4 and 6 expressions and immune infiltrating cells (including B cells, CD8+ T cells, CD4+ T cells, macrophages, neutrophils, and DCs) in KIRC. Interestingly, we observed that SMAD1 expression was positively correlated with the infiltration of five immune cells except B cells (), and the expression levels of SMAD2, 3 and 4 were positively correlated with all infiltrating immune cells evaluated (). In addition, the expression of SMAD6 was positively correlated with the level of infiltrating CD8+ T and CD4+ T cells, while not significantly correlated with the levels of infiltrating B cells, macrophages, neutrophils and DCs (). These data suggested that the expressions of SMAD1, 2, 3, 4 and 6 were closely linked to the level of infiltration immune cells.

Figure 5

The expression of SMAD 1/2/3/4/6 was correlated with immune infiltration levels in patients with KIRC (TIMER). The correlation between the abundance of immune cells and the expression of SMAD1 (A), SMAD2 (B), SMAD3 (C), SMAD4 (D), SMAD6 (E) in KIRC. P<0.05 was considered significant. SMAD, small mother against decapentaplegic protein; KIRC, kidney renal clear cell carcinoma.

Correlation analysis between SMAD1/2/3/4/6 expressions and immune markers

To unravel the effects of the expressions of SMAD1, 2, 3, 4 and 6 on immune infiltration, we further analyzed the relationships between their expressions and multiple specific markers of innate and adaptive immune cells including monocyte, TAM, M1 macrophage, M2 macrophage, DCs, Th1, Th2, Treg, and natural killer cell in KIRC (). We found that SMAD1 and 4 expressions showed a significant correlation with the genetic markers of these infiltrating immune cells, and the expression of SMAD2 was correlated with most immune markers of infiltrating immune cells except for natural killer cells. In addition, the expression of SMAD3 had a significant correlation with immune markers of M1 macrophage, DCs and Th1, and the expression of SMAD6 was remarkably linked to the genetic markers of DCs, Treg and natural killer cells in KIRC. Particularly, we found that the expressions of SMAD1, 2 and 4 were strongly correlated with the levels of specific markers of monocytes, TAM, M1 and M2 macrophages in KIRC ( and ), which has also been further verified in the GEPIA database (). This suggested that SMAD2, 3, and 5 may be involved in regulating the polarity of macrophages in KIRC.

Table 1

Correlation analysis between SMAD1/2/3/4/6 and related markers of immune cells in TIMER

Description	Gene marker	SMAD1			SMAD2			SMAD3			SMAD4			SMAD6
		Cor	P		Cor	P		Cor	P		Cor	P		Cor	P
Monocyte	CD86	0.099	*		0.212	***		0.035	0.458		0.145	**		−0.101	*
	CD115 (CSF1R)	0.197	***		0.242	***		0.086	0.066		0.231	***		0.045	0.332
TAM	IL10	0.172	***		0.232	***		0.059	0.210		0.173	***		−0.035	0.454
	VSIG4	0.109	0.389		0.197	***		−0.026	0.577		0.115	*		−0.077	0.099
	FCGR2A	0.293	***		0.338	***		0.086	0.067		0.269	***		−0.017	0.716
M1 Macrophage	INOS (NOS2)	0.430	***		0.420	***		0.258	***		0.503	***		0.383	***
	PTGS2	0.333	***		0.334	***		0.122	**		0.238	***		0.077	0.098
M2 Macrophage	CD163	0.304	***		0.410	***		0.131	**		0.324	***		−0.062	0.185
	MS4A4A	0.283	***		0.350	***		0.052	0.264		0.281	***		−0.015	0.749
	MRC1	0.315	*		0.619	***		0.334	***		0.617	***		0.212	***
DCs	ITGAX	−0.068	0.142		−0.024	0.609		0.138	**		−0.009	0.854		−0.123	**
	CD1C	0.249	***		0.235	***		0.217	***		0.308	***		0.183	***
	NRP1	0.652	***		0.651	***		0.378	***		0.660	***		0.397	***
	THBD	0.391	***		0.382	***		0.158	***		0.370	***		0.353	***
Neutrophils	CCR7	0.082	0.078		0.094	*		0.085	0.070		0.042	0.370		−0.024	0.606
	ITGAM	0.139	**		0.237	***		0.136	**		0.212	***		−0.013	0.785
	CD59	0.629	***		0.648	***		0.393	***		0.666	***		0.227	***
Th1	STAT4	0.099	*		0.094	*		0.104	*		0.114	*		0.063	0.175
	TBX21	0.155	***		0.077	0.097		0.184	***		0.147	**		0.308	***
	CD4	0.170	***		0.213	***		0.115	*		0.173	***		0.026	0.573
Th2	CXCR4	0.329	***		0.170	***		0.081	0.081		0.169	***		0.015	0.756
	CCR4	0.193	***		0.343	***		0.264	***		0.340	***		0.029	0.535
	CCR8	0.026	0.585		0.144	**		0.121	**		0.055	0.240		−0.145	**
Treg	FOXP3	−0.155	***		−0.113	*		−0.035	0.451		−0.201	***		−0.180	***
	STAT5B	0.531	***		0.684	***		0.567	***		0.784	***		0.363	***
	TGFB1	0.331	***		0.254	***		0.179	***		0.198	***		0.192	***
Natural killer cell	KIR2DL1	0.135	**		0.057	0.219		0.039	0.403		0.098	*		0.218	***
	KIR2DL3	0.115	*		0.088	0.060		0.085	0.069		0.128	**		0.201	***
	KIR3DL1	0.147	**		0.098	*		0.076	0.102		0.145	**		0.231	***
	KIR3DL2	0.092	*		−0.022	0.632		−0.001	0.986		0.017	0.710		0.139	**

P<0.05 was considered statistically significant (*, P<0.05; **, P<0.01; ***, P<0.001). TIMER, Tumor Immune Estimation Resource; TAM, tumor-associated macrophage; SMAD, small mother against decapentaplegic protein; KIRC, kidney renal clear cell carcinoma.

Figure 6

Expressions of SMAD1/2/4 correlated with macrophage polarization in KIRC (TIMER). The correlation between the expression of SMAD1/2/4 and monocytes (gene markers: CD86 and CD115), TAM (gene markers: IL10, VSIG4 and FCGR2A), M1 macrophage (gene markers: INOS and PTGS2) and M2 macrophage (gene markers: CD163, MS4A4A and MRC1) infiltration levels was assessed. Scatterplots of correlations between SMAD1 expression and monocytes (A), TAM (B), M1 macrophage (C) and M2 macrophage (D) in KIRC. Scatterplots of correlations between SMAD2 expression and monocytes (E), TAM (F), M1 macrophage (G) and M2 macrophage (H) in KIRC. Scatterplots of correlations between SMAD4 expression and monocytes (I), TAM (J), M1 macrophage (K) and M2 macrophage (L) in KIRC. P<0.05 was considered significant. SMAD, small mother against decapentaplegic protein; KIRC, kidney renal clear cell carcinoma.

Table 2

Correlation analysis between SMAD1/2/4 and relate genes and markers of monocyte and macrophages in GEPIA

Description	Gene marker	SMAD1			SMAD2			SMAD4
		R	P		R	P		R	P
Monocyte	CD86	0.260	***		0.300	***		0.300	***
	CD115 (CSF1R)	0.360	***		0.340	***		0.360	***
TAM	IL10	0.320	***		0.320	***		0.340	***
	VSIG4	0.310	***		0.270	***		0.250	***
	FCGR2A	0.430	***		0.430	***		0.420	***
M1 macrophage	INOS (NOS2)	0.500	***		0.470	***		0.490	***
	PTGS2	0.350	***		0.320	***		0.280	***
M2 macrophage	CD163	0.320	***		0.290	***		0.260	***
	MS4A4A	0.400	***		0.400	***		0.410	***
	MRC1	0.620	***		0.660	***		0.670	***

P<0.05 was considered statistically significant (***, P<0.001). GEPIA, Gene Expression Profiling Interactive Analysis; TAM, tumor-associated macrophage; SMAD, small mother against decapentaplegic protein.

Functional enrichment analysis of SMAD1/2/3/4/6 in patients with KIRC

To better understand the regulatory functions of SMAD1, 2, 3, 4 and 6 in KIRC, the GSEA method was employed to explore their regulatory mechanisms. All samples were separated into the low group and the high group based on the mean value. The bubble diagrams () showed the signal pathways that were significantly enriched in the low expression groups of SMAD1, 2, 3, 4, and 6. We found that immune-related signaling pathways, including the primary immunodeficiency, antigen processing and presentation, natural killer cell-mediated cytotoxicity, and T cell receptor signaling pathway, were significantly enriched and negatively correlated with SMAD1, 2, 3, 4 and 6. Also, SMAD1, 2, 3, 4 and 6 were closely related to cytokine-cytokine receptor interaction, chemokine signaling pathway, NOD-like receptor signaling pathway, and Toll-like receptor signaling pathway. Furthermore, our findings showed that SMAD3, 4 and 6 were closely associated with the well-known pathways such as cell cycle regulation, p53 signaling pathway and JAK/STAT signaling pathway. These results implied that SMAD1, 2, 3, 4 and 6 had the potential to be therapeutic targets.

Figure 7

Functional enrichment analysis of SMAD1/2/3/4/6 in patients with KIRC. The bubble diagram displayed the GSEA results in phenotype low group of SMAD1 (A), SMAD2 (B), SMAD3 (C), SMAD4 (D), SMAD6 (E). The nominal P value <0.05 and false discovery rate q<0.25 were used to screen significantly enriched gene sets. SMAD, small mother against decapentaplegic protein; KIRC, kidney renal clear cell carcinoma; GSEA, Gene Set Enrichment Analysis.

Discussion

Genetic studies have indicated that the TGF- β/SMAD signaling pathway is involved in the regulation of tissue development and homeostasis and its dysregulation is a hallmark of many human diseases, including cancer (25). SMAD, as a family of TGF-β signal transduction molecules, has been reported to play important roles in cell proliferation, differentiation, migration, survival and apoptosis (5). Emerging evidence has proved that SMAD family genes are abnormally expressed in various cancers, including KIRC, and may be used as prognostic biomarkers or targets of chemical drugs (26). Among these SMADs, some SMAD family genes have been well studied in various cancers. For example, SMAD2 has been confirmed to be associated with tumor progression and development in patients with breast cancer, gastric cancer, esophageal squamous cell carcinoma and colorectal cancer (27-30). SMAD4 played a vital role in the tumorigenesis of diverse human cancers such as squamous cell carcinoma of the esophagus, head and neck squamous cell carcinoma and pancreatic ductal adenocarcinoma, etc. (31-33). SMAD6 has been found to be a prognostic marker for patients with esophageal cancer, oral squamous cell carcinoma, and bladder cancer (13,34,35). However, the prognostic value and immune infiltration relevance of SMAD family members, especially at the mRNA level, are not yet fully elucidated in KIRC.

In our study, we first explored the expression pattern of SMAD family genes in KIRC patients. We determined that the expressions of SMAD1, 2, 3, 4 and 6 were significantly reduced whereas the expression of SMAD9 was significantly increased in KIRC tissues than in normal control tissues. Further, we found that the expressions of SMAD1, 2, 3, 4 and 6 decreased considerably with the progression of KIRC tumors. Aberrant DNA methylation is almost universally found in cancer including KIRC, which can lead to silenced gene expression (36). Thus, we also explored the promoter methylation levels of these differentially expressed SMAD genes. And we found that the promoters of SMAD1, 2, 3, 4, and 6 presented higher methylation in KIRC tissues. Further, the correlation analysis displayed that the expressions of all differentially expressed SMADs were significantly and negatively correlated with their methylation levels. These results implied that aberrant methylation of these SMADs may be a primary reason for their abnormal gene expression.

Park et al. evaluated the influence of SMAD2, 3, 4 nuclear expressions on the prognosis of KIRC patients at the protein level, and found that low nuclear expressions of SMAD3 and SMAD4, but not SMAD2, were independent predictors of poor prognosis for KIRC patients (26). Consistently, we found that KIRC patients with low total expressions of the SMAD3 and SMAD4 also had a poor prognosis. On the contrary, we found that the low total expression of SMAD2 was also associated with shorter OS and DFS of KIRC patients. This may be due to the subcellular localization change of SMAD2 protein when it functions (37), so the total expression level of SMAD2 may not be correlated with their nuclear protein level. In addition, we found that the low total expressions of SMAD1 and SMAD6 were correlated with poor prognosis of KIRC patients, indicating that SMAD1 and SMAD6 had potential as prognostic markers.

Next, we evaluated the genetic mutations of SMAD family genes and found that there were missense mutations in SMAD2 and 4, and truncation mutations in SMAD2 and 3. Further, analyses of gene co-expression and PPI network demonstrated that there were a strong expression correlation and a close functional link among these differentially expressed SMAD members.

Immune infiltration of the tumor microenvironment is a crucial determining factor for immune response and is closely related to clinical prognosis in cancer including KIRC (15). Thus, we further evaluated the correlation between the differential expression of these SMADs and immune cell infiltration. In this study, we found that the expression levels of SMAD2, 3, 4 were significantly and positively related to the levels of infiltrating B cells, CD4+ T cells, CD8+ T cells, macrophage, neutrophils and dendritic cells. The expression of SMAD1 has a significant positive correlation with the infiltration of CD4+ T cells, CD8+ T cells, macrophages, neutrophils and dendritic cells. And the expression level of SMAD6 was remarkably positively correlated with the infiltrating levels of CD4+ T cells and CD8+ T cells. In addition, the correlation between SMAD1, 2, 3, 4, 6 expressions and gene expressions of immune cell markers further implicated their roles in regulating tumor immunology in KIRC. Especially, SMAD1, 2, 4 showed evident correlations with monocyte, TAM, M1/M2 macrophage, indicating that SMAD1, 2, 4 may have a vital role in regulating the polarity of macrophages. In summary, these findings indicated that SMAD1, 2, 3, 4, 6 may play an important role in the regulation and recruitment of immune infiltrating cells in KIRC.

Finally, we focused on the function of these SMADs using the GSEA method. As expected, we found that the functions of SMAD1, 2, 3, 4, 6 were relevant to immune-related pathways such as natural killer cell-mediated cytotoxicity, primary immunodeficiency, T cell receptor signaling pathway, etc. They were also closely related to cytokine-cytokine receptor interaction, chemokine signaling pathway, NOD-like receptor signaling pathway, and Toll-like receptor signaling pathway. In addition, SMAD3, 4 and 6 were found to be closely associated with cell cycle regulation, p53 signaling pathway and JAK/STAT signaling pathway. These data indicated that SMAD1, 2, 3, 4, 6 can be potential therapeutic targets in KIRC.

Conclusions

We comprehensively and systematically analyzed the expression pattern, prognostic value, and immune infiltrating relevance of SMAD family genes in KIRC, and further explored their potential functions. Our results indicate that SMAD1, 2, 3, 4, 6 are potential biomarkers to predict the prognosis and immune cell infiltration of KIRC patients. We hope our results will provide novel insights to identify prognostic biomarkers, assisting the prediction of the survival of KIRC patients, and helping clinicians choose effective treatment targets for KIRC patients. However, there are certain limitations to this study. Transcription levels analysis cannot completely reflect the functional changes of genes. Thus, further in vivo and in vitro investigations are needed to verify our findings and promote a better understanding of SMADs in KIRC.