Proteogenomic analyses of PitNET specimens

To obtain the proteogenomic landscape of PitNET, whole-exome sequencing (WES), transcriptomics, proteomics, and phosphoproteomics datasets were collected from 200 fresh-frozen tumors and 7 APGs as controls, based on pathological criteria (see Materials and methods). Clinicopathological features, including TF lineage, clinicopathological subtypes, surgery invasion status, patient gender, tumor diameter, and KNOSP grade are summarized in Supplementary information, Table S1. Figure 1a illustrates the sample distribution across the three TF lineages and NULL, which were further divided into 10 clinicopathological subtypes: PIT1 lineage (n = 101, including 21 GH, 23 PRL, 15 TSH, 22 silent PIT1, and 20 plurihormonal), TPIT lineage (n = 46, including 21 ACTH, and 25 silent TPIT), SF1 (n = 31, including 12 GN and 19 silent SF1), and NULL (n = 22) (Supplementary information, Fig. S1b).

Fig. 1: Proteogenomic landscape of PitNETs.

a Top panel, pie charts of clinical indicators. Bottom panel, sample numbers and multi-omics datasets of the cohort. b Genomic profile and associated clinical features of patients with PitNETs. SMGs in this dataset identified by MutSigCV and OncodriveCLUST (q value < 0.1) are shown. Right panel, percentage of samples affected. Top panel, number of mutations per sample. Middle panel, distribution of significant mutations across sequenced samples, color coded by mutation type. Bottom panel, percentage of somatic base changes per sample. c Comparison of the TMB of our PitNETs cohort and 33 cancer types in TCGA studies. d Boxplot showing the VAF of the top 20 SMGs. e Bar plot showing the genes with significantly different mutation frequencies based on Fisher’s exact test by clinicopathological subtype (Fisher’s exact test, P value < 0.01). The numbers listed on the right side of the barplot represented the mutation frequencies in the indicated clinicopathological subtype tumors. The numbers listed on the left side of the barplot represented the mutation frequencies in the rest tumors. f, g Arm-level and focal-level amplifications and deletions. GISTIC analysis was performed to determine significant regions and genes included in the recurrent CNAs identified in patients with PitNETs. h PCA analysis of proteomics data from 200 PitNETs and 7 APGs based on clinicopathological subtypes.

WES was conducted on 200 tumor tissues and paired peripheral blood mononuclear cells (PBMCs) to identify possible genetic variants in the cancer genome. RNA sequencing (RNA-seq) was carried out for 194 tumors and 7 APGs. A mass spectrometry (MS)-based label-free quantitative (LFQ) method was used to characterize the proteomes of the 200 tumors and 7 APGs. A Fe-NTA-enrichment-based strategy was employed for phosphoproteomics profiling of 194 tumors and 7 APGs (Fig. 1a; Supplementary information, Fig. S1b).

Overview of the proteogenomic landscape of PitNET

WES data analysis revealed 7333 mutated genes, including 11,092 non-silent point mutations and 419 small insertions or deletions (indels) (Fig. 1b; Supplementary information, Table S1). In the 200 patients, we observed several significantly mutated genes (SMGs, q < 0.1) associated with PitNET functions, including GNAS (11%), KRT76 (8%), TCHH (6%), ZMIZ2 (6%), DGKZ (6%), KRTAP9-2 (5%), and USP8 (4%) (Fig. 1b). Examination of the proportion of somatic base changes revealed that PitNET patients carried a high proportion of C > A transitions compared to the other five substitution types (Fig. 1b). A comparison with previous studies using The Cancer Genome Atlas (TCGA)23 indicated that, in this study, the tumor mutation burden (TMB) in the PitNET cohort remained a lower-middle level among the 33 cancer types (Fig. 1c).

We next compared the variant allele frequencies (VAFs) of the SMGs and found that USP8 mutation (median: 0.45) ranked first, followed by GNAS mutation (median: 0.36) (Fig. 1d). All nucleotide variants in these two genes detected in PitNET patients were previously reported.15,19 Enrichment analysis using Fisher’s exact test to identify mutations associated with clinicopathological subtype showed that USP8 mutations were enriched in the ACTH subtype (P = 3.12e‒8), while GNAS mutations were enriched in the plurihormonal subtype (plurihormonal subtype vs remaining samples, P = 0.002) and GH subtype (GH subtype vs remaining samples, P = 0.00059) (Fig. 1e). Somatic CNA analysis identified arm-level amplifications (Chr 5, 7, 8, 9, 12q, 14q, 19, 20, 21p) and deletions (Chr 1p, 2q, 11) (Fig. 1f). Focal peaks included amplifications of 5p15.33, 7p14.1, 12p13.2, 13q12.12, 14q11.2 and 16p13.13 and deletions of 1q21.3, 11q11, 17q24.3 and 22q11.23, among others (see Materials and methods; Fig. 1g; Supplementary information, Table S2).

Our transcriptomics, proteomics, and phosphoproteomics datasets exhibited a unimodal distribution and passed the quality control (QC) procedure (Supplementary information, Fig. S1c). RNA-seq identified 18,397 genes with fragments per kilobase of transcript per million fragments mapped (FPKM) values over 1, providing an opportunity to explore the relationship between transcriptome and proteome. For proteomics analysis, whole-cell extracts of human embryonic kidney-derived HEK293T cells were used as controls for quality. Quantitative MS analysis of HEK293T cells confirmed the robustness and consistency of the MS data, indicated by a high Spearman’s correlation coefficient of 0.91 among the proteomes of QC samples (Supplementary information, Fig. S1d). Moreover, the dataset used in this study provided in-depth coverage of the human proteome. A total of 10,011 proteins (with ≥ 2 unique peptides per protein) were identified in the 200 tumors and 7 APGs, while a total of 29,219 phosphosites were detected, corresponding to 5483 phosphoproteins. Among them, 6160 proteins and 9905 phosphosites from 3276 phosphoproteins were selected for downstream analysis based on their presence in more than 50% of cases of at least one clinicopathological tumor subtype.

Principal component analysis (PCA) of our multi-omics data revealed a significant separation between the PIT1 lineage (GH, PRL, TSH, silent PIT1, and plurihormonal subtypes) and the APG, ACTH, silent TPIT, and other PitNETs at the proteomics level, whereas the PIT1 lineage showed higher similarity to the APG group at the transcriptomics and phosphoproteomics levels (Fig. 1h; Supplementary information, Fig. S1e, f). Transcriptomics and proteomics data further indicated that lineage-specific TFs (PIT1, TPIT, and SF1) and hormone-related genes (GH1, PRL, TSHB, POMC, LHB, and FSHB) were expressed in specific clinicopathological subtypes (Supplementary information, Fig. S1g). These cumulative results thus provide a multi-omics landscape to improve our understanding of the molecular mechanisms of PitNETs.

Impact of genomic alterations on the transcriptome, proteome, and phosphoproteome

Correlation analysis of the paired transcriptomics and proteomics datasets showed that 92.95% of 6115 mRNA-protein pairs were positively correlated in tumor samples. Genes with strong correlations were enriched in several pathways related to neuronal system, epithelial-mesenchymal transition (EMT), and hormone metabolic process, which indicate that these pathways are overrepresented in PitNETs (Fig. 2a). In addition, the global mRNA-protein correlation was moderate with sample-wise median spearman correlations of 0.45 and the correlation of each clinicopathological subtypes ranged from 0.42 to 0.46 (Supplementary information, Fig. S2a, b), which were consistent with previous reports.21,24

Fig. 2: Impact of CNAs on the transcriptome, proteome and phosphoproteome of PitNETs.

a Gene-wise mRNA-protein Spearman’s correlation in tumors. Red, pathways involving positively correlated genes; blue, pathways involving negatively correlated genes (Spearman’s correlation, FDR < 0.05). b The correlation of CNAs to mRNA (left) or protein abundance (right), with significant positive correlations in red and negative correlations in green (Spearman’s correlation, FDR < 0.05). Genes were sorted by chromosomal location on the x- and y-axes. c Cascading effects of CNAs and the overlap between cis events via the transcriptome and proteome analyses (Spearman’s correlation, FDR < 0.05). d Prioritized cis effect CNA drivers were used for pathway enrichment analysis in ConsensusPathDB. e Venn diagram showing the CAGs with significant CNA cis effects via multi-omics data analyses (Spearman’s correlation, FDR < 0.1). f Genes with cascading copy number cis regulation of their cognate mRNA, protein, and phosphoprotein levels. Shapes indicate the cis effects across the indicated datasets.

We examined the regulatory effects of 23,109 somatic CNAs on mRNA, protein, and phosphoprotein abundances of genes at the same loci (cis effects) and genes at other loci in the genome (trans effects) (see Materials and methods; Fig. 2b; Supplementary information, Fig. S2c and Table S2). We observed cis effects for 6545 and 684 CNAs affecting mRNAs and proteins, respectively. Among them, 502 significant cis effect events overlapped (Spearman’s correlation, FDR < 0.05) (Fig. 2c; Supplementary information, Table S2); these 502 genes were enriched in pathways related to post-translational protein modification, NAD+ biosynthesis, and metabolism (Fig. 2d). We then assessed how CNA events influenced mRNA, protein and phosphoprotein abundances of cancer-associated genes (CAGs) via either cis or trans effects, focusing on alterations in 593 previously described genes (Supplementary information, Table S2).25 We found that CNAs have cis effects on both mRNA and protein abundances of 15 CAGs, while 6 CAGs showed significant overlapped CNA cis effects (FDR < 0.1) at the mRNA and phosphoprotein levels (Fig. 2e). Figure 2f shows the annotations of these 21 CAGs. The cis or trans effects of these 21 genetic alterations were also comprehensively investigated (Supplementary information, Fig. S2d, e). In particular, we observed that GNAS copy number had cis effects on GNAS, and trans effects on EEF2, ELL, and RAB8A at the mRNA and protein levels (Supplementary information, Fig. S2d).

Impact of GNAS mutation and GNAS copy number gain in the PIT1 lineage

In our cohort, GNAS, enriched in the PIT1 lineage, was the most frequently mutated gene, harboring two mutation hotspots, R186C/G/L/H and Q212L (Figs. 1b, 3a). GNAS copy number gain, as a CAG with cis effect (Fig. 2f; Supplementary information, Fig. S2d), had particularly strong impacts on the PIT1 lineage (Spearman’s correlation in 101 PIT1 lineage vs all 200 PitNETs: R = 0.38, P = 0.0001, vs R = 0.19, P = 0.0071 at the mRNA level; R = 0.41, P = 2.29e‒05 vs R = 0.21, P = 0.0024 at the protein level) (Fig. 3b). Samples with either GNAS mutations (VAF > 5%) or GNAS copy number gain were significantly enriched in the PIT1 lineage (Fig. 3c). Thus, by integrating WES data, we could further divide the PIT1 lineage into three subgroups, including wild-type (WT), GNAS copy number gain, and GNAS mutant. Compared with the WT group, samples carrying GNAS mutation showed reduced GNAS protein levels (Wilcoxon rank-sum test, P = 7.6e‒7), while those with GNAS copy number gain showed increased GNAS at both mRNA and protein levels (Wilcoxon rank-sum test, P = 0.014 and P = 0.021) (Fig. 3d; Supplementary information, Fig. S3a, b).

Fig. 3: Impact of GNAS mutation and GNAS copy number gain in the PIT1 lineage.

a Lollipop plot and boxplot showing the position and tumor VAF of the GNAS mutation in the PIT1 lineage. b Spearman’s correlation of chromosome 20q and the copy number, mRNA expression and protein abundance of GNAS in all PitNET samples and PIT1 lineage samples. Spearman’s correlation, *P < 0.05, **P < 0.01, ***P < 0.001. c Distribution of GNAS altered samples in different categories among the PIT1 lineage and other lineages (Fisher’s exact test, **P < 0.01, ***P < 0.001). d Heatmap visualizing multi-omics profiles of the levels of GNAS copy number, mRNA expression and protein abundance. e Volcano plots displaying the differentially expressed proteins in GNAS mutant and GNAS WT patients after applying a two-fold change in expression with P < 0.05 (Wilcoxon rank-sum test). Proteins significantly enriched in the GNAS mutant and GNAS WT patients are represented as red/blue-filled dots. f Pathways enriched for the differentially expressed mRNAs and proteins. Pathways that were significantly upregulated/downregulated in the GNAS mutants are represented as red/blue-filled dots. g Heatmap of multi-omics features of GH secretion-related genes. The pathway diagram on the right depicts how the features included in the heatmap regulate GH synthesis, secretion and activity. Red boxes indicate upregulated genes and blue boxes indicate downregulated genes. Green rectangles indicate kinases and orange circles indicate phosphorylated proteins. Bar chart next to the heatmap shows the fold changes of GNAS mutant/WT (*P < 0.05, **P < 0.01, ***P < 0.001). h GSEA plots for proliferation-related pathways based on the rank of GNAS copy number-mRNA (bottom) or protein (upper) abundance correlations. i Boxplots showing the difference of MGPS and tumor volume between WT and GNAS copy number gain group. The significance was calculated by Wilcoxon test. j Heatmap of multi-omics features of proliferation-related genes. The pathway diagram on the left depicts how the features included in the heatmap regulate cell cycle S-phase and DNA biosynthesis. Red boxes indicate upregulated genes and blue boxes indicate downregulated genes. Green rectangles indicate kinases and orange circles indicate phosphorylated protein. Bar chart next to the heatmap shows the Spearman’s correlation coefficient between GNAS copy number and proliferation-related genes (*P < 0.05, **P < 0.01, ***P < 0.001). k, l Bar plots showing the proportion of CDK6 high H-score cells between GNAS high H-score group and GNAS low H-score group in all PitNETs and PIT1 lineage PitNETs. The significance was calculated by Fisher’s exact test.

GNAS mutations in PitNET patients have been linked to a gain of function in G protein-coupled receptor (GPCR) signaling pathways,26 although the specific downstream impacts remain unknown. Compared with the WT group, we found that genes involved in growth hormone (GH) synthesis, secretion, and action pathways (e.g., GH1 and GH2) were upregulated at the mRNA, protein, and phosphorylation levels (q < 0.05) in the GNAS mutant group (Fig. 3e, f; Supplementary information, Fig. S3c, d and Table S3). More specifically, we identified the upregulation and phosphorylation (ADCY5_S666 and ADCY6_S576) of adenylate cyclase (AC), along with the expression of PKAs in the GNAS mutant PIT1 lineage tumors. We further focused on the components in the GNAS-PKA downstream pathways. Combining with the known mechanism, we speculated that PKAs might promote the phosphorylation of CREBBP and subsequent accumulation of the CREB complex (CREB1 and ATF2), ultimately leading to hypersecretion of GH through PIT1 activation. Likewise, we also infer that the hypersecretion of GH in these samples might affect the levels of SHC (i.e., SHC1 and SHC3) and STAT5B, as well as AKT1 and PKC protein activities (i.e., PRKCA and PRKCB), which are known to promote cell growth and metabolism.

Amplification of 20q has been reported in PitNET,27 while the cis and trans effects of 20q amplification and GNAS CNA (located at 20q) remain unclear. Gene set enrichment analysis (GSEA) of transcriptomics/proteomics datasets by Spearman’s correlation showed upregulation in proliferation-related pathways, such as cell cycle and DNA replication pathways, in patients with GNAS copy number gain (FDR < 0.05) (Fig. 3h; Supplementary information, Fig. S3e, f and Table S3). Furthermore, the GNAS copy number gain group had higher multigene proliferation score (MGPS) and clinical tumor volume as compared with the WT group (Wilcoxon test, P < 0.05) (Fig. 3i). However, the correlation was non-significant in tests with GH PitNETs alone.28 To determine the proliferation characteristics of the PIT1 lineage driven by GNAS copy number gain, we systematically characterized the signal cascade related to cell cycle and DNA synthesis. Among cell cycle-related molecules, PRKDC and CDK6 were the top two proteins positively correlated with GNAS copy number (Fig. 3j; Supplementary information, Fig. S3g). Chemical inhibition or knockdown of GNAS has been shown to decrease the expression of cyclin proteins such as cyclin D, which is closely related to CDK6.29,30 This combined evidence suggested that CDK6 could contribute to the enhanced proliferation rate of PIT1 PitNETs as a result of GNAS copy number gain. It was also noteworthy that Rb mRNA level was positively correlated with GNAS copy number, mRNA and protein levels (Spearman’s correlation: GNAS copy number, R = 0.26, *P < 0.05; GNAS mRNA, R = 0.30, **P < 0.01; GNAS protein, R = 0.26, **P < 0.01) (Supplementary information, Fig. S3h). In addition, Rb phosphorylation levels at the RB1_S37 site were significantly correlated with GNAS at the mRNA and protein levels (Spearman’s correlation: GNAS mRNA, R = 0.23, *P < 0.05; GNAS protein, R = 0.33, ***P < 0.001). Finally, our data showed that E2F and CDK2 were downregulated, which might lead to the upregulation of ORC family members. The upregulation of the ORC family, RFC family, and PCNA in patients with GNAS copy number gain likely led to the elevated DNA biosynthesis and the enhanced tumor cell proliferation (Fig. 3j).

To further confirm the impacts of GNAS copy number gain, we performed immunohistochemistry (IHC) for GNAS and CDK6 and calculated IHC staining scores (H-scores). H-scores of GNAS and CDK6 were divided into high and low H-score groups based on the median score, respectively. As expected, the proportion of CDK6 high H-score cells was greater in the GNAS high H-score group than in the GNAS low H-score group in both the PIT1 lineage PitNETs and all PitNETs (Fisher’s exact test: PIT1 lineage PitNETs, P = 0.026; all PitNETs, P = 0.031) (Fig. 3k, l; Supplementary information, Fig. S3i).

In conclusion, these findings illustrate the diverse impacts of genomic events in the GNAS gene, such as mutations that drive hormone hypersecretion. Moreover, the finding that GNAS copy number gain can markedly enhance tumor cell proliferation implied that an inhibitor therapy targeting CDK6 may be effective for PIT1 lineage patients harboring GNAS copy number gain.

Multi-omics classification of PitNETs

To comprehensively explore the phenotypic and genotypic PitNET diversity in this cohort, classification by consensus clustering31 was performed with the combined transcriptomics, proteomics, and phosphoproteomics data. This analysis identified seven proteomic (Supplementary information, Fig. S4a), five transcriptomic (Supplementary information, Fig. S4b), and seven phosphoproteomic (Supplementary information, Fig. S4c) clusters among the PitNETs (Supplementary information, Table S4), which were subsequently named according to their similarities to clinicopathological subtypes and predominant pathway associations.

At the protein level, the seven proteomic clusters included GHenrich, EMTPRO, PRLenrich, TSH/silent PIT1enrich, ACTHenrich, silent TPITenrich, and SF1/NULLenrich (Fig. 4a). Pathway enrichment analysis (see Materials and methods; Supplementary information, Table S4) showed that the Hedgehog signaling pathway was differentially upregulated in GHenrich, and MYC targets v1 was upregulated in PRLenrich (Fig. 4a). TSH and silent PIT1 were co-clustered and enriched for pathways such as interferon-ɑ response, and antigen processing and presentation. In addition, the SF1 lineage and NULL PitNETs clustered together, forming the SF1/NULLenrich cluster, which showed upregulation in metabolism-related pathways, including fatty acid metabolism and the citrate cycle. Moreover, males were more prevalent (78.3%), average age was higher (> 60, 39.1%), and tumor diameter was larger (≥ 40 mm, 23.9%) in the SF1/NULLenrich cluster compared to other clusters (Fig. 4b). Notably, the TPIT lineage was divided into two smaller clusters, ACTHenrich and silent TPITenrich, at the protein level, which was consistent with clinicopathological subtypes (Supplementary information, Fig. S4d, e). The ACTHenrich cluster was enriched for USP8 mutations and both ACTHenrich and silent TPITenrich had an extremely high proportion of females (90.5% and 90.9%, respectively) (Fig. 4b). In addition, ACTHenrich, silent TPITenrich, and SF1/NULLenrich clusters were all associated with higher MGPS (Kruskal-Wallis test, P = 7.4e‒06) (see Materials and methods; Supplementary information, Fig. S4f), which was aligned well with the upregulation of proliferation and energy metabolism pathways in these three clusters.

Fig. 4: Molecular subtypes of PitNETs based on proteogenomic analysis and association studies.

a Heatmap illustrating the characterization of seven proteomic clusters. Each column represents a patient sample and rows indicate proteins. The color of each cell shows the z score of the protein in that sample. PitNET classification, hormone secretion status, invasion status, clinical features, and mutation status annotations are shown above the heatmap. The chi-square test was used to evaluate the association of proteomic clusters with the 9 variables on the heatmap (*P < 0.05, **P < 0.01, ***P < 0.001). Single-sample Gene Set Enrichment Analysis (ssGSEA) based on proteomics data was also applied to identify the dominant pathway signatures in each proteomic cluster. b Summary of the variables with significant differences among the seven proteomic clusters. The percentage represents the proportion of the population. c Sankey diagram depicting the result of integrative multi-omics analysis, showing the flow of cluster assignments across multiple classification of PitNETs. d Boxplots depicting the distribution of stromal scores inferred by ESTIMATE based on the RNA data (left) and protein data (right) among tumors of the seven proteomic clusters. Kruskal-Wallis test was used to test if any of the differences among the subgroups were statistically significant. The Wilcoxon rank-sum test was used to estimate the difference between two subgroups, *P < 0.05, **P < 0.01, ***P < 0.001. e Boxplot depicting the distribution of stromal scores based on stromal feature percentage_ML among tumors of the seven proteomic clusters. Kruskal-Wallis test was used to test whether any of the differences among the subgroups were statistically significant. Wilcoxon rank-sum test was used to estimate the difference between two subgroups, *P < 0.05, **P < 0.01, ***P < 0.001. f Representative IF staining of pan-cytokeratin (panCK) and fibronectin1 (FN1) in EMTPRO and non-EMTPRO clusters. Scale bar, 50 μm. g Volcano plot showing differential mRNA expression of TFs between EMTRNA and Hormone clusters (the horizontal axis is log2(fold change), and the vertical axis is –log10 FDR). The upregulated TFs in EMTRNA are highlighted in red and EMT-TFs are highlighted in green. h Correlation heatmaps showing the correlation among the mRNA expression of five EMT-TFs in EMTRNA and EMTPRO clusters. Spearman’s correlation, **P < 0.01, ***P < 0.001. i Taking POU1F1 as the positive control, heatmap showing the molecules significantly differentially expressed between EMTRNA and Hormone clusters at the mRNA, protein, and TF activity levels, including EMT-TFs and EMT-related markers. j IHC staining validated the correlation between EMT-TFs and tumor invasion. Scatterplots showing the correlation of H-scores of TWIST1 and ZEB2 with tumor volume (Spearman’s correlation). The boxplots show the association of H-scores of TWIST1 and ZEB2 with surgery invasion status (Wilcoxon rank-sum test). k Summary of the multi-omics classification of the PIT1 lineage.

To further characterize the proteogenomic classification of PitNETs, we performed integrative analysis of the ten clinicopathological subtypes, three TF lineages and NULL, five transcriptomic clusters, seven proteomic clusters, and seven phosphoproteomic clusters for PitNETs. Interestingly, the ACTHenrich, silent TPITenrich, and SF1/NULLenrich clusters identified using proteomics data were highly consistent with clusters identified using transcriptomics and phosphoproteomics data (Fig. 4c). Furthermore, a cluster of PitNETs was also identified with clear EMT characteristics at the transcriptomics, proteomics and phosphoproteomics levels.

An invasive cluster characterized by EMT was identified within the PIT1 lineage

At the protein level, we found that hemostasis-related and EMT-related molecules,32,33,34 including GP1BB, FGB, MMP8, FN1, and ITGB3, were highly expressed in EMTPRO, compared with other proteomic clusters (Supplementary information, Table S4). The pathways related to EMT, TNFA signaling via NFκB, and cytokine‒cytokine receptor interaction were all upregulated in the EMTPRO cluster (Fig. 4a), which covered eight of the ten clinicopathological subtypes, excluding ACTH and silent SF1 (Supplementary information, Fig. S4d). Strikingly, EMTPRO showed strong invasiveness (Fig. 4b), with a high level of KNOSP (grade = 4, 37.9%), surgery invasion (69%), and tumor diameter (≥ 40 mm, 17.2%).

Given the non-negligible role of EMT in cancer metastasis,35,36 we next used the ESTIMATE algorithm37 to deconvolute the contribution of stromal cells in the tumors based on transcriptomics data (Stromal score_RNA) and proteomics data (Stromal score_Protein). EMTPRO showed overall higher stromal scores in both Stromal score_RNA (Kruskal-Wallis test, P = 9.9e‒06) and Stromal score_Protein (Kruskal-Wallis test, P = 3.6e‒04) (Fig. 4d). Hematoxylin and eosin (HE) staining was also processed to evaluate the proportion of tumor cells that featured stromal morphology with quantification by QuPath bioimage analysis (Stromal feature percentage_ML, see Materials and methods) and confirmed that the EMTPRO cluster had the highest proportion of cells with a stromal phenotype (Kruskal-Wallis test, P = 0.04) among all the proteomic clusters (Fig. 4e; Supplementary information, Fig. S4g). To further investigate the EMT status of tumor cells in each of the seven proteomic clusters, immunofluorescence (IF) co-staining was performed to detect the epithelial marker, pan-cytokeratin (panCK), and mesenchymal marker, fibronectin1 (FN1), in a large subset of tumors. The IF results showed significantly higher percentage of areas with co-staining of panCK and FN1 in EMTPRO cluster than in other clusters (Fig. 4f; Supplementary information, Fig. S4h). All these pieces of evidence supported that EMTPRO cluster was characterized by tumor cells with EMT status.

EMT-inducing transcription factors (EMT-TFs)38 are those confirmed as key drivers of the EMT-phenotype. We compared TF mRNA expression level between the two transcriptomic clusters in the PIT1 lineage (Fig. 4g), which confirmed that five EMT-TFs were significantly upregulated in EMTRNA including: SNAI1 (FC = 1.97, FDR = 0.000074), SNAI2 (FC = 6.63, FDR = 2.3e‒06), ZEB2 (FC = 3.04, FDR = 2.1e‒06), TWIST1 (FC = 5.3, FDR = 5.93e‒05), and TWIST2 (FC = 17, FDR = 2.02e‒05). The remaining PIT1 lineage cluster was designated Hormone due to upregulation of hormone secretion proteins. Furthermore, the mRNA expression patterns of the five EMT-TFs were significantly positively correlated in both the EMTRNA (Spearman’s correlation, R = 0.55‒0.74, P < 0.01) and EMTPRO (Spearman’s correlation, R = 0.70‒0.76, P < 0.01) clusters (Fig. 4h). Notably, the levels of transcriptional activity of the five EMT-TFs were higher in the EMTRNA cluster compared with those in the Hormone cluster (Fig. 4i). Similarly, EMT-related molecules35,39 including CDH2, VIM, CD44, FN1, and ITGB1 (mesenchymal markers) were upregulated in the EMTRNA cluster, while CDH1, KRT8, and KRT18 (epithelial markers) were downregulated in the transcriptomics and proteomics datasets (Fig. 4i) further confirming the EMT status of this cluster.

In addition, IHC staining of TWIST1 and ZEB2 in our cohort verified that EMT-TFs were activated in EMTRNA (Supplementary information, Fig. S4i). There were significant positive correlations between the H-scores of TWIST1 and ZEB2 and their corresponding tumor volumes (Spearman’s correlation: TWIST1, R = 0.33, P = 0.037; ZEB2, R = 0.34, P = 0.034), and the high H-scores of TWIST1 and ZEB2 were associated with surgery invasion (Wilcoxon rank-sum test, TWIST1: P = 0.0083; ZEB2: P = 0.025) in the PIT1 lineage (Fig. 4j).

In summary, integrated proteogenomic characterization of PitNETs identified a previously unrecognized, highly invasive cluster defined by EMT in PitNETs, primarily containing PIT1 lineage tumors (Fig. 4k).

Proteogenomics data revealed three modes of EGFR activation in the TPIT lineage

EGFR is associated with a variety of human cancers, including head and neck squamous cell carcinoma and lung adenocarcinoma,40,41,42 and has been proposed as a therapeutic target in ACTH PitNETs.11,43 Here, we found that the levels of EGFR mRNA expression, protein abundance, and phosphorylation modifications were higher in the TPIT lineage than in other tumors (Supplementary information, Fig. S5a‒c), which expanded the previous perception that EGFR was highly expressed in ACTH PitNETs.11 Subsequent analysis of EGFR-related pathways defined three groups that showed diverse mechanisms of EGFR activation, including ACTH tumors with USP8 mutation (ACTH_USP8 mutant), ACTH tumors without USP8 mutation (ACTH_ USP8 WT), and silent TPIT tumors (Supplementary information, Fig. S5b, Table S5). In the ACTH_USP8 mutant group, our data supported the known mechanism that USP8 gain-of-function mutations rescue EGFR from ubiquitination, leading to the enhanced EGFR activity, and further promoting POMC biosynthesis (Supplementary information, Fig. S5c, d).15,16 In the ACTH_ USP8 WT group, the average mRNA expression of EGFR ligands (i.e., AREG, TGFA, EGF, BTC, EPGN, HBEGF, and NRG4) was significantly higher than in other groups (Kruskal-Wallis test, P = 0.009) (Supplementary information, Fig. S5e) and positively correlated with peptide hormone biosynthesis (Spearman’s correlation, R = 0.4, P = 0.0079) and serum ACTH level (Spearman’s correlation, R = 0.34, P = 0.026) (Supplementary information, Fig. S5f, g). In silent TPIT tumors, EGFR T693 phosphorylation showed a significant enrichment (Kruskal-Wallis test, P = 0.0087), and EGFR downstream pathways or components, including PI3K-AKT-mTOR (Spearman’s correlation, R = 0.38, P = 0.0093), MAPK (Spearman’s correlation, R = 0.33, P = 0.023) and the cell cycle pathways were also enriched (Spearman’s correlation, R = 0.19, P = 0.2) (Supplementary information, Fig. S5h‒j), suggesting that EGFR T693 phosphorylation may lead to activation of these pathways. Furthermore, IHC staining of EGFR revealed that its expression was higher in the TPIT lineage than in non-TPIT lineages, and a higher positive staining rate of EGFR T693 phosphorylation was found in silent TPIT tumors as compared with ACTH tumors (Supplementary information, Fig. S5k). Based on these results, we summarized the potential therapeutic options for each of the three modes (Supplementary information, Fig. S5l).

In addition to the finding of the effects of EGFR on POMC biosynthesis, we also explo