Genetic and Epigenetic Characterization of Growth Hormone–Secreting Pituitary Tumors

Somatic SNV and indel background of Gsp⁻ somatotropinomas

The somatic landscape of SNVs and SVs of the tumors ST2-ST12 was reported previously (8). The tumors had an average of 2.3 coding region SNVs per tumor, with GNAS being the only recurrently mutated gene. Here, we examined the somatic SNVs and SVs of six additional Gsp⁻ tumors: ST13, ST16, ST18-20, and ST22. Supplementary Table S1 gives all the coding region (missense, premature stop codon, frameshift) variants that passed the somatic filtering. These additional Gsp⁻ tumors had an excess of somatic variants, in total 92 coding region SNVs and on average 15.3 SNVs per tumor, simply because rare germline variants may have passed the population-based filtering. The majority of the somatic variants (47 SNVs) arise from ST22 due to patient's Italian ancestry being underrepresented among the population-based controls. Among the 21 tumors studied the only recurrently mutated gene was GNAS. No focal deletions or complex chromosomal rearrangements were observed with exception of ST3, a previously reported Gsp⁺ tumor (Supplementary Fig. S1; ref. 8).

SCNA and expression profiling differences in CA- and CS-tumors

Data from 21 somatotropinoma normal/tumor pairs was used for SCNA analysis. Results of eleven normal/tumor pairs (ST2-ST12) were from our earlier work (8). Analyses revealed that 12 Gsp⁻ tumors contained frequent and recurrent (≥4 tumors share the event) chromosomal deletions (chr 1, 4, 6, 14, 15, 16, 18, 22). Also gains (chr 5, 7, 9, 19, 20) of entire chromosomes were detected, although with considerably lower frequency (Supplementary Fig. S1). Copy-neutral LOH or homozygous deletions were not observed. All detected chromosomal gains were duplications of a single chromosome. Gsp⁺ tumors contained limited amount of aneuploidy, mostly gains of single chromosomes. An exception was a tumor ST4 with genetic instability (GI%) 22% (Table 1; Supplementary Fig. S1). In addition, two Gsp⁻ tumors (ST5 and ST6) were totally aneuploidy-free. GI% did not correlate with clinical variables (Supplementary Table S2). To test whether the recurrent aneuploidy is statistically significant, the permutation test was employed (8). The test indicated that the observed recurrent aneuploidy indicates selective advantage during tumorigenesis rather than random copy number alteration (permutation test P < 10⁻⁴; Supplementary Fig. S2).

An unsupervised hierarchical clustering analysis of gene expression demonstrated that aneuploidy was the major determinant of somatotropinoma subtypes. Tumors clustered according to the aneuploidy rate as follows (Fig. 1A; principal components analysis in Supplementary Fig. S3). Twelve Gsp⁻ tumors with aneuploidy, from now on called CA (chromosome aneuploidy)- tumors, clustered together. Accordingly, chromosome stable (CS)-tumors with limited amount or no aneuploidy (7 Gsp⁺ and 2 Gsp⁻) formed their own distinct group. Because CA- and CS-tumors had their own expression signatures (Fig. 1A) expression profile comparison was performed between CA- and CS-tumors. We identified 881 differentially expressed transcripts (FDR < 5%, |FC| > 2; Supplementary Table S3). Integration of expression and SCNA data revealed that 69.8% of differentially expressed transcripts locate at chromosomes with recurrent aneuploidy (≥4 tumors with shared aneuploidy) from which 81.6% positively correlate with the chromosomal copy number change.

• Download figure
• Open in new tab
• Download powerpoint

Figure 1.

A, Gene expression of somatotropinomas clustered according to the aneuploidy rate. Left, result of unsupervised hierarchical clustering of expression array data from 21 somatotropinomas; middle, patients' age at diagnosis, gender, and Gsp mutation status; right, somatic copy number aberrations. Supplementary: Supplementary Fig. S1 shows all somatic chromosomal aberrations in more detail. Of note, the chr 1p of ST3 contains chromothripsis event (8). B, Enriched pathways result from differentially expressed genes (CA- vs. CS-tumors). Supplementary Table S4 shows the complete list of significantly (P < 0.05) enriched expression pathways and genes. The blue horizontal bars denote the association P values for each pathway on a logarithmic scale (dashed vertical line at P = 0.05). The ratio between the number of query genes found and total number of genes in a pathway is shown in orange. C, The enriched G₁–S signaling pathway. The colored molecules identified as differentially regulated in CA- versus CS-tumor groups. FCs and P values are listed in Table 3. Orange label, upregulated; green label, downregulated.

To understand the biological relevance of differentially expressed genes in CA-somatotropinomas, pathway analysis was performed. The top canonical pathway emerged from the tumor subtype-specific expression was the PKA pathway (Fig. 1B; Supplementary Table S4). In this cAMP-mediated signaling pathway a majority of molecules were upregulated in CS-tumor group indicating activated PKA signaling (Fig. 2; Table 2). The major inhibitors of cAMP levels (19), phosphodiesterases (PDE) were upregulated in CS-tumor group and seven protein tyrosine phosphatase (PTP) receptors were differentially expressed between tumor groups (Table 2).

• Download figure
• Open in new tab
• Download powerpoint

Figure 2.

The enriched PKA signaling pathway. The colored molecules identified as differentially expressed genes or gene groups in CA- versus CS-tumors. FCs and P values are listed in Table 2. Orange label, upregulated; green label, downregulated.

The expression pathway analysis also highlighted several pathways associated with cell-cycle regulation and explicitly in the retinoblastoma 1/E2F transcription factor (RB1/E2F)–mediated G₁–S-checkpoint transition. Among the enriched pathways were molecular mechanisms of cancer, chronic myeloid leukemia, glioma, glioblastoma multiforme, and cell cycle: G₁–S checkpoint regulation (Fig. 1B; Supplementary Table S4). The genes showing expression differences between tumor groups accumulated in the RB1/E2F–mediated cell-cycle regulation and the subsequent G₁–S-checkpoint transition indicating difference in mitotic progression (Fig. 1C). E2F4 transcription factor, a regulator of cell cycle (20), was downregulated in CA-tumor group. By contrast, histone deacetylase 5 (HDAC5) and RB transcriptional corepressor like 1 (RBL1/p107), both repressors of E2F family members, were upregulated in CA-tumor group when compared to CS-tumors. Also, the antimitogenic growth factor TGFβ2/SMAD3 signaling was downregulated in CA-tumor group (Fig. 1C; Table 3).

IHC

Because RB1 is the major component of the complex regulating G₁–S phase transition, protein levels of total- and phosphorylated-RB1 were semiquantitatively assessed in tumors using IHC. In all tumors >90% of nuclei showed positive total- and phospho-RB1 staining. Because fractions of stained cells were comparable between tumor groups, staining intensities were compared. In both CA- and CS-tumor groups was seen weak to moderate nuclear immunoreactivity of total-RB1 (CA- vs. CS tumors; 1.67 vs. 1.78, P = 0.6, Student t test). In addition, some occasional cytoplasmic total-RB1 staining was detected in both tumor groups. Phospho-RB1 showed weak to moderate nuclear staining in both tumor groups (CA- and CS- tumors; 1.59 vs. 1.67, P = 0.71; Supplementary Table S5; Supplementary Fig. S4).

To test pathway analysis findings of the PKA signaling (Fig. 2) at the protein level, PTPRD and PTPRK IHCs were performed. For this, we used tumor material from four CA- (ST10, ST12, ST18, ST19) and three CS-tumors (ST15, ST17, ST21). PTPRD showed negative immunostaining in CA-tumors (Gsp⁻) and moderate immunoreactivity in CS-adenomas (Gsp⁺; Supplementary Fig. 5A). PTPRK immunostaining was negative or weak in CA-tumors and moderate in CS-tumors (Supplementary Fig. 5B). In immunopositive tumors, both PTPRD and PTPRK localized mainly in cytoplasm and >90% of cells gave positive immunostaining.

To investigate the role of PTTG1 in the formation of pituitary neoplasia and proliferation rates in these tumors, PTTG1 and Ki-67 immunostainings were performed. Ki-67 and PTTG1 proteins were expressed in all tumors. Ki-67 gave a clear nuclear immunostaining. There was no significant difference in Ki-67 scores between tumor groups (x = 2.7 ± 1.6 vs. 1.7 ± 1.2; P = 0.27; Table 1). PTTG1 immunopositive cells showed both nuclear and cytoplasmic localization of the protein. Most of the cells showed predominant cytoplasmic staining, although there were also cells with predominant nuclear localization of PTTG1 in both tumor groups (Supplementary Fig. S6). All stained cells were scored and number of immunopositive cells were significantly more abundant in CA-tumors compared with CS-adenomas (x = 1.8 ± 1 vs. 0.9 ± 0.5; P = 0.016; Wilcoxon–Mann–Whitney; Table 1). PTTG1 RNA expression levels, emerging mainly from nonproliferating cells, were comparable between CA- and CS-tumor groups (FC = 1.05, P = 0.733). As has been seen earlier (8, 21), the PTTG1 protein levels correlated with Ki-67 scores (r = 0.62, P = 0.002; Pearson coefficient).

DNA methylation

We surveyed DNA methylation and DMRs in 21 somatotropinomas by targeted bisulfite sequencing (Supplementary Table S6). An unsupervised clustering of 50,000 CpG regions with the largest variance between tumors showed, that apart from ST7 and ST21, the tumors clustered according to the Gsp mutation status (Fig. 3A; principal components analysis in Supplementary Fig. S7). The genome-wide distributions of CpG methylation levels for different genomic contexts revealed that methylation levels of promoter regions were comparable across tumors (Fig. 3A). The rest of the annotated regions showed hypomethylation of Gsp⁺ tumors (see below for the supervised analysis of Gsp⁺ and Gsp⁻). Eighteen percent of the differentially expressed (CA- vs. CS) transcripts (119/670 coding transcript clusters; Supplementary Table S3) were correlated with DNA methylation.

• Download figure
• Open in new tab
• Download powerpoint

Figure 3.

A, CpG methylation of somatotropinomas clustered according to the Gsp mutation status. Left, unsupervised hierarchical clustering of CpG methylation data from 21 somatotropinoma samples (based on 50,000 regions with the largest variance between tumors); middle, patient's age at diagnosis, gender, and Gsp mutation status; right, distribution of estimated CpG methylation levels for each genomic context: gene body (from +2 kb relative to the TSS to the end of gene), promoter (−1 kb to +2 kb relative to the TSS), CpG island/shore, enhancer, DNase (DNaseI Hypersensitivity Cluster), and TFBS (Transcription factor cluster; see Materials and Methods for details). Methylation levels are quantified by value ranging from zero (unmethylated) to one (fully methylated). B, Pathway analyses from cis-eQTMs. Supplementary Table S10 shows the list of significantly (P < 0.05) enriched cis-eQTM pathways and genes. The blue horizontal bars denote the association P values for each pathway on a logarithmic scale (dashed vertical line at P = 0.05). The ratio between the number of query genes found and total number of genes in a pathway is shown in orange. C, The enriched CREB signaling in neurons pathway. The colored molecules identified as differentially expressed Gsp⁺ versus Gsp⁻ tumors. eQTM and FDR values are listed in Table 4. Orange label, upregulated; green label, downregulated.

DNA methylation patterns are maintained and regulated by DNA methyltransferases (DNMT), including DNMT1, DNMT3A, and DNMT3B (22). We examined tumor-specific associations of DNMT expressions and median CpG methylation and found that CpG methylation rates correlated with the DNMT1 expression (Spearman rank correlation 0.49; P = 0.025; Supplementary Fig. S8), while DNMT3A (0.02, P = 0.929) and 3B (0.24, P = 0.294) did not show correlation.

DNA methylation associates also with replication timing (17). Supplementary Figure S9 shows an overview of CpG methylation at different quartiles of replication timing. Majority of the tumors displayed the expected hypomethylation of late replicating regions (Spearman ρ < 0 and P < 0.019; Supplementary Table S7). The Gsp⁻ tumors ST10, ST11, ST16, and ST20 had an outstanding, positive correlation to replication time (Spearman ρ >0 and P < 0.018), which suggests methylation maintenance also at late replicating regions. No clinical associations were found to explain the methylation maintenance difference in these four tumors (Table 1).

Because the Gsp mutation status was the major factor behind the DNA methylation rates and profiles across tumors, DMRs were determined between Gsp⁺ versus Gsp⁻ somatotropinomas. Altogether, we found 1,369 DMRs out of which 1,339 (97.8%) were hypomethylated in Gsp⁺ tumors: see Supplementary Table S8 for a complete list of all 1,369 regions' genomic coordinates and annotation of nearby (-20 Kbp upstream; 2 Kbp downstream) genes (1,560 gene annotations). Supplementary Table S8 is sorted by the absolute value of the test statistic, where negative effect direction denotes hypometylation among the Gsp⁺ tumors compared to the Gsp⁻ tumors. Both the outstanding number of DMRs and the enrichment of Gsp⁺ tumors toward hypomethylation can likely be attributed to the genome-wide CpG methylation characteristics between the tumor types (see the unsupervised analysis and Fig. 3A). The DMRs did not enrich among the aneuploidy chromosomes; 51% of the DMRs reside at recurrent aneuploidy, while the expected proportion was 55% based on the distribution of the TBS regions.

In addition to the DMR analysis, we also examined the genome-wide methylation profiles between Gsp⁺ and Gsp⁻ tumors in different genomic contexts. Supplementary Figure S10 displays the CpG methylation levels in the context of promoters, enhancers, CpG islands and CCCTC-binding factor (CTCF) sites. The only systematic difference between the Gsp⁺ and Gsp⁻ tumors at these genomic contexts was the overall hypomethylation in Gsp⁺ tumors.

To assess the DNA methylation–expression association, the CpG methylation data were integrated with matched transcriptomes by eQTM analysis (23, 24). Consistent with previous studies, these associations account for only a small fraction of the assayed CpG sites and expressed genes (24). In a supervised comparison, the DMRs between Gsp⁺ and Gsp⁻ tumors resulted in a total of 400 DMR regions with significant (FDR < 0.05) expression to methylation eQTMs pertaining 155 genes (Supplementary Table S9; see the Supplementary Fig. S11 for an unsupervised analysis of eQTMs). Altogether, 392 DMR regions (98%) were hypomethylated in Gsp⁺ tumors and from these 333 (84.9%) showed positive eQTM⁺ association (hypomethylated in Gsp⁺ together with gene underexpression in Gsp⁺). This enrichment can be attributed mostly to the tumor subtype-specific differences in CpG methylation characteristics at gene body regions (Fig. 3A).

To identify biological functions associated with the methylation difference between Gsp⁺ and Gsp⁻ somatotropinomas, we performed pathway analysis from the eQTM gene list. The majority of emerged pathways (22/31) were associated with inhibitory Gα protein (Gα_i) – and/or voltage-gated calcium channel (CaCn) signaling (Fig. 3B; Supplementary Table S10). Both Gα_i and Ca²⁺ signaling are previously connected to pituitary neoplasia (25, 26). The inhibitory G protein subunit alpha I-2 (GNAI2; Gα_i-2) and G protein subunit beta 1 (GNB1, Gβ₁) were upregulated in Gsp⁺ tumors (Fig. 3C; Table 4). The CaCn members, calcium voltage-gated channel subunit alpha 1A and 1E (CACNA1A and CACNA1E) and calcium voltage-gated channel auxiliary subunit gamma 2 (CACNG2), were downregulated and correlated with hypomethylation (Supplementary Table S8).