Skip to main content

Single-cell transcriptomics reveal distinct immune-infiltrating phenotypes and macrophage–tumor interaction axes among different lineages of pituitary neuroendocrine tumors



Pituitary neuroendocrine tumors (PitNETs) are common gland neoplasms demonstrating distinctive transcription factors. Although the role of immune cells in PitNETs has been widely recognized, the precise immunological environment and its control over tumor cells are poorly understood.


The heterogeneity, spatial distribution, and clinical significance of macrophages in PitNETs were analyzed using single-cell RNA sequencing (scRNA-seq), bulk RNA-seq, spatial transcriptomics, immunohistochemistry, and multiplexed quantitative immunofluorescence (QIF). Cell viability, cell apoptosis assays, and in vivo subcutaneous xenograft experiments have confirmed that INHBA-ACVR1B influences the process of tumor cell apoptosis.


The present study evaluated scRNA-seq data from 23 PitNET samples categorized into 3 primary lineages. The objective was to explore the diversity of tumors and the composition of immune cells across these lineages. Analyzed data from scRNA-seq and 365 bulk RNA sequencing samples conducted in-house revealed the presence of three unique subtypes of tumor immune microenvironment (TIME) in PitNETs. These subtypes were characterized by varying levels of immune infiltration, ranging from low to intermediate to high. In addition, the NR5A1 lineage is primarily associated with the subtype characterized by limited infiltration of immune cells. Tumor-associated macrophages (TAMs) expressing CX3CR1+, C1Q+, and GPNMB+ showed enhanced contact with tumor cells expressing NR5A1 + , TBX19+, and POU1F1+, respectively. This emphasizes the distinct interaction axes between TAMs and tumor cells based on their lineage. Moreover, the connection between CX3CR1+ macrophages and tumor cells via INHBA-ACVR1B regulates tumor cell apoptosis.


In summary, the different subtypes of TIME and the interaction between TAM and tumor cells offer valuable insights into the control of TIME that affects the development of PitNET. These findings can be utilized as prospective targets for therapeutic interventions.


Pituitary neuroendocrine tumors (PitNETs) are a frequently occurring form of intracranial brain tumor. The estimated prevalence rate of clinically significant PitNETs is approximately 1 case per 1000 individuals [1]. About 2% of cases display aggressive behavior, while up to 15% show invasion, with a few cases of craniospinal or distant metastases occurring in 0.1–0.2% of cases. PitNETs can be either functional (i.e., producing hormones and causing hormonal imbalances) or non-functional, leading to neurological symptoms based on size and location. Although significant therapeutic progress in strategies, such as medication and surgery, has been made, many PitNET patients have poor clinical outcomes. The 2022 World Health Organization classification reported the role of key transcription factors in the differentiation of distinct types of pituitary cells, indicating the complexity and heterogeneity of PitNETs [2]. Revolutionary advances in single-cell RNA sequencing (scRNA-seq) and bulk-RNA sequencing technology allow a more comprehensive understanding of the origin and differentiation status of PitNET cells [3,4,5,6,7,8,9]. In addition to deciphering the tumor cells, efforts have been made to characterize the tumor immune microenvironment (TIME), improving the clinical benefit of immunotherapy, especially anti-PD-1/L1 and anti-CTLA4 antibodies [10, 11]; however, many patients still have poor outcomes [12, 13]. Therefore, a deeper insight into PitNET TIME and identifying potential targets is crucial.

TIME is a well-recognized factor in the development and progression of PitNETs [14,15,16]. Macrophages and T lymphocytes are the most prevalent infiltrating immune cells in PitNETs [17, 18]. Also, other immune cell subtypes such as B lymphocytes, FOXP3+ cells, neutrophils, NK cells, or mast cells may be present less frequently [15]. However, the specific function and heterogeneity of these immune cell subtypes in PitNETs remain poorly understood, which can help identify new therapeutic targets. Some studies have endeavored to characterize the immune microenvironment of PitNETs using immunohistochemistry (IHC) [16], flow cytometry (FC) [19], and RNA-seq (including both bulk RNA-seq and single-cell RNA-seq) [3, 5,6,7, 9]. Macrophages, a highly heterogeneous group of immune cells, have many functions, including phagocytosis, antigen presentation, cytokine production, and tissue remodeling [20]. Growing evidence has suggested that immune cells, particularly macrophages, play an essential role in the development and progression of pituitary adenomas [21, 22]. Macrophages are present in both normal pituitary and PitNETs, and their types and populations vary depending on the subtype of PitNETs [23,24,25]. M2 macrophages are prevalent in the pituitary gland and may play a role in tumorigenesis [26]. The number of macrophages correlates with the size and invasiveness of PitNETs; higher macrophage content is seen in sparsely granulated somatotrophinomas and null-cell PitNETs compared to densely granulated somatotrophinomas or corticotrophinomas [22]. M2 macrophages are also more prevalent in non-functioning PitNETs (NF-PitNETs) with cavernous sinus invasion than in non-invasive NF-PitNETs. A higher M2/M1 macrophage ratio was associated with larger and more proliferative NF-PitNETs [25,26,27]. Moreover, the M2 macrophage-conditioned medium increased the proliferation, invasion, and migration of primary NF-PitNET cells [19, 27, 28]. However, the heterogeneity of macrophages in various PitNETs and their subpopulation classification and communication with tumor cells remain unexplored.

Here, we performed integrated analysis on scRNA-seq and bulk-RNA seq data of 3 tumor lineages from 23 adult humans to establish the transcriptional landscape of immune cells in PitNETs. The pituitary tumors were classified into three immune infiltration states based on immune infiltration, and five subtypes of related macrophages were identified. IHC by TMA (tissue microarray) staining and FC confirmed that POU1F1 lineages had more immune cells, while NR5A1 lineages had lower immune cell infiltration. CX3CR1+, LYVE1+, C1Q+, and GPNMB+ macrophages in PitNETs were identified by multiplex IHC (mIHC). Furthermore, INHBA secreted by CX3CR1+ macrophages partially induced cell apoptosis. In contrast, blocking the INHBA receptor ACRV1B reversed this effect. Overall, our findings provide a comprehensive atlas of PitNETs TIME subtypes, particularly from the aspect of macrophages, and highlight the role of CX3CR1+ macrophages in modulating transcriptional status and intercellular communication patterns.


Clinical samples

Patient recruitments three cohorts of patients were recruited: 23 scRNA-seq, 365 bulk RNA-seq PitNETs, 128 Immunohistochemistry (IHC) of TMA, and 45 flow cytometry for analysis and validation. All patients received surgery at the Department of Neurosurgery at Ruijin Hospital, an affiliate of Shanghai Jiao Tong University School of Medicine. The scRNA-seq and flow cytometry PitNETs cohort underwent surgery between 2020 and 2023, while the bulk RNA-seq and IHC PitNETs cohorts were collected between 2017 and 2022. Both cohorts only included subjects without previous malignancies. The clinical data of all patients were retrospectively obtained from medical records (Additional file 2: Table S1-5). The patients are assigned numbers matching the GroupsIndex column in Supplementary Table S1-5. Written informed consent was obtained from all patients, and the Ruijin Hospital ethical committee approved the study.

Tissue dissociation and primary tumor cell extraction

Tissues were transported in DMEM (Gibco, Cat. no. 11875–093) with 1 mM protease inhibitor (Solarbio, Cat. no. P6730) on ice to preserve viability, washed two to three times with phosphate-buffered saline (PBS; Hyclone, Cat. no. SH30256.01), then minced on ice. We used dissociation enzyme cocktail 1 mg/ml Type VIII Collagenase (Sigma-Aldrich, Cat. no. C2139), 2 mg/ml Dispase II (Sigma-Aldrich, Cat. no.4942078001), 1 mg/ml Trypsin Inhibitor (Sigma-Aldrich, Cat. no. T6522), and 1 unit/ml DNase I (NEB, Cat. no. M0303S) dissolved in serum-free DMEM to digest the tissues. Neoplastic tissues were dissociated at 37 °C with a shaking speed of 50 r.p.m for about 40 min. We repeatedly collected the dissociated cells at intervals of 20 min to increase cells yield and viability. Cell suspensions were filtered using a 40-μm nylon cell strainer (Falcon, Cat. no. 352340) and red blood cells (RBC) were removed by RBC lysis buffer (Invitrogen, Cat. no. 1966634) with 1 unit/ml DNase I. Dissociated cells were washed with PBS containing 0.04% bovine serum albumin (BSA; Sigma-Aldrich, Cat. no. B2064) with step-by-step descending centrifuging speed and increasing time. Cells were stained with 0.4% Trypan blue (Invitrogen, Cat. no. T10282) to check the viability and then cultured in DMEM containing 10% FBS and 1% antibiotic mixture.

RNA sequencing alignment

To quantify gene expression in the transcriptome, the RNA sequencing raw FASTQ files were aligned to the human reference genome GRCh38 (release 40). The human reference genome and its annotation file were acquired from the GENCODE database, accessible at Salmon (v1.8.1) [29] was used to generate the count and transcripts per kilobase of exon model per million mapped reads (TPM) matrix. The transcript counts were combined using DESeq2 [30] and then transformed into fragments per kilobase million (FPKM) to normalize the gene length, using the TPM matrix to evaluate the gene expression level. To identify differentially expressed genes, we used the limma [31] package with a significance level set at an adjusted P-value < 0.05 and a log2-transformed fold change (log2FC) > 0.58 (FC > 1.5). To identify differentially expressed genes between more than two conditions, we used analysis of variance (ANOVA) to calculate the P-values.

Functional enrichment analysis

Gene Set Enrichment Analysis (GSEA) [32] was employed with the pre-ranked algorithm and the R package clusterProfiler [33] to perform functional enrichment analysis. Genes were ranked based on their log2FC. The gene sets from the Molecular Signatures Database (MSigDB, v7.5.1) [34] of the Broad Institute were used. Specifically, we used the HALLMARK gene sets (H) [34], which represent 50 well-defined biological processes, and the Kyoto Encyclopedia of Genes and Genomes (KEGG) [35, 36] database to identify differences induced by various drugs. A significance cutoff of P-value < 0.05 was considered. To visualize the enrichment results, we utilized the enrichplot R package ( Another functional enrichment analysis algorithm, Enrichr ( [37] web service, was conducted to evaluate gene signatures. P-value < 0.05 was set to the significance level.

Single-cell RNA sequencing alignment and generation of gene expression matrix

Raw sequencing reads of the cDNA library was processed through the BD Rhapsody Whole Transcriptome Assay Analysis Pipeline (v1.8), which included filtering by reads quality, annotating reads, annotating molecules, determining putative cells, and generating a single-cell expression matrix. The pipeline also selected the sample origin of every single cell via the sample determination algorithm according to the sequencing reads of the SampleTag library. Among all the output files, a matrix of UMI counts for each gene per cell was used for downstream analysis. Genome Reference Consortium Human Build 38 (GRCh38) was used as a reference for the BD pipeline.

Analysis of single-cell expression data and cell type annotation

Gene expression matrices were imported for downstream analysis using the R Seurat package [38].

Single-cell quality control (Additional file 2: Table S6)

In the initial step, we conducted a filtration process to eliminate cells exhibiting low gene expression. The gene expression matrices of individual cells were imported for subsequent analysis using the R Seurat package. Cells with expressed genes < 200, or unique counts > 50,000 or < 500, or expressed mitochondrial RNA > 30% were removed for quality control, resulting in 98,523 cells being retained. Upon analyzing the cell distributions across the 23 samples, the median number of captured cells was 2837. Among the 23 samples, 6 had a cell count exceeding 5000, with 2 samples notably containing a high number of captured cells (> 10,000 cells). To ensure effective integration, we employed random downsampling to limit the maximum number of enrolled cells for each sample to 5000 using the R sample function. Ultimately, a total of 69,539 cells from PitNETs were included for downstream data analysis (Additional File 2: Table S2). We performed the doublets analysis by using the R package DropletUtils (; we identified 108 cells (0.15%) as potential doublets among the total 69,539 cells.

Batch effect adjustment

Following the generation of the Seurat object containing 69,539 cells, we identified the top 3000 highly variable genes and utilized them to conduct principal component analysis (PCA). Batch effect correction across samples was executed using the R package Harmony [39], with the parameter “max.iter.harmony” set to 5. Subsequently, the top 20 harmony coordinates were chosen for clustering analysis and dimensionality reduction.

Clustering and dimensionality reduction

The top 20 harmony coordinates were selected for graph-based unsupervised cell clusters with the resolution was set to 0.8. Dimensionality reduction methods t-Distributed Stochastic Neighbor Embedding (t-SNE) [40] and Uniform Manifold Approximation and Projection (UMAP) [41] was performed for visualization. The t-SNE was selected for final interpretation.

Annotation of major cell populations

After clustering, 30 clusters were calculated. Marker genes for each cluster were calculated using the “FindAllMarkers” function under the following criteria: log2 fold changes (log2FC) > 0.25, min.pct > 0.1, and adjusted P < 0.05. Cells were annotated using both machine-learning-based software SingleR [42] and high expression of canonical markers (i.e., EPCAM for epithelial cells, NCAM1 for neuron cells, VWF for endothelial cells, DCN for fibroblast. PTPRC for immune cells, C1Q for macrophage/Dendritic-cells (MDC), and CD3 for T cells, CD79A for B cells) in each cluster. For the second step of cell annotation, the immune populations were extracted separately and performed the entire workflow of normalization, batch effect correction, unsupervised cell clustering, and cell annotation and identification. The small cell population inner the immune cell populations were finally identified.

Identification of tumor cells in PitNETs

The tumor cells’ identification was described as (i) expression of neuron markers (NCAM1) and canonical markers of PitNETs, such as POU1F1, TBX19, and NR5A1 and (ii) higher copy number variation level. The copy number karyotyping analysis was performed by inferCNV ( and CopyKAT ( [43] package with default parameter.

Differentially expressed genes calculation and transcription factors activities estimation

We employed Seurat's “FindMarkers” feature to detect genes expressed at varying levels between separate cell types. The log fold change (logFC) threshold was set at 0.01, and the minimum percentage of cells expressing the gene (min.pct) was set at 0.01. The significance of DEGs was determined by log2FC > 0.5 and adjusted P-value < 0.05. Protein activities of transcription factors (TFs) were estimated using DoRothEA [44]. A significance criterion of P < 0.05 was employed to identify variations in TF activity across various cell types.

Generation of cell signatures based on scRNA-seq

The top 200 expressed genes (log2FC > 0.58 and adjusted P-value < 0.05) in each cell type were selected as an in-house cell marker database of PitNETs, which was enrolled for analysis on bulk RNA-seq data. The enrichment score to infer cell abundance from bulk RNA-seq was calculated using a single sample gene set enrichment analysis (ssGSEA) algorithm via the R GSVA [45] package.

Trajectory analysis of macrophages

In order to assess the possible dynamic changes in cell state across several subclusters of tumor-associated macrophages (TAM), we employed trajectory analysis using the R package Monocle2 (version 2.24.0) [46]. The DDRTree method was used to reduce dimensionality based on the top 50 DEGs in each TAM subcluster and monocyte. Based on the trajectory results, monocytes were defined as the initiating point of the trajectory. Genes with branch-dependent (CX3CR1+ macro and C1Q+ macro) expressions were identified using the BEAM subprogram, filtered with a q-value < 0.05, and visualized by plot_genes_branched_heatmap function. ScTour and DifussionMap were utilized with the default parameters.

Cell–cell communication analysis

We applied the R package Cellchat to evaluate the interaction weights between tumor cells and TAMs [47]. First, we created a cell chat object with default parameters. The ligand–receptor interaction database we used for analysis was “CellChatDB.human” without additional supplement.

NicheNet was applied to infer the interaction mechanisms between TAMs and malignant cells [48]. We defined the niches of interest for every subcluster. Clustered cells with gene expression over 10% were considered for ligand and receptor interactions. The top 100 ligands and top 1000 targets of differentially expressed genes of “senders” and “receivers” were extracted for paired ligand–receptor activity analysis. When evaluating the regulatory network of TAMs on tumor cells, NR5A1 + tumor cells were considered “receivers”.

In contrast, POU1F1+ tumor cells and TBX19+ tumor cells were used as reference cells to check the regulatory potential of TAMs on tumor cells. The ligand_activity_target_heatmap in Nichenet_output was used to show the regulatory activity of ligands. Activity scores ranged from 0 to 1.

Spatial transcriptomics analysis

The stereo-seq data is preprocessed using the STOmics Analysis Workflow (SAW). The read mapping pipeline uses STAR to align the sequenced reads to the human reference genome (GRCh38). Reads for which mass value (Q = − 10log(err rate)) is less than 15 account for more than 40% of the total number of bases, or N reads with more than 5 bases or containing linker sequences were filtered. Valid CID sequence must completely match chip barcode sequence. The gene quantification pipeline uses Bam2Gem to count the number of reads mapped to each gene. Bin200 was set as the basic unit for further data statistics.

The SpatialFeaturePlot function in Seurat generated spatial feature expression plots (Additional file 2: Table S7). Signature scores were integrated into the “metadata” of the spatial transcriptomics (ST) dataset by calculating the mean expression levels of each gene from the scRNA-seq dataset. The estimation of cellular composition for the spatial transcriptomics spots was performed using the SpatialDecon function, and pie plots visualizing the deconvolution results were generated using the DeconPieplot function in the Cottrazm package (version 0.1.1) [49].

Immunohistochemistry staining of tissue microarray

The tissue microarray (TMA) dataset (n = 128, Additional file 2: Table S3) was used to evaluate the immune cell infiltration. Following deparaffinization and rehydration, heat-induced epitope retrieval (HIER) was performed by submerging the slides in antigen unmasking solution (Solarbio).

After blocking endogenous peroxidase and nonspecific binding sites (0.3% H2O2 and 5% normal goat serum, sequentially), primary antibodies were applied at 4 °C overnight. Slides were incubated with Dako REAL EnVision HRP rabbit/mouse (belong to K5007, DAKO, Glostrup, Denmark) at RT for 20 min, followed by treatment with Dako REAL DAB + CHROMOGEN and Dako REAL substrate buffer (belong to K5007, DAKO, Glostrup, Denmark) to visualize staining signals under light microscopy, finally counterstained using hematoxylin solution. Stained slides were scanned using Ocus (Grundium, Tampere, Finland) and analyzed with Qupath software (see below).

IHC image analysis

Stained slides were scanned using Ocus (Grundium, Tampere, Finland) and analyzed with Qupath software v0.3.0. The built-in stain vector estimator preprocessed images. Cells with shape and stain parameters in each area were identified by build-in cell detection via nucleus stain (hematoxylin). The threshold for positive mean DAB optical density (OD) was determined based on the staining pattern and intensities observed in all photos for each antibody. The percentage of CD45, CD68, CD8, Ki67, and Cleaved Caspase 3 positive cells were quantified using a customized cellular multiplex algorithm. Scripts were produced for the analysis methodology of all the above slide pictures. These scripts were then executed in batches for each set of images and subsequently reviewed by two pathologists who are experts in the field. All quantifications were evaluated blinded to patient clinical information and outcomes.

Multiplex immunohistochemistry and analysis

Formalin-fixed and paraffin-embedded tissue Sects. (3 mm) were de-paraffinized and rehydrated. Next, heat-induced epitope retrieval (HIER) was performed, followed by blocking with 3% hydrogen peroxide in TBST for 10 min and staining with the multiplex mIHC kit (PerkinElmer, NEL861001KT, Shanghai Kelin Institute). Briefly, after the first primary antibody staining, slides were incubated using the HRP-polymer detection system for 10 min, then visualization using Opal TSA working solution (1:100) for another 10 min. Afterward, antigen retrieval was conducted again to prepare the slides for the next antibody. Using this Opal staining method, primary antibodies were applied sequentially. Lastly, slides were counterstained with DAPI (Sigma, 1:1000) for nuclei visualization and subsequently coverslipped using the Hardset mounting media (VectaShield, H-1400).

All tissue sections that underwent multiplex fluorescent staining for each fluorophore were imaged using the Vectra Polaris imaging system (PerkinElmer, Shanghai Kelin Institute) under the appropriate fluorescent filters to produce the spectral library required for multispectral analysis. A whole slide scan of the multiplex tissue sections produced multispectral fluorescent images visualized in Phenochart (PerkinElmer).

Cell culture and reagents

The GH3, AtT20, MMQ, and RC-4BC cell lines were purchased from the American Type Culture Collection (ATCC, Manassas, VA, USA). GH3, AtT20, and MMQ cell line Twere cultured in Ham's F-12 K medium (L450KJ, BasalMedia) supplemented with 2.5% FBS (S615JY, BasalMedia), 15% horse serum (26,050,088, ThermoFisher), and 1% penicillin/streptomycin (C100C5, NCM biotech. RC-4BC cells were cultured in DMEM (L130KJ, BasalMedia) with 10% FBS, 5 ng/ml recombinant rat EGF (ab290070, Abcam), and 1% penicillin/streptomycin. The following inhibitor and recombinant protein were used: recombinant INHBA (C687, Novoprotein), SB-505124 (HY-13521, MedChemExpress), A 83–01 (HY-10432, MedChemExpress), and recombinant human Follistatin (10685-H08H, Sino Biological).

Flow cytometry and cell sorting

For surface staining, cells were resuspended in 50 μl of PBS containing antibody cocktails and stained at room temperature in the dark for 30 min. Antibodies used for flow cytometry are listed (Additional file 2: Table S8). For intracellular staining, cells were fixed and permeabilized by Foxp3 Fixation/Permeabilization kit (Thermo Fisher Scientific) at 4 °C for 45 min and stained with 50 μl of 1 × permeabilization buffer containing antibody cocktails at 4 °C in the dark for 45 min.

For cell sorting, among the immune cells from human pituitary tumor cell suspensions, CD45-positive were divided into T cells (CD45 + CD3 +), B cells (CD45 + CD3-CD19 +), NK cells (CD45 + CD3-CD19-CD56 +), and macrophages (CD45 + CD3-CD19-CD56-CD11b +); among the CD45-negative, non-immune cells were divided into epithelial cells (CD45-EPCAM + CD31-), endothelial cells (CD45-CD31 + EPCAM-), and stromal cells (CD45-CD31-EPCAM-CD90 +). Among CD14-macrophages, CX3CR1 + , C1Q + , and Lyve + were used to classify macrophages, respectively. The gating strategy is shown in Additional file 1: Fig S11.

All flow data were acquired by BD FACSDiva software v8.0.2 and analyzed by FlowJo VX.

TAM and tumor cell co-culture system

We isolated CX3CR1+ TAMs from three independent patients (334, 345, and 356) belonging to the SF1 lineage. The clinical details for these patients are provided in Supplementary table s3. We co-cultured the sorted CX3CR1+ TAMs with matched CD45- tumor cells. Briefly, 1 × 104 TAMs were sorted and seeded onto Transwell polycarbonate filters (0.4 µm pore, 6.5 mm membrane diameter; Corning Incorporated, Corning, NY, USA). Simultaneously, 5 × 105 primary tumor cells were seeded into the lower compartments of Transwell chambers, allowing co-culture with TAMs in the upper compartment for 72 h. The initially identified tumor cells from the lower compartment were then collected for further investigation.

CellTiter-Glo luminescence assay

CellTiter-Glo luminescence assay (Promega, Madison, WI) was used to determine the suppressive effect of the corresponding drug on the control or gene-edited cell line. The cells were distributed at 2000 cells per well on a 96-well plate. The cells were exposed to the appropriate drugs in a 10% FBS-F12K medium at the prescribed doses. Cells were cultured for the indicated time before adding 100 μL of the CellTiter-Glo® luminescence assay reagent in each well. Cells were incubated for an additional 10 min at room temperature to stabilize luminescent signals and transferred to 96-well black plates. Measurements were performed using a luminescence reader (TECAN, Männedorf, Switzerland). Data was analyzed by GraphPad Prism 5 and normalized to the control group. The p-value was calculated using the Student's t-test. Data was generated from at least three independent experiments.

Cell apoptosis

Flow cytometric analysis was conducted using the Annexin V, FITC Apoptosis Detection Kit (AD10, Dojindo) to evaluate cellular apoptosis. Primary tumor cells were harvested and rinsed twice with PBS. Subsequently, the cells were incubated with Annexin V-FITC/PI per the manufacturer’s instructions. The cells were acquired and analyzed using BD FACS Calibur with FlowJo (version 7.6.1) software.

Immunofluorescence and analysis

The cells were derived from samples 375, 381, and 382. Sorted 1 × 104 CX3CR1+ and CX3CR1 macros were centrifuged on a slide by Thermo Scientific Shandon Cytospin. After that, cells were washed once with PBS and then fixed with 4% paraformaldehyde for 5 min at room temperature. Cells were washed thrice with PBS and permeabilized with 0.3% Triton-X for 5 min at room temperature. After three times PBS washing, cells were blocked at room temperature with 5% FBS blocking solution for 30 min. The cells were incubated overnight at 4℃ with anti-INHBA and anti-CD68 primary antibodies in a blocking solution. The cells were then washed for 5 min at room temperature three times with PBS. The cells were then incubated with secondary antibodies for 1 h at room temperature.

Anti-goat FITC for INHBA and anti-rabbit Corelite 594 for CD68 were used. All secondary antibodies were diluted in PBS. After that, the cells were washed at room temperature three times with PBS and stained with Hoechst at 1/1000 (40731ES10, yeasen) for 5 min at room temperature. Finally, the slides were washed three times with PBS and mounted on glass slides with AntiFade Mounting Medium ( G1401-5ML, Servicebio). Image acquisition was performed with an Axio Imager M2 (Carl Zeiss Ltd) and an Apotome. The Zen software piloted 2 (363 oil immersion objective) (Carl Zeiss Ltd). Additional file 2: Table S8 summarizes the antibodies used and how they were diluted.

To quantify the fluorescence intensity, we used ImageJ software with default parameters. The grayscale value assigned to each pixel in the single-channel fluorescence image reflects the corresponding fluorescence intensity. The INHBA mean fluorescence intensity (MFI) was determined using the formula: Mean Fluorescence Intensity (MFI) = Total Fluorescence Intensity in the Region / Area of the Region. The grayscale value and the corresponding area were determined for each region, and the MFI for each photo was calculated. Statistical analysis was performed using a Student’s t-test.

Real-time RT-PCR

Total RNA was extracted from tissue samples and cells using TRIzol reagent (AG21102, Accurate Biology) after washing with PBS. According to the manufacturer’s instructions, cDNA was synthesized from purified RNA using an Evo M-MLV RT Mix Kit with gDNA Clean for qPCR (AG11728, Accurate Biology). SYBR Green PCR Master Mix (AG11718, Accurate Biology) was used for PCR amplification and a real-time PCR machine (ABI-7500, ThermoFisher) was used to quantify the expression of mRNAs. β-actin was used as an endogenous control, and the expression levels were quantified using the 2−ΔΔCt method. All primer sequences are listed in Additional file 2: Table S9, and each primer was detected in triplicate.

Western blotting

Cells and tissues were lysed by RIPA buffer (P0013C, Beyotime Biotechnology) with protease and phosphatase inhibitor cocktail (P002, NCM Biotech), and total protein concentration was measured using a bicinchoninic acid protein assay kit (YSD-500 T, Yoche). Samples were denatured by boiling in 1 × loading buffer and run on a 10% SDS–PAGE gel. Membranes were incubated with primary antibodies overnight at 4 °C, washed three times with TBST, incubated with secondary antibodies for 1 h at room temperature, and developed using Ultra High Sensitive ECL Kit (G2020-25ML, Servicebio). The primary and secondary antibodies used and their dilution are listed in Additional file 2: Table S8. Western blots shown in the accompanying figures are derived from three independent experiments.

Mouse studies

Mouse Studies Athymic nude mice (BALB/cA nu/nu) aged 4 to 5 weeks (Shanghai Jiao Tong University School of Medicine, Shanghai, China) were housed in sterile cages under laminar airflow hoods in a specific pathogen-free room at 22–25° with a 12-h light and 12-h dark schedule and fed with autoclaved chow and water ad libitum.

In vivo subcutaneous xenograft

Female BALB/c nu/nu mice were purchased at 4 weeks of age and AtT20 cells were harvested and washed in PBS, resuspended with Hank's balanced salt solution (HBSS) about 1 × 106 (per side/100 μL) concentrations were injected subcutaneously on the suitable sites per mouse. Treatment was started when the tumor sizes reached approximately 100 mm3; mice were randomized into different groups. An intraperitoneal injection of Activin A (5 mg/kg in 100 μl PBS) or PBS was administered every other day. Tumor volumes (per group) were measured with a digital caliper and calculated as length × width2 × 0.52. Relative tumor volume (RTV) was calculated by (RTV = TVt / TV0, where TV0 is the tumor volume measured when starting drug treatment). The anti-tumor effect of drug treatment was calculated by drug-treated/control (T/C) ratio (T/C = RTVtreated / RTVcontrol × 100%). The mouse viscera were analyzed by HE staining to evaluate the safety of the inhibitor. The growth (Ki-67) and hormone (Pomc) of each tumor were measured by immunofluorescence (IF). The HE staining and IF analyses were performed, and the individuals performing the experiments were unaware of the samples being analyzed.

Statistical analysis

All statistical analyses were performed using R 4.2.1. All heatmaps were generated by the R package pheatmap ( The P-values were two-sided. P-values of less than 0.05 were considered statistically significant. For adjustment of P-values, P-values were adjusted using the Benjamini and Hochberg method (alias false discovery rate, FDR).


Single-cell transcriptional landscape of PitNETs

In order to analyze the distribution of tumor cells and the cells in the tumor microenvironment (TME), we performed single-cell RNA sequencing (scRNA-seq) on the tumor tissues of 23 patients with PitNET. These tumors were classified into 11 PIT1 tumors, 2 silent TPIT tumors, and 10 SF1 tumors based on the expression levels of three transcription factors (Fig. 1A). After stringent quality control, we obtained 69,539 cells for further analysis. The batch effect of the scRNA-seq data across 23 samples was corrected with the Harmony algorithm [39, 50]. Subsequently, a total of 11 prominent cell types were discerned. Canonical markers were used to identify each cell type, and their visualization was performed using t-distributed stochastic neighbor embedding (t-SNE) (Fig. 1B–D) [39]. Despite the tumor cells which were annotated based on lineage-specific transcription factors, namely POU1F1+, TBX19+, and NR5A1+ tumor cells, every cell type appeared across the three lineages (Fig. 1D, Additional File 1: Fig. S1A-C). Besides lineage-specific transcription factors, tumor cell clusters shared expression of epithelial marker gene EPCAM and neuron marker gene NCAM1 (Fig. 1E). The transmembrane protein tyrosine phosphatase gene PTPRC (CD45) is expressed in immune cells. We identified five subtypes of immune cells, including macrophages and DCs expressing C1QA and MS4A7, neutrophils expressing S100A8 and FCGR3B, B cells expressing MS4A1 and CD79A, T/NK cells expressing CD3E and NKG7, and mast cells expressing KIT and MS4A2. Two subtypes of stromal cells were identified: endothelial cells expressing VWF and PLVAP and fibroblasts expressing COL1A2 and FN1 (Fig. 1F). To further validate our findings, flow cytometry analysis was conducted on 45 PitNET samples. The results illustrated that macrophages exhibited the highest expression among CD45+ cells, while epithelial cells dominated the expression among CD45 cells (Additional file 1: Fig. S1D).

Fig. 1
figure 1

Single-cell analysis reveals the transcriptomic landscape in PitNETs. A Graphic overview of the study design and workflow. B Uniform Manifold Approximation and Projection (UMAP) map visualizes 69,539 cells, with color-coding representing individual samples. C UMAP overview color-coded by the three lineages observed. D UMAP overview color-coded by the identified cell types. E Violin plots depict marker gene expression levels across three types of tumor cells and immune cells. F Dot plot illustrating the specific marker genes associated with the 11 major cell types. G Heatmap displaying the chromosomal map of single-cell large-scale copy number variations (CNVs) inferred through scRNA-seq. Immune cells were used as a reference, and amplifications (red) or deletions (blue) were inferred by averaging expression over 100-gene stretches on each chromosome

Tumor heterogeneity can be characterized based on genetic and copy number variations (CNVs). Most PitNET are benign tumors whose subtypes depend on lineage-specific transcription factors [2]. The presence of CNVs in benign tumors and their impact on tumor control remains undisclosed. In order to examine the differences in CNVs among various lineages of PitNET, we measured the CNVs using the inferCNV method [50]. We found that the three lineages of PitNET had more frequent CNV events in tumor cells, suggesting that PitNETs also have many CNV events. PIT1 lineage harbored the highest number of CNVs (Fig. 1G).

We screened differentially expressed genes (DEGs) in three lineages of tumor cells using the FindMarkers function and performed functional enrichment analysis. The results reveal the most significantly up-regulated enriched signal pathways (Additional file 1: Fig. 2A-C). Notably, the POU1F1+ tumor cells exhibited the upregulation of genes involved in intracellular protein trafficking and extracellular protein secretion. The up-regulated DEGs in TBX19+ tumor cells were related to apoptotic signaling pathways and vesicle-mediated transport. NR5A1+ tumor cells showed the upregulation of genes associated with response to extracellular stimuli and regulation of neuronal death.

Fig. 2
figure 2

Tumor cells in PitNETs can be divided into three lineages based on their gene expression profile. A Gene expression patterns and clinical features of three distinct lineages in PitNETs. The columns represent PitNET patients, and the rows indicate gene expression levels or clinical features. B Bubble diagram illustrating the Gene Ontology (GO) biological processes (BP) associated with the three distinct tumor lineages. C Heatmap displaying the activation level of TF and their expression in the three PitNET lineages, as determined by scRNA-seq and bulk RNA-seq methods

Distinct function modules of PitNETs subtypes

Our analysis of gene expression patterns at single-cell levels showed that the tumor cells derived from the PIT1, TPIT, and SF1 lineages of PitNETs have distinct tumor heterogeneity (Fig. 1C–D, and G). In addition, we previously classified six groups (G1 to G6) for PitNET patients based on bulk RNA-seq (n = 180), demonstrating the heterogeneity of PitNETs [51]. To further explore the heterogeneity of PitNET from the tissue and single-cell level, we integrated scRNA-seq with our in-house extended bulk RNA-seq (n = 365) to distinguish the PIT1, TPIT, and SF1 lineages of PitNETs, corresponding to G1–G3, G4–5, and G6, respectively (Fig. 2A). The classification shed light on a more specific clinical treatment. In the PIT1 lineage, the characterized genes were PRL, PRLR, DLK1, ARHGAP36, PLCXD3, GH1, GHRHR, and the transcription-factor-coding POU1F1 (Additional file 1: Fig. S3A). The TPIT lineage had TBX19, POMC, GZMK, SST, SOX3, AVPR1B, CARTPT, and CALB1 genes as feature genes (Additional file 1: Fig. S3B) [6, 51, 52]. Many of the characteristic genes in both PIT1 and TPIT lineages regulate the synthesis and secretion of hormones. On the other hand, the highly expressed genes in the SF1 lineage were NR5A1, GATA2, IDH1, EGR2-4, and NSG2 (Additional file 1: Fig. S3C). These genes primarily participate in tumor development, differentiation, proliferation, and growth. More importantly, the characteristic transcription factor expression levels were coherent between scRNA-seq and bulk RNA-seq data (Fig. 2B). Furthermore, Gene Ontology (GO) enrichment analysis of DEGs for the enriched signal pathway identified key differences among the three lineages of PitNETs cells (Fig. 2C). The PIT1 lineage tumor cells were more involved in glycoprotein metabolic and biosynthetic processes; the TPIT lineage was involved in the transport and exocytosis of ions; and SF1 lineage was involved in gland development and synapse-related process (Fig. 2C). Collectively, both scRNA-seq and our expanded bulk RNA-seq data validated three distinct function modules of PitNETs.

Three distinct subtypes of TIME associated with PitNETs lineages

In addition to specific tumor cell composition across PitNETs lineages, we also found significant differences among the microenvironment. The highest abundance of stromal cells was found in SF1, and the most immune cell infiltration was seen in PIT1 compared to the other two subtypes (Additional file 1: Fig. S3D). However, the specific biological composition and functions of immune cells in the microenvironment of PitNET remain uncertain [2, 53]. Here, we drafted an immune landscape of PitNETs to explore the functional roles of TIME in different lineages. Based on distinct characteristic marker genes, we divided the 12,079 cells into 10 groups from cell clusters that were annotated to be immune cells, including macrophages, monocytes, dendritic cells (pDC and cDC), neutrophils, mast cells, B cells, T cells (CD4+T and CD8+T), and NK cells (Fig. 3A-B, Additional file 1: Fig. S4A-B). Notably, macrophages had C1QA, C1QB, and CD163 as feature genes, while highly expressed S100A10, LYZ, and VCAN were observed in the monocyte cluster. Although each immune cell subset was present in most samples, we observed significant differences in their proportions in different lineages (Additional file 1: Fig. S4C-D). CD4+ T cells, CD8+ T cells, and NK cells were enriched in PIT1; mast cells were increased in SF1; neutrophils were enriched in both PIT1 and SF1, while macrophages showed the highest infiltration in TPIT.

Fig. 3
figure 3

Profiling immune cell populations in PitNET. A Uniform Manifold Approximation and Projection (UMAP) map displaying 12,079 immune cells, with color-coding indicating individual cell types. B Dot plot illustrating the specific marker genes associated with the 10 major immune cell types. C–F The unsupervised clustering of PitNET RNA-seq cohorts based on the abundance of immune cells revealed three distinct TIME subtypes. The plot shows the abundance of immune cells in each subtype (C). Expression levels of immune checkpoint genes in the PitNET TIME subtypes (D). The pathways are enriched based on the DEGs (E, F). P-values were calculated using ANOVA and adjusted using FDR. G The composition of immune classification in samples from groups G1 to G6. H Bar plot indicating the proportion of positive cells stained for CD45, CD68, and CD8 in PitNET Tissue Microarray (TMA) using IHC with G1 (n = 6) and G6 (n = 52). P-values were calculated using the Wilcoxon rank-sum test. Results are presented as mean ± standard error of the mean in bar plots. I–J. Bar plot displaying the proportion of CD68/CD45 positive cells stained in PitNETs TMA using IHC with PIT1 (n = 21) and SF1 (n = 45) (I) and flow cytometry with PIT1 (n = 12) and SF1 (n = 19) (J). P-values were calculated using the Wilcoxon rank-sum test. Results are presented as mean ± standard error of the mean in bar plots

To further investigate the TIME differences among PitNET lineages, we performed unsupervised clustering of PitNETs RNA-seq cohorts based on immune cell infiltration, identifying three distinct TIME subtypes (Additional file 2: Table S10), namely the “Immune Low” (IL), “Immune Intermediate” (II), and “Immune High” (IH) groups (Fig. 3C). The degree of immune cells infiltration increases progressively from group IL to IH (Additional file 1: Fig. S5A). Immune checkpoints are essential in immune responses and are critical as immunotherapy targets [54]. Either inhibitory (e.g., LAG3, HAVCR2) or stimulatory checkpoints (e.g., TNFSF9, TNFNR5A18) rose steadily from group IL group to IH, suggesting a co-regulation of stimulatory and inhibitory immune signaling pathways to maintain TIME equilibrium among different PitNET lineages. Furthermore, single-sample Gene Set Enrichment Analysis (ssGSEA) highlighted distinct signaling pathways enriched in the IH group (e.g., NFKB, PI3K-AKT-mTOR, NOTCH, WNT, P53), as well as cellular processes (e.g., EMT, hypoxia, and apoptosis). Moreover, the IH group exhibited immune-related enrichment of IFN-γ response, IFN-α response, and macrophage infiltration, indicating higher immune activity and immune cell infiltration in these tumors. These enriched pathways regulate cellular functions, including growth, metabolism, survival, differentiation, epithelial-mesenchymal transition, and apoptosis.

Furthermore, we found that the IH group accounted for 82.76% of G1 samples, while 55.42% of G6 samples belonged to the IL group (Fig. 3G). However, the three PitNET lineages showed different proportions of immune infiltration (Additional file 1: Fig. S5B). Further, we conducted TMA-IHC analysis, validating that the G1 subtype has significantly higher infiltration of immune cells, including CD45+ cells, CD68+ cells, and CD8+ T cells, compared to the G6 subtype (Fig. 3H). Intriguingly, compared to the G1 subtype, the SF1 lineage (G6), primarily characterized as “IL” tumors, displayed increased infiltration of TAMs (Fig. 3I–J). These findings indicate that PitNETs can be categorized into three TIME subtypes. The IH group was predominantly observed in the G1 subtype, while the IL group was primarily observed in the G6 subtype with a higher infiltration of macrophages.

Distinct subtypes of tumor-associated macrophage in different PitNETs lineages

Our scRNA-seq data showed macrophages were the major immune cell population in PitNET tissue (Fig. 1D), and distinctive infiltration was seen in different PitNET lineages (Fig. 3I–J). To explore the functions of macrophages, we clustered the 5708 macrophages from 23 patients and identified 5 subclusters, including C1Q+, ISG15+, LYVE1+, GPNMB+, and CX3CR1+macropahges (Fig. 4A and Additional file 1: Fig. S6A-B). These macrophage subtypes had diverse distributions across the three PitNET lineages (Fig. 4C–D and Additional file 1: Fig. S6C-D); C1Q+ macrophages were mainly enriched in TPIT lineage tumors; GPNMB+ macrophages were enriched in PIT1 lineage tumors, and CX3CR1+ macrophages were enriched in SF1 lineage tumors. The distribution of macrophage subtypes in different PitNET lineages was validated using in-house bulk RNA-seq data based on a signature score of the top macrophage subtype-specifically expressed genes (Fig. 4E, Additional file 2: Table S11). We further examined the localization of C1Q+, GPNMB+, and CX3CR1+ macrophages in TPIT, PIT1, and SF1 lineages, respectively. Our multiple immunofluorescence labeling analysis provided strong evidence of the co-localization of CX3CR1+ macrophages and CD68 in SF1 lineage (Fig. 4F). Similarly, GPNMB+ macrophages were mainly co-localized with CD68 in PIT1 lineage (Fig. 4G). In contrast, C1Q+ macrophages were mostly co-localized with CD68 in the TPIT lineage (Fig. 4H).

Fig. 4
figure 4

Three distinct subtypes of tumor-associated macrophage in three different PitNETs lineages. A UMAP map depicting macrophage cells color-coded by subtypes. B Dot plot displaying the specific marker genes for the five subtypes of macrophages. C The percentage of the five subtypes of macrophages observed in the three lineages. P-values are calculated using ANOVA. D–E. Violin plots show marker genes (D) expression levels and cell abundance (E) for the five subtype macrophages of three lineages. P-values were calculated using the Wilcoxon rank-sum test. F–H. Immunofluorescence (IF) image demonstrating the cell distribution of CX3CR1+ (F), GPNMB+ (G), and C1Q+ (H) macrophages in the three lineage tumor cells. I. Flow cytometry data presents the proportion of CX3CR1+ macrophages in the population across the three lineages. P-values were calculated using the Wilcoxon rank-sum test. Results are presented as mean ± standard error of the mean in bar plots. J Trajectory plot of monocytes and the five subtypes of macrophages in a two-dimensional state-space inferred by Monocle 2 analysis (up) and the transition trajectories along pseudotime (down). K Three-phase distribution of C1Q+ macrophages and CX3CR1+ macrophages along pseudotime. Heatmap depicting the dynamic expression changes of genes and related pathways of C1Q+ macrophages and CX3CR1+ macrophages along pseudotime

Furthermore, trajectory analysis using three different methods witnessed a developmental path from monocytes to LYVE+ and GPNMB+ macrophages and finally to CX3CR1+ and C1Q+ macrophages (Fig. 4J, and Additional file 1: Fig. S7). Gene Ontology (GO) enrichment analysis on the genes enriched in each state revealed the different functions of macrophage subtypes. Furthermore, we performed KEGG pathway enrichment analysis to uncover distinct functions of five macrophage subtypes (Additional file 1: Fig. S8). The ribosomal signaling pathway supporting protein synthesis was enriched in TPIT-lineage-associated C1Q+ macrophages and thus might help to maintain elevated hormone secretion. The oxidative phosphorylation pathway providing energy and materials was enriched in POU1F1-lineage-associated GPNMB+ macrophages. The pathway of negative regulation for transforming growth factor beta (TGF-β) was enriched in SF1-lineage-associated CX3CR1+ macrophages and helped control the tumor growth. C1Q+ macrophages (characterized by the C1QA/B/C gene) showed a transitional phenotype from the initiate state to the immune-stimulated state, involving the upregulation of HLA genes and stimulation of the interferon signal pathway (Fig. 4K). The two clusters of genes exhibiting high expressions in CX3CR1+ macrophages play a critical role in immune response and cytokine–cytokine receptor interaction.

Different regulation of macrophage-tumor cell interaction axis in PitNETs

To understand the biology of macrophage–tumor cell interactions, we performed the CellChat analysis along the macrophage–tumor cell axis in different PitNET lineages. The results demonstrated that CX3CR1+, C1Q+, and GPNMB+ macrophages exhibited the highest interaction with NR5A1+, TBX19+, and POU1F1+ tumor cells, respectively (Fig. 5A). This highlighted distinct ligand–receptor pairs at each macrophage–tumor cell axis (Fig. 5B). To explore the mediators and downstream targets of the macrophage–tumor cell axis, we performed NichenetR analysis. We found that CX3CR1+ macrophages regulated NR5A1+ tumor cells through regulation of signaling receptor pathway and TGF-β stimulus (Fig. 5D). The target genes and cluster-enriched pathways of C1Q+ and GPNMB+ macrophages were associated with protein maturation, response to corticosteroid and JAK-STAT, T cell proliferation and activation, and response to cAMP functions (Additional file 1: Fig. S9A-D). Additionally, we investigated the functional profiles of the three macrophage clusters. CX3CR1+ macrophages exhibited higher activity in pro-inflammatory responses and angiogenesis. In contrast, C1Q+ macrophages demonstrated a high level of phagocytosis (Fig. 5E). Conversely, GPNMB+ macrophages were involved in angiogenesis and phagocytosis (Fig. 5E). Furthermore, CX3CR1+ macrophages displayed a pro-inflammatory profile, while GPNMB+ macrophages represented to be anti-inflammatory. C1Q+ macrophages exhibited minimal regulation of inflammation. Gene expression analysis within each macrophage cluster demonstrated that CX3CR1+ macrophages showed the upregulation of pro-inflammatory genes, including IL1B, TNF, CXCL9, and CXCL10 (Fig. 5F). This suggested that CX3CR1+ macrophages contributed to the pro-inflammatory TIME of PitNETs. Furthermore, we observed a downregulation of CX3CR1+ macrophages in cases of MIB1-high and cavernous invasion PitNETs (Fig. 5G). Together, these findings suggested a potential association of reduced numbers of CX3CR1+ macrophages with the increased proliferation and aggressive characteristics of SF1 lineage PitNETs.

Fig. 5
figure 5

CX3CR1+ macrophages interact with NR5A1+ tumor cells. A Interaction weights between TAM subpopulations and the three lineages of tumor cells. B Ligands and receptor pairs are involved in the interaction between TAM subpopulations and the three lineages of tumor cells. C Heatmap displaying the expression levels of ligands highly expressed in CX3CR1+ macrophages (left) and the expression levels of corresponding target genes on NR5A1+ tumor cells (right). D Signal pathways enriched from the ligand–receptor interaction between CX3CR1+ macrophages and NR5A1+ tumor cells. E Evaluation of the anti/pro-inflammatory and angiogenesis/phagocytosis functions of the three subtypes of macrophages. F Heatmap representing the mean expression levels of anti-tumor genes in the three subtypes of macrophages. G Violin plots illustrating the cell abundance of CX3CR1+ macrophages in samples categorized as low (< 3) or high (≥ 3) MIB1 percentage (on the left), as well as in samples with or without invasion. P-values were calculated using the Wilcoxon rank-sum test

The INHBA-ACVR1B axis promotes tumor cell apoptosis

To investigate the regulatory role of CX3CR1+ macrophages on NR5A1+ tumor cells, we co-cultured PitNETs isolated CX3CR1+ or CX3CR1 macrophages with autologous CD45 tumor cells. We found that CX3CR1+ macrophages enhanced the expression of the corresponding target tumor inhibitory genes, such as BTG2, EGR2, ERG3, NR0B1, and SDC4, on NR5A1+ tumor cells (Fig. 6A). INHBA-ACVR1B was found to be the top-ranked ligand-receptor pair among CX3CR1+ macrophages and NR5A1+ tumor cells interaction. Notably, CX3CR1+ macrophages enriched NR5A1+ tumor cells and INHBA-ACVR1B were locationalized closely in the spatial transcriptome of SF1 lineage tissue (Fig. 6B–D and Additional file 1: Fig. S10A). Furthermore, we found that INHBA protein treatment up-regulated three tumor inhibitory genes, namely EGR2, ERG3, and NR0B1 (Fig. 6E). Also, immunofluorescence analysis on flow-sorted CX3CR1+ macrophages validated the expression of INHBA in the CX3CR1+ macrophages population (Fig. 6F–G).

Fig. 6
figure 6

CX3CR1+ macrophages inhibit the growth of NR5A1 + tumor cells through INHBA/activin A—ACVR1B axis. A Histogram plot revealed the gene expression alteration of the tumor cells co-cultured with CX3CR1+ macrophages (n = 3). B–D HE images with pie charts in each spot colored by annotation showing the cellular composition in the ST sample (B). The spatial feature plot shows the location and proportions of macrophages in ST spots using deconvolution results (C). Spatial feature plots show the expression of INHBA and its receptor ACVR1B in the macrophage-containing spots (D). E Histogram plot revealed the gene expression alteration of the tumor cells co-cultured with activin A. F–G The representative image of immunofluorescence revealing the specific expression of INHBA in CX3CR1+ and CX3CR1 macrophages (F) and the statistics of INHBA mean fluorescence intensity (MFI) in CX3CR1+ and CX3CR1 macrophages in three patients (G). H–I The cell viability (H) and apoptosis (I) in the four NR5A1+ primary cell treated activin A (10 ng/ml) with or without SB-505124 (1 μM) or follistatin (100 ng/ml). The histogram shows the cells' apoptosis in different conditions for the primary cell (I). J Cell viability assessment in the AtT20 cell line treated with Activin A (10 ng/ml) in the presence or absence of SB-505124 (1 μM), A 83–01 (1 μM), or follistatin (100 ng/ml). K Cell viability analysis in AtT20 cell line with shEV-GFP or shAcvr1b treated with Activin A (10 ng/ml). L–M. Apoptosis assessment in the AtT20 cell line treated with Activin A (10 ng/ml) in the presence or absence of SB-505124 (1 μM), A 83–01 (1 μM), or follistatin (100 ng/ml). The histogram illustrates the apoptosis levels of cells under various conditions in the AtT20 cell line (M). N–P Tumor growth curve (N). The xenografts harboring mice were injected intraperitoneally with activin A release preparations, 5 mg/kg weight, every other day with PBS. Representative images of H&E and IF staining (Ki-67 and Cleaved-caspase 3) of resected AtT20 tumors from an experimental study between different groups (O). Scale bar as indicated. Statistical analysis of Ki67 and Cleaved-caspase 3 percentages in control and treated Att20 tumors (P). Q The staining revealed the co-localization of CD68 and CX3CR1 in xenografts harboring mice. Arrow indicated the CD68 and CX3CR1 positive staining. Scale bar as indicated. *P < 0.05, **P < 0.01, ***P < 0.001, ****P < 0.0001. P-values were calculated using the Wilcoxon rank-sum test. Results are presented as mean ± standard error of the mean in bar plots

To confirm the role of INHBA on the CX3CR1macrophages NR5A1+ tumor cells interaction axis, we treated primary cells and PitNET cell lines with INHBA. We found that INHBA inhibited cell viability in SF1 lineage primary cells (Fig. 6H) and induced apoptosis (Fig. 6I). We further used SB-505124, which reversed the inhibitory effect of INHBA (Fig. 6H and I). Similar outcomes were shown in the AtT20 cell line, where INHBA caused apoptosis and decreased cell proliferation. Both effects may be reversed by follistatin and reversible ATP competitors that are selective for ALK4 and ALK5 (SB-505124 and A 83–01) (Fig. 6J, L, and M). To further confirm the dependency of the inhibitory effect on the INHBA receptor Acvr1b, we conducted experiments demonstrating a significant reversal of the inhibitory effect upon downregulating Acvr1b using shRNA (Fig. 6K). However, no inhibitory effect of INHBA was observed in the MMQ and GH3 PitNET cell lines (Additional file 1: Fig. S10B). This can be attributed to the low expression of ACVR1B (Additional file 1: Fig. S10C). Moreover, in vivo drug treatment experiments were subsequently performed to confirm the validity of the above results. Tumor growth of AtT20 xenografts in a mouse subcutaneous model was reduced by continuous therapy with INHBA (activin A) (Fig. 6N). The hematoxylin and eosin (H&E) staining results confirmed that activin A treatment induced similar tumor necrosis in the AtT20 xenografts. However, the Ki-67 and Cleaved-caspase 3 staining results showed decreased proliferation and increased apoptosis in treatment with activin A (Fig. 6O and P). To confirm the expression of CX3CR1+ macrophages in AtT20 xenografts, immunofluorescence (IF) staining results confirmed the co-localization of CD68 and CX3CR1 in the tissue (Fig. 6Q). Consequently, these results propose that CX3CR1+ macrophages may trigger apoptosis in tumor cells through the INHBA-ACVR1B axis.


Different phenotypes of TIME, such as immune-inflamed microenvironment, immune-excluded microenvironment, and immune-desert microenvironment, have been associated with multiple cancer types and immunotherapy efficacy [21, 55]. In particular, many recent studies have revealed that different macrophage subtypes execute specific functions in tumors, and therefore, macrophage-targeted therapy is an attractive immunotherapy strategy. Despite the emerging studies [21, 56] on the components of tumor immune microenvironment (TIME), the TIME of PitNETs is not fully deciphered at the transcriptomic level, limiting the progress of clinical research on PitNETs immunotherapy. This study provides a comprehensive analysis of the immune landscape and macrophage heterogeneity in PitNETs, highlighting the potential roles of specific macrophage populations in tumor progression. We performed single-cell, bulk, and spatial RNA-seq analyses of different PitNET lineage samples to identify tumor heterogeneity and distinguish function modules of PitNET subtypes. Finally, three distinct TIME subtypes associated with PitNETs lineages were identified. TAM subtypes interact with different PitNET lineages; CX3CR1+ TAMs interact with NR5A1+ tumor cells through the INHBA-ACVR1B axis, promoting tumor cell apoptosis.

Our previous study used bulk RNA-seq to classify 180 PitNET patient samples from three lineages into 6 subtypes (G1 to G6), revealing the heterogeneity of the PitNET lineage [51]. Using proteomics, Zhang et al. also identified distinct clusters of PitNET heterogeneity showing different functional associations [7]. Here, we further highlighted the tumor heterogeneity and TIME of three transcription factor(TF)-based lineages in PitNETs from a single-cell perspective. Our results complement previous studies concerning the heterogeneity of different lineages in PitNETs [3, 5, 6]. We identified that different lineages of PitNET not only have distinct characteristics in TF-related transcriptomic profiles and functional modules but also exhibit genetic variation and CNVs. We screened differentially expressed genes (DEGs) in three types of tumor cells. Among these, the POU1F1+ tumor cells demonstrated intracellular protein trafficking and extracellular protein release. This aligned with the clinical manifestation of abnormally high secretion of certain hormones in the PIT1 lineage [5]. A recent study also supported these findings and demonstrated an inverse correlation between differentiation status and recurrence within the PIT1/TPIT lineage compared to the SF1 lineage [57]. In particular, we combined 23 single-cell and 365 in-house bulk RNA-seq to reveal that all 3 PitNET subtypes have distinct immune-infiltrated patterns. The G1 subtype (PIT1 lineage) exhibited the highest immune infiltration, while the G6 subtype (SF1 lineage) showed the lowest immune infiltration. Concisely, the molecular phenotypes of PitNETs can be characterized by the immune infiltration phenotypes, providing a novel insight into the clinical diagnosis of PitNETs. This also suggests that the functional differentiation between G1-3 (PIT1 lineage) and G4-5 (TPIT lineage) subgroups may be controlled by TIME.

The heterogeneity of macrophages endows them with diverse functions, such as phagocytosis, antigen presentation, cytokine production, and tissue remodeling [20]. Using scRNA-seq, previous studies identified multiple TAM subtypes with distinct gene signatures and functions [58, 59]. Earlier studies often determined the proportion of macrophages by employing pan markers such as CD68, CD163, or CD204. These indicators investigated possible relationships between treatment response and tumor progression [60]. M2 macrophage infiltration is the primary immunological landscape that the PitNETs subtypes exhibit [61]. The cell-to-cell interactions of pitNETs and macrophages indicate potential therapeutic adjuvant for PitNET treatment [27, 62]. However, current research does not always point to successful treatment outcomes when using classical immune checkpoints like PD-L1 CTLA4 as a target [55, 63]. Here, we identified five distinct macrophage clusters that were differentially distributed and interacted with tumors of three PitNET lineages. Our research indicates that the traditional method of classifying TAMs into distinct M1 and M2 subtypes may not be fully relevant when studying PitNETs. TAMs in PitNETs exhibit features consistent with M1 and M2 phenotypes, disregarding the simplistic dichotomy commonly used in macrophage classification [64]. However, according to a previous study, most TAMs expressed M2-related genes [61]. CX3CR1+ macrophages displayed the upregulation of pro-inflammatory properties and intense interaction with NR5A1+ tumor cells; C1Q+ macrophages showed increased phagocytosis activity and strong interaction with TBX19+ tumor cells; GPNMB+ macrophage with an intermediate level of angiogenesis and phagocytosis showed preferential interactions with POU1F1+ tumor cells. Concisely, different macrophage subtypes in PitNET have different functions and drive differentiated macrophage–tumor interaction axis. Gaining insight into the unique functions of various macrophage subsets and discerning the disparities between different types of PitNETs enhance the development of therapeutic approaches for PitNETs. CX3CR1+ macrophages were reported to have both anti- or pro-tumor effects in malignant tumors, and therefore, their functions are controversial [65,66,67,68]. Here, we demonstrated that CX3CR1+ macrophages exhibit anti-tumor effects, displaying a pro-inflammatory phenotype.

Furthermore, the expression of anti-tumorigenic genes (ERG2, ERG3, and BTG2) was promoted by the interaction between CX3CR1+ macrophages and NR5A1+ tumor cells through the ligand-receptor INHBA-ACVR1B pair [69, 70]. This was in agreement with earlier research [71] that ACVR1B can restore activin antiproliferative effects in human pituitary tumor cells. Interleukin-1 (IL-1), interferon-γ (IFN-γ), and tumor necrosis factor-α (TNF-α) were among the pro-inflammatory cytokines that were reduced by blocking FKN-CX3CR1 interaction, suggesting that this pathway could be a target for the therapy of rheumatoid arthritis [69]. Treatments for our condition may be viewed in the opposite direction. Our findings suggest a potential inhibitory role of CX3CR1+ macrophages in tumor growth within the SF1 lineage of PitNETs. Clinically, INHBA agonists can promote the INHBA-ACVR1B pathway to eliminate SF1 lineage pituitary tumors or increase the infiltration of CX3CR1+ macrophages in other lineages, helping PitNET treatment.

Our current analysis is limited by the scarcity of secreting TPIT lineage tumors in the scRNAseq dataset, as only two silent TPIT available samples exist. The primary reason for this limitation is the low frequency of secreting TPIT-related cases, which makes it challenging to gather a large enough sample size for thorough research. Due to the infrequency of secreting TPIT lineage tumors, it is crucial to use caution when interpreting results related to this specific subtype. Furthermore, limitations in T cells analysis have been identified, characterized by inadequate data quantity and inconclusive grouping. To overcome these problems, we suggest enhancing the presence of T cells in sequencing data by enriching CD45 in future investigations. This strategic improvement will enable more accurate categorization of immune cells and thorough analysis of the functions of T cells, therefore enhancing our understanding of the immune response in the tumors under investigation.

In summary, integrated analysis of bulk, single-cell, and spatial transcriptomic profiles on PitNET offered a comprehensive perspective on the diverse cellular compositions of PitNETs (Fig. 7). The identification of immune-infiltrating phenotypes, TAM subtypes, and distinct interaction axes of TAM-tumor cells provided valuable resources for future in-depth studies. Mainly, CX3CR1+ macrophages may hold significant promise in regulating tumor angiogenesis and inflammation. Here, our findings enable the categorization of molecular phenotypes of PitNETs according to the infiltration of immune cells and imply the possibility of treating PitNETs with macrophage subtype-based therapies. Additional investigation is required to elucidate the functions of different immune cells during the TIME of PitNETs.

Fig. 7
figure 7

Working model of immune-infiltrating phenotypes and macrophage–tumor interaction axes in diverse lineages of pituitary neuroendocrine tumors


We conducted an in-depth exploration of the immune microenvironment within PitNETs, with a specific focus on the interactions between macrophages and tumor cells. Utilizing single-cell RNA sequencing analysis, our study unveiled three distinct subtypes of the TIME, each characterized by varying degrees of immune cell infiltration and macrophage activity. These findings underscore the significance of considering the tumor microenvironment in the therapeutic approach to PitNETs, particularly highlighting potential strategies targeting TAMs. Moreover, our research identified specific axes of interaction between macrophages and tumor cells, which may serve as novel targets for therapeutic intervention. The identification of these interaction axes not only enhances our understanding of the complexity of the immune microenvironment in PitNETs but also provides valuable insights for future clinical interventions. These discoveries hold promise for stimulating further research aimed at uncovering the potential of immunotherapy in treating PitNETs.

Availability of data and materials

The raw sequence data reported in this paper have been deposited in the Genome Sequence Archive in National Genomics Data Center, China National Center for Bioinformation/Beijing Institute of Genomics, Chinese Academy of Sciences (GSA-Human: HRA005096) that are publicly accessible at [70]. Code used in the paper is available on GitHub ( [71] or FigureShare ( [72].


  1. Melmed S, Kaiser UB, Lopes MB, Bertherat J, Syro LV, Raverot G, Reincke M, Johannsson G, Beckers A, Fleseriu M, et al. Clinical biology of the pituitary adenoma. Endocr Rev. 2022;43:1003–37.

    Article  PubMed  PubMed Central  Google Scholar 

  2. Asa SL, Mete O, Perry A, Osamura RY. Overview of the 2022 WHO classification of pituitary tumors. Endocr Pathol. 2022;33:6–26.

    Article  CAS  PubMed  Google Scholar 

  3. Zhang Z, Zamojski M, Smith GR, Willis TL, Yianni V, Mendelev N, Pincas H, Seenarine N, Amper MAS, Vasoya M, et al. Single nucleus transcriptome and chromatin accessibility of postmortem human pituitaries reveal diverse stem cell regulatory mechanisms. Cell Rep. 2022;38: 110467.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  4. Svensson V, Vento-Tormo R, Teichmann SA. Exponential scaling of single-cell RNA-seq in the past decade. Nat Protoc. 2018;13:599–604.

    Article  CAS  PubMed  Google Scholar 

  5. Zhang Q, Yao B, Long X, Chen Z, He M, Wu Y, Qiao N, Ma Z, Ye Z, Zhang Y, et al. Single-cell sequencing identifies differentiation-related markers for molecular classification and recurrence prediction of PitNET. Cell Reports Med. 2023;4: 100934.

    Article  CAS  Google Scholar 

  6. Cui Y, Li C, Jiang Z, Zhang S, Li Q, Liu X, Zhou Y, Li R, Wei L, Li L, et al. Single-cell transcriptome and genome analyses of pituitary neuroendocrine tumors. Neuro Oncol. 2021;23:1859–71.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  7. Zhang F, Zhang Q, Zhu J, Yao B, Ma C, Qiao N, He S, Ye Z, Wang Y, Han R, et al. Integrated proteogenomic characterization across major histological types of pituitary neuroendocrine tumors. Cell Res. 2022;32:1047–67.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  8. Zhang D, Hugo W, Bergsneider M, Wang MB, Kim W, Vinters HV, Heaney AP. Single-cell RNA sequencing in silent corticotroph tumors confirms impaired POMC processing and provides new insights into their invasive behavior. Eur J Endocrinol. 2022;187:49–64.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  9. Lyu L, Jiang Y, Ma W, Li H, Liu X, Li L, Shen A, Yu Y, Jiang S, Li H, et al. Single-cell sequencing of PIT1-positive pituitary adenoma highlights the pro-tumour microenvironment mediated by IFN-γ-induced tumour-associated fibroblasts remodelling. Br J Cancer. 2023;128:1117–33.

    Article  PubMed  PubMed Central  Google Scholar 

  10. Wang P-F, Wang T-J, Yang Y-K, Yao K, Li Z, Li YM, Yan C-X. The expression profile of PD-L1 and CD8+ lymphocyte in pituitary adenomas indicating for immunotherapy. J Neurooncol. 2018;139:89–95.

    Article  PubMed  Google Scholar 

  11. Xi Z, Jones PS, Mikamoto M, Jiang X, Faje AT, Nie C, Labelle KE, Zhou Y, Miller KK, Soberman RJ, Zhang X. The upregulation of molecules related to tumor immune escape in human pituitary adenomas. Front Endocrinol. 2021;12: 726448.

    Article  Google Scholar 

  12. Raverot G, Ilie MD, Lasolle H, Amodru V, Trouillas J, Castinetti F, Brue T. Aggressive pituitary tumours and pituitary carcinomas. Nat Rev Endocrinol. 2021;17:671–84.

    Article  PubMed  Google Scholar 

  13. Ilie MD, Vasiljevic A, Jouanneau E, Raverot G. Immunotherapy in aggressive pituitary tumors and carcinomas: a systematic review. Endocr Relat Cancer. 2022;29:415–26.

    Article  CAS  PubMed  Google Scholar 

  14. Principe M, Chanal M, Ilie MD, Ziverec A, Vasiljevic A, Jouanneau E, Hennino A, Raverot G, Bertolino P. Immune landscape of pituitary tumors reveals association between macrophages and gonadotroph tumor invasion. J Clin Endocrinol Metab. 2020;105:520.

    Article  Google Scholar 

  15. Marques P, Barry S, Carlsen E, Collier D, Ronaldson A, Awad S, Dorward N, Grieve J, Mendoza N, Muquit S, et al. Chemokines modulate the tumour microenvironment in pituitary neuroendocrine tumours. Acta Neuropathol Commun. 2019;7:172.

    Article  PubMed  PubMed Central  Google Scholar 

  16. Mei Y, Bi WL, Agolia J, Hu C, Giantini Larsen AM, Meredith DM, Al Abdulmohsen S, Bale T, Dunn GP, Abedalthagafi M, Dunn IF. Immune profiling of pituitary tumors reveals variations in immune infiltration and checkpoint molecule expression. Pituitary. 2021;24:359–73.

    Article  CAS  PubMed  Google Scholar 

  17. Lin S, Han C, Lou X, Wu ZB. Hydroxychloroquine overcomes cabergoline resistance in a patient with Lactotroph Pituitary neuroendocrine tumor: a case report. Front Endocrinol (Lausanne). 2022;13: 955100.

    Article  PubMed  Google Scholar 

  18. Heshmati HM, Kujas M, Casanova S, Wollan PC, Racadot J, Van Effenterre R, Derome PJ, Turpin G. Prevalence of lymphocytic infiltrate in 1400 pituitary adenomas. Endocr J. 1998;45:357–61.

    Article  CAS  PubMed  Google Scholar 

  19. Yagnik G, Rutowski MJ, Shah SS, Aghi MK. Stratifying nonfunctional pituitary adenomas into two groups distinguished by macrophage subtypes. Oncotarget. 2019;10:2212–23.

    Article  PubMed  PubMed Central  Google Scholar 

  20. Mantovani A, Allavena P, Marchesi F, Garlanda C. Macrophages as tools and targets in cancer therapy. Nat Rev Drug Discov. 2022;21:1–22.

    Article  Google Scholar 

  21. Ilie MD, Vasiljevic A, Bertolino P, Raverot G. Biological and therapeutic implications of the tumor microenvironment in pituitary adenomas. Endocr Rev. 2023;44:297–311.

    Article  PubMed  Google Scholar 

  22. Lu JQ, Adam B, Jack AS, Lam A, Broad RW, Chik CL. Immune cell infiltrates in pituitary adenomas: more macrophages in larger adenomas and more T cells in growth hormone adenomas. Endocr Pathol. 2015;26:263–72.

    Article  CAS  PubMed  Google Scholar 

  23. Han C, Lin S, Lu X, Xue L, Wu ZB. Tumor-associated macrophages: new horizons for pituitary adenoma researches. Front Endocrinol (Lausanne). 2021;12: 785050.

    Article  PubMed  Google Scholar 

  24. Marques P, Barry S, Carlsen E, Collier D, Ronaldson A, Grieve J, Dorward N, Mendoza N, Nair R, Muquit S, et al. The expression of neural cell adhesion molecule and the microenvironment of pituitary neuroendocrine tumours. J Neuroendocrinol. 2021;33: e13052.

    Article  CAS  PubMed  Google Scholar 

  25. Luo M, Tang R, Wang H. Tumor immune microenvironment in pituitary neuroendocrine tumors (PitNETs): increased M2 macrophage infiltration and PD-L1 expression in PIT1-lineage subset. J Neurooncol. 2023;163:663.

    Article  CAS  PubMed  Google Scholar 

  26. Fujiwara K, Yatabe M, Tofrizal A, Jindatip D, Yashiro T, Nagai R. Identification of M2 macrophages in anterior pituitary glands of normal rats and rats with estrogen-induced prolactinoma. Cell Tissue Res. 2017;368:371–8.

    Article  CAS  PubMed  Google Scholar 

  27. Zhang A, Xu Y, Xu H, Ren J, Meng T, Ni Y, Zhu Q, Zhang W-B, Pan Y-B, Jin J, et al. Lactate-induced M2 polarization of tumor-associated macrophages promotes the invasion of pituitary adenoma by secreting CCL17. Theranostics. 2021;11:3839–52.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  28. Sato M, Tamura R, Tamura H, Mase T, Kosugi K, Morimoto Y, Yoshida K, Toda M. Analysis of tumor angiogenesis and immune microenvironment in non-functional pituitary endocrine tumors. J Clin Med. 2019;8:695.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  29. Patro R, Duggal G, Love MI, Irizarry RA, Kingsford C. Salmon provides fast and bias-aware quantification of transcript expression. Nat Methods. 2017;14:417–9.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  30. Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15:550.

    Article  PubMed  PubMed Central  Google Scholar 

  31. Ritchie ME, Phipson B, Wu D, Hu Y, Law CW, Shi W, Smyth GK. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. 2015;43: e47.

    Article  PubMed  PubMed Central  Google Scholar 

  32. Subramanian A, Tamayo P, Mootha VK, Mukherjee S, Ebert BL, Gillette MA, Paulovich A, Pomeroy SL, Golub TR, Lander ES, Mesirov JP. Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc Natl Acad Sci U S A. 2005;102:15545–50.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  33. Yu G, Wang LG, Han Y, He QY. clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS. 2012;16:284–7.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  34. Liberzon A, Birger C, Thorvaldsdottir H, Ghandi M, Mesirov JP, Tamayo P. The molecular signatures database (MSigDB) hallmark gene set collection. Cell Syst. 2015;1:417–25.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  35. Kanehisa M, Goto S. KEGG: kyoto encyclopedia of genes and genomes. Nucleic Acids Res. 2000;28:27–30.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  36. Kanehisa M, Sato Y, Kawashima M, Furumichi M, Tanabe M. KEGG as a reference resource for gene and protein annotation. Nucleic Acids Res. 2016;44:D457–462.

    Article  CAS  PubMed  Google Scholar 

  37. Kuleshov MV, Jones MR, Rouillard AD, Fernandez NF, Duan Q, Wang Z, Koplev S, Jenkins SL, Jagodnik KM, Lachmann A, et al. Enrichr: a comprehensive gene set enrichment analysis web server 2016 update. Nucleic Acids Res. 2016;44:W90–97.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  38. Butler A, Hoffman P, Smibert P, Papalexi E, Satija R. Integrating single-cell transcriptomic data across different conditions, technologies, and species. Nat Biotechnol. 2018;36:411–20.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  39. Korsunsky I, Millard N, Fan J, Slowikowski K, Zhang F, Wei K, Baglaenko Y, Brenner M, Loh PR, Raychaudhuri S. Fast, sensitive and accurate integration of single-cell data with Harmony. Nat Methods. 2019;16:1289–96.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  40. Amir E-AD, Davis KL, Tadmor MD, Simonds EF, Levine JH, Bendall SC, Shenfeld DK, Krishnaswamy S, Nolan GP, Pe’er D. viSNE enables visualization of high dimensional single-cell data and reveals phenotypic heterogeneity of leukemia. Nat Biotechnol. 2013;31:545–52.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  41. Becht E, McInnes L, Healy J, Dutertre C-A, Kwok IWH, Ng LG, Ginhoux F, Newell EW. Dimensionality reduction for visualizing single-cell data using UMAP. Nat Biotechnol. 2018;37:38–44.

    Article  Google Scholar 

  42. Aran D, Looney AP, Liu L, Wu E, Fong V, Hsu A, Chak S, Naikawadi RP, Wolters PJ, Abate AR, et al. Reference-based analysis of lung single-cell sequencing reveals a transitional profibrotic macrophage. Nat Immunol. 2019;20:163–72.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  43. Gao R, Bai S, Henderson YC, Lin Y, Schalck A, Yan Y, Kumar T, Hu M, Sei E, Davis A, et al. Delineating copy number and clonal substructure in human tumors from single-cell transcriptomes. Nat Biotechnol. 2021;39:599–608.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  44. Garcia-Alonso L, Holland CH, Ibrahim MM, Turei D, Saez-Rodriguez J. Benchmark and integration of resources for the estimation of human transcription factor activities. Genome Res. 2019;29:1363–75.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  45. Hanzelmann S, Castelo R, Guinney J. GSVA: gene set variation analysis for microarray and RNA-seq data. BMC Bioinformatics. 2013;14:7.

    Article  PubMed  PubMed Central  Google Scholar 

  46. Qiu X, Hill A, Packer J, Lin D, Ma YA, Trapnell C. Single-cell mRNA quantification and differential analysis with Census. Nat Methods. 2017;14:309–15.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  47. Jin S, Guerrero-Juarez CF, Zhang L, Chang I, Ramos R, Kuan CH, Myung P, Plikus MV, Nie Q. Inference and analysis of cell-cell communication using cell chat. Nat Commun. 2021;12:1088.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  48. Browaeys R, Saelens W, Saeys Y. NicheNet: modeling intercellular communication by linking ligands to target genes. Nat Methods. 2020;17:159–62.

    Article  CAS  PubMed  Google Scholar 

  49. Xun Z, Ding X, Zhang Y, Zhang B, Lai S, Zou D, Zheng J, Chen G, Su B, Han L, Ye Y. Reconstruction of the tumor spatial microenvironment along the malignant-boundary-nonmalignant axis. Nat Commun. 2023;14:933.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  50. Patel AP, Tirosh I, Trombetta JJ, Shalek AK, Gillespie SM, Wakimoto H, Cahill DP, Nahed BV, Curry WT, Martuza RL, et al. Single-cell RNA-seq highlights intratumoral heterogeneity in primary glioblastoma. Science (New York, NY). 2014;344:1396–401.

    Article  CAS  Google Scholar 

  51. Yao H, Xie W, Dai Y, Liu Y, Gu W, Li J, Wu L, Xie J, Rui W, Ren B, et al. TRIM65 determines the fate of a novel subtype of pituitary neuroendocrine tumors via ubiquitination and degradation of TPIT. Neuro Oncol. 2022;24:1286–97.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  52. Zhang S, Cui Y, Ma X, Yong J, Yan L, Yang M, Ren J, Tang F, Wen L, Qiao J. Single-cell transcriptomics identifies divergent developmental lineage trajectories during human pituitary development. Nat Commun. 2020;11:5275.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  53. Wang Z, Guo X, Gao L, Deng K, Lian W, Bao X, Feng M, Duan L, Zhu H, Xing B. The immune profile of pituitary adenomas and a novel immune classification for predicting immunotherapy responsiveness. J Clin Endocrinol Metab. 2020;105:e3207–3223.

    Article  PubMed  PubMed Central  Google Scholar 

  54. Sharma P, Goswami S, Raychaudhuri D, Siddiqui BA, Singh P, Nagarajan A, Liu J, Subudhi SK, Poon C, Gant KL, et al. Immune checkpoint therapy-current perspectives and future directions. Cell. 2023;186:1652–69.

    Article  CAS  PubMed  Google Scholar 

  55. Ilie MD, De Alcubierre D, Carretti AL, Jouanneau E, Raverot G. Therapeutic targeting of the pituitary tumor microenvironment. Pharmacol Ther. 2023;250: 108506.

    Article  CAS  PubMed  Google Scholar 

  56. Marques P, Korbonits M. Tumour microenvironment and pituitary tumour behaviour. J Endocrinol Invest. 2023;46:1047–63.

    Article  CAS  PubMed  Google Scholar 

  57. Zhang Q, Yao B, Long X, Chen Z, He M, Wu Y, Qiao N, Ma Z, Ye Z, Zhang Y, et al. Single-cell sequencing identifies differentiation-related markers for molecular classification and recurrence prediction of PitNET. Cell Rep Med. 2023;4: 100934.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  58. Cheng S, Li Z, Gao R, Xing B, Gao Y, Yang Y, Qin S, Zhang L, Ouyang H, Du P, et al. A pan-cancer single-cell transcriptional atlas of tumor infiltrating myeloid cells. Cell. 2021;184:792.

    Article  CAS  PubMed  Google Scholar 

  59. Ma R-Y, Black A, Qian B-Z. Macrophage diversity in cancer revisited in the era of single-cell omics. Trends Immunol. 2022;43:546–63.

    Article  CAS  PubMed  Google Scholar 

  60. Chiloiro S, Giampietro A, Gessi M, Lauretti L, Mattogno PP, Cerroni L, Carlino A, De Alessandris QG, Olivi A, Rindi G, et al. CD68+ and CD8+ immune cells are associated with the growth pattern of somatotroph tumors and response to first generation somatostatin analogs. J Neuroendocrinol. 2023;35: e13263.

    Article  CAS  PubMed  Google Scholar 

  61. Yeung JT, Vesely MD, Miyagishima DF. In silico analysis of the immunological landscape of pituitary adenomas. J Neurooncol. 2020;147:595–8.

    Article  PubMed  PubMed Central  Google Scholar 

  62. Matsuzaki H, Komohara Y, Yano H, Fujiwara Y, Kai K, Yamada R, Yoshii D, Uekawa K, Shinojima N, Mikami Y, Mukasa A. Macrophage colony-stimulating factor potentially induces recruitment and maturation of macrophages in recurrent pituitary neuroendocrine tumors. Microbiol Immunol. 2023;67:90–8.

    Article  CAS  PubMed  Google Scholar 

  63. Zhou W, Zhang C, Zhang D, Peng J, Ma S, Wang X, Guan X, Li P, Li D, Jia G, Jia W. Comprehensive analysis of the immunological landscape of pituitary adenomas: implications of immunotherapy for pituitary adenomas. J Neurooncol. 2020;149:473–87.

    Article  CAS  PubMed  Google Scholar 

  64. Shi Y, Zhang Q, Bi H, Lu M, Tan Y, Zou D, Ge L, Chen Z, Liu C, Ci W, Ma L. Decoding the multicellular ecosystem of vena caval tumor thrombus in clear cell renal cell carcinoma by single-cell RNA sequencing. Genome Biol. 2022;23:87.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  65. Lee M, Lee Y, Song J, Lee J, Chang S-Y. Tissue-specific role of CX3CR1 expressing immune cells and their relationships with human disease. Immune Netw. 2018;18:e5.

    Article  PubMed  PubMed Central  Google Scholar 

  66. Bassez A, Vos H, Van Dyck L, Floris G, Arijs I, Desmedt C, Boeckx B, Vanden Bempt M, Nevelsteen I, Lambein K, et al. A single-cell map of intratumoral changes during anti-PD1 treatment of patients with breast cancer. Nat Med. 2021;27:820–32.

    Article  CAS  PubMed  Google Scholar 

  67. Chen X, Jiang M, Li H, Wang Y, Shen H, Li X, Zhang Y, Wu J, Yu Z, Chen G. CX3CL1/CX3CR1 axis attenuates early brain injury via promoting the delivery of exosomal microRNA-124 from neuron to microglia after subarachnoid hemorrhage. J Neuroinflammation. 2020;17:209.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  68. Otobe S, Hisamoto T, Miyagaki T, Morimura S, Suga H, Sugaya M, Sato S. CX3CR1 deficiency attenuates DNFB-induced contact hypersensitivity through skewed polarization towards M2 phenotype in macrophages. Int J Mol Sci. 2020;21:7401.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  69. Imai T, Yasuda N. Therapeutic intervention of inflammatory/immune diseases by inhibition of the fractalkine (CX3CL1)-CX3CR1 pathway. Inflamm Regen. 2016;36:9.

    Article  PubMed  PubMed Central  Google Scholar 

  70. Lin S. Single-cell transcriptomics reveal distinct immune-infiltrating phenotypes and macrophage-tumor interaction axes among different lineages of pituitary neuroendocrine tumors, Genome Sequence Archive,

  71. Dai Y, Lin S, Wu Z. Source code for PitNETs single-cell project. GitHub; 2024. (2024)

  72. Dai Y, Lin S, Wu Z. Source code for PitNETs single-cell project. Figshare; 2024

Download references


Not applicable.


This work was supported by the National Natural Science Foundation of China under grant number 82141114, 81972339 (Z.B.W.), 82373131(S.L), 82172605 (L.X.), 82200153 (Y.D.), and 32100730 (S.L), Fundamental Research Funds for the Central Universities No.YG2023ZD06 (Z.B.W.); health care leader of Shanghai Municipal Health Commission No.2022LJ006 (Z.B.W.); National Research Center for Translational Medicine under grant number NRCTM (SH)2023–15 (Z.B.W.), Shanghai Municipal Science and Technology Commission under Grant Number 18XD1403400 (Z.B.W.), QiusuoFunds for Distinguished Young Scholars under grant number 20JC1419500 (L.X.); and Interdisciplinary Program of Shanghai Jiao Tong University project number YG2022QN008 (Y.D.).

Author information

Authors and Affiliations



Z.B.W., S.L., X.L., and Y.Y. conceived and supervised the project. S.L. and Y.D. designed and performed the research. S.L., C.H., B.Z., N.H., Y.L., Y.W., M.L., Y.B., and Y.C. collected the PitNETs samples. S.L., Y.D., T.H., C.H., L.Z., R.W., and J.L. performed the data analyses. S.L., C.H., N.H., and Y.L. performed experiments. S.L., Y.D., Y.Y., X.L.,C.H., N.H., Y.L., L.Z., R.W., J.L., H.Y., H.T., Y.Z., X.W.,S.L., and H.F. interpreted the results. X.L., S.L., Y.D., and Y.Y. wrote paper with input from all the other authors. All authors read and approved the final manuscript.

Corresponding authors

Correspondence to Lei Shen, Youqiong Ye, Li Xue or Zhe Bao Wu.

Ethics declarations

Ethics approval and consent to participate

All procedures performed with the use of samples obtained from patients were approved by Ruijin Hospital ethical committee approved the study. Our institution evaluated and approved animals experiments performed under the Guideline for Animal Care at Shanghai Jiao Tong University, School of Medicine (IACUC: B-2021–004). All the patients signed informed consent to publish. The research conformed to the principles of the Helsinki Declaration.

Consent for publication

Not applicable.

Competing interests

The authors declare no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Additional file 1.

Figure S1-S11. All supplementary figures.

Additional file 2. Table S1-5

. Clinical information of Bulk-RNA seq (Table S1), scRNA-seq (Table S2), TMA (Table S3), FC (Table S4) and clinical information for experiments (Table S5). Table S6. Quality control and filtration. Table S7. Spatial feature expression plots of ACVR1B coexpressing with INHBA and in adjacent spot. Table S8. List of antibodies and chemical dyes. Table S9. qPCR Oligonucleotide Sequences. Table S10. Cell type-specific markers in the UMAP representation of PitNETs. Table S11. Specific markers of macrophages.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit The Creative Commons Public Domain Dedication waiver ( applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Lin, S., Dai, Y., Han, C. et al. Single-cell transcriptomics reveal distinct immune-infiltrating phenotypes and macrophage–tumor interaction axes among different lineages of pituitary neuroendocrine tumors. Genome Med 16, 60 (2024).

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: