Gene expression profiling in human corticotroph tumours reveals distinct, neuroendocrine profiles.

Adrenocorticotrophic hormone (ACTH)-secreting pituitary adenomas give rise to a severe endocrinological disorder, comprising Cushing's disease, with multifaceted clinical presentation and treatment outcomes. Experimental studies suggest that the disease variability is inherent to the pituitary tumour, thus indicating the need for further studies into tumour biology. The present study evaluated transcriptome expression pattern in a large series of ACTH-secreting pituitary adenoma specimens in order to identify molecular signatures of these tumours. Gene expression profiling of formalin-fixed, paraffin-embedded specimens from 40 human ACTH-secreting pituitary adenomas revealed the significant expression of genes involved in protein biosynthesis and ribosomal function, in keeping with the neuroendocrine cell profile. Unsupervised cluster analysis identified 3 distinct gene profile clusters and several genes were uniquely overexpressed in a given cluster, accounting for different molecular signatures. Of note, gene expression profiles were associated with clinical features, such as the age and size of the tumour. Altogether, the findings of the present study show that corticotroph tumours are characterised by a neuroendocrine gene expression profile and present subgroup-specific molecular features.

INTRODUCTION

Cushing's disease, that is adrenocorticotrophic hormone (ACTH)‐secreting pituitary adenoma, is a rare and complex endocrine disease. Hypersecretion of ACTH by the tumour gives rise to excess cortisol secretion by the adrenal glands, which, in turn, promotes a variety of clinical signs encompassing osteoporosis, muscle atrophy, hypertension and diabetes. These patients carry a considerable disease burden, and mortality, if untreated, is high.1 Studies into the pathophysiology of corticotroph tumours have shed light on the involvement of individual factors, such as epidermal growth factor (EGF),2 ubiquitin‐specific peptidase 8 (USP8),3, 4 Brahma‐related gene 1 (Brg‐1),5 although little is known about other, possible alterations in these tumours.

Of note, an increasing body of evidence indicates that ACTH‐secreting adenomas differ markedly among themselves regarding their responses to the main modulators of corticotroph secretion, such as corticotrophin‐releasing factor and steroids,6, 7 as well as ancillary modulators, such as vasopressin analogues, gonadotrophin‐releasing hormone.8, 9 Differences in in vitro ACTH secretory patterns have also been observed in response to dopamine or somatostatin receptor agonists,10, 11 which obviously translates into different responsiveness to these drugs. Furthermore, somatic mutations and individual polymorphisms have been linked to differences in tumour size and secretory pattern.12, 13

It is therefore apparent that wide‐range evaluation of gene expression patterns in corticotroph adenomas is needed to understand which features characterise tumoral corticotrophs and underlie the phenotypical diversity and thus, possibly, pave the way for a personalised clinical strategy. The present study evaluated gene expression profiles in a large series of human ACTH‐secreting adenomas. We used RNA extracted from formalin‐fixed, paraffin‐embedded archival samples collected during transsphenoidal surgery, thus overcoming limitations as a result of the rarity of the disorder.

MATERIALS AND METHODS

Sample collection

Forty formalin‐fixed, paraffin‐embedded adenomatous samples were retrieved and sections were microdissected under a surgical microscope. Specimens were collected from the centre of the biopsy to minimise RNA degradation as a result of nucleic acid oxidation. Immunohistochemistry and neuropathological examination confirmed the diagnosis and the absence of normal pituitary tissue on sections. MIB‐1 staining was <3%, attesting to typical corticotroph adenoma (grade I). Overall, specimens had been stored 2‐15 years. Clinical data from patient presentation, surgical outcome and follow‐up were collected. The entire series comprised 9 men and 31 women, aged 17–69 years, 29 microadenomas (diameter ≤10 mm) and 11 macroadenomas (diameter >10 mm), 6 surgical failures and the remainder surgical remissions. Mean ± SD follow‐after surgery up was 140.3 ± 13.6 months (range 2‐15 years). The study was approved by the Ethical Committee of the Istituto Auxologico Italiano.

RNA formalin‐fixed, paraffin‐embedded sample extraction

Four 20‐μm thick sections from each specimen were extracted using Recover All Total Nucleic Acid Isolation Kit (Life Technologies, Carlsbad, CA, USA). Tissue blocks were deparaffinised in xylene, RNA‐bound proteins removed by protease, and RNA captured on glass‐filters and eluted. Recovery of RNA averaged 1.5‐2 μg per sample. To maximise recovery of adequate transcripts and reduce RNA damage through oxydation, slices were cut from the center of the paraffin block.

RNA quality control

RNA (200 ng) was reverse‐transcribed (Superscript‐Vilo cDNA synthesis kit; Life Technologies) and amplified by real‐time polymerase chain reaction (PCR) for ribosomal protein L13A (RPL13A),14 using Taqman probe Hs03043885_g, (Applied Biosystems, Foster City, CA, USA). All 40 samples yielded the expected 81‐bp transcript at <30 cycles, attesting to adequate quality RNA.

Whole genome‐DASL HT

The ILLUMINA Whole Genome DASL High‐Throughput assay (WG‐DASL‐HT; Illumina, San Diego, CA, USA) for RNA extracted from formalin‐fixed, paraffin‐embedded tissues was used. The Human HT_12 v4 Bead Chip evaluates 29 285 probes (20 815 genes) and is highly accurate with respect to estimating differences in gene expression by head‐to‐head comparison with quantitative real‐time PCR.15 Each sample carried 300 ng of RNA. Fluorescence was imaged and data captured into HiScan, a high‐resolution laser imager (Illumina). Microarray data have been deposited with Gene Expression Omnibus (accession number GSE93825).

USP8 sequencing

Approximately 100 ng of RNA was reverse‐transcribed (Superscript VILO cDNA synthesis kit; Life Technologies) with the oligonucleotide primers: 5′‐CTTGACCCAATCACTGGAAC‐3′ (forward); 5′‐TTACTGTTGGCTTCCTCTTCTC‐3′ (reverse)13 for amplification of USP8 exon 14, the most frequent site of mutations reported so far.13 Touch‐down PCR was performed using GO TAQ DNA polymerase (Promega, Madison, WI, USA) at 64–57°C annealing. PCR products were purified by ExoProStar Illustra enzyme (Ge Healthcare, Chicago, IL, USA) and Sanger sequencing performed using the ABI PRISM Big DYE Terminator V3.1 cycle sequencing kit (Applied Biosystems) on ABI PRISM 3500 analyser.

Data analysis

Genome Studio software (Illumina) was used to detect genes significantly expressed in all specimens (ie, genes with detection P < 0.05 based on fluorescence intensity [average signal]), thus establishing a common expression profile. Subsequently, we attempted to identify different gene expression profiles in human corticotroph adenomas using unsupervised cluster analysis. This approach looks for similarities among gene expression profiles and identifies subgroups of samples on the basis of similarity between expression profiles.16, 17 Unsupervised methods for clustering among disease samples do not require prior knowledge on samples themselves or control tissues. Unsupervised cluster analysis was performed with hClust package (RStudio Inc., Boston, MA, USA) using the Ward minimum variance method18 and patients organised in accordance with their gene expression profile. Dissimilarity between clusters was measured by Euclidean distance.

Comparisons of clinical variables between clusters (ANOVA followed by Fisher's least significant difference or chi‐squared test, as appropriate) were performed with statview, version 4.5 (Abacus Concepts, Berkeley, CA, USA). P < 0.05 was considered statistically significant. Quantitative variables are given as the mean ± SEM.

Differential gene expression analysis

Once clusters had been established on the basis of similarity in gene expression profiles, differential gene expression analysis enabled the identification of genes that are expressed to a greater or lesser extent in each cluster. Clusters were analysed for differential expression analysis by rank invariant normalisation (Genome Studio; Illumina). Only genes with Benjamini and Hochberg P ≤ 0.05 were considered and Diff Scores were calculated based on P value transformation (ie, Diff Score ≥13 [up‐regulated gene] and ≤−13 [down‐regulated gene] P = 0.05) according to the difference between average signals in reference cluster vs the other 2 clusters. A Volcano plot was used to illustrate differential expression because it visualises the magnitude of over‐ and underexpression of a given gene plotted vs the significance of the comparison. This double‐filtering approach allows significantly differentially expressed genes to be easily identified.19

Real‐time gene expression analysis

RNA (100 ng) was reverse‐transcribed (Superscript‐Vilo cDNA synthesis kit; Life Technologies) and quantitative real‐time PCR was performed on a 7900 HT sequence Detection System (Applied Biosystems), using Platinum Quantitative PCR Supermix‐UDG with premixed ROX. Taqman assay (Applied Biosystems) was used for detection of the following genes: BCAS2 (probe Hs00903014_g1), SCG5 (probe Hs00161638_m1), IGFBP5 (probe Hs00181213_m1) and RPLP0 (probe Hs99999902_m1). Basal expression data () were calculated and normalised to RPLP0; differences between clusters were established by a Mann‐Whitney test.

Functional annotation and Gene Ontology

Two approaches were used to extract biological features from gene lists: david, version 6.720 was used to annotate and classify significant genes and perform functional annotation clustering. The minimum value of the enrichment score for significant clusters was ≤1.3. Clusters were annotated to Gene Ontology (GO) project, including Molecular Function (MF), Cellular Component (CC), Biological Process (BP), Kyoto Encyclopaedia of Genes and Genomes (KEGG) and Protein Information Resource (SP_PIR). Furthermore, cytoscape with Cluego plug‐in21 was used to identify molecular interaction networks. The biological pathways included were: GO project, KEGG and Inter Pro for protein domains.22 K score = 0.4 and Benjamini and Hochberg P value ≤0.001 were used for cytoscape analysis of significantly expressed genes. Of note, the group leading term is the most significant and genes can be included in several terms.22 The size of nodes reflects statistical significance of terms, whereas the degree of connectivity between terms (ie, edge) is calculated using Kappa statistics.23

RESULTS

Clustering analysis

Figure 1 shows unsupervised clustering and heat map profiling of significantly expressed genes in all 40 corticotroph adenoma samples. Of 20 815 genes, 1259 genes were expressed at P ≤ 0.05 and formed 5 distinct gene groups with decreasing expression. Of note, genes known to be associated with the corticotroph phenotype (eg, POMC, TBX19 [ie, Tpit]) were highly expressed (group 2). Conversely, the expression of genes associated with other anterior pituitary cells (for example, GH, PRL, PIT1, LHB, FSHB and CGA) was not significant across specimens, thus confirming absent contamination by normal anterior pituitary tissue.24, 25, 26 Highest enrichment scores were annotated to protein biosynthesis, ribosomes and RNA processing by GO, KEGG and SP‐PIR (see Supporting information, Table S1). In detail, the highest enrichment score was annotated to protein biosynthesis, followed by ribosomal function in group 1. The second‐ranked cluster revealed significant enrichment in ribosomal function, RNA processing, membrane transport and lumen, ATP metabolic processes and cytoskeleton organisation. Enrichment scores were lower in groups 3 and 4, mainly annotated to vesicle and lumen function, protein transport and localisation. Enrichment in group 5 was the lowest, with lumen and protease the only annotated functions (see Supporting information, Table S1).

Figure 1

Unsupervised clustering of 1259 significantly expressed genes. Colour coding of gene expression values is shown in the upper left corner. USP8 mutation status is indicated to the left of the heatmap: white squares identify wild‐type sequence, whereas coloured squares indicate USP8 variants: red, variant c.2159C>G; green, variant c.2152T>C; black, variant c.2155_2157delTCC; yellow, variant c.2157_2171delCCCAGATATAACCCA. Height of nodes represents dissimilarity between clusters as measured by Euclidean distance

cytoscape analysis revealed that significantly expressed genes formed consistent networks related primarily to protein targeting to the endoplasmic reticulum (ER), cellular macromolecule metabolic processes, RNA and polyRNA binding (Figure 2). Other networks were related to cellular transport of the protein to the ER and Golgi system and to membrane activity, such as vesicle forming and coating. Networks related to post‐transcriptional regulation of gene expression and the energy derivation by oxidation of organic compounds were also identified. Altogether, functional annotation by both approaches identified major enrichment of pathways related to protein synthesis and intracellular transport. Of note, this analysis identified genes that characterise the expression profile of corticotroph adenomas per se, and not genes uniquely expressed in tumoral corticotrophs.

Figure 2

cytoscape analysis showing networks formed by principal significant terms. The number of nodes is proportional to the number of the term‐forming genes; nodes are colour‐coded according to functional annotations and can have more than 1 colour because genes may belong to more than 1 term. The 2 major terms are shown in bold. COPI, coat protein I; ER, endoplasmic reticulum

Clinical data analysis

Gene expression profiles led to the grouping of patients in 3 separate clusters (ie, A, B and C). Analyses of differences in clinical variables among clusters revealed a higher proportion of macroadenomas and extrasellar tumours in cluster A (Table 1); indeed, 91% of macroadenomas and 92% of extrasellar tumours from the entire series fell in cluster A, with the remainder in cluster B. With respect to a quantitative comparison of tumour size, adenomas in cluster A were bigger than those in clusters B and C, although this failed to reach statistical significance (F = 1.88, not significant). Patients in cluster A were also older (F = 3.258, P < 0.05) (Table 1), although age and tumour diameter were not associated with each other (r = 0.01, not significant). No differences regarding hormonal values, pathological findings and surgical outcomes were detected among clusters (Table 1).

Table 1

Hormonal and clinical variables in individual patient clusters

	Cluster A	Cluster B	Cluster C	Significance
Demographical features
Sex distribution (male/female)	5/15	2/14	2/2	NS
Age (years)	45.1 ± 3.01	36.7 ± 3.18	30.3 ± 3.74	P < 0.05
Hormonal parameters
Urinary free cortisol (μg 24 h‐1)	547.4 ± 153.4	690.9 ± 120.8	474.7 ± 179.2	NS
Plasma ACTH (pg mL‐1)	89.2 ± 13.9	63.1 ± 6.15	45.3 ± 3.22	NS
Morning serum cortisol (μg dL‐1)	23.5 ± 1.28	19.6 ± 1.60	18.5 ± 2.20	NS
Midnight serum cortisol (μg dL‐1)	20.5 ± 1.48	19.8 ± 2.52	18.2 ± 6.00	NS
Cortisol after 1 mg of dexamethasone (μg dL‐1)	16.5 ± 1.98	15.5 ± 2.83	7.95 ± 4.55	NS
Cortisol suppression after 8 mg dexamethasone (% baseline)	34.9 ± 9.44	43.0 ± 13.51	11.6 ± 1.81	NS
ACTH peak after corticotrophin‐releasing hormone (% baseline)	341.9 ± 93.8	271.5 ± 79.4	392.1 ± 57.4	NS
Macroadenoma (%)	50%	6.2%	0%	P < 0.005
Extrasellar extension (%)	55%	6.2%	0%	P < 0.005
Pathology findings
Tumour diameter (mm)	10.1 ± 1.52	7.0 ± 0.64	6.7 ± 0.25	NS
ACTH immunoreactive cells (%)	86.0 ± 2.48	84.1 ± 3.54	93.7 ± 1.25	NS
Crooke cells (absent/slight/moderate)	45%/45%/10%	69%/25%/6%	75%/25%/0%	NS
Cellular atypia (absent/slight)	55%/45%	75%/25%	100%/0%	NS
Surgical outcomes
Immediate remission (%)	85%	87.5%	100%	NS
Relapse (%)	6%	0%	0%	NS
Length on steroid replacement therapy (months)	19.3 ± 5.9	16.9 ± 3.5	23.2 ± 12.1	NS

Data are the mean ± SEM or percentage. ACTH, adrenocorticotrophic hormone; NS, not significant.

Five patients presented USP8 variants: 2 patients showed missense variant rs672601311 (c.2159C>G, p.720Pro>Arg), 1 patient showed missense variant rs672601307 (c.2152T>C, p.718Ser>Pro), 1 patient showed the 3‐bp deletion variant rs672601306 (c.2155_2157delTCC, p.Ser718del) and 1 patient showed a novel, 15‐bp deletion (c.2157_2171delCCCAGATATAACCCA, p.Pro720_Gln724del). Gene expression profiles in specimens bearing USP8 variants did not group together, except for the 2 specimens with the same missense variant, which clustered next to each other (Figure 1).

Differentially expressed genes

Comparisons of gene profiles from the 3 different clusters allowed identification of over‐ or under‐expressed genes between clusters. A Volcano plot (Figure 3) shows genes up‐regulated with Diff Score ≥13 and genes down‐regulated with Diff Score ≤−13 between cluster A vs clusters B and C (Figure 3A), cluster B vs clusters A and C (Figure 3B), and cluster C vs clusters A and B (Figure 3C).

Figure 3

Volcano plot. Genes up‐ and down‐regulated in the 3 clusters. Effect (average signal, AVG) is shown on the x‐axis and significance (Diff Score) is shown on the y‐axis. Up‐regulated genes appear to the right and down‐regulated genes appear to the left of the x‐axis. (A) Comparison between cluster A and clusters B and C. (B) Comparison between cluster A and clusters A and C. (C) Comparison between cluster C and clusters A and B. Significant genes are highlighted in red (Diff Score >13) and selected genes are identified by name

Of 1259 significantly expressed genes, 4 genes were uniquely overexpressed in cluster A, 313 genes were overexpressed only in cluster B and 29 genes overexpressed only in cluster C (Table 2). Among the 4 genes uniquely overexpressed in cluster A, 3 genes were associated with tumoral processes (ie, BCAS2 [breast carcinoma amplified sequence 2] and RPS6KA6 [ribosomal protein S6 kinase 6] with breast cancer27, 28 and RBX1 [ringbox 1] with gastric and lung cancer). With regard to genes overexpressed in cluster B compared to the other 2 clusters, we observed overexpression of the POMC gene itself, as well as of SCG5, which encodes for secretogranin 5 (also known as neuroendocrine protein 7B2) associated with prohormone convertase 2 activation.29 Furthermore, genes involved in neuroendocrine tumours (eg, NOTCH3, LGALS3, PITX1 and NGRN) were also overexpressed in cluster B compared to clusters A and C. Additional genes of interest to corticotroph adenoma pathophysiology were RXRB and TIMP1 in view of their links to retinoic acid sensitivity30, 31 and matrix metalloproteinase, respectively.32, 33 Given the number of genes overexpressed in this cluster, we could perform functional analysis and detected significant enrichment for 123 genes. Highest enrichment scores were relative to ribosomal function and translation (enrichment 4.89), ribosome assembly (enrichment 1.97), protein biosynthesis (enrichment 1.67) and mRNA processing (enrichment 1.60), thus comprising items recurring with features identified for functional analysis of all corticotroph adenomas (see Supporting information, Table S1). The remaining 190 genes did not fall into specific functional annotation groups. Gene expression profiles in 4 patients differed markedly from all others and clustered separately (cluster C). Amongst genes uniquely overexpressed in this cluster (Table 2), IGFBP5 (insulin‐like growth factor binding protein 5) and genes associated with intracellular calcium flux (ie, CALCOCO1, CAMTA2) appeared to be considerably overexpressed as measured by fold‐change vs other clusters and as of interest to tumoral corticotroph pathophysiology.

Table 2

Genes up‐regulated in individual clusters

Symbol	Diff Score	Definition
Cluster A
BCAS2	22,57	Homo sapiens breast carcinoma amplified sequence 2, mRNA
RPS6KA6	15,18	Homo sapiens ribosomal protein S6 kinase, 90 kDa, polypeptide 6, mRNA
RBX1	13,5	Homo sapiens ringbox 1, mRNA
FAM155A	13,2	Homo sapiens family with sequence similarity 155, member A, mRNA
Cluster B (first 20 genes with highest Diff Scores listed here, the entire list is given in the Supporting information, Table S2)
RAB33B	64,25	Homo sapiens RAB33B, member RAS oncogene family, mRNA
RAX2	64,25	Homo sapiens retina and anterior neural fold homeobox 2, mRNA
EIF2AK4	54,75	Homo sapiens eukaryotic translation initiation factor 2 alpha kinase 4, mRNA
RBM25	51,91	Homo sapiens RNA binding motif protein 25, mRNA
ZNF358	51,62	Homo sapiens zinc finger protein 358, mRNA
ACTR1A	49,43	Homo sapiens ARP1 actin‐related protein 1 homolog A, mRNA
SCG5	49,18	Homo sapiens secretogranin V (7B2 protein), mRNA
NT5C2	48,13	Homo sapiens 5′‐nucleotidase, cytosolic II, mRNA
BTBD2	47,99	Homo sapiens BTB (POZ) domain containing 2, mRNA
TNPO2	46,04	Homo sapiens transportin 2, mRNA
UBXN2A	45,06	Homo sapiens UBX domain protein 2A, mRNA
SH3YL1	43,87	Homo sapiens SH3 and SYLF domain containing 1 (SH3YL1), mRNA
HNRNPUL2	43,84	Homo sapiens heterogeneous nuclear ribonucleoprotein U‐like 2, mRNA
HARS2	43,51	Homo sapiens histidyl‐tRNA synthetase 2, mitochondrial, mRNA
DARS	43,26	Homo sapiens aspartyl‐tRNA synthetase, mRNA
RNGTT	43,15	Homo sapiens RNA guanylyltransferase and 5′‐phosphatase, mRNA
NDUFS8	42,47	Homo sapiens NADH dehydrogenase (ubiquinone) Fe‐S protein 8, 23 kDa (NADH‐coenzyme Q reductase), mRNA
THBS3	42,47	Homo sapiens thrombospondin 3, mRNA
FOXO4	41,59	Homo sapiens forkhead box O4 (FOXO4), mRNA
GCN1	41,08	Homo sapiens eIF2 alpha kinase activator homolog, also known as general control of amino‐acid synthesis 1‐like 1 (GCN1L1), mRNA
Cluster C
SNX29P1	15,12	Homo sapiens sorting nexin 29 pseudogene 1, also known as RUN domain containing 2B, mRNA (RUNDC2B)
IGFBP5	14,78	Homo sapiens insulin‐like growth factor binding protein 5, mRNA
PLD2	14,78	Homo sapiens phospholipase D2, mRNA
SNORA54	14,34	Homo sapiens small nucleolar RNA, H/ACA box 54, small nucleolar RNA
SNORD1141	14,30	Homo sapiens small nucleolar RNA, C/D box 114‐1, small nucleolar RNA
ACIN1	14,28	Homo sapiens apoptotic chromatin condensation inducer 1, mRNA
SNORA10	14,18	Homo sapiens small nucleolar RNA, H/ACA box 10, small nucleolar RNA
SNORD104	14,05	Homo sapiens small nucleolar RNA, C/D box 104, small nucleolar RNA
SNORD41	13,86	Homo sapiens small nucleolar RNA, C/D box 41, small nuclear RNA
SNORA73B	13,74	Homo sapiens RNA, U105A small nucleolar, small nucleolar RNA a.k.a RNU105A
CALCOCO1	13,71	Homo sapiens calcium binding and coiled‐coil domain 1, mRNA
SCAF1	13,69	Homo sapiens SR‐related CTD‐associated factor 1, mRNA
SNORA64	13,65	Homo sapiens small nucleolar RNA, H/ACA box 64, small nucleolar RNA
CAMTA2	13,49	Homo sapiens calmodulin binding transcription activator 2, mRNA
SNORA73B	13,45	Homo sapiens small nucleolar RNA, H/ACA box 73B, small nucleolar RNA
ARSG	13,42	Homo sapiens arylsulfatase G, mRNA
SNORA32	13,36	Homo sapiens small nucleolar RNA, H/ACA box 32, small nucleolar RNA
RNU5A	13,31	Homo sapiens RNA, U5A small nuclear, small nuclear RNA
SNORD12	13,30	Homo sapiens small nucleolar RNA, C/D box 12, small nucleolar RNA
RNU1G2	13,24	Homo sapiens RNA, U1G2 small nuclear, small nuclear RNA
PPID	13,23	Homo sapiens peptidylprolyl isomerase D, mRNA
SNORD3D	13,19	Homo sapiens small nucleolar RNA, C/D box 3D, small nucleolar RNA
SNORD33	13,16	Homo sapiens small nucleolar RNA, C/D box 33, small nucleolar RNA
SNORD68	13,14	Homo sapiens small nucleolar RNA, C/D box 68, small nucleolar RNA
SNORD110	13,11	Homo sapiens small nucleolar RNA, C/D box 110, small nucleolar RNA
SCARNA3	13,05	Homo sapiens small Cajal body‐specific RNA 3, guide RNA
KMT5A	13,05	Homo sapiens lysine methyltransferase 5A, also known as SET domain containing 8 (SETD8), mRNA
SNORD66	13,02	Homo sapiens small nucleolar RNA, C/D box 66, small nucleolar RNA
SNORA71A	13,01	Homo sapiens small nucleolar RNA, H/ACA box 71A, small nucleolar RNA

Real‐time PCR confirmed overexpression of BCAS2 in cluster A, SCG5 in cluster B and IGFBP5 in cluster C (Figure 4), in support of the differential expression observed with respect to microarray data.

Figure 4

Quantitative expression of selected genes in different clusters. (A) Expression of BCAS2 in cluster A vs other clusters. (B) Expression of SCG5 in cluster B vs other clusters. (C) Expression of IGFBP5 in cluster C vs other clusters

DISCUSSION

ACTH‐secreting pituitary tumours are a rare and poorly understood disease. The studies performed so far in adenomatous specimens have indicated the involvement of individual genes (ie, mutations in USP8)3, 4 and the variable expression of factors that mediate sensitivity to glucocorticoids, such as HSP90, BRG‐1 and HSD11B2.5, 34, 35 Furthermore, cellular models derived from rat or dog tumoral corticotrophs called EGF36 and testicular orphan nuclear receptor 4 (TR4)37 into play. Studies in zebrafish also linked the development of pro‐opiomelanocortin (POMC)‐cell tumours to pituitary tumour‐transforming gene (PTTG) overexpression.38 Altogether, several mechanisms appear to be involved in corticotroph tumorigenesis, although a clear‐cut picture of the changes occurring in corticotroph tumours is as yet missing.

The present study had 2 main aims: identifying gene expression patterns characteristic to tumoral corticotrophs and evaluating differences in among corticotroph adenoma subgroups. Over recent decades, it has become increasingly evident that ACTH‐secreting adenomas differ considerably among themselves, in terms of secretory parameters and responses to medical therapy. Thus, studies investigating the variability of these tumours are necessary to identify the features associated with tumour subtypes and to aid in the development of targeted therapeutical approaches.

Heatmap profiling revealed the significant expression of 1259 of 29 285 probes in corticotroph adenoma specimens. Obviously, genes associated with the corticotroph phenotype (ie, POMC, TBX19) were among the highly expressed genes39, 40 Other genes associated with POMC‐processing (ie, proprotein convertase inhibitor [PCSKIN], secretogranin [SCG5], cathepsin L1 [CTSL1]), with drug action on corticotrophs (ie, retinoid X receptors [RXRB], GABA receptors [GABARAPL2, GABRB3]), or with steroid feedback (ie, 11β‐hydroxysteroid dehydrogenase [HSD11BIL]), were also significantly expressed. As regards the individual genes noted above, the expression of HDAC2, SMARCA4 (ie, formerly BRG‐1) TR4, FGFR4 and EGF or EGFR was not included among the significantly expressed genes in the corticotroph adenoma set; PPTG1‐interacting protein (PTTG11P) was significantly expressed, although PTTG itself was not. Likewise, several members of the heat shock protein family A (HSP70) (eg, HSPA8, HSPA9) were detected, although there were none related to heat shock protein alpha family (HSP90). Furthermore, given that the response to vasopressin, somatostatin and dopamine receptor agonists in Cushing's disease has been linked to specific receptor expression,8, 10, 11 we sought evidence for the significant expression of these receptors, however none proved significant.

Several genes detected at normal human pituitary expression profiling41 were also found in tumoral corticotrophs, such as PUM1 (pumilio homolog 1), TPT1 (tumour protein translationally controlled 1) and CDK2AP1 (cyclin‐dependent kinase 2‐associated protein 1), as were some genes detected in other pituitary tumours,42, 43, 44 such as AKT2 (v‐akt murine thymoma viral oncogene homolog 2), NOTCH3 (Notch homolog 3) and PITX1 (paired‐like homeodomain transcription factor 1). Our data confirmed some of the findings obtained in 7 corticotroph macroadenomas with the UniGEM microarray probe on some 7000 genes25 (ie, overexpression of HSPH1, SH3BGRL2 and RGS2, with the latter being of particular interest because it increases during corticotrophin‐releasing hormone incubation in murine tumoral corticotrophs).45 Moreover, we confirmed the results obtained in 12 corticotroph adenomas reporting the expression of several genes (eg, SEZ6L, CALY, CD200, TCF7L2, NISCH, RUNDC3A and NGFRAP1).46 We could also confirm increased expression of cold‐inducible RNA binding protein (CIRP), recently proven to stimulate murine corticotroph tumour growth47 and KRT8 (ie, keratin 8 or CAM 5.2), which, together with significant expression of LGALS3, represents a immunophenotype marker for corticotroph tumours.48, 49, 50 Furthermore, our results corroborate those obtained in 9 adenomas by RNA sequencing on 32 186 protein‐coding genes,51 which revealed the significant expression of a variety of genes, some with clear links to neuroendocrine or pituitary tumoral phenotype (eg, CGA, FILIP1, PITX1, IGFBP5) and others with as yet unknown links to corticotroph tumour development (eg, FBOXO31, FZD7, GPX3, MYH6). Our data also matched the results obtained in ectopic ACTH‐secreting tumours52 (ie, significant expression of MAX dimerisation protein 4 [MXD4], neugrin [NGRN] and SHC transforming protein 1 [SHC1]).

Most recent studies have revealed that up to 60% of human pituitary ACTH‐secreting adenomas may harbour somatic mutations in the ubiquitin‐specific protease 8 (ie, USP8) gene3, 4 and, indeed, 5 specimens in our series presented USP8 variants. The detected mutations were either missense or deletion variants and fell within the 14‐3‐3 binding motif of USP8 suggesting enhanced USP8 activity, as most frequently occurs.3, 4 Interestingly, the 2 specimens carrying the same variant clustered next to each other, indicating highly similar gene expression profiles, whereas the 3 remaining variants were associated with differing expression profiles and fell into separate clusters. It appears that, although variants appear to lead to a defect in USP8 activity, the resulting gene expression profile may diverge; indeed, catalytic activity of USP8 mutants reportedly differs.3

As regards enrichment analysis, pathways involved in protein biosynthesis, macromolecule metabolic processes and RNA processing appeared to be the most significant, in keeping with a neuroendocrine, secretory phenotype. Pathways associated with cell death and DNA repair ranked amongst the low enrichment scores, as expected given the benign nature of these tumours. Of note, there was little overlap with pathways involved in prolactinomas53 and nonsecreting pituitary adenomas.54

Lastly, differential expression analysis was able to identify genes uniquely expressed in each of the 3 clusters formed by unsupervised clustering analysis. Interestingly, although POMC is obviously significantly expressed in all specimens, it proved to be overexpressed together with the POMC‐processing co‐adjuvator SCG5 in specimens grouping in cluster B compared to those in clusters A and C. Increased expression of genes involved in other pituitary tumours or ACTH‐secreting neuroendocrine tumours, most notably NOTCH3, PITX1, LGALS3 and NGRN,38, 43, 44, 52 was also observed in cluster B. Furthermore, adenomas in this cluster overexpressed the retinoid X receptor beta gene (RXRB), thus confirming possibly greater sensitivity to RXR agonists30, 31 in selected adenomas. TIMP1 (tissue metallopeptidase inhibitor 1) was also overexpressed in cluster B adenomas, which is an interesting result given the recent reports on matrix metalloproteinase‐932 and TIMP1 itself33 in corticotroph adenomas.

Four adenomas presented a markedly different gene expression phenotype and clustered separately (cluster C). IGFBP5 was among the genes highly overexpressed in this cluster, which may prove to be of particular interest given its role in both carcinogenesis and normal cell growth.55 Indeed, IGFBP5 is known to be expressed in the pituitary itself, independently from growth hormone‐producing cells,56 and, furthermore, IGFBP5 expression and protein synthesis was increased in a subset of nonfunctioning pituitary adenomas.57 Lastly, overexpression of BCAS2 (breast carcinoma amplified sequence 2) and RPS6KA6 (ribosomal protein S6 kinase 6), both associated with the tumour suppressor p53,27, 28 as well as RBX1 (ringbox 1), a gene linked to cullin‐modulated cell cycle progression58, was detected in cluster A; p53 itself has been shown to accumulate in a subset of pituitary ACTH‐secreting lesions, in particular invasive adenomas.59 Accordingly, we observed that patients in cluster A presented the highest prevalence of macroadenomas and invasive tumours and both p53 and cullin‐RING pathways represent attractive therapeutic targets.58 A recent study aimed at identifying the genetic signature of invasive macroadenomas reported an increased expression of CCND2 (cyclin 2) and ZNF676 (zinc‐finger protein 676) in 3 invasive adenomas compared to non‐invasive corticotroph adenomas.46 Not surprisingly, neither gene was significantly expressed in our series, which comprised 12 invasive adenomas, because our study aim differed and required an unforced approach to gene expression analysis (ie, unsupervised clustering without a priori classification); additional studies into this as well as other group‐derived signatures will allow the above results to be extended and confirmed.

In conclusion, the present study describes the gene expression profile of human pituitary ACTH‐secreting adenomas and has identified genes associated with distinct expression clusters. Overall, it appears that these tumours present a neuroendocrine cell profile but, at the same time, there are clearly distinct gene expression patterns in individual subgroups. This evidence provides the basis for future studies into the molecular pathophysiology of these adenomas and, possibly, paves the way to target therapy.

CONFLICT OF INTERESTS

The authors declare that they have no conflicts of interest that could be perceived as prejudicing the impartiality of the research reported.

Click here for additional data file.

Click here for additional data file.