Transcriptome Analysis Identifies an Attenuated Local Immune Response in Invasive Nonfunctioning Pituitary Adenomas.

BACKGROUND: Invasive nonfunctioning pituitary adenomas (NFPAs) remain challenging due to their high complication rate and poor prognosis. We aimed to identify the distinctive molecular signatures of invasive NFPAs, compared with noninvasive NFPAs, using gene expression profiling by RNA sequencing.
METHODS: We obtained frozen fresh tissue samples from 14 patients with NFPAs who underwent primary transsphenoidal surgery. Three non-invasive and 11 invasive NFPAs were used for RNA sequencing. The bioinformatics analysis included differential gene expression, gene ontology analysis, and pathway analysis.
RESULTS: A total of 700 genes were differentially expressed (59 up-regulated and 641 down-regulated genes) between invasive and non-invasive NFPAs (false discovery rate <0.1, and |fold change| ≥2). Using the down-regulated genes in invasive NFPAs, gene ontology enrichment analyses and pathway analyses demonstrated that the local immune response was attenuated and that transforming growth factor-β (TGF-β) RII-initiated TGF-β signaling was down-regulated in invasive NFPAs. The overexpression of claudin-9 (CLDN9) and the down-regulation of insulin-like growth factor-binding protein 5 (IGFBP5), death-associated protein kinase 1 (DAPK1), and tissue inhibitor of metalloproteinase-3 (TIMP3) may be related with invasiveness in NFPAs.
CONCLUSION: Invasive NFPAs harbor different gene expression profiles relative to noninvasive NFPAs. In particular, local suppression of the immune response and TGF-β signaling can make PAs prone to invasiveness.

INTRODUCTION

Nonfunctioning pituitary adenomas (NFPAs) are differentiated from functioning pituitary adenomas (PAs) in that they do not secrete pituitary hormones that cause clinical symptoms. Hence, subjects with NFPAs present with visual impairment due to the large size of these tumors and their location abutting the optic chiasm. Although most PAs are slow-growing and benign intracranial tumors, about 25% to 55% of NFPAs are invasive and infiltrate into neighboring tissues, such as the cavernous sinuses laterally and the sphenoid sinus or bone inferiorly [12]. In invasive PAs, total surgical resection is difficult to achieve owing to the risk of vascular or nerve damage. Moreover, invasive PAs show a higher recurrence rate caused by the tumor remnants, which require additional surgery or radiation, posing a further risk of complications. The definition of invasiveness for PAs has been controversial, and the concept of invasiveness has been confused with aggressiveness and rapid growth. Trouillas et al. [3] proposed a new histopathological classification system distinguishing “invasive” and “proliferative” tumors. In this framework, invasion was defined as histological and/or radiological signs of cavernous or sphenoid sinus invasion [3].

There have been several efforts to establish molecular mechanisms of invasiveness in NFPAs. The overexpression of fibroblast growth factor receptor 4 (FGFR4), matrix metalloproteases (MMPs) and pituitary tumor transforming gene (PTTG) has been established to be related with invasiveness [4]. Hypoxia-inducible factor-1α expression in response to hypoxia or apoplexy may promote angiogenesis via overexpression of vascular endothelial growth factor (VEGF), thereby allowing tumors to acquire invasive abilities [5]. Nonetheless, the whole picture is not clear since the molecular pathogenesis of invasiveness is highly complex, encompassing multiple genes and pathways.

In this study, we aimed to investigate whether invasive NFPAs harbor a distinctive gene expression profile compared with noninvasive NFPAs and to elucidate new mechanisms underlying the invasiveness of NFPAs through whole RNA sequencing.

METHODS

Subjects

This study was approved by the Institutional Review Board of Seoul National University Hospital (IRB no. 1503-040-654). Written informed consent was obtained from all study subjects. We included 14 subjects with clinical NFPAs who underwent primary transsphenoidal surgery and available fresh frozen tumor tissues. Detailed information on the study protocol has been published previously [2]. The tumor tissues were placed in plastic tubes immediately after surgical resection and stored at −80℃ in liquid nitrogen. We obtained magnetic resonance (MR) images of the sella turcica with T1-and T2-weighted spin echo before and after gadolinium-based contrast medium administration in all patients. The maximal tumor size and invasiveness were measured according to preoperative MR imaging according to the Knosp classification system [6]. Invasive PAs were defined by the presence of cavernous sinus or sphenoid sinus invasion on coronal T1-weighted contrast imaging. Based on these criteria, 11 of the 14 subjects were classified as having invasive NFPAs.

Immunohistochemical staining for Ki-67 (MIB-1, Dako, Glostrup, Denmark) was carried out using the Ventana Benchmark XT automated staining system (Ventana Medical Systems, Tucson, AZ, USA). Morphometric analysis of Ki-67 immunostaining was performed in all PAs. The Ki-67 index was calculated as the mean ratio of the number of tumor cells with positive nuclear Ki-67 staining to the total number of tumor cells in the three areas with the highest positive nuclear staining cell ratio.

RNA extraction from tissues

Total RNA was extracted using an RNeasy mini kit from Qiagen (Qiagen, Venlo, the Netherlands), according to the manufacturer's recommendations. RNA quality was evaluated using an Agilent 2100 Bioanalyzer with the Agilent RNA 6000 Nano Kit (Agilent, Santa Clara, CA, USA).

Total RNA sequencing

We used 100 ng of total RNA from all subjects to prepare sequencing libraries by using the TruSeq stranded total RNA sample preparation kit (Illumina, San Diego, CA, USA) which combines RiboZero rRNA depletion with a strand-specific method similar to the deoxyuridine diphosphate method. The quality of these cDNA libraries was evaluated with the Agilent 2100 BioAnalyzer (Agilent). They were quantified with the KAPA library quantification kit (Kapa Biosystems, Wilmington, MA, USA) according to the manufacturer's library quantification protocol. Following cluster amplification of denatured templates, paired-end sequencing progressed (2×100 bp) using the Illumina HiSeq 2500.

RNA-seq data analysis

Sequencing reads were mapped to the human genome (hg19) using STAR [7]. Transcript quantification from the mapped reads and the newly assembled transcripts was performed using HTseq [8]. Fragments per kilobase per million sequenced reads (FPKM)-normalized gene expression were used for the further analysis. Analyses of FPKM data were conducted using iDEP 0.90, using the R package limma [9]. Genes with FPKM values >0.5 in all samples were selected. Genes with a |fold change (FC)| ≥2 and false discovery rate (FDR) for log2-transformed (FPKM+1) <0.1 were considered to be differentially expressed. Gene ontology enrichment analysis and pathway analysis using Reactome 2016 were performed using enrichR [10].

Statistical analysis

Data are shown as mean±standard deviation or number (%). We performed the chi-square test for categorical variables and the Student t test for continuous variables. A P value <0.05 was defined as indicating statistical significance.

RESULTS

Table 1 showed the clinical, radiological, and pathological characteristics of study subjects with noninvasive and invasive NFPAs. Age was not related with invasiveness, but female sex was more prevalent among the patients with invasive NFPAs. The tumor size was not significantly different between the two groups. The proportion of clivus invasion, regrowth after gross total resection, and a Ki-67 proliferation index >3% was higher in subjects with invasive NFPAs than in those with non-invasive NFPAs.

Table 1

Clinical, Radiological, and Pathological Characteristics of Study Subjects with Noninvasive and Invasive NFPAs

Variable	Subjects with non-invasive NFPAs (n=3)	Subjects with invasive NFPAs (n=11)	P value
Age, yr	33.3±14.6	41.8±10.1	0.503
Female sex	1 (33.3)	5 (45.5)	0.015
Maximal tumor size, cm	30.2±5.4	34.5±9.1	0.412
Clivus invasion	0	5 (45.5)	0.011
Regrowth after GTR	0	4 (36.4)	0.012
Ki-67 index >3%	0	5 (45.5)	0.011

Values are expressed as mean±standard deviation or number (%).

NFPA, nonfunctioning pituitary adenoma; GTR, gross total resection.

Transcriptome analysis was performed in frozen tissues of three non-invasive and 11 invasive NFPAs. As shown in Fig. 1, hierarchical clustering and principal component analysis indicated the presence of a substantial difference between the two groups. Using the limma package, we identified 700 significantly differentially expressed genes between non-invasive and invasive NFPAs: 59 up-regulated and 641 down-regulated genes (FDR <0.1 and |FC| ≥2). A plot of magnitude and abundance and a scatter plot show the differences in the gene expression profile between the two groups (Fig. 2). Table 2 presents the top 20 differentially expressed genes in invasive NFPAs compared with noninvasive PAs. In Fig. 3, the log-transformed FPKM expression of genes is shown. Immune system-related genes, such as immunoglobulin kappa constant (IGKC), complement C1s (C1S), complement C1r (C1R), interferon induced transmembrane protein 1 (IFITM1), and transforming growth factor-β (TGF-β) signaling-related genes, such as TGFRB2 and TGFB, were down-regulated in invasive NFPAs. Due to the low number of up-regulated genes in invasive NFPAs, further analysis was performed using only down-regulated genes in invasive NFPAs.

Fig. 1

(A) Heatmap showing gene expression profiles by hierarchical clustering. (B) Principal component analysis (PCA) showing the distinction between noninvasive and invasive pituitary adenomas.

Fig. 2

(A) A magnitude and abundance (MA) plot and (B) a scatter plot show the different gene expression profiles between the two groups.

Table 2

The Top 20 Differentially Expressed Genes in Invasive PAs Compared with Noninvasive PAs (|FC|≥2, FDR <0.1)

Gene	Gene description	Log2FC	FDR
Up-regulated genes
ZCCHC12	Zinc finger CCHC-type containing 12	2.81	0.047
NEFM	Neurofilament, medium polypeptide	2.54	0.073
SCD	Stearoyl-CoA desaturase	2.42	0.055
DCX	Doublecortin	2.24	0.026
EPHA7	EPH receptor A7	2.17	0.026
SLC24A2	Solute carrier family 24 member 2	2.12	0.040
FRMPD1	FERM and PDZ domain containing 1	2.12	0.036
SLC8A2	Solute carrier family 8 member A2	2.01	0.059
RNF157	Ring finger protein 157	1.98	0.034
NRG1	Neuregulin 1	1.96	0.001
SLC32A1	Solute carrier family 32 member 1	1.88	0.056
CLDN9	Claudin-9	1.75	0.041
LRRN2	Leucine rich repeat neuronal 2	1.66	0.017
SCG2	Secretogranin II	1.65	0.041
NDRG4	N-myc downstream-regulated gene 4 protein	1.62	0.034
IL5RA	Interleukin 5 receptor subunit alpha	1.62	0.014
FGF12	Fibroblast growth factor 12	1.62	0.031
NLRP1	NLR family pyrin domain containing 1	1.55	0.089
EEF1A2	Eukaryotic translation elongation factor 1 alpha 2	1.54	0.020
LANCL3	NDRG family member 4 LanC like 3	1.53	0.049
Down-regulated genes
MGP	Matrix Gla protein	−5.22	<0.001
DCN	Decorin	−5.13	<0.001
IGKC	Immunoglobulin kappa constant	−4.81	<0.001
C1S	Complement C1s	−4.58	<0.001
C1R	Complement C1r	−4.44	<0.001
TGM2	Transglutaminase 2	−4.41	<0.001
COL3A1	Collagen type III alpha 1 chain	−4.17	0.012
IGFBP6	Insulin like growth factor binding protein 6	−3.97	<0.001
AGR2	Anterior gradient 2, protein disulfide isomerase family member	−3.94	<0.001
SERPING1	Serpin family G member 1	−3.86	<0.001
BGN	Biglycan	−3.81	0.001
CYBRD1	Cytochrome b reductase 1	−3.80	<0.001
IGFBP4	Insulin like growth factor binding protein 4	−3.66	<0.001
AEBP1	AE binding protein 1	−3.66	<0.001
IFITM1	Interferon induced transmembrane protein 1	−3.65	0.001
FAT1	FAT atypical cadherin 1	−3.60	0.001
CFH	Complement factor H	−3.57	<0.001
COL6A2	Collagen type VI alpha 2 chain	−3.53	0.004
ANXA1	Annexin A1	−3.47	0.001
TAGLN	Transgelin	−3.43	<0.001

PA, pituitary adenoma; FC, fold change; FDR, false discovery rate.

Fig. 3

(A–F) Expression of immune response-related genes (immunoglobulin kappa constant [IGKC], complement C1s [C1S], complement C1r [C1R], interferon induced transmembrane protein 1 [IFITM1]) and transforming growth factor-β (TGF-β) signaling-related genes (TGFBR, TGFB1). PA, pituitary adenoma; FPKM, fragments per kilobase per million.

Gene ontology enrichment analysis was performed using the down-regulated genes in invasive NFPAs using enrichR (Table 3). The biological processes of extracellular matrix (ECM) organization, cytokine-mediated signaling, and neutrophil activation involved in the immune response were down-regulated in invasive NFPAs. In terms of molecular function, collagen binding and integrin binding were down-regulated in invasive NFPAs. Within the category of cellular components, focal adhesion and the integral component of the plasma membrane were down-regulated in invasive NFPAs.

Table 3

Enriched GO Terms of the Down-Regulated Genes in Invasive NFPAs

	GO term ID	GO term	Overlap	Adjusted P value	Combined score
Biological processes	GO:0030198	Extracellular matrix organization	46/230	5.08E-20	330.5
	GO:0019221	Cytokine-mediated signaling pathway	57/634	5.53E-09	75.3
	GO:0030334	Regulation of cell migration	38/317	4.64E-09	99.6
	GO:0042127	Regulation of cell proliferation	60/741	5.50E-08	60.3
	GO:0002283	Neutrophil activation involved in immune response	46/484	5.12E-08	70.3
Molecular functions	GO:0005518	Collagen binding	16/53	5.85E-09	245.0
	GO:0005178	Integrin binding	20/95	1.01E-08	162.7
	GO:0031995	Insulin-like growth factor II binding	5/7	2.54E-04	317.1
	GO:0003779	Actin binding	24/255	7.57E-04	37.7
	GO:0019199	Transmembrane receptor protein kinase activity	11/63	0.001048	67.0
Cellular components	GO:0005925	Focal adhesion	51/357	1.04E-16	191.2
	GO:0005887	Integral component of plasma membrane	96/1,464	2.55E-09	51.6
	GO:0005788	Endoplasmic reticulum lumen	30/271	5.61E-07	67.0
	GO:0045121	Membrane raft	19/120	9.92E-07	91.6
	GO:0070821	Tertiary granule membrane	15/74	9.28E-07	116.3

GO, Gene Ontology; NFPA, nonfunctioning pituitary adenoma.

Pathway analysis of the down-regulated genes in invasive NFPAs was conducted using the Reactome 2016 package (Table 4). ECM organization, immune system, integrin cell surface interactions, cytokine signaling in the immune system, innate immune system, assembly of collagen fibrils, and ECM-proteoglycans were down-regulated in invasive NFPAs.

Table 4

Enriched Pathway Analysis of the Down-Regulated Genes in Invasive PAs Using Reactome 2016

Reactome term ID	Term	Overlap	Adjusted P value	Combined score
R-HSA-1474244	Extracellular matrix organization	49/283	4.00E-19	268.5
R-HSA-168256	Immune system	108/1547	5.66E-12	70.9
R-HSA-216083	Integrin cell surface interactions	17/67	1.24E-08	193.5
R-HSA-109582	Hemostasis	48/552	1.53E-07	58.7
R-HSA-1280215	Cytokine signaling in immune system	50/620	6.88E-07	50.1
R-HSA-168249	Innate immune system	58/807	2.08E-06	41.8
R-HSA-2022090	Assembly of collagen fibrils and other multimeric structures	13/54	2.41E-06	137.6
R-HSA-3000178	ECM proteoglycans	13/55	2.68E-06	133.4
R-HSA-3000157	Laminin interactions	9/23	3.13E-06	217.4
R-HSA-445355	Smooth muscle contraction	10/33	7.75E-6	158.8

PA, pituitary adenoma; ECM, extracellular matrix.

DISCUSSION

In the present study, transcriptome analysis using RNA sequencing revealed that invasive NFPAs expressed different molecular signatures compared with noninvasive NFPAs. In particular, among 700 differentially expressed genes, immune system-related genes were down-regulated in invasive NFPAs.

The down-regulated pathways involved in the immune system, integrin cell surface interactions, cytokine signaling, and innate immune system were the distinctive features of invasive NFPAs. These features of immunosuppression have also been reported in other studies. Richardson et al. showed that silent subtype III PAs may be aggressive due to suppression of the local immune response [11]. The genes arginase-2 (ARG2) and semaphorin 3A (SEMA3A), which are related to T-cell regulation and immunosuppression, were found to be highly expressed in silent subtype III PAs [11]. Through the direct data integration of all published PA-related microarray datasets, Yang et al. [12] discovered that 66 immune-related genes were down-regulated in PAs. They did not investigate the relationship between immune-related genes and invasiveness, but they proposed that immune-related genes may play a role in PA development. Mei et al. [13] directly explored the expression of programmed death ligand 1 in NFPAs, but failed to show any difference in expression according to invasiveness. Yang and Li [5] also suggested that the IP3 pathway and VEGF-induced immune escape might play a role in invasiveness of PAs. Invasive PAs can evade immune surveillance through suppression of the immune response and invade parasellar structures.

We also demonstrated that TGF-β signaling was related with the invasiveness of NFPAs. TGFB1 (TGF-β1) expression was lower in invasive PAs compared with noninvasive PAs in our study (FDR=6.16E-02, log2|FC|=−1.65) (Supplemental Table S1). Zhenye et al. [14] also reported the down-regulation of Smad3, phospho-Smad3, and TGF-β1 expression in invasive PAs [14]. In addition, we demonstrated that the expression of TGFBR2 (TGF-β RII), not TGFBR1, was lower in invasive NFPAs than noninvasive NFPAs (FDR=2.72E-03, log2|FC|=−2.35) (Supplemental Table S1). This finding is compatible with a previous report [15]. TGF-β RII can recruit and phosphorylate TGF-β RI to activate downstream Smad-dependent signaling. Thus, TGF-β RII may play a more pivotal role than TGF-β RI in the development of invasive PAs. Overall, the down-regulation of TGF-β signaling may be involved in the development of invasive NFPAs.

We compared our RNA sequencing analysis results with previous microarray data. Galland et al. [16] suggested that insulin-like growth factor-binding protein 5 (IGFBP5), myosin-Va (MYO5A), FMS-like tyrosine kinase 3 (FLT3), and nuclear factor, erythroid 2 like 1 (NFE2L1) were overexpressed in invasive NFPAs, and that MYO5A was a useful marker of tumor invasiveness. However, there was no significant difference in MYO5A, FLT3, and NFE2L1 expression between the two groups. Instead, IGFBP5 was down-regulated in invasive NFPAs (FDR=0.002, log2|FC|=−3.28) (Supplemental Table S1). IGFBP5 is known to exert anti-cancer activity by inhibiting angiogenesis, which is compatible with our results [17]. Furthermore, drawing on different transcriptome analysis methods, in our study, invasive NFPAs had a higher Ki-67 index than noninvasive NFPAs, whereas in the study of Galland et al. [16], the Ki-67 index was not different between the two groups.

As another microarray-applied study, de Araujo et al. [18] suggested that cyclin D2 (CCND2) and zinc finger protein 676 (ZNF676) were overexpressed, whereas death-associated protein kinase 1 (DAPK1) and tissue inhibitor of metalloproteinase-2 (TIMP2) were down-regulated in invasive corticotropinomas. However, expression of the genes CCND2, ZNF676, and TIMP2 was not significantly different between invasive and noninvasive PAs in our study. Only the gene DAPK1 was down-regulated in invasive PAs (P=0.015, log2|FC|=−1.55, data not shown). DAPK1 codes for death-associated protein kinase 1, which acts as a tumor suppressor, and hypermethylation of the DAPK1 promoter is associated with greater invasiveness in head and neck cancers [19]. Among the inhibitors of the matrix metalloproteinases, tissue inhibitor of metalloproteinase-3 (TIMP3) was down-regulated in invasive NFPAs in our study (FDR=9.39E-03, log2|FC|=−3.00) (Supplemental Table S1). TIMP3 expression in PAs was related with PA fibrosis [20], which may restrain cell proliferation within tumors. Although the subtype of PAs analyzed by de Araujo et al. [18] study was corticotropinomas, invasiveness-related markers can be applied to NFPAs.

Cao et al. [21] also identified dysregulation of leukocyte transendothelial migration and cell adhesion molecules as invasion-related pathways in NFPAs using microarray analyses. The invasion-related genes caludin-7 (CLDN7), contactin associated protein like 2 (CNTNAP2), integrin subunit alpha 6 (ITGA6), junctional adhesion molecule 3 (JAM3), protein tyrosine phosphatase receptor type C (PTPRC), and catenin alpha 1 (CTNNA1) identified by Cao et al. [21] were not validated in our study. Instead, claudin-9 (CLDN9), which is a member of the claudin family, was overexpressed in invasive NFPAs (FDR=4.13E-02, log2|FC|=1.75) (Supplemental Table S1). Higher expression of CLDN9 was also found in invasive pituitary oncocytomas [22]. Claudin-9, which is encoded by CLDN9, may interact with matrix metalloproteinases, weaken the vascular endothelium, and increase paracellular permeability, which results in invasion [22]. The gene JAM2, junctional adhesion molecule 2, was under-expressed in invasive NFPAs (FDR=4.56E-03, log2|FC|=−1.70) (Supplemental Table S1). Down-regulation of JAM2 was correlated with disease progression and metastasis of colorectal cancer [23]. Thus, the overexpression of CLDN9 and the down-regulation of JAM2 may be potential markers of invasive NFPAs.

Transcriptomic research using RNA sequencing to discover markers of invasive NFPAs has scarcely been conducted. Most transcriptome studies of invasive NFPAs have performed microarray analyses. Compared with microarray analyses, RNA-sequencing generates big data regarding the gene expression profile, thereby enabling the discovery of new molecular mechanisms underlying invasiveness in NFPAs. Based on the attenuated local immune response in invasive NFPAs, cancer immunotherapy such as immune checkpoint inhibitors may be another therapeutic option for invasive and unresectable NFPAs.

Several limitations should be mentioned. The sample size was small, especially for noninvasive NFPAs. The unbalanced distribution of the sample size may have led to selection bias. Immunochemistry data were not included. We defined invasive NFPAs in terms of cavernous sinus or sphenoid sinus invasion. However, we could not exclude the possibility of intrinsic invasiveness in tumors that had not yet invaded the cavernous sinus. In addition, we did not validate the significant differentially expressed genes with real-time quantitative reverse transcription polymerase chain reaction or functional studies.

Taken together, invasive NFPAs have different gene expression profiles relative to noninvasive NFPAs. Invasive NFPAs can escape immune attack due to the attenuated local immune response, and express down-regulated TGF-β signaling. However, our findings discovered by RNA-sequencing need to be validated in a large cohort study and with experimental studies.