Blog

By Azusa Yamato, Hidekazu Nagano, Yue Gao, Tatsuma Matsuda, Naoko Hashimoto, Akitoshi Nakayama, Kazuyuki Yamagata, Masataka Yokoyama, Yingbo Gong, Xiaoyan Shi, Siti Nurul Zhahara, Takashi Kono, Yuki Taki, Naoto Furuki, Motoi Nishimura, Kentaro Horiguchi, Yasuo Iwadate, Masaki Fukuyo, Bahityar Rahmutulla, Atsushi Kaneda, Yoshinori Hasegawa, Yusuke Kawashima, Osamu Ohara, Tetsuo Ishikawa, Eiryo Kawakami, Yasuhiro Nakamura, Naoko Inoshita, Shozo Yamada, Noriaki Fukuhara, Hiroshi Nishioka, and Tomoaki Tanaka

Excerpt from the article published in Communications Biology 5, 1304, 27 November 2022, DOI: https://doi.org/10.1038/s42003-022-04272-1

Editor’s Highlights

• Pituitary adenomas/pituitary neuroendocrine tumors (pituitary adenomas/PitNETs) account for ~15% of all primary intracranial neoplasms.
• 57.0% of patients had stimulatory Gs protein mutations (GNAS-MT).
• Sigma nonopioid intracellular receptor 1 (SIGMAR1) is involved in GPCR activity and was among the 458 identified differentially expressed proteins.
• SIGMAR1 was significantly increased in the NFPAs/non-functioning PitNETs group compared with the GNAS-MT group.
• SIGMAR1 plays an important role in the cellular functions of various tissues associated with the endocrine, immune, and nervous systems, in addition to cancer cells.
• Although SSTR5 and SIGMAR1 are not associated with GH change by octreotide loading test and tumor volume change by preoperative SSA treatment, these molecules can regulate GH secretion in response to GPCR signals.
• SIGMAR1 molecules were important factors underlying the clinical and biochemical features involved in the responsiveness of GHomas/somatotroph PitNETs to medical treatment.

Abstract

The clinical characteristics of growth hormone (GH)-producing pituitary adenomas/somatotroph pituitary neuroendocrine tumors (GHomas/somatotroph PitNETs) vary across patients. In this study, we aimed to integrate the genetic alterations, protein expression profiles, transcriptomes, and clinical characteristics of GHomas/somatotroph PitNETs to identify molecules associated with acromegaly characteristics. Targeted capture sequencing and copy number analysis of 36 genes and nontargeted proteomics analysis were performed on fresh-frozen samples from 121 sporadic GHomas/somatotroph PitNETs. Targeted capture sequencing revealed GNAS as the only driver gene, as previously reported. Classification by consensus clustering using both RNA sequencing and proteomics revealed many similarities between the proteome and the transcriptome. Gene ontology analysis was performed for differentially expressed proteins between wild-type and mutant GNAS samples identified by nontargeted proteomics and involved in G protein–coupled receptor (GPCR) pathways. The results suggested that GNAS mutations impact endocrinological features in acromegaly through GPCR pathway induction. ATP2A2 and ARID5B correlated with the GH change rate in the octreotide loading test, and WWC3, SERINC1, and ZFAND3 correlated with the tumor volume change rate after somatostatin analog treatment. These results identified a biological connection between GNAS mutations and the clinical and biochemical characteristics of acromegaly, revealing molecules associated with acromegaly that may affect medical treatment efficacy.

Introduction

Pituitary adenomas/pituitary neuroendocrine tumors (pituitary adenomas/PitNETs) account for ~15% of all primary intracranial neoplasms. Among them, growth hormone (GH)-producing pituitary adenomas/somatotroph pituitary neuroendocrine tumors (GHomas/somatotroph PitNETs) are the second most common hormone-producing adenomas after prolactin-producing adenomas^1,2. Acromegaly, due to GHomas/somatotroph PitNETs, presents with a wide variety of symptoms, such as the enlargement of the distal extremities and tongue, cardiac hypertrophy, osteoarthropathy, metabolic disorders, and malignancy. These complications can lead to a high mortality rate if GH overproduction remains uncontrolled. The first-line treatment for GHomas/somatotroph PitNETs is surgical tumor resection. Medical treatments, including somatostatin analog (SSA) or GH receptor antagonists, may be used in cases of noncurative resection or for preoperative tumor volume reduction. SSAs have efficacy in suppressing the secretion of GH and shrinking tumors when they bind to the somatostatin receptor on the surface of the tumor cell membrane³. Disappointingly, some patients do not achieve remission, regardless of the treatment strategy. Therefore, the elucidation of the underlying biological mechanisms may identify targets with more effective treatment effects.

Somatic mutations in the α-subunit of the stimulatory Gs protein (GNAS) gene, resulting in the constitutive activation of cAMP, are observed in 30–50% of patients with acromegaly⁴. Although whole-exome analyses of GHomas/somatotroph PitNETs revealed that GNAS was the only driver gene of GHomas/somatotroph PitNETs^1,2, previous reports have identified other causative gene mutations. Germline mutations in aryl hydrocarbon receptor-interacting protein (AIP), multiple endocrine neoplasia type 1 (MEN1), and cyclin-dependent kinase inhibitor 1B (CDKN1B) have been reported in a small percentage of young patients with acromegaly^5,6,7,8,9 and germline or somatic microduplications of the chromosome region Xq 26.3, which contains the gene encoding the G protein–coupled receptor 101 (GPR101), have been implicated in younger-onset gigantism¹⁰. Previous reports examining genotype–phenotype relationships in patients with acromegaly showed that patients with GNASmutations (GNAS-MT) presented with smaller tumors and had better SSA responsiveness than those without GNAS-MT^11,12, although the relationships between other gene mutations and clinical features remain unclear. Several reports have presented correlations between somatostatin receptor 2 (SSTR2) mRNA expression levels and the efficacy of SSAs^3,13,14. However, these findings do not sufficiently explain the various clinical characteristics of acromegaly.

Accumulating evidence shows the involvement of genetic events other than single mutations or small indels in pituitary adenomas/PitNETs. Several studies have identified high levels of genomic instability and a high frequency of copy number alterations (CNAs) in pituitary adenomas/PitNETs^{15,16,17,18,19,20}, and genomic CNAs may be especially frequent in GHomas/somatotroph PitNETs^2,21. An association between DNA damage and cAMP activation has also been indicated¹⁹. Epigenetic alterations in GHomas/somatotroph PitNETs have also been reported^19,22, and multiomics analyses have revealed potential alterations in gene expression patterns in GH-producing adenomas. Salmon et al. reported an important relationship between DNA methylation patterns and gene expression profiles in pituitary adenomas/PitNETs²², and Neou et al. presented molecular profiles of pituitary adenomas/PitNETs based on the performance of various omics analyses²³.

Although recent studies of pituitary adenomas/PitNETs have aimed to clarify the underlying tumorigenesis mechanisms and disease etiologies, the association between genetic alterations and clinical characteristics remains unclear, and the factors that affect responsiveness to treatment are unknown. In addition, few genetic studies have been limited to GHomas/somatotroph PitNETs. Here, we performed gene alteration analysis and proteomics analysis on 121 GHomas/somatotroph PitNETs and integrated these results with the clinical characteristics of acromegaly. We attempted to identify key players involved in shaping the clinical features of acromegaly, especially those related to treatment efficacy. Our study revealed the importance ofGNAS-MT in terms of clinical and biochemical characteristics and identified molecules that may be involved in the responsiveness to medical treatment.

Results

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% of patients presented with a high Knosp grade (≥3). All 121 patients underwent transsphenoidal surgery (TSS), and 67 patients (55.4%) 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% 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).

	All patients (n = 121)	GNAS wild-type (n = 52)	GNAS mutant (n = 69)
Clinical characteristics				p value
Age (years)	48.0 (39.0–59.5)	46.5 (40.0–60.5)	49.0 (38.5–59.5)	0.763
Sex, n (male/female)	51/70	20/32	31/38	0.577
Basal GH(ng/mL)	15.7 (7.2–34.1)	13.4 (4.6–31.3)	16.9 (8.3–48.5)	0.061
IGF-1 (ng/mL)	630.0 (473.0–817.5)	602.0 (448.0–829.3)	237.5 (505.0–817.5)	0.745
IGF-1 SDS	7.0 ± 2.4	7.2 ± 2.6	6.9 ± 2.3	0.580
PRL (ng/mL)	13.4 (8.1–29.1)	10.3 (7.4–21.0)	20.2 (9.8–37.3)	0.005
GH change by octreotide test (%)	−87.7 (−93.9–−61.3)	−83.4 (−92.0–−57.1)	−89.6 (−95.0–−69.6)	0.063
GH change by bromocriptine test (%)	−65.6 (−83.8–−23.7)	−39.3 (−68.0–−7.3)	−79.1 (−90.2–−55.4)	<0.001
Tumor volume (mm3)	1633.1 (651.9–4412.4)	2323.7 (1020.3–4990.0)	1135.4 (406.2–3493.1)	0.023
Knosp grade 0–2/3–4, n	93/28	33/19	60/9	0.004
preoperative therapy
preoperative SSA treatment, n (%)	67 (55.4%)	33 (63.5%)	34 (49.3%)
GH change by preoperative SSA (%)	−78.9 (−90.8–−48.1)	−55.5 (−80.8–−41.7)	−83.2 (−92.7–−64.2)	0.007
Tumor volume change by preoperative SSA (%)	−23.1 ± 21.3	−19.9 ± 18.6	−26.2 ± 23.5	0.237
postoperative profile
Nadir GH after OGTT (ng/mL)	0.42 (0.27–0.84)	0.50 (0.26–1.0)	0.39 (0.27–0.67)	0.514
IGF-1 (ng/mL)	189.0 (137.0–227.0)	186.5 (153.0–243.2)	190.0 (130.0–220.0)	0.237
IGF-1 SDS	0.7 ± 1.6	0.9 ± 1.7	0.5 ± 1.5	0.119
postopertive treatment, n (%)	21 (17.4%)	9 (17.3%)	12 (17.4%)	0.990

We observed various activating GNAS-MT (57.0% 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% of the genes showed a positive mRNA-protein correlation, only 65.2% 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).

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).

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)

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).

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.

Discussion

In this study, we integrated clinical characteristics, genetic analysis, and biochemical analysis methods, including IHC and proteomics, to assess GHomas/somatotroph PitNETs in a large cohort. Using transomics classification, we showed that mRNA and protein levels were modestly correlated in pituitary adenomas/PitNETs, similar to reports for human colon and rectal cancer studies²⁴. TCS analysis revealed that GNAS was the only driver gene of GHomas/somatotroph PitNETs, as previously reported, and the nontargeted proteomics analysis revealed that GNAS was an important player in shaping various acromegaly characteristics.

The present analysis revealed that 57.0% of patients had GNAS-MT. The prevalence of GNASin GHomas/somatotroph PitNETs has a wide range, 4–59%, according to previous reports²⁶. The clinical characteristics of patients with GNAS-MT were consistent with those in previous reports^11,12. However, the clinical data obtained for each group presented an extremely wide range, indicating that many patients showed exceptional characteristics in some respects. Actually, among reports with low prevalence, there are reports of significant characteristics with GNAS-MT¹² or no characteristic differences with or without GNAS-MT²⁶. Next, for recurrent genes aside from GNAS in acromegaly, we identified 18 recurrent genes using Neou M et al.‘s method. There are previous reports that the only recurrent gene is GNAS^1,27, while Neou M et al. identified ~30 recurrent genes²³. This finding suggests that the clinical characteristics of acromegaly are not reflective of a single genetic mutation but are the result of several genetic events and gene expression changes, consistent with previous reports regarding the involvement of CNAs or epigenetic alterations²⁷. The investigation of associations between these genetic events and the clinical phenotypes of acromegaly will be explored in future studies. Although GNAS-MT cannot explain all of the clinical characteristics, proteomics analysis showed that GNAS-MT were related to many differentially expressed proteins. Notably, GO analysis revealed the impacts of GNAS-MT on the GPCR pathway. GPCRs play important roles in GHomas/somatotroph PitNETs as signal transducers of GH-releasing hormone, leading to the production and secretion of GH, and as SSTRs, which regulate GH secretion and cell proliferation. GNAS mutations (GNAS-MT) have been proposed to trigger different acromegaly clinical characteristics through the altered expression of components in these important pathways. In particular, the induction of the GPCR pathway may contribute to differences in the GH change rate by octreotide loading test which is known to correlate with GH change rate following SSA treatment, based on the comparison of clinical phenotypes observed in our study. However, SSTR2, which is the typical GPCR expressed in GHomas/somatotroph PitNETs, was not identified as a differentially expressed molecule involved in the GPCR pathway. Conventional SSAs, such as octreotide long-acting release (LAR) or lanreotide, bind primarily with SSTR2, and some previous reports have identified differences in SSTR2 mRNA expression in patients with and without GNAS-MT^28,29. In addition, neither the SSTR2 IHC score nor the protein expression level derived from the proteomics analysis correlated with the GH change rate following octreotide test or with the tumor volume change rate following SSA treatment, although some reports showed that the SSTR2 mRNA expression level was positively correlated with the response to SSA in GHomas/somatotroph PitNETs^3,13,14. However, our analysis was based on protein expression levels rather than mRNA expression levels, and SSTRs may be downregulated in patients preoperatively treated with SSA, although we did not identify differences in SSTR2 protein expression levels between patients with and without preoperative SSA treatment (Supplementary Fig. 5).

SIGMAR1 plays an important role in the cellular functions of various tissues associated with the endocrine, immune, and nervous systems, in addition to cancer cells. SIGMAR1 primarily functions in physiological and pathophysiological processes of the CNS, such as pain, memory, neurodegenerative diseases, stroke, and addiction. SIGMAR1 is activated in response to tissue injury and during disease development to promote cell survival. However, interactions between SIGMAR1 and ion channels may change cellular behaviors in response to the microenvironment, leading to the unexpected consequence of tumor development³⁰. Although SSTR5 and SIGMAR1 are not associated with GH change by octreotide loading test and tumor volume change by preoperative SSA treatment, these molecules can regulate GH secretion in response to GPCR signals.

ATP2A2 and ARID5B were identified as being correlated with the GH change rate in response to octreotide treatment, and WWC3, SERINC1, and ZFAND3 were identified as being correlated with the tumor volume change rate in response to SSA treatment. All of these proteins demonstrated significant differences in expression between the GNAS-WT and GNAS-MT groups. ATP2A2 is an intracellular pump located in the sarcoplasmic or endoplasmic reticulum. Interestingly, a physical association with ATP2A2 has been described as a universal feature among GPCRs, which may be required for the efficient folding or membrane integration of GPCRs³¹. ARID5B forms a histone H3K9Me2 demethylase complex with PHD finger protein 2 and regulates the transcription of target genes³². Recently, several reports have indicated that chromatin modification genes are dysregulated in pituitary adenomas/PitNETs; for example, pituitary adenomas/PitNETs exhibit increased acetylation of H3K9³³. WWC3 is an upstream regulator of the Hippo signaling pathway, which controls cell proliferation and organ growth³⁴. SERINC1 is a membrane protein whose expression is restricted to the CNS³⁵. ZFAND3 is a zinc finger protein involved in nucleic acid recognition, transcriptional activation, protein folding, and assembly and causes tumor invasion in glioblastoma³⁶. Therefore, these molecules may be related to GH secretion or tumor development. In present study, pathological expression pattern was consistent with protein expression level derived from nontargeted proteomics analysis in SIGMAR1, ATP2A2 and ARID5B. Among them, ATP2A2 and ARID5B expression levels were significantly correlated with the GH change rate after the octreotide loading test, which correlated with the effect of SSA treatment. Currently, in clinical practice, SSTR2 immunostaining for surgical specimens of GHomas/somatotroph PitNETs is often performed as a predictor of the effect of additional postoperative SSA therapy. On the other hand, some reports have shown no correlation between SSTR2 pathological expression and SSA efficacy, just as we have shown above in protein expression data derived from nontargeted proteomics analysis. Our results suggested that ATP2A2 and/or ARID5B immunostaining may be a predictor of the efficacy of SSA treatment, either independently or in combination with SSTR2 immunostaining in noncurative resection. Further investigations into how these molecules are involved in the development of GHomas/somatotroph PitNETs or how they impact treatment efficacy are needed.

This study has two major limitations. First, we focused on only 36 genes by TCS; however, many other genes are likely involved in the tumorigenesis of GHomas/somatotroph PitNETs. Accordingly, we also examined the CNAs of only 36 genes. Second, we were unable to evaluate germline genomic information.

In summary, we integrated gene alteration data, proteomics data, and clinical information from a large cohort of patients with GHomas/somatotroph PitNETs. GNAS mutational status could not provide a complete explanation of the clinical characteristics of acromegaly; however, proteomics analysis showed that GNAS-MT influenced the expression of many proteins, including those involved in the GPCR pathway. We confirmed that these molecules were important factors underlying the clinical and biochemical features involved in the responsiveness of GHomas/somatotroph PitNETs to medical treatment.