Macrophage M2 Co-expression Factors Correlate With the Immune Microenvironment and Predict Outcome of Renal Clear Cell Carcinoma

Purpose: In the tumor microenvironment, the functional differences among various tumor-associated macrophages (TAM) are not completely clear. Tumor-associated macrophages are thought to promote the progression of cancer. This article focuses on exploring M2 macrophage-related factors and behaviors of renal clear cell carcinoma. Method: We obtained renal clear cell carcinoma data from TCGA-KIRC-FPKM, GSE8050, GSE12606, GSE14762, and GSE3689. We used the “Cibersort” algorithm to calculate type M2 macrophage proportions among 22 types of immune cells. M2 macrophage-related co-expression module genes were selected using weighted gene co-expression network analysis (WGCNA). A renal clear cell carcinoma prognosis risk score was built based on M2 macrophage-related factors. The ROC curve and Kaplan–Meier analysis were performed to evacuate the risk score in various subgroups. The Pearson test was used to calculate correlations among M2 macrophage-related genes, clinical phenotype, immune phenotype, and tumor mutation burden (TMB). We measured differences in co-expression of genes at the protein level in clear renal cell carcinoma tissues. Results: There were six M2 macrophage co-expressed genes (F13A1, FUCA1, SDCBP, VSIG4, HLA-E, TAP2) related to infiltration of M2 macrophages; these were enriched in neutrophil activation and involved in immune responses, antigen processing, and presentation of exogenous peptide antigen via MHC class I. M2-related factor frequencies were robust biomarkers for predicting the renal clear cell carcinoma patient clinical phenotype and immune microenvironment. The Cox regression model, built based on M2 macrophage-related factors, showed a close prognostic correlation (AUC = 0.78). The M2 macrophage-related prognosis model also performed well in various subgroups. Using western blotting, we found that VSIG4 protein expression levels were higher in clear renal cell carcinoma tissues than in normal tissues. Conclusion: These co-expressed genes were most related to the M2 macrophage phenotype. They correlated with the immune microenvironment and predicted outcomes of renal clear cell carcinoma. These co-expressed genes and the biological processes associated with them might provide the basis for new strategies to intervene via chemotaxis of M2 macrophages.


INTRODUCTION
Renal clear cell carcinoma (RCC) accounts for 80-90% of all renal cell carcinomas; clear cell carcinoma is not sensitive to chemotherapy and radiotherapy (Hsieh et al., 2017). For this reason, radical surgery has become the main treatment method. In clinical practice, although radical nephrectomy can benefit mostly patients, 30% of patients experience distant metastases after surgery (Motzer et al., 2013). Although we have adopted various treatment strategies for these patients with poor status, the long-term outcomes are not ideal (Linehan and Ricketts, 2019). With the development of immunotherapy in recent years, there have been studies showing that immunotherapy can benefit patients with renal clear cell cancer (Chowdhury and Drake, 2020;Díaz-Montero et al., 2020;Wang C. et al., 2020).
Renal clear cell carcinoma is characterized by many new tumor antigen peptides and high mutation burden; it is relatively sensitive to immunotherapies such as targeting PD1 and PD-L1 . Immune regulation plays a crucial role in the renal clear cell carcinoma microenvironment. This process includes immune checkpoints [mainly programmed cell death 1 (PD-1) and programmed cell death 1 ligand 1 (PD-L1)], as well as regulatory T cells, the original source of suppressor cell tumor-associated macrophages, and type 2 innate and adaptive lymphocytes (Xu W. et al., 2020). Macrophages in the primary or secondary tumor tissues are called tumorassociated macrophages (TAMs); these are the largest number of macrophages in the tumor stroma (Herberman et al., 1979). In recent years, clinical and experimental evidence has shown that macrophages promote the progression and metastasis of solid tumors, and this is somewhat different from our previous understanding (Pollard, 2004;Karnevi et al., 2014). Tumorassociated macrophages are divided into two types, M1 and M2 (Herberman et al., 1979;DeNardo and Ruffell, 2019). The biological effects of the two types are exact opposites. As tumors progress, increasing numbers of M2 macrophages appear, resulting in a weaker antigen presentation effect. For this reason, targeting macrophages has become a new therapeutic strategy (DeNardo and Ruffell, 2019). M1 type macrophages, namely, classically activated macrophages, highly express IL-12 and IL-23 that enhance antitumor effects (Lawrence and Natoli, 2011 contrast, M2 type macrophages, namely, alternatively activated macrophages, promote tumor formation and development (Cervantes-Villagrana et al., 2020). The mechanism of this polarization of macrophages is not clear. This article focuses on exploring the M2 macrophage-related genes in renal clear cell cancer, and constructing co-expression networks of M2 macrophages using the WGCNA method. The results of this paper revealed the underlying interaction mechanisms of M2 macrophage co-expressing factors and explained the role of M2 macrophages in the immune microenvironment from the perspective of bioinformatics.

Macrophage M2, Tumor Purity, and Tumor Mutation Burden Evaluation
We downloaded The Cancer Genome Atlas TCGA-KIRC FPKM data (http://cancergenome.nih.gov/) containing 539 renal clear cell cancer tissue samples and 72 normal tissues. GSE8050 (Weinzierl et al., 2008), GSE12606 (Stickel et al., 2009), GSE14762 (Wang et al., 2009), andGSE36895 (Peña-Llopis et al., 2012) were also downloaded from the GEO (http://www.ncbi.nlm.nih. gov/geo/) database. The Robust Multi-Array Average (RMA) algorithm of the "sva" (Leek et al., 2012) package was used to remove batch effects among the four GEO cohorts. The TCGA cohort was used to select M2-related genes. Four GEO cohorts were combined using "sva" packages and to verify the results. The Cell type Identification By Estimating Relative Subsets Of RNA Transcripts (CIBERSORT) is a deconvolution algorithm based on a gene expression profile that characterizes the cell composition of complex tissues, quantifies immune cells, and accurately estimates the immune components of tumor samples. It expands the potential of the genomic database, showing the pattern of Renal Clear Cell Carcinoma with comprehensive immune cells. We calculated macrophage M2 cell proportions based on the LM22 matrix using the CIBERSORT (Chen et al., 2018) algorithm, Cibersort was used as an obvious method to evaluate the significance of infiltration of immune cells in the samples. The assessment results of some samples were not statistically significant, and we used P < 0.05 to screen the samples. The Estimation of Stromal and Immune cells in Malignant Tumor tissues using Expression data (ESTIMATE) is a method that infers the fraction of stromal and immune cells using gene expression signatures (Yoshihara et al., 2013). Using the ESTIMATE package, we calculated tumor purity in each renal clear cell cancer sample. TMB (tumor mutation burden) per megabyte is calculated by dividing the total number of mutations by the size of the target coding region Yang et al., 2020).

Macrophage M2 Co-expression Network Conduction
Weighted gene co-expression network analysis (WGCNA) is a system biology approach that converts co-expression correlations FIGURE 1 | Flowchart of the experimental design. We first calculated immune infiltration to determine the content of M2 macrophages in the immune microenvironment of RCC. Then, we constructed a co-expression network related to M2 macrophages of RCC and analyzed the enriched pathways in this network. We then calculated the survival analysis of these co-expressed genes. We constructed a COX regression prognostic model associated with the co-expression genes of M2 macrophages in RCC and performed a subgroup analysis of this model. We analyzed the relationship between key genes in the model and tumor purity and CD8 + T cells. Finally, we also verified the feasibility of the model with 4 GEO datasets and conducted western blotting experiments on VSIG4.
into connection weights or topology overlap values (Langfelder and Horvath, 2008). We used this method to determine proportions of co-expressed genes in the M2 macrophage. The expression patterns are similar for genes with the same biological process and biological function (Jiang et al., 2017). We built a scale-free topology network, set the soft threshold at 5, R square = 0.89, and set the number of genes in the minimum module at 30. The M2 macrophage cell proportion was considered for phenotype files in WGCNA. In this manner, a cluster of M2 macrophage cell proportion-related genes with similar function were identified in the same module. The factors with M2 macrophage correlation >0.4 in the most relevant modules were determined.

M2 Macrophage-Related Module Analysis
The genes were selected using |correlation coefficient| > 0.4. The Database for Annotation, Visualization and Integrated Discovery (DAVID, v6.8) is an open-source database that performs function enrichment (Huang et al., 2007). We used the Kyoto Encyclopedia of Genes and Genomes (KEGG) (https://www. genome.jp/kegg/) (Kanehisa et al., 2017) and Gene Ontology (GO) (http://geneontology.org/) analysis (Ashburner et al., 2000) to identify the biological function in each co-expression module. In this way, we identified the biological processes associated with M2-type macrophage proportion.

M2 Macrophage-Related Genes Analysis
To verify the correlation between these factors and the clinical phenotype, we measured the overall survival from clear cell carcinoma as the prognostic indicator. Survival analysis was performed to evaluate the prognostic value of these co-expressed factors in M2 macrophages. Subsequently, a Cox regression hazard model was built based on the M2 macrophage-related genes. Next, we generated a model validation of clinical subgroups, which was based on age, gender, tumor metastasis, tumor stage, tumor purity, and degree of tumor mutation FIGURE 2 | (A) A hierarchical clustering tree was built using the dynamic hybrid cutting method, where each leaf on the tree represents a gene, and each branch represents a co-expression module; 21 co-expression models were generated. Frontiers in Genetics | www.frontiersin.org burden. In different subgroups, we evaluated the predictive abilities of M2 macrophage-related prognostic models. Finally, we calculated tumor purity in TCGA samples and explored the correlations between macrophage-related factors and tumor purity.

HPA
To verify the protein expression levels of candidate genes in melanoma and normal tissues, the human protein atlas (HPA, https://www.proteinatlas.org/) database was used to demonstrate differences in co-expressed genes at the protein level (Uhlén et al., 2015).

Western Blotting
Thirty clear renal cell carcinoma tissue samples were obtained from patients who underwent Nephrectomy at the First Affiliated Hospital of China Medical University. This study was authorized by the Ethics Committee of the First Affiliated Hospital of China Medical University. All patients signed informed consent. Protein exaction and western blotting were conducted as described previously (Pripp, 2018). An antibody against VSIG4 was purchased from Sigma-Aldrich.

Statistical Methods
Pearson correlation coefficients measure the strength of the linear relationship between two variables. The correlation coefficients are −1 to +1, respectively, indicating negative correlation and positive correlation, respectively, while 0 indicates no correlation . The key factors in the model score, tumor purity, tumor mutation burden, M2 macrophages, and CD8 + T lymphocytes were assessed using this test.

M2 Macrophages, Tumor Purity, and Tumor Mutation Burden
The results of our methodology are explained in Figure 1.
We summed up the following clinical data composed by M2 macrophages, tumor mutation burden, and clinical following survival data. M2, and M1, and M2/M1 macrophages were inputted as phenotype files to WGCNA. The detailed information is displayed in Supplementary Table 1.

M2 Macrophages Co-expression Network Conduction
We performed WGCNA analysis with TCGA-KIRC. A hierarchical clustering tree was built using the dynamic hybrid cutting method, where each leaf on the tree represents a gene, and each branch represents a co-expression module; 21 co-expression models were generated (Figure 2A). The correlation coefficients between each phenotype and co-expression module of TCGA are shown in Figure 2B. The results showed that the purple module had the strongest negatively correlation with M2 macrophage cell proportion in the TCGA-KIRC cohort (Cor = −0.45; P = 4e −15 ) and had the strongest correlation with CD8 + T cell proportion in the TCGA-KIRC cohort (Cor = 0.73; P = 6e −47 ) ( Figure 2B). Based on these findings, we have supplemented the scatter plots of the correlation between the factors in the purple module (Figures 2C-E). The horizontal axis is the correlation between the gene and the module, which is used to measure the relationship between the gene and the co-expression module, and the vertical axis is the correlation between the gene and the Frontiers in Genetics | www.frontiersin.org macrophage. By drawing the scatter diagram above, we screened out genes that are related to both M2 macrophages and the co-expression module.

M2 Related Genes Function Analysis
Twenty-four M2 macrophage negatively co-expressing genes were identified with coefficient <-0.4 in the TCGA-KIRC purple module. The gene significance for M2 macrophage-related genes in the purple module is shown in Table 1. Top 20 M2 macrophage cell proportion positively co-expressing genes were identified in the TCGA-KIRC pink module. The 24 M2 macrophage negatively co-expressing genes were most significantly enriched in the antigen processing and presentation of exogenous peptide antigen via MHC class I, which suggested a declining effect on the tumor antigen peptide process ( Figure 3A). The 20 M2 macrophage negatively co-expressing genes were most significantly enriched in neutrophil activation involved in immune responses ( Figure 3B).

M2 Related Genes Prognosis Analysis
To analyze their influence on overall survival, we performed survival analysis.   Frontiers in Genetics | www.frontiersin.org were prognosis-protective factors for clear renal cell carcinoma (Figure 4).

Immune Environment Correlation
Significant associations between M2 frequency and the genes involved in the risk signature are indicated in Figure 6, and the highest correlation of MFSD1 was 0.49 ( Figure 6A); the correlation of LCK was the lowest at −0.47 ( Figure 6B). TAP2, PSME2, HLA-E, and LCK were negatively related to M2 macrophage proportions. We then analyzed the correlations with CD8 + T cell and tumor mutation burden of these four genes. TAP2 (P < 0.001; Cor = 0.60), PSME2 (P < 0.001; Cor = 0.52), HLA -E (P < 0.001; Cor = 0.69), and LCK (P < 0.001; Cor = 0.69) ( Figure 7A) positively related to CD8 + T cell and negatively correlated with tumor purity (Figure 7B). This result suggested that M2 macrophages were negatively related to antigen processing.

HPA
The prognostic value and immune phenotype correlation were determined for these M2 macrophage-related genes. We compared the various expression levels of these genes between normal and tumor tissues. HPA001804 is an antibody against F13A1, which showed higher intensity in tumor tissue than in normal tissue. HPA056371 is an antibody against FUCA1, which showed higher intensity in the normal tissue than in tumor tissue. CAB012245 is an antibody against SDCBP, which showed a higher intensity in tumor tissue than in normal tissue. HPA003903 is an antibody against VSIG4, which showed higher intensity in tumor tissue than in normal tissue. HPA031454 is an antibody against HLA-E, which showed lower intensity in tumor tissue than in normal tissue. HPA001312 is an antibody against TAP2, which showed lower intensity in tumor tissue than in normal tissue. The protein levels  of these M2 macrophage genes were similar to the results of prognosis analysis at the transcription level (Figure 8). Subsequently, the M2 correlations for VSIG4, FUCA1, F13A1, SDCBP, HLA-E, and TAP2 were verified in the four GEO datasets (Supplementary Figure 1).

VSIG4
VSIG4 is thought to positively correlate with M2 macrophages; therefore, we conducted a combined analysis of VSIG4 and M2type macrophages. Combining VSIG4 elevated the predictive accuracy of M2 macrophages even more than either of them alone; the hazard of the "high VSIG4 expression + high M2 macrophage" group showed more survival risk than the other group (Kaplan-Meier analysis, low VSIG4 expression + low M2 macrophage; HR = 1.458; Figure 9A). Subsequently, we compared VSIG4 protein expression levels between normal renal tissues and clear renal cell carcinoma and found that VSIG4 protein expression levels in tumors were higher than the normal tissues ( Figure 9B). Then, various tumor infiltration deconvolution methods were applied; we found that VSIG4 was one of the most commonly associated M2 macrophage biomarkers ( Figure 9C).

DISCUSSION
In the tumor microenvironment, the chemotactic effects of the functional differences between the types of tumor-associated macrophages are not completely clear. The biological cytological role of M2/M1 macrophages in tumor tissues still needs to be explored. The present study is based on a bioinformatics algorithm to determine some of the M2 macrophage coexpression networks. Through the analysis of various modules, we tried to explain the biological function of co-expressed genes with M2 macrophages and related pathway changes from the perspective of bioinformatics. Our data processing and analysis processes are shown in the flowchart (Figure 1). F13A1, FCUA1, HLA-E, VSIG4, SDCBP, and TAP2 were the most common co-expressed genes in M2 macrophages. In terms FIGURE 8 | From the HPA database to verify protein expression-level differences of these candidate genes. Of these, F13A1, SDCBP, and VSIG4, and corresponding immunohistochemical samples, the degree of renal clear cell carcinoma tissue staining is higher than in normal kidney tissue. In FUCA1, HLA -E, and TAP2, and corresponding immunohistochemical samples, the degree renal clear cell carcinoma tissue staining is lower than in normal kidney tissue. These M2 macrophage gene protein levels at the transcription level were similar to those of the prognostic analysis.
of function enrichment, the 24 negatively co-expressed genes in M2 macrophages were most significantly enriched in antigen processing and presentation of exogenous peptide antigen via MHC class I. The 16 negatively co-expressing genes in M2 macrophages were most significantly enriched in neutrophil activation involved in immune response. M1 macrophages tend to adopt a Th1 response gene expression pattern and can secrete various cytokines that present MHC II and B7 molecules so as to present antigen efficiently (Herberman et al., 1979). This mechanism resists pathogen invasion, monitors tumor pathological changes, and generates Th1 immune responses in macrophages. By contrast, M2 macrophages have poor tumor antigen processing ability.
F13A1 encodes the coagulation factor XIII A subunit which has a catalytic function. In a human stem cell study, mRNA transcription expressed by F13A1 increased as myeloid progenitors differentiated into macrophages and erythroblasts (De Paoli et al., 2015). The protein encoded by FCUA1 is a lysosomal enzyme involved in the degradation of fucose-containing glycoproteins and glycolipids. Downregulation of FUCA1 enhances autophagy and inhibits macrophage infiltration so as to inhibit tumor growth (Xu L. et al., 2020). VSIG4 is a transmembrane receptor of the immunoglobulin superfamily that is specifically expressed in macrophages and mature dendritic cells. It is a newly discovered B7 family-related macrophage protein that inhibits T cell activation and has a potential role in cancer (Kim et al., 2016). VSIG4 negatively regulates macrophage activation by reprogramming mitochondrial pyruvate metabolism . HLA-E belongs to the HLA class I heavy chain paralogs. This class I molecule is a heterodimer consisting of a heavy chain and a light chain (beta-2 microglobulin). The heavy chain is anchored in the membrane. HLA-E binds a restricted subset of peptides derived from the leader peptides of other class I molecules. HlA-E is a non-classical HLA-I molecule that is best known for its role in protecting natural killer cells. Camilli et al. found that HLA-E was significantly increased during the differentiation of monocytes and macrophages (Camilli et al., 2016). The expression of HLA-E is related to the poor clinical results of anti-PD-1 immunotherapy. From the surface of M2 tumor-associated macrophages (TAMs), HLA-E antigen binds to the receptor CD94/NKG2A, which inhibits the expression of NK cell subpopulations and activated cytotoxic T lymphocytes, protecting cells from being destroyed (Marchesi et al., 2013). Epithelial-derived cancer cells, tumor macrophages, and CD141 + traditional dendritic cells promote the enrichment of HLA-E in carcinomas. CD8 + tumor-infiltrating T lymphocytes with high PD-1 content are prevented from surviving in the tumor microenvironment by the interaction of enriched HLA-E and CD94/NKG2A inhibition (Abd Hamid et al., 2019).
This study has some limitations, including lack of cross-validation of multicenter data. There is also lack of experimental verification of M2 macrophage biomarkers in renal clear cell cancer. We found that using the coexpression method of network-building, we can explicitly identify biomarkers, demonstrating the correctness of the logic based on bioinformatics.
In conclusion, we found that F13A1, FCUA1, HLA-E, VSIG4, SDCBP, and TAP2 were biomarkers of M2-type macrophages using a co-expression network of infiltrated immune cells, and we proposed six candidate-related factors. The biomarkers and related processes of M2 macrophages in the tumor microenvironment were explained from the perspective of bioinformatics, providing a strategy to explore the polarization of macrophages.

DATA AVAILABILITY STATEMENT
Publicly available datasets were analyzed in this study. This data can be found here: The TCGA-BLCA dataset used in this study were obtained from TCGA database (https://cancergenome.nih. gov/). GEO datasets used in this study were obtained from GEO database (https://www.ncbi.nlm.nih.gov/geo/).

ETHICS STATEMENT
The studies involving human participants were reviewed and approved by The First Affiliated Hospital of China Medical University. The patients/participants provided their written informed consent to participate in this study.

AUTHOR CONTRIBUTIONS
YW, KY, and JB designed the study. YW, KY, JLin, JLi, and JB analyzed and wrote the article.