Integrated transcriptional analysis reveals macrophage heterogeneity and tumor-macrophage interactions in the progression of pancreatic ductal adenocarcinoma

Pancreatic ductal adenocarcinoma (PDAC) is a devastating metastatic disease for which better therapies are urgently needed. It highlights the need for an improved understanding of the biological features underlying the progression of PDACs. Macrophages signicantly enhance tumorigenesis and development of pancreatic cancer. However, cellular reprogramming and phenotypic changes of macrophages during PDAC progression remain poorly understood, as well as the molecular mechanism underlying tumor-macrophage interactions.


Conclusion
Together, our ndings deciphered the macrophage heterogeneity and detected a tumor environmentalspeci c gene expression program in macrophages of PDAC, which uncover a potential mechanism laying the tumor-promoting role of TAM. The ligand-receptor interactions identi ed would be potential biomarkers for anti-cancer targeted therapy.

Background
Pancreatic cancer, mostly pancreatic ductal adenocarcinoma (PDAC), remains one of the most lethal cancers with a 5-year survival rate of less than 8% [1]. The intractable characteristics of PDAC are the result of comprehensive factors like demographic factors, genetic factors, medical background conditions and environmental factors, which relates to the accumulation of acquired drive mutations (K-Ras, HER2/neu, TP53), and as well as abnormalities in growth factors and their receptors, resulting in constitutive activation of downstream signaling pathway and the transcription factors [2][3][4]. As regards to the limited survival bene ts of gemcitabine, the rst-line regimen for pancreatic cancer, a better understanding of the mechanisms that contribute to uncontrolled proliferation is urgently needed.
Intriguingly, PDACs have a complex morphology environment where cancer cells are surrounded by a cellular milieu containing diverse cell types such as stellate cells, broblasts, endothelial cells and myeloid cells [5]. Myeloid compartment, including mostly tumor-associated macrophages (TAMs) and dendritic cells (DCs) as revealed in single-cell analyses of multiple cancers, contribute to all aspects of tumor progression [6][7][8]. PDACs could elicit a particularly severe T-cell exhaustion signature and immunosuppressive niche with TAMs [9], which are a heterogeneous cell type that can contribute to gemcitabine resistance by reduced apoptosis of cancer cells and promote malignancy by angiogenesis, immunosuppression and extracellular matrix (ECM) remodeling [10][11][12]. As the mechanism implicated in the communication between TAMs and pancreatic tumor cells during tumorigenesis is complicated, there is need for more knowledge about the spectrum of cell-cell interactions occurring in TME of PDACs and understanding how these interactions in uence the outcome.
Single-cell RNA sequencing (scRNA-seq) enables the characterization of cell diversity and heterogeneous phenotypic states of macrophage compartments in more detail, which has been applied recently for primary pancreatic cancers and mouse PDAC models [13][14][15]. However, the functional heterogeneity of different subpopulations, cellular reprogramming and phenotypic trajectories of macrophages during PDAC progression remain poorly characterized. Developmental trajectory construction in colon cancer has suggested the origin of TAMs from tissue-resident macrophages and recruited monocytes that subsequently polarized into macrophages [7]. Researches have highlighted that tissue-resident macrophages could be originated from different resources with diverse phenotypes or functions [16,17]. Tissue environment itself performs environmental programming of tissue-speci c macrophage, and can regulate the gene expression program regardless of origin [18]. In PDACs, tissue-resident macrophages could expand during tumor progression and contribute to growth of PDACs [19]. While the understandings of its reprogramming trajectory and molecular phenotype are still limited.
Although macrophages divisions based on their polarization states could be into three types: un-activated macrophages (M0 macrophages), the classically activated type 1 (M1 macrophages), and the alternatively activated type 2 [20]. Recent research on the genetically engineered mouse model of PDACs, as well as human primary gastric cancer revealed no such classi cation of clusters in the single-cell level [14,21]. The in vitro programming of M0 macrophage to M1 is thought to enhance the cytotoxicity against tumor cells by pro-in ammatory factors. However, by comparing the early stage and late stage of PDAC mouse model, it was found that PDAC progression is accompanied by an increase of in ammatory features in macrophages [14,22]. It raised questions about M1/M2 classi cation and suggested pro-and anti-in ammatory features might not be dependently exited in the in-vivo system. As novel dichotomous functional phenotypes of TAMs have been revealed across different cancer types with a difference of susceptibility to targeted therapy [7,15,23], extended functional classi cation of macrophage was required to be determined in PDACs.
Here, we utilized scRNA-seq and TCGA RNA-seq datasets to perform computational analysis, through which we constructed a novel cell-cell interaction network to de ne key cell communication involved in regulating malignancy of PDACs. Besides, we found that macrophages in PDACs is highly heterogeneous and featured with different gene expression pro le and regulatory networks (regulons). We further identi ed that the tissue-resident macrophage and in ammatory monocyte are the potential sources of TAMs with changes of a list of novel gene expression during the differentiation trajectories. Our study uncovered cellular reprogramming programs and cell-cell communication features that speci c to TAMs and thus might offer insights on inferring therapeutic targets.

Generation of the Primary PDAC Datasets
The Cancer Genome Atlas (TCGA) Pancreatic Cancer project was composed of gene expression data and relevant clinical information (n=182) which were downloaded through the UCSC data portal ( https://xenabrowser.net ). The raw count data obtained from the RNA-sequencing was normalized concerning library size via the R Deseq2 package, then it was performed log transformation and changed into same comparison level. After excluding samples that were not a histological type of primary ductal cancer, we reserved a total of 147 cases of primary samples for further survival analysis.

Single-cell Data Processing
The published PDAC single-cell dataset [13] containing 24 PDACs and 11 normal pancreas tissues was deposited in the Genome Sequence Archive with the accession number of PRJCA001063. We retrieved the raw data matrix and transformed it into a SingleCellExperiment object using the Seurat R package. After removing cells with poor quality (<200 genes/cell, <10 cells/ gene and >10% mitochondrial genes), raw matrix was processed using "Sctransform" methods, which contains an encapsulated function of "NomalizeData", "ScaleData", and "FindVariableFeatures" for automatically high-variance gene identi cation. Then we generated the ElbowPlot and determined the optimal dimensionality -number of principal components (PCs). As high resolution often leads to an increasing number of clusters, the "FindClusters" function for applying modularity optimization was adjusted to 0.15. The UMAP (Uniform Manifold Approximation and Projection), a dimension reduction technique was used to visualize the subclusters. To assign cell identity to each cluster, we consulted previously published literature and the CellMarker dataset (http://biocc.hrbmu.edu.cn/CellMarker/) for reference. Following characterizing the identity of myeloid cells/macrophages based on the markers AIF1 and CD68, the same work ow was used to identify sub-clusters of macrophage. To perform over-representation (ORA) analysis, we identi ed differentially expressed features (cluster biomarkers) for each cluster against other cells using the "FindMarkers" function in which parameters of min.pct (minimum percent of cells expressing the gene) and log2fold-change are set to 0.1 and 0.5, respectively. Pathway enrichment analysis was based on the items of the Reactome pathway using ClusterPro ler R package.
Correlative cell-cell interactions inferred by combined scRNA-seq and TCGA-PDAC datasets Based on the TCGA RNA-seq pro les, we estimated the relative abundance of each cell subtype by the average expression of the cell type-speci c genes de ned above. Next, we conducted a Pearson correlation analysis between the expression of each gene and the relative abundance of each cell type to identify genes that may infer the co-occurring cells with these particular cell types. Since the expression levels of speci c signature/marker genes in a given cell subtype were obviously in high correlation with the abundance of this speci c cell subtype, we identi ed these genes based on the single-cell dataset according to the criteria of average expression > 1 and cell frequency of expression > 20%. We then ltered these genes and transformed their correlation coe cient to 0. After that, the top 20 highly correlated non-self-expressed genes were selected from the ranked correlation matrix by the correlation coe cient. To identify which cell subtypes could contribute to the highly correlated non-self-expressed genes, we performed a gene set enrichment analysis by calculating the mean value of z-score transformed expression of these genes for a speci c cell subtype. At last, the correlated cell types were identi ed as the z-score transformed enrichment score > 1.28. If the correlated two types were identi ed mutually, the maximum of the enrichment scores was chosen.

Single-cell regulatory network inference and clustering (SCENIC)
For single-cell regulons (i.e., transcription factors and their target genes) inference in the different subclusters of macrophage, the SCENIC (https://aertslab.org/) vignette was used for reference [24]. The regulon inference was performed in the following three steps. (1). co-expression modules between transcription factors and the potential target were identi ed based on the gene-expression matrix through GENIE3 (R package) (2). Next, to remove indirect targets, these modules were pruned by cis-regulatory motif discovery (cisTarget), leaving 94 regulons to be analyzed in the next step. (3) the activity score of each regulon at cellular resolution was computed using the AUcell algorithm. Finally, the cellular activity patterns of these regulons were used for nonlinear projection and clustering. The overexpressed activity scores of regulons for each sub-cluster were visualized using a heatmap.

Cell-cell communication inference
To systematically study the interactions between macrophage and tumor cells, we predicted cell communication via CellphoneDB (www.cellphonedb.org), a publicly available repository of ligandreceptor interacting pairs and multi-subunit protein complex [25]. It considered the expression level of ligands and receptors with each cell type, generate a null distribution of each ligand-receptor pair, and determined the likelihood of cell-speci city of a given receptor-ligand complex based on the P-value. Those with the least number of signi cant P-values are chosen to be biologically relevant. After ltering in a criterion of P < 0.05 and secreted partners, a total of 94 ligand-receptor pairs were identi ed as shown in Supplementary Table 3.

Trajectory analysis
The monocle2 package in R was used to determine the developmental trajectory of macrophages. Following the monocle protocol, we passed the normalized data of indicated sub-clusters into monocle2 to establish a monocle object and selected top 2,000 most signi cantly (q <0.1) differentially expressed genes (DEGs) as the set of ordering genes sorted by the q value. For isolating the root-to-branch speci c gene expression patterns, the monocle BEAM statistical test was utilized. For visualizing the fatedependant gene expression patterns, "plot_genes_branched_heatmap" function was used. For visualization of change tendency for a single gene, "plot_genes_branched_pseudotime" function was adopted.

Statistical analysis
In survival analysis, the means of the average expression level of interacting molecule 1 in tumor cell and interacting molecule 2 in macrophage are calculated using the square root: mean = (molecule1 * molecule 2) 1/2 . The Kaplan-Meier survival curve was built to show the differential survival among the two groups, and the log-rank test was performed to use statistical tests for comparing the survival distributions of two groups (P<0.05). All statistical analyses were performed using R software (3.6.2).

Myeloid and tumor cells comprise the core of a predicted cell-cell interaction network in PDAC
To interrogate global cell-cell interactions in PDACs, we performed the computational analysis by combining scRNA-seq and TCGA-PDAC RNA-seq datasets. The scRNA-seq pro les were obtained from published datasets, which were processed to identi ed different cell clusters and cell type-speci c signatures as previously de ned [13]. For each cluster, we identi ed gene signatures that were expressed in greater than 25% of cells as well with a log2 fold-change greater than 0.5 (Supplementary Table 1). In addition to identifying the abundance of a particular cluster in each TCGA sample using its gene signature, we generated non-self highly correlated genes and matched this gene set back to the scRNAseq dataset, through which co-occurring enriched cell type was identi ed to build the correlative networks (Supplementary Figure 1a). With the analyses, interactions were mainly identi ed between malignant ductal cells, broblasts, endothelial cells and myeloid cells, re ecting the crosstalk of them as previous ndings de ned [10,26,27]. Among the interactions, Myeloid cells and tumor cells harbor the most number of partners. Besides, they act as the core of the interaction network accompanying the highest enrichment score, indicative of the strong connection between them (Supplementary Figure 1a). In contrast, we observed no such interaction in normal pancreas tissue, suggesting it is a biologically meaningful interaction. Then we examined the relative population contribution in tumor and normal

Identi cation ofheterogeneity for myeloid cells
The myeloid cells compartment is comprised of a heterogeneous population of cells [28,29]. By combining myeloid cells from tumor and normal tissue, the 5408 cells were further clustered into 8 subsets (Fig. 1a). Of these, cluster 5, 6 and 7 belong to myeloid DC lineage as shown by low expression of CD68, which was not identi ed in a previous analysis [13]. Cluster 5 corresponds to plasmacytoid dendritic cell (pDC) as shown by DERL3, CXCR3 and IGJ, cluster 6 represents type-2 conventional DC (cDC2) as shown by CCR7, IDO1 and DAPP1, and cluster 7 was identi ed as cDC1 cells with high expressions of CD1C, CD1E (Fig. 1b, Supplementary Figure 2a). The leakage of these DC clusters into macrophages might be owing to the low resolution (granularity 0.15) of the above main-class clustering.
The remaining 5 clusters expressing CD68 and CD163 were designated as macrophages, albeit they showed rather a poor separation of clusters (Fig. 1b, Supplementary Figure 2a), suggesting that they represent diverse cell states on a graded scale rather than separate entities. It was noticed that 91.9% of the cells in cluster 2 came from normal tissue, indicating cluster 2 might represent an entitle of tissueresident macrophage (Fig. 1a, Supplementary Figure 2b). However, we surprisingly found 47.3% of the macrophage population in normal tissue was from cluster 0 and cluster 1 (Supplementary Figure 2b), which suggested the involvement of monocyte-derived macrophages in normal tissues. To further understand the molecular features of the ve sub-clusters of macrophages mentioned above, we conducted a DEG analysis and identi ed cell makers for each cluster (Supplementary Table 2). Cluster 4 expresses S100A8, S100A9, VCAN, and FCN1 which are associated with in ammatory monocytes [29,30],  Figure 2c). Intriguingly, different tumor samples consist of a group of macrophages with high heterogeneity, especially the ratio between cluster 1 and cluster 0, two large populations of macrophages among patients (Fig. 1c). It interested us whether the two populations were representative of the different polarized states of macrophage. As there was no recognizable cluster representative of M1 and M2 polarization macrophage, we examined the expression of marker genes representative of M1 (e.g. CXCL9, IL1β, CCL5) and M2 (e.g. MRC1, CCL18, CCL23, CD163) states across the macrophage clusters. However, the expression of M1/M2 markers shows no evident differences among the sub-clusters (Fig. 1d).

t-SNE atlas reveals gene set enrichment state across different sub-clusters of macrophages
To delineate the functional heterogeneity of sub-clusters of macrophages, and to see which cell cluster principally contribute to the malignant progression of tumor cells, we performed gene set enrichment analysis using MsigDB hallmark gene terms and some immune-related terms, which aggregate redundant biological terms and represent well-de ned biological states or process. Gene-set variation analysis (GSVA) revealed that cluster 0 displayed upregulation of pathways associated with complement cascade, antigen processing and presentation, PD1 signaling. While cluster 1 was characterized by high expression of cancer-related terms: angiogenesis, in ammatory response, glycolysis, cholesterol homeostasis and epithelial-mesenchymal transition, indicating it is more aligned with the property of TAMs. In correlation with its high proliferating property, cluster 3 expressed the highest level of apoptosis, G2M checkpoint and mitotic spindle (Fig. 2a-c). In addition to GSVA analysis, GRO analysis based on the Reactome pathway database further revealed that pathways related to MHC-II antigen presentation, CD28 family and TCR signaling were predominantly enriched in cluster 0. In contrast, cluster 1 was more activated and enriched with chemokines and extracellular matrix pathways, while cluster 3 was featured with signaling pathways related to cell cycle and mitotic process. Besides, cluster 4, monocyte-like macrophage, was associated with a gene expression pattern towards in ammation and innate immune response (Supplementary Figure 4). Thus, these data suggested the functional heterogeneity of sub-clusters of macrophages and indicated cluster 1 likely represented TAMs.

Identi cation of master regulatory networks (regulons) across different sub-clusters of macrophages
The transcriptional state of each cluster emerges from an underlying gene regulatory networks. We thus applied SCENIC to assess which regulatory networks determine the cell state for each sub-cluster of macrophage. By projecting the cell identities on the t-SNE clustering atlas based on the regulons' activities, we found each cluster of macrophages located separately from other clusters, showing differences of transcriptional states (Fig. 3b). A distinct set of regulatory genes, including EZH2, SOX9, HES4, ATF4, IRF1, IRF2 et al., were revealed to de ning the various macrophage populations (Fig. 3a). Consistent with their in ammatory phenotype and functions, SCENIC revealed that genes regulated by in ammatory transcription factors, IRF1, CEBPB, and STAT3, were upregulated in a subset of monocytelike cells. We also identi ed EZH2, best known for its function in macrophage activation [32,33], as a candidate transcription factor driving the proliferating status of cluster 3 macrophages. Besides, we identi ed STAT1 speci cally in cluster 1 macrophages, and upregulation of SOX9 and HES4 seemed responsible for the transcriptional phenotype of tissue-resident macrophages (Fig. 3c).

Pseudotime trajectory reconstruction of macrophage differentiation course
To delineate monocyte-to-macrophage differentiation, the macrophage subpopulation of PDAC tissue was arranged to build single-cell pseudotime ordering. According to a developmental course, the trajectory constitutes one decision point and two termini corresponding to two distinct cell fates ( Fig. 4a  and 4b). By projecting the cell identity on the trajectory plot, we found cluster 4 was mainly distributed in the root state, verifying its monocyte-like phenotype (Fig. 4c). Besides, cluster 2 was mainly in the location of the root, which might suggest the differentiation of a subset of tissue-resident macrophages into TAMs. Besides, the cluster 1 randomly appeared in the whole developing course of the pseudotime path, while cluster 3 and cluster 2 cells populated mainly on the two termini of the course constituting the pool of mature differentiated macrophage (Fig. 4c), suggesting they are in a polarized position. Next, we assessed the expression pattern of genes regulated during differentiation in cells at fates 1 and 2 of the trajectory. Pro ling of gene expression dynamics along the trajectory revealed a reduction of known monocyte markers (S100A9, S100A8) and upregulation of MMP9, CCL2, CTSL et al. (Fig. 4d). Pseudotime kinetics of SPP1 expression showed gradually upregulated from the root to both fates, while STMN1 was increasingly elevated from the root to fate 2. Besides, CTSS expression was upregulated in the late stage of fate 2 differentiation and downregulated early in early fate 1 differentiation. On the contrary, S100A9 expression slightly decreased at cells from the root to fate 1, but markedly increased in the late phase of pseudotime branch via fate 2 (Fig. 4e).

Inference of tumor-macrophage interactions via CellphoneDB
Having  Figure 5c). By contrast, we found cluster 3 express high levels of chemokine receptor CX3CR1, the receptor for CX3CL1, which might account for the in ltration and recruitment of macrophage with the stimulation of tumor cells [34].
We subsequently sought to assess the degree of interactions for each cluster-cluster pair in the TME. As a result, the presented correlation matrix revealed active macrophages-tumor cell interactions (Supplementary Figure 5a). After ltering with a stringent standard, the master ligand-receptor pairs embedded in the communication between tumor cells and macrophages were illustrated in Figure 5A-B. Previous work has de ned the secretion of CSF1 by tumor cells and its interaction with the CSF1R on the macrophages [35]. Besides, Osteopontin (SPP1) -CD44 signaling, which presented in the list with remarkable signi cance, has been shown to lay an important role in cell migration, tumor growth and progression [36]. This result also raised the previous unknown pairs, such as CCL3L1-DPP4, SPP1-PTG, LGALS9-MET, TNF-DAG1 et al, which might serve critical functions. (Fig. 5a).
Having de ned the master ligand-receptor pairs, we investigate the value of these ligand-receptor pairs in predicting patients' survival. We choose the 21 pairs in which ligand was representative of macrophagesorigin to observe the tumor-promoting effect of each ligand. Through analyzing RNA sequencing datasets of PDCA from the TCGA portal, we found 5 pairs (HBEGF-CD44, HBEGF-EGFR, LGALS9-CD44, LGALS9-MET, GRN-EGFR) were correlated with poor outcome of patients (Fig. 5b, Supplementary Figure 5d). Among them, MET and EGFR were found highly expressed in the malignant ductal cells compared with normal ductal cells by projecting the expression levels of these molecules onto the UMAP clustering.
Correspondingly, LGALS9, HBEGF and GRN were con rmed overexpressed in macrophage cluster (Supplementary Figure 5e). We conclude that over-expressions of LGALS9, HBEGF and GRN in macrophages may help to promote PDAC progression, thus predict adverse clinical outcomes in patients.

Discussion
The 'Bulk' transcriptome sequencing is a commonly used modality to see gene expression changes of speci c phenotype and pathway activity. However, it does not resolve problems of observing gene expression changes of speci c cells. Consequently, the discrepancy of gene expression pro ling between neoplastic cells and stromal cells would obscure the true activity identity for speci c molecular signaling. Herein, we employed scRNA-seq methods to provide more insights into phenotypes and gene signatures of macrophage in PDAC.
Our data revealed 8 molecular subsets of myeloid cells that consist of previously ignoring DCs subsets, which perhaps was owing to the different clustering setting for resolution and high expression level of AIF in DCs [13]. The other 5 macrophage subsets identi ed demonstrated distinct phenotypes and functionalities. Among them, cluster 2 was de ned as tissue-resident macrophages, which was not identi ed before [13]. It is generally accepted tissue-resident macrophage plays crucial roles in pancreatic tissue restoration and regeneration [37]. The recent nding indicated that the tissue-resident macrophages become detrimental by promoting pancreatic cancer progression. [19] During tumorigenesis, we found tissue-resident macrophages present signi cant phenotypic change and cellular reprogramming, as shown by the remarkable reduction of them in the tumor comparing with normal tissue ( Figure A2A). Our data revealed one molecular subtype of macrophages conferred with a monocyte-like phenotype for the following reasons. Firstly, this cluster highly expressed CCR2, an in ammatory monocyte marker [30]. Besides, through pseudotime trajectory analysis, we found this cluster was present in the starting point of differentiation trajectory. Moreover, cell-cell interaction inference indicated this cluster has the weakest relationship with tumor cells, supporting it is in a condition of immaturity. Besides, the other 3 clusters: cluster 1 expresses numerous chemokine-and in ammation-associated genes, cluster 3 expresses proliferation markers and enriched in the cell cycle pathway, and cluster 0 is rich in MHC-II-associated genes. We dissected the cell communication frequencies for each cluster and found that cluster 3 and 2 were in more close contact with tumor cells. Also, these two clusters showed enrichment of cytokines-related pathways or in rapid proliferation stage, which is consistent with the characteristics of TAMs [38]. According to the pseudotime cell trajectory, we identi ed both in ammatory monocytes and tissue-resident macrophages as the potential sources of TAMs by trajectory analysis. TAMs compose essential compartment of the cancer-immune microenvironment and play critical roles in the regulation of tumor progression [39,40]. An in-depth understanding of the ontogeny of TAMs could facilitate therapeutic intervention.
As regards to the polarization states of macrophages, Wang et al. revealed that pancreatic cancer cells activate macrophages to the M2 phenotype in a HIF1/2a-dependent manner, which then promotes the migration, invasion, and EMT of pancreatic cancer cells [41]. While other ndings indicated as PDA progresses, an in ammatory feature with macrophage is substantially increased [14]. In our investigation, the ve sub-cluster showed no evident difference between M1 and M2 markers. Consistent with previous studies, neither of these populations are representative of M1 or M2 cells [8,42]. According to our hypothesis, the M1 is more likely a concept indicating that macrophage performs a pro-in ammatory function in particular settings rather than a speci c cell cluster. It was observed that M1-TAMs secrete TNF that induces activation of the NF-κB signaling pathway to drive acinar-to-ductal metaplasia in the beginning stages of PDAC [43]. While as tumor invasion, the M1 macrophages secrete MMPs, like MMP-9, to degrade the extracellular matrix [44]. In other settings such as proliferation and vascularization, the M2-like phenotype seizes advantage [45]. We hypothesis that the polarization state of each cell frequently changes over the course of tumorigenesis as well as the number of cells presenting a speci c polarization state, for which more investigations are required.
Beyond characterizing the cellular composition of PDACs, it is crucial to understand how the different cellular components interact with one another to give rise to the biology of tumorigenesis. Numerous studies have begun to utilize scRNA-seq data to characterize cell-cell interaction [46,47]. However, limited to the small sample size, inter-tumoral heterogeneity tends to skew the results. Here, we have leveraged the advantages of scRNA-seq and large bulk RNA-seq samples to deduce a cell-cell interaction network, through which we identi ed myeloid lineage and tumor cells as key mediators of cellular cross-talk in the tumor microenvironment. Investigation of cell-cell interactions via CellphoneDB database suggested a cell-cell interaction pro le which may be activated and work in mutual supporting processes between macrophage and tumor cells, therefore contributing to GBM progression. Depend on the CSF1-CSF1R interaction, the crosstalk between tumor and macrophage shifts the balance towards an immunosuppressive environment. Previous work by others demonstrated that CSF1R blockade effectively decreases the number of TAMs, reduced immune suppression and elevated immune response [35]. SPP1, also known as Osteopontin, was identi ed as secreted factors in our interaction pairs, showing an association with macrophage activation and polarization [48,49]. It has been substantiated that TAM secreted SPP1 to promote the survival of glioblastoma cell and angiogenesis [50].
Second, we identi ed a correlation between the ligand-receptor pair and poor prognosis of PDAC patients. Previous work highlights the critical role of galectin-9 in regulating immune response. Galectin-9 enhances iTregs function and stability through CD44, which strengthens the iTreg cell phenotype [51]. However, lack of evidence focused on the function of paracrine galectin-9 on tumor cells. Our data demonstrate that galectin-9-CD44/MET generates an autocrine loop underlying the interactions with a unique value in predicting the prognosis of patients. Besides, Galectin-9 has been reported to been highly expressed in human PDAC comparing with the normal pancreas and is capable of polarizing macrophages toward a pro-tumorigenesis M2 phenotype. It is also available to be detected in the blood serum, thus showed great value in diagnosing PDAC [52,53]. A study on pancreatic cancer identi ed expression and cleavage of HB-EGF were required for proteasome inhibitor-induced EGFR activation. Ray et al. demonstrated the macrophage expression of HB-EGF links chronic pancreatitis and K-Ras initiated PDAC [54]. We veri ed macrophages are the primary source of HB-EGF production and rstly illustrated that HB-EGF and EGFR signaling in pancreatic tumor-stroma interactions may explain the limited e cacy of therapy and poor prognosis in pancreatic cancer patients. Another ligand-receptor pair for concern in prognostic determination was granulin (GRN)-EGFR. Recent work demonstrated that macrophage-derived granulin facilitates pancreatic cancer metastasis and resistance to αPD-1 therapy [55,56]. Our data may provide the potential mechanism behind it.

Conclusion
In brief, our study de nes a phenotypically diverse subpopulation of macrophage in TME of PADC through computational analysis. Our results improve understandings of the macrophage's related impact on the tumor, and interactome analysis demonstrated the underlying mechanism of intercellular communication between macrophage and cancer cells. Regarding limited immunotherapy options in hand, tumor-macrophage interactome represents potential treatment targets to inhibit cancer proliferation, overcome the immunosuppressive microenvironment and aid in the implementation of immunotherapy approaches in PDAC. Authors' contributions Y.K.D. was responsible for the study designing and study supervision; G.J. was responsible for raw data processing and data clearing. Y.K.D. and Z.X were responsible for the development of methodology, analysis, and interpretation of data. G.J. and Z.X. reviewed, revised and approved of the manuscript. All authors read and approved the nal manuscript.

Funding
This work was supported by grants from the National Natural Science Foundation of China (Nos. 81602606 to G.J.).
Availability of data and materials