NOD2 deficiency confers a pro‐tumorigenic macrophage phenotype to promote lung adenocarcinoma progression

Abstract Nucleotide‐binding and oligomerization domain‐containing protein 2 (NOD2) was a member of the NOD‐like receptor family and played an important role in the innate immune response. Dysregulated NOD2 had been reported to contribute to tumorigenesis and progression. Here, we investigated that decreased NOD2 expressions could affect the phenotypic polarization of tumour‐associated macrophages and thus lead to the poor prognosis of lung adenocarcinoma patients. We clustered the patients by the single‐sample gene set enrichment analysis of tumour microenvironment and 13 prognostic differentially expressed immune‐related genes (PDEIRGs) were obtained based on prognostic analyses. After multiple assessments on the 13 PDEIRGs, NOD2 was considered to be the central immune gene and had a strong effect on suppressing tumour progression. Decreased NOD2 expression could be induced by cancer cells and lead to the phenotypic polarization of macrophages from protective M1 phenotype to pro‐tumorigenic M2 subtype which might be attributed to the down‐regulating of NF‐κB signalling pathway. This study draw attention to the role of inhibited innate immune function mediated by depletion of NOD2 in the TME. Our work also points to a potential strategy of NOD2‐mediated TAM‐targeted immunotherapy.

induce an immunosuppressive effect in particular through TGFβ pathway. Therefore, the two phenotypes may have opposite effects on the prognoses of cancer patients. 8 M1 macrophages function in immune surveillance, whereas M2 macrophages, in the TME, are closely related to bad clinical prognosis in many kinds of human cancers. 9 Dissection of the roles of TAMs in tumour progression can pave the way to emerging TAM-targeted therapeutic strategies. 10 Here, we performed an immunological analysis of the TME (immunoscore) by clustering the immune components of LUAD samples in The Cancer Genome Atlas (TCGA) database and identified that NOD2 might be a suppressor for LUAD. NOD2 is a cytosolic receptor belonging to the NOD-like receptor (NLR) family that initiates innate immune in response to bacterial peptidoglycan (PGN)-conserved motifs. It has been reported that macrophages could exhibit a phenotypic polarization by up-regulating NLR expressions to adapt to their surrounding microenvironment. Our results revealed downregulated NOD2 in macrophages induced by lung cancer cells could impel the phenotypic conversion of TAMs from the protective M1 phenotype to the pro-tumorigenic M2 subtype. These findings provided a clue to develop an optional TAM-targeted immunotherapeutic strategy for LUAD treatment.

| Differentially expressed genes analysis
The differentially expressed genes (DEGs) were obtained from a comparison between HIC and LIC using the R package limma. 11 The cut-off values were false discovery rate (FDR) < 0.05 and |log2 fold change (FC)| > 1. Differentially expressed IRGs (DEIRGs) were obtained using the intersection between DEGs of HIC plus LIC and IRGs from ImmPort database.

| Function and pathway enrichment analyses
The analyses of function and pathway enrichments of DEGs and DEIRGs were performed by R package cluster profiler 12 with a strict cut-off of FDR < 0.05. Gene Ontology (GO) terms and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways of NOD2 were obtained by GSEA, 13 and the cut-off value was FDR < 0.05 based on 10,000 permutations.

| The risk score calculation
We cited the regression coefficients (coef) of the multivariate Cox hazard regression analysis and the FPKM (fragments per kilobase million) values of 13 PDEIRGs to create risk a score formula for each patient: Here, coef (i) represented the regression coefficient of the ith selected PDEIRGs and FPKM (i) represented gene expression value of the ith selected PDEIRGs.

| Immunohistochemistry
Immunohistochemistry was performed following the protocol pro-     Table S1. The relative levels of target genes were normalized to GAPDH by 2 ΔΔCt method.

| Western blot
Protein extraction and Western blot were performed as we previously reported. 14

| Statistical analyses
For the comparisons of two groups, unpaired Student t-test was used for the variables of normal distribution, and Mann-Whitney U test was used to analyse the variables of non-normal distribution. Non-parametric test, Kruskal-Wallis test, parametric test and one-way ANOVA were used for the comparisons of three or more groups. 15 The Pearson correlation coefficient test was used to estimate the rank correlation among the different variables. 16 The timedependent receiver operating characteristic (ROC) analysis was used to evaluate the accuracy of the prognostic model and calculate area under the curve (AUC). 17,18 p Values, two-sided, of less than 0.05 were considered statistically significant.

| Tumour immune microenvironments of LUAD patients from TCGA database were evaluated by ssGSEA and three immune clusters were established
To evaluated the TIME of LUAD patients in TCGA cohort, we quantified the enrichment levels of 29 immune-associated gene sets (Table S2) including different immune cell types, functions, and signalling pathways by ssGSEA and divided the TCGA patients into three immune clusters according to the ssGSEA scores. 19,20 The three immune clusters were defined as high-immunity cluster (HIC), middle-immunity cluster (MIC) and low-immunity cluster (LIC). Obvious enrichment differences were shown among three clusters ( Figure 1A). Additionally, we applied the index of TME scores to confirm the results of ssGSEA clustering. Being consist with ssGSEA scores, HIC patients had the highest ESTIMATE scores including stromal scores and immune scores, whereas LIC patients had the highest tumour purity scores ( Figure  Moreover, we calculated the expressions of three sets of marker genes including immunotherapy target gene, PD-L1 ( Figure 1D), genes for innate immune ( Figure 1E) and genes for specific immune ( Figure 1F). Box plots were plotted to show that the expressions of all the marker genes ascended from LIC to HIC apparently.
In general, TIME of LUAD patients could be clearly distinguished by ssGSEA. There were the most significant differences between HIC and LIC. LIC patients not only have fewer immune cells proportions, but also have lower levels of immune gene expressions.

| 13 PDEIRGs between HIC and LIC were obtained through Cox regression hazards analyses
To explore the immune genes that were most closely related to the prognosis of LUAD patients, we filtered all the DEGs between HIC and LIC because of the most remarkable discrepancy between the two groups. To begin with, an obvious different tendency of total transcriptomic gene expressions could be observed in the heat map Next, to identify DEIRGs, we utilized 1811 IRGs obtained from the ImmPort Database (Table S3) to intersect with total 1572 DEGs.
As a result shown in the venn diagram ( Figure 2E), 21 332 DEIRGs were gained totally. Compared with genes in the LIC, 309 genes were up-regulated and 23 genes were down-regulated ( Figure 2F and G). Again, we performed GO and KEGG functional enrichment analyses of 332 DEIRGs ( Figure 2H and I). We found that there was a high functional enrichment of cytokine production and cytokinecytokine receptor interaction. Here, we inferred that cytokine production and the interaction between cytokines and immune cells might play a key role in TIME.
To investigate the relationship between expressions of 332 DEIRGs and prognosis of LUAD patients, we involved 332 DEIRGs into Cox regression hazards analyses. Firstly, we applied a univariate Cox hazard regression analysis to screen all the DEIRGs of patients in HIC and LIC. Thirty genes with significant p value were established as candidate PDEIRGs ( Figure 3A). Then, we utilized the Lasso regression to eliminate five PDEIRGs that might over-fit the model (λ = 25, Figure 3B and C). Next, the reserved 25 candidate genes were included in the multivariate Cox hazard regression analysis for further assessment. Finally, 13 PDEIRGs were established including genes associated with antigen processing and presentation such as SLC10A2, THBS1, ERAP2 and genes taking part in anti-microbial such as S100P, NOD2 and LCN15 ( Figure 3D).
Furthermore, we used the 13 PDEIRGs to build up an IPM. The risk scores of each patient were calculated according to our formula mentioned before. All the patients were divided into the high-risk score group and the low-risk score group ( Figure S1A). High-risk patients had much shorter survival times ( Figure S1B). In addition, the expressions of 13 PDEIRGs in two groups were displayed by the heat map ( Figure S1C). The Kaplan-Meier curve was used to show that the significantly higher survival rates in the low-risk group (p = 2.116E-07, Figure S1D). After that, the time-dependent ROC curve with the AUC value of 0.709 was plotted to confirm the predictive accuracy of the model for LUAD patients ( Figure S1E).  Figure 4D). Likewise, lower NOD2 gene expression was found in low-risk score cluster which was thought to have poor prognosis ( Figure 4E). After that, since the multivariate Cox analysis had shown that TNM stage was an independent prognostic factor, we es-

| The expressions of NOD2 in macrophages were decreased after co-culturing with LUAD cells
Since our previous results suggested that cytokine production and the interaction between cytokines and immune cells might play a key role in TIME ( Figure 2H and I), we made further efforts to discuss the relationships between NOD2 and immune cells. A total of 422 TCGA LUAD patients were divided into two groups according to median expression level of NOD2. The proportions of 22 kinds of immune cells of patients in the two groups were statistically analysed.
As shown in Figure 5A, decreased NOD2 expressions were found in patients with low monocytes and macrophage fractions (p < 0.05).
Furthermore, we analyse the correlations between NOD2 expressions and CIBERSORT algorithm results of TAMs. It was proved that NOD2 expressions were positively related to TAMs ( Figure 5B-E).
Moreover, the results from immunohistochemistry showed that there was the same immune-positive signal for NOD2 and CD68, a macrophage marker molecule, in tumour tissues ( Figure 5F and G).
There was an obvious decline in NOD2 expressions in primary lung cancer tissues than those of normal lung tissues. In order to further confirm these results, we used the single-cell RNA sequencing (scRNA-Seq) data of 58 LUAD patients from GEO database to make a conjoint analysis. 22 The NOD2 expression values were calculated in the samples of 11 distant normal lungs, 11 primary tumours and 10 metastatic brain tissues. The decreased total and average NOD2 expressions in TAMs were found in both malignant tissues than in normal tissues ( Figure 6A and B). To explore whether the depletion of NOD2 in TAMs was induced by cancer cells, we co-cultured either HBEpiC or LUAD cells with THP-1 ( Figure 6C). After 24 h and 48 h co-culture, lower NOD2 expressions were found in THP-1 cocultured with cancer cells (Figure 6D-F). These results suggested that the expressions of NOD2 in monocytes and macrophages were decreased when these monocytes and macrophages were recruited into the TME.

| Tumour cells spurred the conversion of the protective M1 phenotype to pro-tumorigenic M2 subtype by down-regulating NOD2 expression in TAMs
It had been well reported that macrophages could be differentiated into M1 or M2 phenotype under microenvironment stimulus. 23 M1 macrophages had a pro-inflammatory and antitumour effect, whereas M2 macrophages could secrete cytokines and chemokines to enhance proliferative and pro-metastatic effect on the tumour cells. 24 TCGA data showed that there were fewer M1 and more M2 proportions in advanced stage patients than early-stage patients ( Figure 7A), which was consistent with the changeable tendency of NOD2 expressions ( Figure 4F). To further explain whether there was To further confirm these findings, we first used Q-PCR to examine the mRNA levels of phenotypic markers in THP-1 cells within co-culture system ( Figure 6C). After co-cultured with cancer cells, THP-1 cells expressed less M1 markers and more M2 markers ( Figure S2). Next, after 48 h of transfection with NOD2 siRNAs, the expressions of NOD2 in THP-1 cells were significantly knocked down ( Figure S3). Subsequently, LPS, a M1 stimulus, or IL-4, a M2 stimulus, was applied for the cellular co-culture system, respectively ( Figure 7G). Compared to the control group, knockdown of NOD2 could inhibit the expression of M1 markers under LPS stimulation. Conversely, the expression of M2 markers was obviously increased after addition of IL-4 to the cellular co-culture system ( Figure 7H-K). Taken together, our data suggested that the reduction of NOD2 in macrophages might be involved in the phenotypic conversion from the protective M1 phenotype to protumorigenic M2 subtype.
In addition, to answer the possible molecular mechanism of phenotypic conversion mediated by NOD2, we performed the GSEA using TCGA data. As shown in Figure S4A, NF-κB signalling pathway was significantly enriched (p < 0.05). Since it had been reported that NF-κB pathway could be activated through CARD-CARD domain interactions between NOD2 and RIPK2, we again analysed the previous scRNA-Seq data of TAMs to evaluate the expression values and correlations of RIPK2. Consequently, a sharp decline of RIPK2 expressions in the low NOD2 group and a strong correlation was obtained and shown in Figure S4B and C. Further, we used NF-κB pathway antibody kit to evaluate the activation of the NF-κB pathway in LPS-stimulated NOD2-silencing THP-1. After NOD2 expression was knocked down by SiRNA2 and SiRNA3, the phosphorylation of NF-κB pathway inhibitors, IKKα/β (Ser176/180) and IκBα (Ser32), was significantly decreased. Moreover, the activation of transcription factor p65/RelA (Ser536) was obviously impeded ( Figure S4D). It appeared that descended NOD2 in TAMs impeded the formation of the active NOD2-RIPK2 complex and then down-regulated the NF-κB signalling pathway by which means, at least partly, macrophages converted from M1 phenotype to M2 subtype.

| D ISCUSS I ON
Serving as a critical role in cancer immune escape leading to tumour growth and aggressiveness, 25 the TME signature was a robust biomarker for predicting survivals of tumour patients and providing effective immunotherapy strategies. 26 into an activated functional status. 30 On the other, the cancer cells in the TME could reconstruct IRGs expression patterns by mimicking normal cells to avoid the destruction which was called cancer immunoediting. 31,32 In this study, we were the first to create the IPM based on the ssGSEA clustering and identified a tumour inhibitor, NOD2 was a member of the NLR family and played an important role in both innate and adaptive immune response, apoptosis, autophagy and reactive oxygen species generation. 33 The unbalanced level of NOD signalling was associated with diseases by breaking immune homeostasis. 34 It had been previously shown that a lack of NOD2 signalling could increase Crohn's disease susceptibility using mouse models. 35 Recently, NOD2 had been connected to cancer development and treatment. 36 48 Serving as one of intracellular NTRs, NOD2 contained many similarities to TLRs, membrane-bound pattern recognition receptors (PRPs), in the processing of recognizing pathogen-associated motifs and inducing inflammatory signalling cascades. 49  Prospective studies should be further conducted and the exact roles and mechanisms of NOD2 on TAMs in the TME should be deeply investigated to provide benefit strategies of TAM-targeted immunotherapy for LUAD patients.
In summary, we clustered the TIME of LUAD patients by ssGSEA and performed an immune gene prognosis analysis. Our data supported the notion that decreased NOD2 leads to an inhibited innate immune status by mediating TAMs phenotypic conversion in the TME which contributed to tumour progression and poor prognosis of LUAD patients. Notably, this study provides a potential strategy to boost TAM-targeted immunotherapies in order to bring clinical interest for LUAD patients.

ACK N OWLED G EM ENTS
This study was supported by the National Natural Science Foundation of China (81672881) and Department of Science and Technology of Liaoning Province (20180530020).

CO N FLI C T O F I NTE R E S T
The authors have no conflicts of interests and declare no competing financial interests.

AUTH O R CO NTR I B UTI O N S
Yibei Wang: Formal analysis (equal); Software (equal); Visualization

DATA AVA I L A B I L I T Y S TAT E M E N T
The data sets used and/or analysed during the current study are available from the corresponding author on reasonable request.