Proteogenomic landscape and clinical characterization of GH-producing pituitary adenomas/somatotroph pituitary neuroendocrine tumors

Mutational landscape assessment by targeted capture sequencing in GHomas/somatotroph PitNETs

A summary of the clinical data for all 121 GHomas/somatotroph PitNETs patients is shown in Table 1. The mean basal GH level was high (15.7 ng/mL, IQR; 7.2–34.1), and the average IGF-1 SD score was 7.0 ± 2.4 (SD). The mean tumor volume was 1633.1 mm³ (IQR; 651.9–4412.4), and 23.1{6d6906d986cb38e604952ede6d65f3d49470e23f1a526661621333fa74363c48} of patients presented with a high Knosp grade (≥3). All 121 patients underwent transsphenoidal surgery (TSS), and 67 patients (55.4{6d6906d986cb38e604952ede6d65f3d49470e23f1a526661621333fa74363c48}) received preoperative SSA treatment. We evaluated biological remission by postoperative oral glucose tolerance test (OGTT) or normalization of IGF-1 SDS, resulting in 21 patients requiring additional therapy. The clinical data for all patients are provided in Supplementary Data 1. We performed targeted capture sequencing (TCS) of 36 genes (Supplementary Data 2) in all tumor tissue samples. The average sequencing depth was 1838.3 (range 29×–11973×), and 99.7{6d6906d986cb38e604952ede6d65f3d49470e23f1a526661621333fa74363c48} of the target regions were covered at least 50×. We detected a total of 83 mutations in 30 genes, with a median of 3.0 somatic coding region variants per tumor (range, 0–9). Among these mutations, 37 mutations in 18 genes were defined as recurrent (Fig. 1a).

Table 1 Clinical characteristics of GH-producing pituitary adenoma/ somatotroph pituitary neuroendocrine tumor patients and comparison of those with wild-type versus mutant GNAS.

**Fig. 1: The landscape of gene mutations and CNAs in GHomas/somatotroph PitNETs.**

We observed various activating GNAS-MT (57.0{6d6906d986cb38e604952ede6d65f3d49470e23f1a526661621333fa74363c48} of patients) at known hot spots in 69 patients in our cohort; p.Arg201Cys was identified in 45 patients, p.Arg201His in 5 patients, p.Arg201Ser in 3 patients, p.Gln227Leu in 13 patients, p.Gln227Arg in 2 patients, and p.Gln227Glu in 1 patient. One patient carried both p.Arg201Cys and p.Gln870Leu mutations. We also detected 3 unreported GNAS-MT (p.Gly49Arg, p.Ser111Asn, and p.Ala249Asp). The number of patients was one. We mapped the mutations to the crystal structure of the Gsα protein (PDB: 6AU6) and found that p.Gly49 is located in proximity to the GTP binding domain (Supplementary Fig. 1a, b).

We detected several gene mutations that have previously been reported in GHomas/somatotroph PitNETs [AIP, GPR101, somatostatin receptor 5 (SSTR5), arrestin beta 1/2(ARRB1/2)] and associated with the cAMP pathway [protein kinase cAMP-activated catalytic subunit alpha (PRKACA), protein kinase cAMP-activated catalytic subunit beta (PRKACB), and protein kinase cAMP-dependent type I regulatory subunit alpha (PRKAR1A)]. AIP mutations were observed in 8 patients, and a GPR101 mutation was detected in 1 patient. SSTR5 was mutated in 3 patients. We detected 2 patients with ARRB1 mutations and 2 patients carrying ARRB2 mutations. PRKACA mutation was detected in 1 patient, and PRKAR1A and PRKACB mutations were observed in 2 patients each. Most of these mutations were not recurrent in our study, and no hotspots have previously been reported for these genes.

Copy number analysis revealed the loss of SDHx

We calculated copy numbers from the TCS data using CNVkit, as described in the Methods, and detected 21 cases with CNAs (Fig. 1a). Fourteen patients showed copy number gains, and 10 patients presented copy number losses. Seven patients harbored CNAs in more than 2 genes (range, 1–5). Notably, among the 14 patients with identified copy number gains, 7 patients harbored gains in SSTR5, 5 patients harbored gains in GPR101, and one patient presented CNAs in both genes. Among patients with identified copy number losses, 6 patients harbored losses in PRKACB, and 4 patients presented losses in succinate dehydrogenase complex (SDHx) genes, including SDHB, SDHC, and SDHD.

Clinical characteristics of patients with GNAS mutations

Although we identified 37 recurrent mutations in this study, we focused specifically on GNAS-MT. GNAS mutations (GNAS-MT) are well-known driver mutations of GHomas/somatotroph PitNETs, and more than half of the patients in our study harbored these mutations. A comparison between patients with GNAS-MT and those without mutations (GNAS-WT) is shown in Table 1 and Fig. 1b, c. The GNAS-MT group presented significantly better responsiveness to bromocriptine, consistent with a previous report that GNAS-MT adenomas are associated with higher dopamine receptor 2 mRNA expression than GNAS-WT adenomas²³. A bubble plot analysis indicated that the GNAS-MT group also presented with smaller tumors and lower Knosp grades, although each dataset was associated with an extremely wide range. The GNAS-MT group tended to present with higher basal plasma GH levels and better responsiveness to octreotide loading test than the GNAS-WT group. The GNAS-MT group also showed a significantly better GH change rate following preoperative SSA treatment than the GNAS-WT group. However, the percentage of patients who required postoperative therapy was not significantly different between the two groups.

Consensus clustering-based transomics classification of pituitary adenomas/PitNETs

We performed proteomics and RNA sequencing analyses in pituitary adenomas/PitNETs (45 non-functioning pituitary adenomas/non-functioning pituitary neuroendocrine tumors [NFPAs/non-functioning PitNETs] for comparison, 60 GHomas/somatotroph PitNETs among 121 cases of the TCS cohort). The clinical data for all NFPAs/non-functioning PitNETs patients are provided in Supplementary Data 3. First, to clarify the transomics analysis, we examined the correlation between RNA and protein expression from RNA sequencing and proteomics data. Figure 2a shows a scatter plot of the correlation coefficient of a significantly positive correlation gene, with a mean Spearman’s correlation coefficient of 0.476 (Fig. 2a). This is similar to the results of the correlation analysis of RNA sequencing and proteomics in human rectal colon cancer²⁴ and may represent tumor characteristics in pituitary tumors. Next, the whole gene was analyzed. Although 89.1{6d6906d986cb38e604952ede6d65f3d49470e23f1a526661621333fa74363c48} of the genes showed a positive mRNA-protein correlation, only 65.2{6d6906d986cb38e604952ede6d65f3d49470e23f1a526661621333fa74363c48} were significantly correlated (Fig. 2b). The average Spearman’s correlation between mRNA and protein variations was 0.330. There were uncorrelated genes and negatively correlated genes. These results suggest that there are networks that could only be identified by transomics. To test whether the concordance between protein and mRNA variation was related to the biological function of the gene product, we performed Kyoto Encyclopedia of Genes and Genomes enrichment analysis. Genes involved in several hypothalamic-pituitary-adrenal axis signaling pathways showed concordant mRNA and protein variations (Fig. 2c). These results suggest that the transomics data also accurately represented pituitary tumor characteristics. Consensus clustering subtyping was applied to RNA sequencing datasets, proteomics datasets, and the combination of these two datasets (defined as the transomics data set), exploring between 2 and 8 K-means clusters. Consensus cumulative distribution functions (CDF) and delta area (change in the CDF area) plots were generated for each dataset to determine the optimal K value (Supplementary Fig. 2). The Sankey diagram depicts the flow of cluster assignments for NFPAs/non-functioning PitNETs and GHomas/somatotroph PitNETs across data types. To clarify the significance of performing transomics analyses, we examined the classification by consensus clustering using both RNA sequencing and proteomics data. The appropriate cluster sizes were determined to be K = 5 for RNA sequencing data, K = 4 for proteomics data, and K = 4 for the transomics combination (Fig. 2d). The results showed that transomics analysis better reflects protein expression classification compared to RNA analysis.

RNA sequencing classification may differ from proteomics classification because the expression of RNA and protein molecules is not always correlated (Fig. 2a, b). However, similar expression patterns were identified for pituitary-specific transcription factors in the classification of NFPAs/non-functioning PitNETs and GHomas/somatotroph PitNETs (Fig. 2e, f). Although some differences were observed in the classification between the proteomics and transomics analyses, many similarities were also identified (Fig. 2d). Transomics analysis suggests that there may be molecules that exhibit variation in only protein expression that reflect the characteristics of pituitary adenomas/PitNETS to the best of our knowledge.

To explore the intrinsic cohort structure using the full complement of single-omics data, clustering was performed for mRNA and protein using non-negative matrix factorization (NMF)²⁵, which yielded between 2 and 4 clusters for the single-omics analyses (Fig. 2g). A GO analysis of a group of molecules with significantly higher expression in Cluster 3 by proteomics revealed the terms GH synthesis, secretion and action (Supplementary Data 4). Transcriptomic data were only able to distinguish between NFPAs/non-functioning PitNETs and GHomas/somatotroph PitNETs, whereas proteomics data were able to distinguish among GHomas/somatotroph PitNETs with and without GNAS mutation. Therefore, we focused on both proteomics and genomics approaches for the characterization of GHomas/somatotroph PitNETs.

Proteogenomic characterization of nonfunctioning and GHomas/somatotroph PitNETs

The samples in which the proteins could not be identified by proteomics analysis were excluded. The protein expression of key transcription factors was compared between NFPAs/non-functioning PitNETs and GHomas/somatotroph PitNETs using uniform manifold approximation and projection (UMAP) analysis (Supplementary Fig. 3). The results showed that NFPAs/non-functioning PitNETs and GHomas/somatotroph PitNETs could be divided into 2 distinct groups, likely because NFPAs/non-functioning PitNETs and GHomas/somatotroph PitNETs have distinctly different protein expression characteristics. However, among GHomas/somatotroph PitNETs, no clear differences in protein expression characteristics were identified between cases with and without GNAS-MT. The expression of key transcription factors that contributed to pituitary differentiation visualized using a feature plot, revealed that nuclear receptor subfamily 5 Group A member 1 (NR5A1) and GATA-binding protein 3 (GATA3) are expressed at high levels in most NFPAs/non-functioning PitNETs samples, with some GHomas/somatotroph PitNETs also demonstrating high expression levels. In contrast, POU class 1 homeobox 1 (POU1F1), a key regulatory factor, was highly expressed in GHomas/somatotroph PitNETs, whereas POU1F1 and GH were not expressed in NFPAs/non-functioning PitNETs. SSTR5 was expressed heterogeneously in most GHomas/somatotroph PitNETs, with some expression observed in NFPAs/non-functioning PitNETs (Fig. 3a).

**Fig. 3: Nontargeted proteomics analysis of pituitary adenomas/PitNETs.**

The Brunner–Munzel test was used to quantitatively compare protein expression levels between NFPAs/non-functioning PitNETs and GH-producing adenomas and between GNAS-WT and GNAS-MT GH-producing adenomas. A total of 7300 differentially expressed molecules were identified between NFPAs/non-functioning PitNETs and GHomas/somatotroph PitNETs, and 714 differentially expressed molecules were identified between GNAS-WT and GNAS-MT GH-producing adenomas. To evaluate the protein expression profiles of the 458 molecules that were identified as significantly differentially expressed in both comparisons (Fig. 3b, c), gene ontology (GO) analysis was performed. GO analysis (http://geneontology.org) revealed that the GNAS mutation influenced several binding functions and GPCR pathways (Supplementary Data 5). Considering these results collectively, we hypothesize that GNAS-MT induce GPCR pathways in GHomas/somatotroph PitNETs, influencing endocrinological characteristics, such as the response to drug treatments.

Correlation between the protein expression profile and clinical characteristics of Acromegaly

We first evaluated the expression of some characteristic proteins by pathological analysis. In our study, the immunohistochemistry (IHC) score for SSTR2, which is known to serve as the primary SSA receptor in GHomas/somatotroph PitNETs, did not differ between the GNAS-WT and GNAS-MT groups (Fig. 4a). Similarly, the IHC scores for GH and the CAM5.2 cytokeratin IHC pattern did not differ between the two groups (Fig. 4b, c).

**Fig. 4: Protein expression levels in NFPAs/non-functioning PitNETs, *GNAS*-WT GHomas/somatotroph PitNETs, and *GNAS*-MT GHomas/somatotroph PitNETs.**

Based on the results of the nontargeted proteomics analysis, we compared protein expression levels between three groups: NFPAs/non-functioning PitNETs, GNAS-WT GHomas/somatotroph PitNETs, and GNAS-MT GHomas/somatotroph PitNETs. The protein expression level of T-box transcription factor 19 (TBX19, also known as TPIT) did not differ among the 3 groups, whereas the expression level of POU1F1 was significantly higher in the GNAS-WT and GNAS-MT groups than in the NFPAs/non-functioning PitNETs group (Fig. 4d, e), consistent with the regulation of transcription factors during pituitary development. The protein expression levels of AIP were significantly higher in the NFPAs/non-functioning PitNETs group than in the GHomas/somatotroph PitNETs group (Fig. 4f). Conversely, SSTR2 expression was significantly higher in the GNAS-WT and GNAS-MT groups than in the NFPAs/non-functioning PitNETs group (Fig. 4g).

Next, we focused on 4 molecules: SSTR5, sigma nonopioid intracellular receptor 1 (SIGMAR1), adhesion G protein–coupled receptor V1 (ADGRV1), and sortilin-related VPS10 domain–containing receptor 3 (SORCS3). These molecules are involved in GPCR activity and were among the 458 identified differentially expressed proteins. The expression levels of these molecules were significantly different between the NFPAs/non-functioning PitNETs and GHomas/somatotroph PitNETs groups: SSTR5 expression was higher in the GHomas/somatotroph PitNETs group than in the NFPAs/non-functioning PitNETs group, and there was no difference in expression with or without GNAS-MT (Fig. 4h). The SSTR2/5 protein expression ratio was not significantly different between the NFPAs/non-functioning PitNETs and GHomas/somatotroph PitNETs groups, nor were they expressed with or without GNAS-MT in the GHomas/somatotroph PitNETs group (Fig. 4i). SIGMAR1 expression was lower in the GHomas/somatotroph PitNETs group than in the NFPAs/non-functioning PitNETs group and was lower in the GNAS-MT groups than in the GNAS-WT group (Fig. 4j). ADGRV1 expression was not significantly different between the NFPAs/non-functioning PitNETs and GHomas/somatotroph PitNETs groups, nor was it different with or without GNAS-MT in the GHomas/somatotroph PitNETs group (Fig. 4k). SORCS3 expression was higher in GHomas/somatotroph PitNETs GNAS-WT than in NFPAs/non-functioning PitNETs and lower in GNAS-MT expression than in GNAS-WT (Fig. 4l). To determine whether these GPCR-related molecules affect the clinical characteristics of acromegaly, we analyzed the correlations between the expression levels of these four proteins and clinical characteristics, including the GH change rate following the octreotide loading test and the tumor volume change rate following SSA treatment. The expression levels of ADGRV1 and SORCS3 were correlated with the GH change rate following the octreotide loading test (Supplementary Fig. 4).

The GH change rate following the octreotide loading test, which has been reported to correlate with the SSTR2 mRNA expression level, did not correlate with the SSTR2 or SSTR5 protein expression levels or with the SSTR2/5 protein expression ratio (Fig. 5a–c). We then performed correlation analyses between all of the protein expression values derived from the nontargeted proteomics analysis and the GH change rate following the octreotide loading test. We detected positive correlations for ATPase sarcoplasmic/endoplasmic reticulum Ca2+ transporting 2 (ATP2A2) and AT-rich interaction domain 5B (ARID5B) with the GH change rate following the octreotide loading test, although no correlation was observed for ATPase sarcoplasmic/endoplasmic reticulum Ca2+ transporting 1 (ATP2A1) (Fig. 5d–f). The protein expression levels of ATP2A2 and ATP2A1 were significantly lower in the GHomas/somatotroph PitNETs group than in the NFPAs/non-functioning PitNETs group, and ARID5B expression was significantly higher in the GHomas/somatotroph PitNETs group. These 3 molecules were expressed at significantly lower levels in the GNAS-MT group than in the GNAS-WT group (Fig. 5g–i)

**Fig. 5: Correlation between protein expression levels and GH change rates in the octreotide loading test.**

To clarify the pathology of functional pituitary adenomas/PitNETS, we performed correlation analyses between the nontargeted proteomics data and the tumor volume change rate following SSA treatment because it is important to consider not only the production and secretion of hormones but also tumorigenesis. SSTR2 expression, SSTR5 expression, and the SSTR2/5 ratio did not correlate with the tumor volume change rate following SSA treatment (Fig. 6a–c). The tumor volume change rate following SSA treatment was correlated with the expression levels of WWC family member 3 (WWC3) and serine incorporator 1 (SERINC1) and tended to correlate with zinc finger AN1-type containing 3 (ZFAND3) expression (Fig. 6d–f). The expression of WWC3 was significantly lower in the GNAS-MT group than in the GNAS-WT group, whereas SERINC1 and ZFAND3 expression levels were significantly higher in the GNAS-MT group than in the GNAS-WT group (Fig. 6g–i).

**Fig. 6: Correlation between protein expression levels and tumor volume change rate in the SSA test.**

These differentially expressed protein candidates from the nontargeted proteomics analysis were evaluated by pathological examination. Sections were stained with SIGMAR1, ATP2A2, ARID5B, WWC3, and SERINC1 (Fig. 7a–e), and IHC scoring was performed as described in the Methods. SIGMAR1 was significantly increased in the NFPAs/non-functioning PitNETs group compared with the GNAS-MT group. (Fig. 7f). ATP2A2 was significantly increased in the NFPAs/non-functioning PitNETs group compared with the GHomas/somatotroph PitNETs groups (Fig. 7g). ARID5B was significantly increased in the GNAS-WT group compared with the GNAS-MT group (Fig. 7h). There were no significant differences in WWC3 and SERINC1 among the three groups. (Fig. 7i, j). These expression patterns were similar to the trends observed using nontargeted proteomics analysis.

**Fig. 7: Immunohistochemistry of differentially expressed candidates in nontargeted proteomics.**

Tags: adenomassomatotroph, characterization, clinical, GHproducing, Landscape, neuroendocrine, pituitary, Proteogenomic, tumors