The Epigenomic Landscape of Pituitary Adenomas Reveals Specific Alterations and Differentiates Among Acromegaly, Cushing's Disease and Endocrine-Inactive Subtypes Free

This article is featured in Highlights of This Issue, p. 4057

Pituitary adenomas account for approximately 15% of all primary brain tumors (1). They can be differentiated on the basis of the presence or absence of hormonal secretion (endocrine-active versus endocrine-inactive), clinical syndrome and manifestations, tumor size, and pathological features. On the basis of WHO criteria, approximately 85% of pituitary tumors are typical adenomas, whereas 11% to 18% are atypical, exhibiting more aggressive growth patterns (2, 3). Of the four most common adenoma subtypes, those associated with acromegaly (GH-secreting), Cushing's disease (ACTH-secreting) and endocrine-inactive adenomas are generally treated with transsphenoidal surgery as first-line therapy and are the focus of this analysis. In contrast, prolactinomas are typically treated with first-line dopamine agonist therapy. Because many pituitary adenomas invade the cavernous sinus or skull base structures, or in the case of some patients with Cushing's disease, the adenoma itself is so small that it cannot be found or completely removed at surgery. As a result, up to 40% of patients with acromegaly, Cushing's disease, or endocrine-inactive adenomas may need additional therapies beyond surgery. These treatments may include hormone-regulating medications or stereotactic radiosurgery to achieve tumor remission. Furthermore, after initial remission, up to 10% to 20% of pituitary adenoma patients may have disease recurrence. In patients with persistent or recurrent acromegaly and Cushing's disease, although each have distinct clinical syndromes, both can result in poorly controlled hypertension, diabetes, cardiac dysfunction, poor quality of life, and shortened lifespan. In patients with persistent or recurrent endocrine-inactive macroadenomas causing mass effect, patients may suffer from pituitary gland failure (hypopituitarism), vision loss, ophthalmoplegia, and/or persistent headaches. Given the significant differences in clinical manifestations in patients with these three adenoma subtypes, the relatively high rates of disease persistence and recurrence, and the associated morbidity due to ongoing hormonal hypersecretion and tumor mass effect, a more precise and comprehensive understanding of the molecular differences among pituitary adenomas is needed to improve tumor-specific management and long-term remission rates.

Recent studies indicate that pituitary adenomas have low mutation burdens (4–6), with approximately 60% of pituitary adenomas demonstrating no detectable genetic mutational background (7). To date, only two recurrent somatic mutations have been identified; one affecting the USP8 gene in ACTH-secreting pituitary adenomas (4, 5, 8) and the other affecting the GNAS gene in GH-secreting pituitary adenomas (9), the latter of which may play a role during the transition from prolactin (PRL)-secreting to GH-secreting pituitary adenomas (10). Likewise, syndromes caused by rare germline variants, including alterations in the genes menin 1 (MEN1) and Aryl hydrocarbon receptor Interacting Protein (AIP) found in pediatric and young patients with sporadic GH-secreting or PRL-secreting pituitary adenomas, account for a small number of familial cases (11). Thus, although recurrent somatic alterations and rare inherited variants have been identified, their frequencies among pituitary adenomas are low, suggesting that the etiology of a large fractions of sporadic pituitary adenomas are driven by factors other than gene mutation(s).

Increasing evidence points to alterations in gene expression as a major contributor to disease states in pituitary adenomas (12, 13). Transcriptome studies have highlighted several genes with altered patterns of expression among pituitary adenomas as compared with normal pituitary glands (13). In Cushing's disease, the proopiomelanocortin (POMC) and T-Box 19 (TBX19) genes were found to be upregulated compared with silent-corticotrophin adenomas and non-functioning pituitary adenomas (12), whereas high expression of EGFR is believed to be related to the hypersecretion of ACTH, a hallmark of Cushing's disease. In addition, recent studies have evaluated the role of epigenetic changes in the differential gene expression program of pituitary adenomas (14, 15), showing that endocrine-inactive pituitary adenomas have higher levels of DNA methylation compared with endocrine-active pituitary adenomas. Although some of these studies have associated epigenetic alterations with clinicopathological features of pituitary adenomas, the role and clinical utility of epigenetic profiling of these highly heterogeneous tumors remain poorly understood.

To date, genomic studies of pituitary adenomas have focused on contrasting the differences between endocrine-active and endocrine-inactive pituitary adenomas (15). However, grouping different endocrine-active pituitary adenomas into a single functional class can obscure important differences that exist among various subtypes of functional pituitary adenomas. To address this limitation, in this study, genomic, transcriptomic and DNA methylation data were generated from the three major subtypes of surgically resectable pituitary adenomas, GH-secreting (acromegaly), ACTH-secreting (Cushing's disease) and endocrine-inactive pituitary adenomas. By contrasting these two functional pituitary adenoma subtypes to endocrine-inactive pituitary adenomas, specific alterations that contribute to the unique biology of each tumor subtype can be uncovered. Furthermore, the effects of the somatostatin analogues (e.g. octreotide and somatuline) on the transcriptome and epigenome of a subset of GH-secreting adenomas is shown.

This retrospective study enrolled adult patients ages 18 to 90 years who had undergone transsphenoidal resection of pituitary adenomas at Providence Health System, Saint John's Health Center and had been consented for tissue and bodily fluid procurement and analysis. This study was conducted in accordance with recognized ethical guidelines (e.g., Declaration of Helsinki, CIOMS, Belmont Report, U.S. Common Rule) and was approved by the Western Institutional Review Board (WIRB) and Saint John's Health Center/John Wayne Cancer Institute Institutional Review Board. All ethnicities and both genders were included. Only patients who were clinically diagnosed and pathologically confirmed to harbor GH-secreting, ACTH-secreting or endocrine-inactive pituitary adenomas (including silent-ACTH adenomas) were included in the study. Specifically, a diagnosis of acromegaly was reached on the basis of clinical manifestations of acromegaly or gigantism as well as age and gender-specific elevated insulin-like growth factor type I (IGF-1) along with unsuppressed GH levels (16). Similarly, a diagnosis of Cushing's disease was based on clinical signs and symptoms of Cushing's syndrome (which include obesity, diabetes mellitus, hypertension, moon facies, dorsocervical fat pad, supraclavicular fat pad, abdominal striae, and/or easy bruising) and a combination of elevated urine free cortisol (24-hours collection), elevated midnight salivary cortisol levels, elevated or inappropriately normal serum ACTH levels, and lack of suppression of serum cortisol following a dexamethasone suppression test (17). All patients also had MRI-evidence of a pituitary adenoma on high-quality sellar MRI with and without gadolinium. Both GH and ACTH-secreting pituitary adenomas were confirmed by immunohistochemical (IHC) staining. Endocrine-inactive pituitary adenomas were macroadenomas without clinical evidence of hormone secretion. Silent-ACTH tumors stained for ACTH, but lacked clinical signs, symptoms or biochemical data supportive of Cushing's disease.

Patients age <18 or >90 years, patients with prolactin- or thyrotrophic hormone (TSH)-secreting adenomas, or patients who were clinically diagnosed as Cushing's disease or acromegaly without pathological evidence of adenoma, patients with familial pituitary adenomas, including Multiple Endocrine Neoplasm Type1 (MEN1) syndrome, Carney complex, familial isolated pituitary adenoma (FIPA), and McCune-Albright syndrome, were excluded from this study. Patients with radiation therapy before surgery were excluded as well.

Patients' clinical records, imaging studies, hormonal studies, and pathological reports were reviewed. All DICOM data were loaded into OsiriX (Pixmeo) to calculate the adenoma maximal dimension. Follow-up MRI images and clinical data were obtained from the medical records. Follow-up information in some patients was obtained via email or telephone.

Forty-eight patients (57% male, 43% female, mean age 52 ± 11 years, range: 21–70 years) were enrolled in the study between August 2008 and September 2014. There were 18 patients with endocrine-inactive adenomas, 17 with GH-secreting adenomas (acromegaly), and 13 with ACTH-secreting adenomas (Cushing's disease, including three silent-corticotroph adenomas, one of which was also a Crooke's cell adenoma). Of the 48 patients, 46 (96%) were typical adenomas and two (4%) were atypical (one endocrine-inactive and one GH-secreting). Ki-67 index was estimated at less than 3% in 33 patients (69%), 3%–6% in 14 patients (29%), and greater than 6% in one patient (2%). The mean maximal adenoma dimension was 19 ± 10 mm (range: 5–44 mm; Supplementary Table S1). Seven patients had prior history of pituitary adenoma resection. Post-operative radiation was not performed on any of these patients.

Complete follow-up was achieved in 42 patients (85%), with a median clinical follow-up of 24 months (range, 2–76 months). Adenoma recurrence rate was 0% at 12 months after surgery. Disease control was achieved in 41 of these patients (98%), according to post-operative MRI and endocrinological evaluations. Only one patient with invasive GH-secreting adenoma underwent repeat endoscopic surgery for further adenoma debulking 10 months after the first surgery due to persistent disease. Of the 17 patients with acromegaly, four (25%) were treated with at least three doses of octreotide or lanreotide before surgery, whereas the other patients received no therapy pre-operatively.

Blood specimens collected in sodium citrate tubes were immediately processed for the isolation of peripheral blood lymphocytes (PBL) with Purescript RBC Lysis Solution (Gentra Systems). Genomic DNA was extracted from PBL using DNAzol BD Reagent (Life Technologies) or the Quick-gDNAMiniPrep Kit (Zymo Research) with RNase A (Qiagen), according to the manufacturer's recommendations. When blood specimens were not available, genomic DNA was isolated from 500 μL of saliva collected in the ORAgene Discover Kit (DNA Genotek, Inc.) using prepIT-L2P reagent (DNA Genotek, Inc.), following the manufacturer's recommendations for DNA purification. DNA isolation from formalin-fixed paraffin-embedded (FFPE) tissue was carried out as previously described (18). Briefly, deparaffinization was performed using standard xylene and ethanol rinse, and DNA isolation was then carried out using the ZR FFPE DNA MiniPrep Kit (Zymo Research).

RNA isolation from FFPE was performed as previously described (19). Briefly, deparaffinization was performed using standard xylene and ethanol rinse, and then digested by proteinase K. RNA was extracted using RNA-Solv Reagent (Omega Bio-Tek) with chloroform, followed by overnight precipitation in isopropanol and the addition of Pellet Paint NF (Novagen).

Frozen tissue was lysed using ZR BashingBead Lysis Tubes with the Xpedition Sample Processor (Zymo Research) in conjunction with the QIAshredder tubes (Qiagen). DNA and RNA isolation from frozen tissue was performed using either the ZR-Duet DNA/RNA MiniPrep Kit or the ZR Tissue & Insect DNA MiniPrep Kit (Zymo Research), according to the manufacturers' protocols.

Quantitation of double-stranded DNA (dsDNA) was performed using the Quant-iT PicoGreen dsDNA Assay (Life Technologies) or using the Qubit dsDNA HS Assay (Life Technologies). Quality of DNA was determined using Agilent 2200 TapeStation Genomic DNA Analysis (Agilent Technologies). Total RNA was quantified using the Quant-iT RiboGreen dsDNA Assay (Life Technologies), as previously described (18). RNA quality was assessed using Agilent 2200 TapeStation RNA Analysis (Agilent Technologies) to determine RNA Integrity Number (RIN) or DV200, following the manufacturer's protocol.

To assess the genetic changes between the three disease groups, whole-exome sequencing (WES) was performed on 39 frozen tissue adenoma specimens and 7 FFPE adenoma specimens and 41 peripheral blood leukocyte (PBL) specimens. Five saliva specimens were also sequenced as a paired germline control. Sequencing libraries were prepared from genomic DNA using the Agilent SureSelect XT All Exon V5 Kit (Agilent Technologies), following the manufacturer's library preparation protocol for 200 ng of DNA input. For FFPE samples, we followed the 200 ng DNA input protocol with the following modifications: (1) input of 1 to 3 μg of DNA and (2) 1:3 dilution of the SureSelect Adaptor Oligo mix was used during the paired-end adapter ligation step. The normalized library pools were sequenced on the Illumina HiSeq 2500 (Illumina Inc.) in rapid mode with 100 bp paired-end reads with a loading concentration of 8.5 pmol/L.

Because of the input requirements of the different DNA and RNA platforms used in this study it was not always feasible to collect both nucleic acid types from the exact same extraction prep in all patients. In cases where nucleic acid yields were insufficient, additional tissue was sampled from the same piece of tumor.

For each sample, 150 ng of total RNA was used for library construction with the TruSeq RNA access library prep kit, following the manufacturer's protocol (Illumina Inc.). The normalized library pools were sequenced on the Illumina HiSeq 2500 in rapid mode with 100 bp paired-end reads with a loading concentration of 8.5 pmol/L.

The methylation profiles of the adenoma specimens were determined using the Infinium HumanMethylation450 BeadChip Array (HM450K; Illumina Inc.) following the manufacturer's protocol. Sodium bisulfite modification (SBM) was performed on 1 μg of DNA using Zymo EZ DNA methylation kit (Zymo Research). After analyzing an aliquot of the resulting SBM-DNA with MethyLight-based quality control to test bisulfite completeness, the SBM-DNA was whole-genome amplified and enzymatically fragmented. For FFPE DNA, the SBM-DNA was treated with the Infinium HD FFPE DNA Restore Kit (Illumina Inc.) before whole-genome amplification. Finally, the amplified and fragmented SBM-DNA was hybridized to the HM450K BeadChip. The chips were scanned with Illumina iScan.

Raw genomic sequence reads were mapped to the 1000 Genomes (b37) build of the human genome reference using BWA-MEM (version 0.7.5a) with default settings (20). The resulting alignments were further processed using the GATK (version 2.8-1) base quality score recalibration, indel realignment, duplicate removal (picardtools, version 1.103) following the GATK Best Practices recommendations (21). The picardtools CollectOxoMetrics tool was used to screen alignments for evidence of excessive oxidative damage.

Somatic single nucleotide variants (SNV) and small indels were identified using the MuTect2 algorithm (22). MuTect2 was run on matched tumor and normal alignments using default parameters along with dbSNP (version 137) and COSMIC (version 67) databases. The Mutect2 internal filters were used to remove suspected germline variants and possible errors due to sequencing artifacts (i.e., Oxo-G error). Only somatic mutation calls assigned the “PASS” filter were included in downstream analyses. The Bioconductor package Maftools (23) was used to summarize mutation data, estimate mutation load in comparison with the TCGA cohorts, and generate lollipop plots of mutated genes. Suspected driver genes were identified using the algorithm OncodriverCLUST (24) as implemented in Maftools. We identified overall potential drive mutations in all pituitary adenoma samples, excluding the three hypermutated individuals, using an FDR threshold of 0.05 and a minimum number of mutations per gene of five. In addition, we identified driver mutations in each of the three pituitary adenoma subtypes by running OncodriveCLUST on each subtype separately using an FDR threshold of 0.05 and setting the minimum number of mutations per gene to three.

Somatic copy-number alterations (SCNA) were identified using CNVkit (version 0.9; ref. 25) to analyze the read coverage and segment changes of the tumor genome against a pooled set of all pituitary adenoma normal exome samples following the authors recommended settings for “Tumor Analysis.” Because we used a pooled reference with mixed sexes, the XY chromosomes were excluded in the analysis. The genomic range of gain, loss and total genomic fraction were calculated from the sum of segment length with the segment log2 ratio of >0.2, <−0.25, and both, respectively. Segments with extremely low log₂ ratios (i.e., below −15) were removed as artifacts. To detect recurrent regions, GISTIC2 (version 6.10; ref. 26) was run using the segment results from CNVkit using default settings, except for the confidence level threshold which was set to 0.99. GISTIC2 was run for six different groupings; All pituitary adenoma samples, only GH-secreting pituitary adenomas, only ACTH-secreting pituitary adenomas, only endocrine-inactive pituitary adenomas, invasive versus noninvasive pituitary adenomas, and recurrent versus non-recurrent pituitary adenomas to identify specific regions common to all pituitary adenomas regardless of subtype, regions specific to a pituitary adenoma subtype, or regions specific to a particular clinical characteristic.

Raw RNA sequencing reads were first checked for overall quality and filtered for adapter contamination using Trimmomatic (version 0.36; ref. 27). The filtered reads were then mapped to the GENCODE comprehensive gene annotation set version 19 using STAR (version 2.4.2a; ref. 28) with default parameters. Read counts for each feature were generated using the “–quantModeGeneCounts” function in STAR.

Significantly differentially expressed genes were identified using the DESeq2 (29) Bioconductor package. We tested for differentially expressed genes between pairs of pituitary adenoma subtypes using a model that accounted for tissue source (i.e., frozen or FFPE) and sex (design ∼ tissue + sex + pituitary adenoma subtype). The list of significantly differentially expressed genes was further filtered using a log fold change threshold of 1 and an FDR threshold of 0.1. The filtered set of genes was used for downstream analyses. Pathway analysis was performed using the Reactome DB Bioconductor package (30) with an adjusted significance level of 0.05.

The ChAMP Bioconductor package was used to process the raw methylation array data starting from the idat files generated by the array platform. Following the recommendations in the ChAMP manual, we filtered the raw data to remove probes that: had a detection P < 0.01, probes with < three beads in 5% or more of the samples, were non-CpG probes, probes associated with a known SNP, probes that are known to bind to multiple genomic locations, and probes located on sex chromosomes. Normalization was performed on the filtered probe set using the beta-mixture quantile normalization (BMIQ) method (31) and ComBat (32) was used to adjust for the effects of array slide and sex. Differentially methylated regions (DMR) were identified using the DMRcate method as described in (33) and regions with an adjusted P < 0.1 were considered significant and included in downstream analysis. We used generalized linear models to test for differences among the pituitary adenoma subtypes for each genomic annotation feature.

To better understand the relationship between methylation state and gene expression, we compared the top 1% most significantly differentially expressed genes for both the GH-secreting and ACTH-secreting pituitary adenomas compared with endocrine-inactive pituitary adenomas to generate “starburst” plots (34). In addition, we tested for the relationship between gene promoter methylation and gene expression in the context of protein interaction networks suing the FEM Bioconductor package (35). The FEM analysis was run with default settings and using the normalized RNA-seq read counts and the batch corrected and normalized methylation beta values generated from the analyses described above. Modules with an adjusted P < 0.05 were considered significant.

Multiplex immunohistochemistry (mIHC) was performed on pituitary adenoma samples. Tissue sections were cut at 4 μm from FFPE blocks. All the sections were deparaffinized and subjected to heat-induced epitope retrieval in citrate buffer pH 9.0 (Biogenex). A six-plex panel IHC was performed for each tissue slide using the following antibodies (Ab): anti-FoxP3 Ab (clone 236A/E7, dilution 1:100, Abcam), anti–PD-L1 Ab (clone E1L3N, dilution 1:250, Cell Signaling Technology), anti-CD8 (clone SP16, dilution 1:50, Spring Bioscience), anti-CD3 AB (clone SP7, dilution 1:50, Spring Bioscience), anti-CD163 Ab (clone MRQ26, Ventana). Antigen-antibody binding was visualized with TSA-Cy5 (PerkinElmer), TSA-Cy3 (PerkinElmer), TSA-FITC (PerkinElmer), TSA-Alexa594 (Life Technologies), and TSA-Cy5.5 (PerkinElmer), respectively. Microwave treatment in citrate buffer pH 6.0 was performed between Ab detection to prevent cross-reactivity. Tissue slides were incubated with DAPI as counterstain and cover slipped with VectaShield mounting media (Vector Labs). Control tissue samples were stained for each different marker. Digital images were captured with the PerkinElmer Vectra platform. Representative tumor areas were scanned at ×20 and selected for analysis with InForm Software (PerkinElmer). Hematoxylin and eosin staining was performed on each sample and reviewed by a pathologist to ensure the quality of the tissue sample.

Using whole-exome sequencing, we characterized the landscape of somatic single-nucleotide variants across the exomes of 48 pituitary adenomas. Among the 48 pituitary adenoma samples, the median number of variants per exome was 14. However, three patient specimens, two GH-secreting pituitary adenomas and one endocrine-inactive pituitary adenoma, with unusually high mutation loads of 11,355, 2,264, and 3,215 (synonymous and non-synonymous) variants per exome were identified. We screened these three hypermutated samples (GH.1, GH.4, and Inactive.21) for somatic mutations and germline variants in DNA mismatch repair (MMR) related genes. We first focused on somatic variants related to MMR, including the genes MSH2, MSH3, MSH4, MSH5, MSH6, MLH1, MLH3, PMS1, PMS2, and PMS2L3. Of the three samples, only Inactive.21 harbored somatic mutations in MMR-related genes, with four missense mutations in the genes MSH5, MSH6, MLH1, and MLH3. We also looked for somatic variants in polymerase genes involved in double strand brake repair. In sample GH.1, we identified several missense mutations in polymerase genes, including POLE4, POLQ, and a frame shifting deletion in POLQ. Likewise, we observed a POLE missense mutation in sample Inactive.21. We did not find any MMR, polymerase, or BRCA somatic or germline variants in sample GH.4. These results suggest that subset of pituitary adenomas (∼6% in the current study) may display a previously unidentified hypermutation phenotype.

We tested for presence of driver mutations across all 45 non-hypermutated pituitary adenomas using the OncodriveCLUST method (see Materials and Methods) and identified only a single candidate in the gene GNAS (P = 2.29e⁻⁰⁸), with no other genes reaching significance. We further stratified mutations by pituitary adenoma subtype and tested for drive mutations. Within GH-secreting pituitary adenomas only GNAS reached significance (P = 1.40e⁻⁰⁸). In ACTH-secreting pituitary adenomas, recurrent mutations were identified in the gene USP8 (P = 2.00e⁻⁰⁵) and the ploy(A)–binding cytoplasmic protein 1 (PABCP1; P = 1.43e⁻⁰³) gene. The two recurrent mutations in the genes GNAS (50%) and USP8 (23%) in GH-secreting and ACTH-secreting samples, respectively, showed pituitary adenoma subtype frequencies in this study similar to previous reports (4, 8, 36). There were no significant recurrent mutations identified in the endocrine-inactive pituitary adenoma subtype.

Whole-exome sequencing data were used to identify SCNAs across all 48 pituitary adenoma samples. Overall, the proportion of the genome that was disrupted by SCNAs ranged from approximately 0% to 87% of the genome. Both highly disrupted and non-disrupted pituitary adenomas where observed among all three pituitary adenoma subtypes (Fig. 1), with GH-secreting pituitary adenomas showing a significantly (P = 0.034) higher portion of the genome disrupted, followed by inactive pituitary adenomas and ACTH-secreting pituitary adenomas. Interestingly, the three corticotroph ACTH-secreting pituitary adenomas (Fig. 1. orange highlight) had a greater fraction of the genome disrupted than did any of the other ACTH-secreting pituitary adenomas.

We then identified chromosome arm level gains and losses among all 48 pituitary adenomas. Across the 48 pituitary adenomas, chromosome arm level losses were observed more frequently than gains, with 12 losses and one gain (P = 0.003). Significant losses were observed on chromosome arms 1p, 1q, 2q, 3q, 4p, 4q, 6q, 8q, 12p, 13q, 18q, and 21p. Only a single gain on chromosome arm 19q was detected among all 48 pituitary adenomas.

To identify pituitary adenoma subtype-specific alterations, we tested for focal level SCNAs across all pituitary adenomas together as well as within each pituitary adenoma subtype independently. In addition, we compared pituitary adenomas that were defined clinically (see Materials and Methods) as invasive versus non-invasive to identify SCNAs related to overall invasiveness. Supplementary Table S2 summarizes the top 10 most significant focal level SCNAs identified across all pituitary adenomas, within a specific pituitary adenoma subtype, and between invasive and non-invasive pituitary adenomas. Across all pituitary adenomas, the most significant loss was found in the region 12q13.2 and the most significant gain in region 12q13.3. Among the three pituitary adenoma subtypes, GH-secreting pituitary adenomas had a subtype-specific loss at 12q21.1 that was not observed in either the overall analysis or in the other two pituitary adenoma subtype analyses and ACTH-secreting pituitary adenomas were found to have a subtype-specific loss at 12p13.31 that is unique to this pituitary adenoma subtype. In addition, invasive pituitary adenomas were found to have a unique SCNA at 3p22.3, while non-invasive pituitary adenomas had a unique loss in the region 12q21.31.

The impacts of epigenetic alterations on the disease biology of pituitary adenomas was investigated using the HM450K array to characterize the genome-wide DNA methylation profiles of 37 pituitary adenomas. Unsupervised hierarchical clustering of the top 5,000 most variable CpG sites revealed three distinct methylation profiles that clearly separated each pituitary adenoma subtype (Fig. 2). Notably, the three silent-corticotroph adenomas were found to cluster within the functional ACTH-secreting pituitary adenomas, highlighting the overall similarity of the genome-wide DNA methylation profiles between silent-corticotroph and ACTH-secreting adenomas and reflecting their more recalcitrant clinical behavior. Overall, GH-secreting samples showed a DNA methylation profile dominated by an abundance of hypomethylated regions, as compared with ACTH-secreting and endocrine-inactive pituitary adenoma samples (Fig. 3A). To check whether the genomic load of SCNAs was influencing the overall methylation profiles, we performed a PCA analysis using all probes in the methylation data set. Again, we see three distinctive clusters of samples that correspond directly to pituitary adenoma subtype and do not observe any pattern with relation to the percentage of the genome altered by SCNAs (Supplementary Fig. S1).

We then compared the genome-wide DNA methylation profiles for specific genomic features across the three pituitary adenoma subtypes. DNA methylation levels were significantly different (P < 2.0e⁻⁴) among the three pituitary adenoma types for CpG sites in the 1^st exon, gene body, and intergenic regions (IGR). However, methylation levels were significantly different (P < 2.0e⁻⁴) only for GH-secreting pituitary adenomas compared and not between ACTH-secreting and endocrine-inactive pituitary adenomas for CpG sites in the promoter region (Fig. 3B). In addition, significant differences (P < 2.0e⁻⁴) were observed among the three pituitary adenoma subtypes in the DNA methylation levels of CpG sites in low-CpG density regions (known as CpG shore, shelves and open sea), but only for GH-secreting pituitary adenomas at high-CpG density regions (known as CpG islands; Fig. 3C).

To evaluate the unique epigenetic changes related to GH-secreting and ACTH-secreting pituitary adenomas (including the three silent-corticotroph pituitary adenomas), we used the set of endocrine-inactive pituitary adenoma samples as a baseline DNA methylation profile. Overall, the number of significantly differentially methylated CpG sites (FDR < 0.05) was greater in GH-secreting pituitary adenomas with 169,487 differentially methylated sites, compared with 55,860 sites in ACTH-secreting specimens. The differentially methylated sites were further classified into differentially methylated regions (DMR) based on the DMRcate method (33). Again, a greater number of DMRs was identified in the GH-secreting samples (4,127) than in the ACTH-secreting specimens (1,307).

To support the translational relevance of the DNA methylation analysis, transcriptomic differences were investigated among pituitary adenomas using whole-transcriptome sequencing of 38 pituitary adenomas. Similar to the pattern observed in the DNA methylation analysis above, an unsupervised hierarchical clustering of the top 500 most variable genes revealed three distinct gene-expression profiles corresponding to pituitary adenoma subtype (Fig. 4). Consistent with the DNA methylation results, there were a greater number of differentially expressed genes in the GH-secreting versus endocrine-inactive comparison than in the ACTH-secreting verses endocrine-inactive comparison, with 4,546 and 3,923 significantly differentially expressed genes at FDR < 0.05, respectively. Consistent with the observation of the methylation profiling, the three silent-corticotroph pituitary adenomas clustered with the ACTH-secreting pituitary adenomas, further reinforcing the overall similarities in the biology of these two groups of pituitary adenoma subtypes. Finally, we performed a PCA analysis using the normalized expression data for all genes to look for patterns of SCNA load among the pituitary adenoma samples. Again, we observed three separate clusters that corresponded to pituitary adenoma subtype and did not reflect the percentage of the genome altered by SCNAs (Supplementary Fig. S2).

Next, the relation between DNA methylation levels and gene expression patterns was investigated among the three pituitary adenoma subtypes for 37 samples that had matched methylation and expression data. Figure 5 shows the resulting “starburst” plots for both GH-secreting and ACTH-secreting pituitary adenomas (Fig. 5A and B). GH-secreting pituitary adenomas showed the greatest number of differentially expressed genes with a predominant skew toward the hypomethylated state, whereas ACTH-secreting showed overall fewer differentially expressed genes with no apparent skew toward methylation state as compared with the endocrine-inactive pituitary adenomas. Supplementary Tables S3 and S4 show the top 50 most significantly differentially expressed genes contained within Fig. 5 for GH-secreting pituitary adenomas and ACTH-secreting pituitary adenomas, respectively.

A closer examination of the starburst plot for GH-secreting compared with the endocrine-inactive samples shows that the elevated expression of Somatostatin Receptor 5 (SSTR5), Growth Hormone 1 (GH1), and Growth Hormone 2 (GH2) in GH-secreting pituitary adenomas results from the hypomethylation of the promoters of these genes. For the ACTH-secreting and endocrine-inactive pituitary adenoma comparison, we observed overexpression of proopiomelanocortin (POMC) influenced by a hypomethylated promoter. POMC is the precursor protein of ACTH, the pituitary hormone that causes increased serum cortisol levels in Cushing's disease.

To better understand how the relationship between epigenetic and transcriptional changes work in concert, we integrated DNA methylation and gene-expression profiles to identify functional epigenetic modules (FEM; ref. 35) in GH-secreting and ACTH-secreting pituitary adenomas. This analysis unraveled interaction networks centered on GH1 significantly affected (FDR < 0.05) in GH-secreting pituitary adenomas and on POMC significantly affected (FDR < 0.05) in ACTH-secreting pituitary adenomas (Fig. 6A and B).

Hypophysitis has emerged as a distinctive endocrine adverse effect of anti-CTLA4 and anti–PD-1 therapies (37). Therefore, we examined the expression of PD-L1 among all three pituitary adenoma subtypes using multiplexed immunohistochemistry (mIHC) staining (38). High levels of PD-L1 expression were found among all three pituitary adenoma subtypes, and although GH-secreting pituitary adenomas showed higher levels of PD-L1 expression, followed by endocrine-inactive adenomas then ACTH-secreting samples, the difference among pituitary adenoma subtypes was not significant (P = 0.113). Similarly, we did not detect any significant differential expression of PD-L1 among the three pituitary adenoma subtypes in the RNA-seq data. Across all three pituitary adenoma subtypes, PD-L1 expression was heterogeneous among the tissue specimens, with endocrine-inactive adenoma samples showing only small pockets of PD-L1 expression (data not shown). The presence of immune infiltrates was observed in all three pituitary adenoma subtypes which coincided with the PD-L1 expression (data not shown).

To investigate the molecular underpinnings of tumor recurrence in pituitary adenomas we identified genes that were differentially expressed in primary tumors taken from patients that later experienced a recurrence. Of the patients in our cohort with RNA-seq data, there were seven recurrent and 23 non-recurrent (see Supplementary Table S1). By contrasting these two groups, we identified 68 genes that were significantly differentially expressed. Pathway enrichment analysis of this set of genes showed significant (P = 3.8e⁻⁰³) enrichment for the Reactome pathway chemokine receptors bind chemokines (R-HSA-380108) containing genes involved in chemokine receptor binding. In particular, the chemokine related genes IL8, CXCR1 and CXCR2 were enriched in this pathway and are known to drive aggressiveness in other tumor types (39).

Finally, the influence of pre-operative somatostatin analogue treatment was assessed on GH-secreting adenomas in patients with acromegaly by comparing the methylation and gene expression profiles of three adenoma specimens treated with octreotide or lanreotide to the untreated samples. In all, 95 genes were found to be differentially expressed between the two groups at a significance level of FDR ≤ 0.1. Interestingly, a significant decrease in expression of marker of proliferation Ki-67 (MKI67) was observed in somatostatin analogue-treated samples indicative of reduced tumor proliferation. An increase in expression was also observed among treated samples of the gene Mucin 1 (MUC1), which is known to play a crucial role in cancer development (40), and increased expression may contribute to treatment efficacy in pituitary adenomas (41). In addition, B-cell surface antigen CD40 (CD40), a member of the tumor necrosis factor receptor (TNF-R) family that modulates B-cell activation and differentiation, was significantly upregulated in the somatostatin analogue-treated pituitary adenomas.

In this study, we show that the three major subtypes of surgically resectable pituitary adenomas display tumor-specific genome-wide patterns of DNA methylation and gene expression and that these profiles can be used to differentiate among each specific pituitary adenoma subtype. We also demonstrate that although pituitary adenomas harbor some recurrent somatic mutations, DNA methylation alterations play a major role in the disease etiology.

Our DNA methylation analysis demonstrates well-separated clustering of endocrine-inactive adenomas, GH-secreting adenomas and ACTH-secreting/silent-corticotroph adenomas. GH-secreting adenomas demonstrated differentially hypomethylated genes, particularly in comparison with ACTH-secreting adenomas, suggesting that genetic targets may be more accessible for the treatment of acromegaly. Interestingly, the silent-corticotroph adenomas clearly and consistently clustered with the ACTH secreting tumors rather than the endocrine-inactive adenomas, which is how they are grouped clinically. Silent-corticotroph adenomas tend to be more invasive and have a higher recurrence rate requiring re-operation and radiotherapy (42). The similar genetic patterns between Cushing's ACTH-secreting adenomas and silent-corticotroph adenomas suggest that similar pharmacotherapeutic strategies may be effective for both of these tumor subtypes.

We also demonstrate that promoter hypomethylation drives overexpression of SSTR5 and GH2 in GH-secreting and POMC in ACTH-secreting pituitary adenomas even in the absence of USP8 mutations. Silent-ACTH adenomas have a higher recurrence rate, similar to ACTH tumors causing Cushing's disease (43). SSTR5 is the predominant somatostatin receptor subtype that is identified in GH-secreting adenoma cells (including the GH2 variant). This finding supports the increased efficacy demonstrated in acromegaly treatment with pasireotide (high affinity to SSTR5) compared with octreotide and other somatostatin analogues with predominant affinity for SSTR2 (44). Interestingly, variable SSTR5 expression or hypomethylation was not identified in the Cushing's disease cohort when compared with the endocrine-inactive group, despite efficacy of pasireotide therapy for these patients. SSTRs are traditional targets for somatostatin targeted therapy. However, these drugs are effective in only about 30% of patients with refractory acromegaly or Cushing's disease (45, 46). Hence, additional targeted therapy is necessary to achieve better biochemical control and disease remission in these patients.

Separately, the overexpression of POMC has been demonstrated to have an indirect effect of USP8 mutations on EGFR signaling (4, 8). However, in our DNA methylation and gene-expression datasets, only 3 of 13 ACTH-secreting samples had a detectable USP8 mutation. When we compared both POMC expression levels and POMC promoter DNA methylation levels between USP8 mutant and wild type, there was no significant difference, suggesting that the DNA methylation state of the POMC promoter drives expression regardless of USP8 mutation status. Therefore, the functional relevance of the low level of USP8 mutation in pituitary adenoma diseases may be questionable. Similarly, we observed no significant difference in the expression of GNAS between GNAS mutated and wildtype GH-secreting adenomas. This further illustrates that the expression of key driver genes may be independent of mutation status and warrants further investigation. In addition, the gene SLIT1 was highly expressed with a hypomethylated gene body. Interestingly, SLIT1 is involved in axon guidance, cell proliferation, cell motility and angiogenesis pathways and its regulation is known to be influenced by DNA methylation in other cancers (47).

Recent advancements in immune check-point inhibition, including anti–PD-1/PD-L1 therapies, have shown strong antitumor activity in several cancers (48). Thus, PD-1/PD-L1 treatments directed toward pituitary adenomas may provide a novel therapeutic strategy for these common diseases. In addition, these results also indicate that T-cell–mediated inflammation could explain the widely reported occurrence of hypophysitis autoimmune responses related to checkpoint immune inhibitor therapies such as anti-CTLA4 and anti-PD1 antibody treatment (49). Furthermore, endocrine-active pituitary adenomas may be more susceptible to immune reactivity in its various disease states and related to autoimmunity activity. The interaction of CD40L with CD40-expressing cancer cells can directly inhibit the growth of various carcinoma cells through cell-cycle blockage and/or apoptotic induction (50), and thus may be indicative of increased apoptosis in somatostatin analogue-treated acromegaly pituitary adenomas. Our study suggests that patients with refractory acromegaly may be responsive to anti–PD-1 or anti–PD-L1 therapy based on both PD-L1 expression and immune cell infiltration data. Prospective clinical studies investigating this possibility will be necessary to assess efficacy.

B.A. Fox reports receiving commercial research grants from and is a consultant/advisory board member for PerkinElmer. G. Barkhoudarian is a consultant/advisory board member for Vascular Technologies Inc. No potential conflicts of interest were disclosed by the other authors.

Conception and design: B.A. Fox, G. Barkhoudarian, D.F. Kelly, D.S.B. Hoon

Development of methodology: N. Nelson, C. Ballesteros-Merino, D.F. Kelly, D.S.B. Hoon

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): X. Wang, S.C. Hsu, N. Nelson, X. Zhang, G. Barkhoudarian, D.S.B. Hoon

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): M.P. Salomon, X. Wang, D.M. Marzese, S.C. Hsu, X. Zhang, C. Matsuba, C. Ballesteros-Merino, B.A. Fox, G. Barkhoudarian, D.S.B. Hoon

Writing, review, and/or revision of the manuscript: M.P. Salomon, X. Wang, D.M. Marzese, S.C. Hsu, N. Nelson, B.A. Fox, G. Barkhoudarian, D.F. Kelly, D.S.B. Hoon

Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): S.C. Hsu, N. Nelson, B.A. Fox, D.S.B. Hoon

Study supervision: B.A. Fox, G. Barkhoudarian, D.F. Kelly, D.S.B. Hoon

Other (pathology assessment): Y. Takasumi

We would like to thank Dr. Ian Hutchinson for help with editing the article. Also, we would like to thank Mike Stegner and Adam Morris for IT support and Annie Heng, Amy Eisenberg, Charlene Welker, Shanthy Chenathukattil, Ling Takeshima, and Garret Cheung for their support. This study was supported in part by a grant from Novartis (CSMS995BUS66T; to D.S.B. Hoon and D.F. Kelly) for the pituitary sequencing studies and the Department of Translational Molecular Medicine at the John Wayne Cancer Institute.

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.