Comprehensive Analysis of the Butyrate-Metabolism-Related Gene Signature in Tumor Microenvironment-Infiltrating Immune Cells in Clear Cell Renal Cell Carcinoma

A wealth of experimental evidence has validated that butyrate is capable of inhibiting tumorigenesis, while the potential role of butyrate metabolism in the tumor immune microenvironment (TIME) has been rarely explored. This study aims to explore the potential of butyrate-metabolism-related genes as prognostic biomarkers and their correlations with immune infiltrates in clear cell renal cell carcinoma (ccRCC) patients. Based on The Cancer Genome Atlas dataset (TCGA; n = 539), a total of 22 differentially expressed genes (DEGs) related with butyrate metabolism in ccRCC and normal samples were identified. Among them, a prognostic signature involving six butyrate-metabolism-related genes was created (Bu-Meta-GPS) in the training set (n = 271) and validation set (n = 268), and risk scores were calculated based on them. ccRCC patients with high-risk scores exhibited an unfavorable prognosis, high immunoscore, upregulated immuno-oncological targets (PD1, PD-L1, CTLA4, and CD19), and distinct immune-cell infiltration than those with low-risk scores. High-risk ccRCC patients without radiotherapy had a better survival rate than radiotherapy-treated patients. The negative regulation of cytokine production and cytokine-mediated signaling pathways was remarkably enriched in ccRCC patients with high-risk scores. A nomogram was then formulated to assess the overall survival (OS) of ccRCC patients. In summary, we illuminated the key role of butyrate metabolism in ccRCC TIME. The developed Bu-Meta-GPS was a sensitive predictive biomarker for the prognosis of ccRCC, which also provided new perspectives in improving immunotherapeutic efficacy.


INTRODUCTION
Renal cell cancer (RCC) has multiple subtypes, and among them, 70% of them belong to ccRCC (Jonasch et al., 2021). Early-stage ccRCC usually presents an acceptable outcome following surgery or ablation. However, one-third of ccRCC patients suffer metastases, which are highly lethal and different from non-metastatic lesions (Jonasch et al., 2014;Jonasch et al., 2021). So far, large-scale clinical trials on ccRCC have been conducted, although its OS has not been remarkably improved by any systemic adjuvant therapy (Ravaud, 2017;Gross-Goupil et al., 2018). Stromal infiltration, immune infiltration, and transformed cells constitute the complex environment of the tumor. Depending on the type of cancer or tumor model, tumor-infiltrating cells are capable of indicating tumor suppression or promotion based on the cancer type (Senbabaoglu et al., 2016). Previous evidence demonstrated that ccRCC is a highly immune-infiltrated tumor (Thompson et al.,2007;Yoshihara et al., 2013). ccRCC is also one of the earliest malignant tumors responsive to immunotherapy in history and the most responsive one currently (Topalian et al., 2012;Herbst et al., 2014). About 1% of ccRCC patients present spontaneous regression, which is attributed to the immune mediation (Janiszewska et al., 2013), while the underlying mechanisms of spontaneous remission, immune infiltration, and immunotherapy response are poorly understood. The high tumor mutational burden has triggered a great breakthrough in the immune checkpoint blocking treatment of tumors (Snyder et al., 2014;Rizvi et al., 2015). It is expected that more tumor mutations will result in more MHC-binding neoantigens that are responsible for mediating tumor immune infiltration and chemotherapy response (Igarashi et al., 2002;Gubin et al., 2014;Snyder et al., 2014;Yadav et al., 2014). Previous data have demonstrated exciting outcomes in metastatic ccRCC patients intervened with immune checkpoint inhibitors (Motzer et al., 2018;Ghatalia et al., 2019), and their potential benefits are being assessed. Therefore, it is of great significance to illustrate ccRCC TIME, which presents a potential prognostic value and assists decision-making of adjuvant therapy.
A growing amount of evidence has shown the close link between human gut microbiota and carcinoma via specific pathogens or a wide range of microbial communities, especially their metabolites (Louis et al., 2014). Butyrate is a typical metabolite with an anti-cancer role against inflammation and uncontrolled malignant behaviors of cancer cells as well as accelerating their apoptosis and differentiation (Andoh et al., 2002;Hinnebusch et al., 2002;Zeng and Briske-Anderson, 2005;Comalada et al., 2006;O'Keefe, 2016). A decreased fecal butyrate level serves as an indicator for cancer risk, progression, and severity (Hu et al., 2016). In addition, intracellular butyrate promotes hyperacetylation of histones by inhibiting histone deacetylase (HDAC) activities in cancer and immune cells, which also regulates proteins and transcription factors involved in key signaling pathways (Wilson et al.,2010;Fung et al., 2012). As a result, butyrate is able to mediate expression and differentiation of various genes and proinflammatory cytokines (Chang et al., 2014). A latest study revealed that butyrate performed anti-tumoral effects (Smith et al., 2013) by regulatory T cells inhibiting local inflammatory responses (Arpaia et al., 2013;Furusawa et al., 2013). Moreover, extracellular butyrate may interact with surface-exposed G proteincoupled receptors (GPRs) in host cells. Serving as tumor suppressors, GPR43 and GPR109A are involved in the anticancer role of butyrate associated with a high fiber intake (Ganapathy et al., 2013). Butyrate is found to downregulate primary miR-17-92a and precursor and mature miR-92a in cancer cells (Hu et al., 2015), which is attributed to the inhibited proliferative and accelerated apoptotic rate via downregulating the oncogene MYC and upregulating CDKN1C (also known as p57), respectively (Hu et al., 2011).
Recent studies have found that butyrate directly enhances the in vitro and in vivo anti-cancer response of CD8 + T cells via activating the IL-12 signaling, indicating the anti-cancer therapy potential of butyrate (He et al., 2021). In the metastatic liver of mice with colorectal cancer, butyrate intervention positively mediates the section of IL-10 and IL-17 via T regulatory cells and natural killer T cells and T helper 17 cells, respectively (Ma et al., 2020). However, it is still not fully understood whether the butyrate could regulate immune-infiltrate and how it modulates immune response in ccRCC.
In this work, the prognostic value of the butyrate metabolism was assessed via bioinformatics in the ccRCC dataset downloaded from TCGA. The risk score was calculated to construct Bu-Meta-GPS based on a prognostic signature involving six identified genes strongly correlated with ccRCC prognosis. The influences of risk scores on immune-oncological targets (PD1, CTLA4, and CD19), immunoscore, and immune cell infiltration were then assessed, thus revealing the effect of the butyrate metabolism on ccRCC TIME. A nomogram was formulated, aiming to assess the prognosis of ccRCC patients with varying clinical characteristics. Our findings are expected to guide ccRCC immunotherapies in clinical practice.

MATERIALS AND METHODS ccRCC Dataset
RNA-seq data of 539 ccRCC patients and 348 butyratemetabolism-related genes were obtained from TCGA (https:// portal.gdc.cancer.gov/) and the Gene Set Enrichment Analysis (GSEA) database (https://www.gsea-msigdb.org/gsea/index.jsp), respectively, and the latter were identified based on available mRNA expression data of the former.

Bioinformatics Analysis
DEGs related with butyrate metabolism between ccRCC and normal tissues were screened by the limma R package, with the false discovery rate controlled via p-value (FDR <0.05). The principal component analysis (PCA) was conducted using the R v4.1.0 package to identify gene-expression patterns between ccRCC and normal tissues. The tumor mutation burden (TMB) and somatic cell copy number alternations (CNAs) were calculated to determine the genomic variance.
Bu-Meta-GPS in the training set was created using the least absolute shrinkage and selection operator (LASSO) regression. The top six DEGs were selected for the following analysis. Risk scores were calculated according to expression levels and weighted coefficients obtained from the LASSO regression algorithm as follows: risk score = sum of coefficients × expression levels of six DEGs. High and low risks of ccRCC in the training and validation sets were classified using the median risk score as the cut-off.
ESTIMATE in R was adopted to estimate the immunoscore, stromal score, and ESTIMATE score-based prognostic signature in ccRCC TIME. Using the CIBERSORTx tool (https://cibersort. stanford.edu/), the infiltration levels of 22 types of immune cells were estimated, and those with CIBERSORT p-values of less than 0.05 were screened out for further analysis. Java-based GSEA was performed to identify the survival difference among ccRCC subtypes. Gmt files (go.v7.4. symbols.gmt) were analyzed as the reference gene set in the MSigDB database. A random permutation of 1,000 was conducted. Pathways with p < 0.05 and a normalized enrichment score (NES) were considered significantly enriched.

Statistical Analysis
Expression levels of butyrate-metabolism-related genes in ccRCC and normal tissues were examined using the Mann-Whitney U-test. Differences between groups and among multiple groups were compared using Student's t-test and one-way ANOVA, respectively. Kaplan-Meier curves were depicted to compare OS between groups. Pearson's correlation test was conducted to assess the correlation between risk scores, clinical features, immuneoncological targets, and immune infiltration levels. The independent prognostic potential of Bu-Meta-GPS was assessed by the univariate and multivariate Cox regression analyses. Receiver operating characteristic (ROC) curves and the area under the curve (AUC) were introduced to assess the prognostic value. A nomogram based on multivariate Cox regression analysis was created to determine the accuracy in predicting the prognosis. A significant difference was set at p < 0.05, which indicated statistical significance. GraphPad Prism 7.0 and R v. 4.1.0 were used to perform statistical analyses.

Expression Levels of Butyrate-Metabolism-Related Genes in ccRCC
A dataset involving 539 ccRCC samples and 79 adjacent normal ones was downloaded from TCGA. Differentially expressed Frontiers in Cell and Developmental Biology | www.frontiersin.org May 2022 | Volume 10 | Article 816024 3 butyrate-metabolism-related genes were identified, and finally, 22 genes were screened out for further explorations ( Figures 1A,B). To characterize the heterogeneity, PCA was performed to analyze mRNA levels of 22 butyrate-metabolism-related genes, which presented a significant difference ( Figure 1C). The aforementioned analysis indicated the potential effect of butyrate-metabolism-related genes on distinguishing normal samples from ccRCC.

Genomic Variance of Butyrate-Metabolism-Related Genes in ccRCC
A landscape of the incidence of somatic mutations of butyratemetabolism-related genes in ccRCC patients from TCGA dataset was obtained. In the available samples, 72/236 (30.51%) ccRCC patients experienced mutations of butyrate-metabolism-related genes (top 20 genes were shown), with a frequency ranging from 1 to 4% (Figure 2A). Among them, the highest number of mutations was detected in the COL6A3 gene (4%). In addition, the missense mutation was the most detected variant classification, and the main variant type was single-nucleotide polymorphism (SNP) (Supplementary Figure S1). Copy number variations (CNVs) were detected in 97/348 (27.87%) genes, and most of them were copy number amplifications ( Figure 2B). Notably, among the 22 DEGs, nine genes (40.91%) had CNVs, including ADH6, C7, STC2, DPP6, PLCB2, CD300C, TYMS, MAPK4, and ALPK2 ( Figure 2C). Thus, it is believed that CNVs may be responsible for DEGs in ccRCC.

Construction and Validation of Bu-Meta-GPS in ccRCC
A total of 539 ccRCC patients in TCGA dataset were randomly assigned to the training set (n = 271) and validation set (n = 268). LASSO regression analysis identified 6/22 DEGs in the training set, which were used to create Bu-Meta-GPS. Risk scores in both sets were then calculated by the following equation: Risk scores = 0.7145 × DLX4 level +0.0439 ×PLCB2 level +0.8463 ×FIRRE level +0.0752 × IL-11 level -0.0178 × ADH6 level +0.0026 × EIF4A1 level.
ccRCC patients were further subgrouped as high-risk and lowrisk groups. The risk score distribution, OS, and Bu-Meta-GPS signature in both training and validation sets are depicted in Figures 3A,B. Low-risk patients had a significantly longer OS than high-risk patients in both sets (p < 0.001, Figures 3C,E). To further evaluate the prognostic accuracy of the Bu-Meta-GPS, 1year, 3-year, and 5-year ROC curves were depicted, and the corresponding AUCs were calculated. It is shown that the 1-, 3-, and 5-year AUC values of the Bu-Meta-GPS in the training set were 0.699, 0.687, and 0.783, and those for the validation set were 0.793, 0.715, and 0.696 , respectively ( Figures 3D,F). These findings indicated that the Bu-Meta-GPS involving butyrate-metabolismrelated genes was a favorable prognostic tool for ccRCC.

Independent Prognostic Value of Bu-Meta-GPS and Construction of the Bu-Meta-GPS-Based Nomogram
The univariate Cox regression analysis illustrated that the risk scores calculated by DEGs related with butyrate metabolism were correlated with OS of ccRCC patients (HR, 4.536; 95%CI, 3.124-6.585; p < 0.001; Figure 4A). Furthermore, Bu-Meta-GPS was validated by multivariate Cox regression analysis as an independent prognostic factor for ccRCC (HR, 2.919; 95%CI, 1.909-4.463; p < 0.001; Figure 4B). The AUC value for the Bu-Meta-GPS was 0.733, which was higher than that for the age (0.556), gender (0.494), and N stage (0.530) ( Figure 4C).
A Bu-Meta-GPS-based nomogram was then created to test the 1-, 2-, and 3-year OS using risk scores and other clinical factors like the age, gender, grade, stage, and T and N stages ( Figure 4D), which were consistent with the ones predicted by the nomogram calibration curves ( Figure 4E).

Bu-Meta-GPS Correlated With Clinical Features and Radiotherapy in ccRCC
We further estimated the influence of risk scores on clinical features of ccRCC. Expression levels of the Bu-Meta-GPS classified by the risk level of the training set were depicted in a heatmap ( Figure 5A). IL-11, PLCB2, EIF4A1, DLX4, and FIRRE were remarkably upregulated in high-risk patients than those of the other group, while ADH6 was downregulated. We further examined the correlation of Bu-Meta-GPS with gender, age, stage, grade, and TNM staging in ccRCC patients, which revealed that the risk scores were positively correlated with tumor grade, stage, and TNM staging (p < 0.05, Figure 5B). Compared to female patients, male patients had high risk scores (p = 0.036, Figure 5B). It is suggested that Bu-Meta-GPS was significantly correlated with clinical features in ccRCC patients. We next investigated the potential influence of risk scores on the prognostic performance in ccRCC patients subjected to radiotherapy. The results showed that high-risk ccRCC patients who did not receive radiotherapy had a better OS ( Figure 5C). Although we did not obtain a significant difference in the survival benefit of radiotherapy in low-risk patients, it remained similar to that of the high-risk group ( Figure 5D). Accordingly, high-risk ccRCC patients without radiotherapy had a better survival rate than radiotherapytreated patients.

Correlation Between Bu-Meta-GPS and Distinct Immune Cell Infiltration in ccRCC
The difference in immune infiltrate levels classified by risk levels in ccRCC patients was examined, thus revealing the role of Bu-Meta-GPS in ccRCC TIME. Significantly higher immunoscore, stromal score, and ESTIMATE score were detected in the former (p < 0.01, Figure 6A). Subsequently, the fraction of 22 immune cell types between the two groups was analyzed. Higher infiltration levels of follicular helper T cells (Tfhs), regulatory T cells (Tregs), and activated mast cells were yielded in the highrisk group, while gamma delta T cells (γδ T cells), resting dendritic cells, and resting mast cells were mainly pronounced in the low-risk group ( Figure 6B). In addition, we calculated the relationship between the immune cell types. A positive correlation was obtained between the number of Tfhs and CD8 + T cells ( Figure 6C).

Correlation Between Immuno-oncological Targets and Bu-Meta-GPS in ccRCC
To assess the correlation between immuno-oncological targets and Bu-Meta-GPS, their expression levels which differed from ccRCC risk were examined. PD1, PD-L1, and CTLA4 were remarkably upregulated in ccRCC tissues than those of controls (p < 0.01; Figure 7A). Compared with those in low-risk patients, PD1, CTLA4, and CD19 were upregulated in the high-risk group (p < 0.01; Figure 7B). In addition, PD1, CTLA4, and CD19 levels were positively correlated with those of PLCB2, EIF4A1, DLX4, and FIRRE. The PD-L1 level was positively correlated with those of PLCB2 and EIF4A1 but negatively correlated with DLX4 ( Figure 7C). GSEA was performed to illustrate the mechanism underlying the differences in the TIME and immuno-oncological targets between groups. It is shown that the cytokine-mediated signaling pathway and negative regulation of cytokine production were remarkably correlated with high-risk ccRCC patients (NES = 1.72 and 1.64; normalized p < 0.020 and 0.034, respectively; Figure 7D). Therefore, the cytokinemediated signaling pathway may be attributed to different TIMEs and immuno-oncological targets classified by ccRCC risk.

DISCUSSION
Previous evidence has shown the anti-cancer effect of butyrate on multiple types of tumor cells via varying mechanisms (O'Keefe, 2016). It induces mitochondrial fusion and apoptosis in colorectal cancer cells via DRP1 (Tailor et al., 2014). Butyrate also inhibits migration of bladder cancer cells via inducing autophagy and apoptosis . Additionally, butyrate triggers apoptosis by overproducing reactive oxygen species (ROS) and causing mitochondrial dysfunction in breast cancer cells (Salimi et al., 2017). Multiple functions of butyrate in cancers indicate an extremely complicated regulatory mechanism. At present, cancer-intrinsic oncogene pathways have been well concerned. To clarify the potential regulatory mechanism of butyrate metabolism in TIME, it is necessary to fully analyze the function of butyrate. So far, the research on butyrate metabolism in ccRCC TIME is scant. We have identified 348 DEGs related with butyrate metabolism and their prognostic values and effects on ccRCC TIME, which were capable of distinguishing tumor tissues from normal tissues. Among them, expression levels of 22 genes were significantly altered, including upregulated MCAM, CXCR4, ALPK2, STC2, PYGL, TYMS, DLX4, FIRRE, PLCB2, CD300C, THBS4, IFI30, and EIF4A1 and downregulated C7, IL-11, MAPK4, DPP6, CNKSR1, LAD1, MT1H, G6PC, and ADH6, in ccRCC samples. Solid evidence suggested that these 22 genes functioned either as oncogenes or tumor-suppressor genes in various types of tumors. It is concluded that dysregulated genes related with butyrate metabolism may be important in ccRCC.
The prognostic value of Bu-Meta-GPS involving six DEGs related with butyrate metabolism in ccRCC patients was assessed and then validated, which effectively stratified the risk level of ccRCC. Notably, high-risk ccRCC patients in the training and validation sets had a worse OS than those of the other group. Bu-Meta-GPS was also identified to be correlated with clinical features of ccRCC. The risk scores calculated by Bu-Meta-GPS were positively correlated with the grade, stage, and TNM in ccRCC patients, which was consistent with the finding of a worse prognosis in high-risk patients. Furthermore, the prognostic role of Bu-Meta-GPS in ccRCC was confirmed. Compared with clinicopathological features, Bu-Meta-GPS presented a superb accuracy in predicting the prognosis of ccRCC. In addition, a nomogram combining Bu-Meta-GPS and clinical features was created and validated as an effective tool to predict the survival of ccRCC. We found that Bu-Meta-GPS was also able to predict radiotherapy sensitivity in ccRCC patients, of whom patients with low-risk scores gained more benefits from radiotherapy. Current studies revealed that radiotherapy caused immunosuppressive effects in the tumor microenvironment due to the attraction of immunosuppressive cells such as regulatory T cells (Tregs) as well as because of the release of immunosuppressive cytokines (TGF-β and IL-10) and chemokines (Shevtsov et al., 2019;Sato et al., 2020). It is suggested that radiotherapy should be carefully considered for high-risk ccRCC patients.
In this study, the influence of Bu-Meta-GPS on ccRCC TIME was demonstrated. Risk scores of high-risk and low-risk ccRCC patients were correlated with distinct immune cell infiltration levels, immunoscore, and expression levels of immunooncological targets. The stromal score and immunoscore were significantly higher in high-risk ccRCC patients, indicating that the high immune/stromal score patients had poorer survival outcomes. Classified by cell types, the infiltration levels of Tfhs, Tregs, and activated mast cells in high-risk ccRCC patients were higher than those in the low-risk group, while those of γδ T cells, dendritic cells, and resting mast cells presented opposite trends. We considered that the immunoscore and immune cell infiltration were the main causes for different survival in ccRCC patients stratified by risk levels. In addition, higher levels of PD-1, CTLA4, and CD19 were detected in highrisk patients, while the PD-L1 level was comparable. PD-L1 was significantly upregulated in ccRCC tissues compared to that in controls. Previous studies have demonstrated vital functions of CTLA4 and PD-L1 in tumorigenesis, tumor immunity, and prognosis (Zheng et al., 2019;Liu et al., 2020). Differentially expressed CD274 and CD19 genes are reported to be associated with the prognosis of ccRCC . Our study also found that risk scores were significantly correlated with expression levels of PD-1, PD-L1, CTLA4, and CD19, suggesting that the four immune checkpoints may be potential immunotherapy targets for ccRCC. Yukihiro et al. reported that butyrate ameliorates the development of colitis via inducing the differentiation of Tregs (Furusawa et al., 2013). Yang et al. (2020) showed that butyrate promotes IL-22 production in human CD4 + T cells. These findings suggested that butyrate metabolism mediated ccRCC TIME. GSEA results demonstrated that the cytokine-mediated signaling pathway was mainly enriched in high-risk ccRCC patients, suggesting that vital genes involved in this pathway may be vital targets for butyrate metabolism. Considering this, butyrate metabolism and cytokine-mediated signaling pathways were responsible for different ccRCC TIMEs between groups. Several limitations in the present study should be noted. First, further verification failed due to a small sample size, which should be further validated in multi-center large-scale trials. Second, the regulatory mechanism of butyrate metabolism in ccRCC TIME requires an in-depth exploration.
In conclusion, we created a prognostic signature involving butyrate-metabolism-related genes and evaluated their influence on ccRCC TIME, which provided prognostic markers and immunotherapy targets for ccRCC.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found at: https://portal.gdc. cancer.gov/.

AUTHOR CONTRIBUTIONS
GT and WY conceived and designed the study. GT and HG carried out the acquisition, analysis, and interpretation of data. GT, HG, and ZD wrote the manuscript. WY supervised the study.  GSEA showed that the cytokine-mediated signaling pathway was differentially enriched in the high-risk group (NES, normalized ES; *p < 0.05; **p < 0.01; ***p < 0.001; ns, not statistically significant).