Transcriptome of GH-producing pituitary neuroendocrine tumours and models are significantly affected by somatostatin analogues

Abstract

Pituitary neuroendocrine tumours (PitNETs) are neoplasms of the pituitary that overproduce hormones or cause unspecific symptoms due to mass effect. Growth hormone overproducing GH-producing PitNETs cause acromegaly leading to connective tissue, metabolic or oncologic disorders. The medical treatment of acromegaly is somatostatin analogues (SSA) in specific cases combined with dopamine agonists (DA), but almost half of patients display partial or full SSA resistance and potential causes of this are unknown. In this study we investigated transcriptomic landscape of GH-producing PitNETs on several levels and functional models—tumour tissue of patients with and without SSA preoperative treatment, tumour derived pituispheres and GH3 cell line incubated with SSA to study effect of medication on gene expression. MGI sequencing platform was used to sequence total RNA from PitNET tissue, pituispheres, mesenchymal stromal stem-like cells (MSC), and GH3 cell cultures, and data were analysed with Salmon—DeSeq2 pipeline. We observed that the GH-producing PitNETs have distinct changes in growth hormone related pathways related to its functional status alongside inner cell signalling, ion transport, cell adhesion and extracellular matrix characteristic patterns. In pituispheres model, treatment regimens (octreotide and cabergoline) affect specific cell proliferation (MKI67) and core functionality pathways (RYR2, COL8A2, HLA-G, ARFGAP1, TGFBR2). In GH3 cells we observed that medication did not have transcriptomic effects similar to preoperative treatment in PitNET tissue or pituisphere model. This study highlights the importance of correct model system selection for cell transcriptomic profiling and data interpretation that could be achieved in future by incorporating NGS methods and detailed cell omics profiling in PitNET model research.

Background

Pituitary neuroendocrine tumours (PitNETs) are common intracranial neoplasms with variable health effects and prognosis with the prevalence of clinically significant tumours 1 per 1000 individuals [1]. Despite the fact that PitNETs are rarely metastatic, they are usually the cause of significant morbidity and mortality [2]. Although currently WHO uses expression of specific transcription factors as basis of PitNET classification (six morphologically distinct types based on transcription factors with three main cell lineages: SF1 (gonadotroph cells), PIT1 (lactotroph, somatotroph, mammosomatotroph, and thyrotroph cells), and TPIT (corticotroph cells) they can be characterised according to the type of synthesised hormone and distinct clinical symptoms [3]. Hormonally active PitNETs tend to overproduce either growth hormone (GH), prolactin (PRL), adrenocorticotropic hormone (ACTH), thyroid-stimulating hormone (TSH) and gonadotropins like luteinizing (LH) and follicle-stimulating (FSH) hormones or some combination of them [4, 5]. The majority (70–75%) of non-functioning PitNETs are gonadotroph tumours as they are immunopositive in 60% of the cases (hormone negativity in 40% of cases) for FSH and LH but do not secrete them in significant quantities [3, 6]. Current treatment strategies for patients with PitNETs are aimed to control biochemical activity and reduce tumour size while maintaining normal pituitary function or entirely remove the tumour via surgery [7]. Due to the fact that cyclic peptide somatostatin regulates pituitary hormone secretion, somatostatin binding receptors (SSTR) have become valuable therapeutic targets [8]. About 90% of GH-producing PitNETs express SSTR2 and SSTR5 which are targeted by somatostatin analogues (SSA), when intervention fails to control tumour activity dopamine agonists (DA) that target dopamine receptor 2 (D2R) are introduced to improve disease management [9]. Pegvisomant (GH receptor antagonist) can be added for acromegaly patients with only partial response to SSAs. As a result, this combination can normalize insulin-like growth factor 1 (IGF-1) levels in patients more effectively than SSA monotherapy. Since SSAs reduce excessive production of GH by PitNET, pegvisomant decreases GH actions in peripheral tissues by blocking the increased production of IGF-I in the liver [10].

Multiple studies have been devoted to investigating the PitNET transcriptome in recent years with the aim to provide a basis for PitNET classification and understanding association between tumour development, its clinical characteristics and transcriptome perturbation [11,12,13,14]. As a result, numerous genes have been identified controlling a vast array of PitNET characteristics [1]. POU1F1 gene has shown to be involved in development of GH, PRL and TSH adenomas [12]. It is also responsible for the regulation of hormone secretion as the expression of POU1F1 correlates with the levels of secreted PRL, TSH and GH [15]. Several coding genes (CLDN9, IGFBP5, DAPK1 and TIMP3) and non-coding genes (LINC00473 and CDKN2BAS) have been shown to be associated with the invasiveness of PitNETs [13, 14, 16]. Expression of MUC16, MACC1 and GRHL2 has been reported to be altered in response to SSA treatment [17]. Non-coding RNA LL21NC02-21A1.1 and protein coding gene NOL6 have been associated with the recurrence of PitNETs [18]. The elaborate pan-genomic study of transcriptomic landscape of PitNET subgroups based on 2017 WHO transcription-factor classification [11]. The report detected that tumours with USP8 and GNAS mutations have distinct transcriptomic profiles compared to tumours without these tumour driver variants, indicating that transcriptomic studies can bring valuable information on PitNET functional aspects [11].

So far, the impact of SSA therapy on the transcriptomic landscape of PitNETs has been investigated in two studies [11, 17]. It has been demonstrated that SSA preoperative therapy downregulated Ki67 levels and upregulated MUC1 and CD40 expression in tumour tissue transcriptomes [11]. Additionally, expression of several tumourigenesis related factors are decreased upon SSA/DA treatment—MUC16, MACC1, and GRHL2 and extracellular matrix related collagen pathways might have implications in PitNET response to medication [17]. Therefore, more studies are necessary to more comprehensively assess the impact of SSA/DA treatment on the transcriptomic landscape of PitNETs.

The most often used model cell lines used in PitNET functional studies are derived from Rattus norvegicus (MMQ, GH3, RC-4B) or Mus musculus (AtT-20) PitNET tumours, and the representation of human pathobiology of PitNETs in these cell line remains questionable. Although there are reports of human PitNET cell lines used in functional studies, widely accepted human PitNET models are not commercially available. Many authors have demonstrated free-floating sphere (pituisphere) formation obtained from the primary pituitary tumour [19,20,21,22,23]. According to reports, PitNETs have detectable levels of cancer stem cells (CSC), which were found in a variety of tumour types and are thought to promote tumour growth and tissue invasion, as well as resistance to therapy. CSCs are considered to be salient players in PitNET development. These cells mainly display expression of stemness markers OCT4, CD133, nestin, SOX2, and CXCR4 and demonstrate self-renewal competence [24]. We have previously demonstrated that pituispheres genetically correspond to the PitNET tissues, showing that pituispheres might be a valuable novel model system for PitNET tumourigenesis and therapy response studies [21]. Mesenchymal stromal stem-like cells (MSC) are thought to represent tumour microenvironment [25]. It has been recorded that MSCs play an essential role in tumour formation and progression via causing epithelial-mesenchymal transition (EMT). This cell type is adherent in standard culture and express surface markers CD73, CD90, and CD105 [21, 25]. We have previously demonstrated that PitNET derived MSCs do not contain somatic variants linked to the tumour [21, 26].

In this study we aimed to characterise transcriptomic patterns of GH-producing PitNETs specifically assessing impact of SSA treatment effects on gene expression on several levels and functional models of GH-producing PitNETs: tumour tissue of patients with and without SSA preoperative treatment, tumour derived pituisphere model and classically used GH3 cell lines treated with SSA. For the first time, we incorporate various complementary tumour models to derive a comprehensive assessment of SSA impact on GH-producing PitNET functionality.

Materials and methods

Study group

Tumour tissue samples were collected from 82 patients who underwent planned resection at Pauls Stradins Clinical University Hospital, Latvia from 2010 till 2021. Patients’ clinical data and medical treatment history is presented in Additional file 1: Table S1. All patients were recruited to the Latvian National biobank—Genome Database of the Latvian population [27]. Broad informed consent for biobank and project-specific consent for research involving the pituitary tumours were obtained from all patients (approved by the Central Medical Ethics Committee of Latvia protocol No. 22.03.07/A7 and 01.29.1/28/renewed prot. No. 01-29.1/5035, respectively).

PitNET tissue samples (after resection) were divided into two parts. One part was submerged in RNAlater Solution (Thermo Fisher Scientific, USA) for DNA/RNA extraction, and another part was immersed in Dulbecco’s Modified Eagle Medium (DMEM) (Thermo Fisher Scientific, USA) containing 1 × penicllin/streptomycin solution (GIBCO, USA) for cell culture development.

Study design

We developed a complex design to investigate various aspects of GH-producing PitNETs transcriptomics (Fig. 1). Comparison of GH-producing and non-functional PitNETs included the largest sample set with 21 and 61 samples, respectively. Within the GH-producing PitNET group there were nine patients with SSA/DA drug therapy and 11 without it. Primary derived cell cultures (MSC and pituitspheres) derived from eight PitNETs were propagated for 72 h incubated with two common PitNET drugs: octreotide and cabergoline. Finally, widely used commercial cell line GH3 was used to evaluate reproducibility of octreotide incubation related DEGs in rat model cell line. Information regarding tested sample groups and their related information is available in Additional file 2: Tables S1–S3.

Fig. 1

Graphic study design displaying aims, comparison groups and sample sizes

Transcriptome sequencing

Total RNA for transcriptome sequencing was extracted using AllPrep DNA/RNA/miRNA universal kit (Qiagen, Germany) according to manufacturer’s instructions from tumour tissue samples stored in RNAlater Tissue Storage Reagent (Sigma-Aldrich, USA). The concentrations of extracted RNA were measured using Qubit 2.0 with Qubit RNA HS kit (Thermo Fisher, USA). The quality of extracted RNA was assessed using Agilent 2100 Bioanalyzer (Agilent Technologies, USA). Following this the RIN and DV200 values were calculated to calculate the RNA input amounts for RNA-seq compatible library preparation. Prior to the library preparation rRNA removal was performed using MGIEasy rRNA Depletion kit (MGI, PRC). After the rRNA depletion reverse transcription, second strand synthesis, and cDNA library preparation for NGS were carried out using MGIEasy RNA Directional Library Prep Set (MGI, PRC). The paired-end libraries were sequenced on DNBSEQ-G400 platform (MGI, PRC) with the aim of 35 million reads per sample at 150 bp read length.

Culturing of PitNET tissue material and pituisphere sequencing

Within 12 h after surgery, PitNET tissue samples were processed for propagation. Tissue material was mechanically sliced into small pieces and washed in DMEM with 1 × Antibiotic–Antimycotic solution (Thermo Fisher Scientific, USA). Enzymatic dissociation using Accutase solution (Thermo Fisher Scientific, USA) was carried out on a rotating platform for 20 min at 37 °C in a humidified environment with 5% CO2. The cells were centrifuged for 5 min at 360 × g after the incubation period to obtain cell pellets. The cell pellet was treated with a red blood cell lysis solution (154 mM NH4Cl, 10 mM KHCO3, 0.1 mM EDTA, pH 7.4) for 10 min to reduce contamination with red blood cells. To remove red blood cell debris, the sample was centrifuged, and the cell pellet was washed twice. The obtained pellet was split into two portions. To obtain PitNET tissue-derived free-floating spheres, cells were grown in DMEM-F12 (Thermo Fisher Scientific, USA), containing 1 × penicillin/streptomycin solution, 20 ng/ml epidermal growth factor (EGF) (Sigma-Aldrich, Germany), 10 ng/ml basic fibroblast growth factor (bFGF) (Sigma-Aldrich, Germany), and 1 × B27 supplement (GIBCO, USA). To obtain MSC, culture cells were grown in DMEM-F12, supplemented with 10% Fetal Bovine Serum (FBS) (Thermo Fisher Scientific, USA), 1% ITS (Corning, USA), and 100 µg/ml primocin (InvivoGen, USA) until confluent, then propagated and passaged 2–6 times. All cell culture cultivations and incubations were performed at 37 °C, 95% air, and 5% CO2. For the gene and protein expression experiments, pituispheres were grown on a 6-well plate and treated with 10 μM octreotide or 10 μM cabergoline for 4, 8, 24, 48, and 72 h.

For pituisphere transcriptome sequencing we generated and amplified 62 cDNA from tubes containing pituispheres using REPLI-G WTA single cell kit according to manufacturer’s instructions. The kit included lysis reagents and gDNA removal reagents therefore no prior RNA extraction from tubes containing pituispheres was required. cDNA was generated only from Poly A transcripts to avoid rRNA reads in sequencing. After cDNA was generated, the cDNA was amplified using multiple displacement amplification to meet the required concentrations for downstream processing. The quality and concentration of final cDNA was tested using Agilent 2100 bioanalyzer and Qubit 2.0. For fragmentation before the library preparation, we optimised a method using Covaris S220 ultrasonicator (Covaris, USA). For ultrasonication 1000 ng of cDNA was used for each sample. The following settings were used for fragmentation: target bp—300 bp, peak incident power—75 W, duty factor—20%, cycles per burst—1000, fragmentation time—45 s. Fragmented cDNA samples were visualised on Agilent 2100 bioanalyzer prior to the library preparation. 62 libraries were prepared with DNBSEQ-G400 compatible library preparation kit—MGIEasy PCR–Free DNA Library Prep Set (MGI, PRC) according to the manufacturer’s instructions. The quality and concentrations of libraries were evaluated using Agilent 2100 Bioanalyzer and Qubit 2.0. We sequenced the pituisphere paired-end libraries on DNBSEQ-G400 sequencer (MGI, PRC) with aim of 35 million reads per sample at 150 bp read length.

GH3 cell line culturing and sequencing

Rat pituitary derived GH3 cell line was obtained from ATCC (American Type Culture Collection, USA). The GH3 cells were maintained in F12-K medium containing 15% horse serum, 2.5% FBS and 1 × penicillin/streptomycin solution. For studying the effects of PitNET drugs on gene expression, the cells were incubated for 4, 8, 24, 48, and 72 h with 10 μM octreotide. The GH3 cells were sequenced as previously described in the context of tissue material sequencing.

Data analysis

Raw sequencing data quality control was performed using FastQC (v0.11.9) and MultiQC (v1.10) software [28]. Afterwards, the paired-end data was trimmed with fastp (v0.23.2) software to retain reads with average base quality of at least 20 (Phred score) and minimum read length 100 base pairs [29]. MGI/BGI sequence ID correction and overlapping read base correction were also enabled to minimise chance of incompatibility with the aligner and possibility of base level sequencing errors. Read count and quality was again inspected using the aforementioned QC software which was followed by ribosomal RNA (rRNA) removal with SortMeRNA (v4.3.4) before final quality and read count evaluation [30]. Reads were further quasi-mapped and quantified using Salmon (v1.6.0) against the GENCODE (v38) Homo sapiens gentrome for tissue and pituisphere samples and against Ensembl (v105) Rattus Norvegicus gentrome for the GH3 cell lines with GC bias, sequence level bias and positional bias correction enabled [31]. Tissue samples with quasi-mapping rate of at least 45% were selected for differential expression analysis. R software (v4.1.1) was employed to summarise gene level counts using the tximeta package (v1.12.4) [32, 33]. Differential expression analysis was performed by DESeq2 (v1.34.0) [34]. For tissue samples, quantified read counts were filtered by frequency setting the count threshold at 10 and the sample frequency threshold as 33% of the smallest comparison group contrast size, accordingly 3 for the SSA/DA therapy comparison and 7 for the PitNET type comparison.

Wald test was used to determine expression differences. The default independent filtering function was replaced with independent hypotheses weighing from the IHW package (v1.22.0) to increase statistical power by adjusting p values based on the mean expression level of each gene which was followed with multiple testing correction using the Benjamini–Hochberg adjustment [35]. Log Fold Change (L2FC) shrinkage algorithm from the apeglm package (v1.16.0) was used to correct L2FC values for genes with low counts and high dispersions [36]. Heatmap of gene level normalised count values was graphed by pheatmap (v1.0.12) package and used for visualisation of the differences between the selected differentially expressed genes (DEG) [37]. To inspect p-value variability across the full range of L2FC values, a volcano graph was drawn using the EnhancedVolcano (v1.12.0) package, where the threshold for p-value was < 0.05 and L2FC threshold was the same as in the shrunken results table filtering step [38]. Box plots for each DEG were constructed with the ggplot2 (v3.3.5) package to inspect changes in normalised count values in the possible candidate genes [39]. To gain insight in DEG involvement with signalling pathways or association with disease, enrichment was performed by using STRING-db’s (v11.5) online tool [40].

Most of the mentioned differential expression analysis steps above were also applied to the Rattus norvegicus model data of GH3 cells and in large part also to data of the pituisphere model experiments, except the Likelihood Ratio Test (LRT) was used, to perform a time-series analysis against a reduced model which allowed us to observe DEG’s with a more concurrent trend across all included time points. Subsequently a Wald test was performed between each incubation period and the unincubated controls. No L2FC shrinkage or IHW p-value correction was applied to the pituisphere LRT test or GH3 cell model results. For the incubated GH3 cell analysis, null hypothesis correction was applied to infer more accurate p-values and subsequent false discovery rate detection using the fdrtool (v1.2.17) [41]. Information about sample groups and related information available in Additional file 2: Tables S1–S3.

Results

Transcriptomic patterns of GH-producing PitNETs

As a result of PitNET tissue transcript quasi-mapping, quantification and count summarization to gene level, expression profiles of 60230 genes were obtained, which after exclusion of low expressed genes (ten counts in at least seven samples) were reduced to 22669.

Differential expression tests between the subtypes of GH-producing PitNETs (N = 21) and non-functioning PitNET (N = 61) resulted in 1595 significant differentially expressed genes (L2FC >  ± 1.5 and p-adjusted < 0.05 thresholds). There was a similar count of upregulated (829 (54%), median L2FC = 2.1, IQR = 0.92) and downregulated DEGs (740 (46%), median L2FC = − 2.01, IQR = 0.85) in the GH-producing PitNET group (Fig. 2, Additional file 1: Table S2, Additional file 3: Fig. S1).

Fig. 2

Volcano plot representing the relationship between differential expression test (Wald test) p-values and shrunken (“apeglm” algorithm) transformed L2FC values. The dashed vertical lines represent absolute L2FC threshold of 1.5 and the horizontal dashed line represents p-value threshold of 0.05. Red points denote genes passing both thresholds, blue points represent genes not passing the L2FC threshold and grey points represent points not passing any of the mentioned thresholds. See full list of differentially expressed genes in Additional file 1: Table S2

Functional enrichment and protein–protein interactions in GH-producing PitNETs

Functional enrichment analysis of 955 (59.8%) out of 1596 genes returned 426 statistically significant results with FDR < 0.05 (Additional file 1: Table S3. and S4.). “GO Process” represented the largest number of detected enrichments with 212 entries (Additional file 3: Fig. S2). This analysis revealed significant and strong enrichment (Median strength = 1.08, IQR = 0.35) for seven pathways associated with growth hormone secretion and signalling (Additional file 1: Table S5., Fig. 3A). Strong enrichment was also detected for several different ion transport channels enriched for RYR1 and RYR3 DEG’s (Additional file 1: Table S6, Fig. 3B).

Fig. 3

A. Reduced STRING-db protein interaction network for growth hormone related significantly altered enrichment terms (Additional file 1: Table S5). B. Reduced STRING-db protein interaction type network for RYR1 and RYR2 associated molecule transport channels (Additional file 1: Table S6). Line thickness indicates confidence level

Other relevant interaction networks specific to GH-producing PitNETs in our data were related to ion channel functionality (“Potassium ion homeostasis”, “Sodium ion homeostasis”, “Voltage-gated calcium channel activity” with enrichment strength 0.8, 0.75, 0.71 respectively), calcium signalling (“Calcium channel activity”, “Calcium transport”, “Calcium channel” with enrichment strength of 0.65, 0.64, 0.64 respectively), neuronal and synaptic regulatory networks (“Regulation of postsynaptic density organisation”, “Neurexins and neuroligins”, “Synaptic membrane adhesion” with enrichment strength 0.89, 0.86, 0.74 respectively) and head structure development factors (“Pituitary gland development” and “Diencephalon development” with enrichment strength of 0.73 and 0.67 respectively).

A total of 722 protein–protein interactions were identified between nodes of genes which is more than twice expected at random sample (320 expected) providing protein–protein enrichment p-value < 1.0e-16. 1444 reciprocal interactions were detected (See also Additional file 3: Fig. S3, S4; Additional file 1: Table S7). 380 (40.2%) genes were excluded from gene enrichment and protein–protein interaction analysis due to missing information in the STRING database.

Previously encountered markers in GH-producing PitNETs

To test whether our data is an overall representation of the genetic alteration landscape of PitNET tissues, we curated a list of 428 genes associated with genetic changes of said tumour tissue types as found in literature [1]. As a result, 55 genes from the curated list were found to be differentially expressed in our data set, 29 of whom were upregulated with median L2FC of 2.31 (IQR = 1.4) and 26 of whom were downregulated with median L2FC of − 2.43 (IQR = 0.78) (Additional file 1: Table S8.). Genes like calcium binding protein 1 (CABP1), insulin like growth factor binding protein 5 (IGFBP5), dual specificity phosphatase 4 (DUSP4), paired like homeodomain 2 (PITX2), and regulator of G protein signalling 16 (RGS16) have previously been suggested as candidate markers for PitNETs [1]. Epidermal growth factor receptor pathway substrate 8 (EPS8) controls various cellular protrusions by regulating actin cytoskeleton dynamics and architecture as well as participates in growth factor activation, therefore promotes cell proliferation and cell survival within tumour. Overexpression of EPS8 has previously been detected in PitNETs compared to normal pituitary [42]. Major histocompatibility complex, class I, G (HLA-G) has also been detected in various types of PitNETs, primarily in lactotroph PitNETs suggesting tumour immuno-surveillance suggesting tumour immune-surveillance issue [1], (Fig. 4, Additional file 1: Table S8).

Fig. 4

Box plot visualisation of seven proposed candidate genes involved in altered pathway signalling matching with information found in literature for somatotropinoma and non-functioning PitNET subtypes. Gene expression levels transformed with variance stabilising transformation

Differential expression in preoperative SSA treatment group

Frequency filtering (ten counts in at least three samples), left 22348 genes available for further analysis. Parametric dispersion model was determined to have the best fit with this sample set. Differential expression test between GH-producing PitNET with SSA/DA therapy (N = 11) and without therapy (N = 9) initially produced 143 DEG’s with 108 (75.52%) upregulated and 35 (24.48%) downregulated genes in the therapy group. After L2FC shrinkage and results table subsetting with p-adjusted threshold of 0.05 and L2FC threshold of 0.58, a total of 95 DEG’s were obtained. 85 of these DEG’s were upregulated in the therapy group with a median L2FC of 2.68 (IQR = 1.67) and 10 were downregulated with a median L2FC of -1.11 (IQR = 0.63) (Fig. 5, Additional file 1: Table S9).

Fig. 5

A. (left) Volcano plot of differential expression results between PitNET SSA therapy and non-therapy groups with ‘apeglm’ transformed log fold change values. Dashed vertical lines represent absolute log2 fold change threshold of 0.58 and horizontal dashed line represents p-value threshold of 0.05. Red dots mark genes passing both thresholds. B. (right) Heatmap visualising gene level and mean normalised expression values of differentially expressed genes across all samples of the treated and untreated groups. Pink bar above the heatmap represents samples treated with SSA therapy, the green bar represents the untreated samples

Functional enrichment and protein–protein interaction results for SSA treated/untreated GH-producing PitNET comparison

By performing functional enrichment analysis, we obtained three terms associated with biological processes (Gene Ontology), eight with a cellular component (Gene Ontology) and six with annotated keywords (UniProt) with an average enrichment strength of 0.58 (SD = 0.28) DEG’s were revealed to be involvement in regulation of ion transport, generation of neurons and nervous system development with a protein–protein enrichment p-value of 0.0005 (Additional file 1: Table S10, S11, Additional file 3: Fig S5, S6). Enrichment was not possible for 16 (16.8%) DEG’s because of missing information or protein names in the string database.

Comparison of GH-producing PitNET SSA treated/untreated differential expression results in an unrelated cohort and a curated list of literature

Comparison of obtained DEG’s in SSA/DA treatment group to a curated list of PitNET relevant genes and two publications: Saksis et. al and Neou et. al. with the response variable description matching that of this publication: PitNET treatment with SSA therapy or lack of it [1, 11, 17]. Both publications had similar sample sizes for the treatment groups, respectively five for Neou et. al. and six for Saksis et. al, while the non-treated groups of acromegaly patients were 15 and six respectively. In total four DEG’s were found to be overlapping between mentioned sources of data. No overlap was detected between Neou et. al. and the current results. None of the overlapping DEG’s had a differing change of expression: all four were upregulated in the therapy group and downregulated in the untreated group. Notably, our discovery DE results for these four genes had smaller standard error values pointing to more evenly matched contrast group samples (Table 1, Fig. 6).

Table 1 DEG’s between SSA treated and untreated GH-producing PitNET groups replicated in an independent sample set and relevant literature

Fig. 6

Box plot graphs of four DEG’s from this study overlapping with DEG’s from Saksis et. al. 2021 of the same design and with a curated list of PitNET associated genes. Gene expression values transformed with variance stabilising transformation (Table 1).

Finally, to test whether and how effectively the SSA/DA therapy affects the expression of genes upregulated in the GH secreting acromegaly patients, we compared the differential expression results from both groups. Six genes were found to be overlapping (Table 2.). All of the matching differentially expressed genes had an increased expression in the SSA/DA therapy contrast, with the median positive L2FC difference of 3.49 (IQR = 2.79). Most of the observed genes are associated with cell membrane and cell adhesion.

Table 2 Six differentially expressed genes overlapping between SSA/DA comparison and GH/NF comparison groups

Transcriptomic patterns of pituisphere model

To detect more stable alterations in the pituisphere expression profile we first used time-series differential expression analysis with Likelihood Ratio Test between the control (non-incubated) and incubated with octreotide or cabergoline. groups. Thresholds (L2FC ± 0.58 and p-adjusted < 0.05) revealed 1941 differentially expressed genes for the octreotide and control comparison (Fig. 7A). Majority (1237, (64%) DEG’s were downregulated with a median L2FC of -2.35 (IQR = 2.33) while 703 (36%) were upregulated with a median L2FC of 1.92 (IQR = 1.96) in the incubated pituisphere group (Additional file 1: Table S12). For the cabergoline incubated and non-incubated time-series contrast, a list of 1991 differentially expressed genes was obtained (L2FC ± 0.58 and p-adjusted < 0.05 thresholds). 1186 (59%) DEG’s were downregulated with a median L2FC of − 3.89 (IQR = 3.97) and 805 (41%) were upregulated with a median L2FC of 3.74 (IQR = 4.08) in the incubated PitNET sphere group. (Fig. 7B, Additional file 1: Table S13).

Fig. 7

Volcano plot of LRT differential expression results between control and with (A) octreotide and (B) cabergoline incubated pituispheres. Dashed vertical lines represent absolute log fold change threshold of 0.58 and horizontal dashed line represents p values threshold of 0.05. Red dots represent genes passing both thresholds and grey dots represent dots not passing any of the thresholds

To investigate the more minute perturbations of pituisphere expression profile in response to the incubation length with octreotide or cabergoline, we performed differential expression analysis at each time point (4 h, 8 h, 24 h, 48 h and 72 h) using the Wald test and identified differentially expressed genes after subsetting the data with L2FC ± 0.58 and p-adjusted < 0.05 thresholds. For the octreotide contrasts, most DEGs were found in the PitNET sphere group incubated for 48 h. The majority of DEGs at each time point were downregulated except for the group treated for 72 h (54% of DEGs were upregulated) (Table 3, Additional file 3: Fig. S7–S11.; Additional file 1: Table S15–S19). When testing the incubation length contrasts for cabergoline, we observed that the largest number of DEGs was found in the PitNET sphere group incubated for 72 h. The majority of DEGs at each time point were downregulated except for the group treated for 4 h (65% of DEGs were upregulated) (Additional file 3: Fig. S12–S16; Additional file 1: Table S20–S24) (Table 3).

Table 3 Number of differentially expressed genes at each time point treated with octreotide or cabergoline

When comparing genes with altered expression from the time-series tests in both octreotide and cabergoline incubated sample groups, we observed that ryanodine receptor 2 (RYR2), marker of proliferation KI-67 (MKI67), and collagen type VIII alpha 2 chain (COL8A2) were downregulated (Fig. 8).

Fig. 8

Differentially expressed time stable downregulated genes in the pituispheres incubated with octreotide (A), and in the pituispheres incubated with cabergoline B. Gene expression values normalized with variance stabilizing transformation

For the three mentioned DEG`s, 59 time-series enrichments were detected in the octreotide incubated group (Additional file 1: Table S25) from the 163 observed in total (Additional file 1: Table S26), and 91 time-series enrichments were detected in the cabergoline incubated group (Additional file 1: Table S27) from the 263 observed in total (Additional file 1: Table S28). We also observed that Major Histocompatibility Complex, Class I, G (HLA-G) and Solute Carrier Family 2 Member 1 (SLC2A1) genes were downregulated in the presence of octreotide incubating for 4 h and 8 h and for 24 h and 48 h, respectively. Incubating pituispheres with cabergoline for 4 h, 24 h, 48 h, and 72 h (except 8 h) we observed ADP Ribosylation Factor GTPase Activating Protein 1 (ARFGAP1) gene downregulation. When incubating pituispheres with cabergoline for 8 h and 24 h, expression alterations of Transforming Growth Factor Beta Receptor 2 (TGFBR2) and Epidermal Growth Factor Receptor Pathway Substrate 8 (EPS8) were noticed, these pathways have been elaborately studied in relation to PitNETs [43,44,45,46,47,48,49,50,51,52]. The significance of cadherins in PitNETs was already made clear in the early 2000s [53, 54]. Alterations in the ability of cells to adhere and interact with neighboring cells and extracellular matrix are proven to influence tumour progression and strongly correlate with its aggressiveness [53,54,55]. All previously mentioned genes were alluded to in the context of PitNETs and are considered as candidate markers for PitNET tumourigenesis.

Comparison of transcriptome expression profile similarity between matched, untreated samples of PitNET tissue, pituispheres and mesenchymal stromal stem-like cells

To evaluate whether MSC could be used as a potential model for SSA therapy effects on GH-producing PitNETs for future studies, we calculated the Euclidean distances on variance stabilising transformation transformed gene level counts for five PitNET samples involved in the tissue, pituisphere and MSC experiments, resulting in 15 total samples. Euclidean distances were calculated for 500 most variable genes amongst all three sample group origins after filtering genes for frequency to increase gained information. Afterwards the PCA method was applied to said counts to determine the similarity by transcriptional profile distances of mentioned sample groups. Between the compared groups, transcription profile for MSC samples was more dissimilar to PitNET tissue samples, with an explained variance for the first principal component of 49.17%. On the other hand, pituisphere samples overlap PitNET tissue samples almost completely suggesting that they are more similar than MSC samples. Nevertheless, some variation can be observed in the pituisphere group mostly along the second principal component, which explains 20.36% of the variation in the data. Genes with the largest loading scores across principal component 1 are associated with cell extracellular matrix, ion binding, organogenesis and tissue differentiation, while genes with the largest effect on differences between pituisphere and normal tissue are from the lncRNA class (Fig. 9).

Fig. 9

PCA bi-plot for the Euclidean distances of 500 most variable genes amongst PitNET tissue, pituisphere and MSC sample groups. The ellipses represent a 95% confidence level of the multivariate t distribution for each group. Colour indicates biological samples from the same PitNET patient from all three sample origins. Ten gene symbols with the largest loading scores across principal components 1 and 2 are also shown (ENSEMBL notation used for genes with no HGNC symbol). MSC—mesenchymal stromal stem-like cells

Multiple time point evaluation of octreotide incubated GH3 cells

To investigate whether the effects of octreotide therapy could be observed and replicated in an animal model which also expresses GH1 hormone similar to GH-producing PitNETs, we performed a differential expression test (Wald test) between GH3 cells incubated with octreotide at different time points (4, 8, 24) against unincubated GH3 cells. The test yielded 39 DEG’s (p-adjusted < 0.05) amongst all tested time points, 26 (66.6%) of whom were upregulated and 13 (33.4%) of whom were downregulated (Additional file 3: Figs. S17–S19, Additional file 1: Table S29–S31). At least 60% of DEG’s from each time point comparison were upregulated, gradually increasing along with the incubation time. The number of genes with altered expression also increased along longer incubation periods, indicating that the effects of octreotide incubation not only lasts but, furthermore, increases for at least 24 h (Table 4).

Table 4 Summarization of DE results for multiple time point comparisons of octreotide treated vs untreated GH3 cells (p-adjusted < 0.05)

Next, we checked for Homo sapiens orthologous genes and performed enrichment, and protein–protein interaction analysis for said genes of both organisms to determine pathways and molecular functions potentially affected by detected differentially expressed genes. Only the incubation period of 24 h had a statistically significant number of enrichment terms and protein–protein interactions in Rattus norvegicus, while there were no significant interactions in the human orthologues in any of the time points (Table 5, Additional file 3: Fig. S21–S22; Additional file 1: Table S32–S33).

Table 5 Results of functional enrichment and protein–protein interaction for differentially expressed genes in the octreotide incubated GH3 cells in comparison to control in both Rattus norvegicus and Homo sapiens

Comparative pathway analysis

Finally, we performed a pathway analysis by determining which significantly affected pathways were common in each of the characterised groups and subsequently searching for matching DEG’s within the matched pathways as possible time stable drivers for changes associated with PitNET therapy. First, we compared functional enrichment results from the PitNET tissue therapy sample group with the cabergoline incubated PitNET pituispheres. From the 17 enriched pathways of tissue SSA therapy comparison and 939 pathways of the combined cabergoline incubation time point comparison results (Additional file 1: Table S34–S38), nine matching pathways were detected, six of which were categorised as ‘GO Process’ and three as ‘GO Component’. The number of overlapping pathways slightly increased with the incubation time, suggesting that long lasting alterations associated with SSA therapy are established after 24–72 h. Strength for the overlapping pathways was higher in the tissue therapy group (mean 0.52 ± 0.13) compared to the combined cabergoline group (0.2 ± 0.05), indicating that while the long-term effects can be determined after 24–72 h, the number of genes with altered expression increases after application of prolonged therapy (Table 6).

Table 6 Overlapping significantly enriched pathways (FDR < 0.05) and genes (L2FC >  ± 0.58, p-adjusted < 0.05) between PitNET tissue therapy and Cabergoline multi-time point incubation pituisphere therapy results

Within the overlapping pathways of tissue therapy and cabergoline incubation contrasts, three distinct genes were found to be driving described changes—RGS4 (Regulator of G protein signalling 4), EEF2K (Eukaryotic elongation factor 2 kinase) and KCNH5 (Potassium voltage-gated channel subfamily H member 5). RGS4 is consistent in its expression upregulation across tissue and pituisphere samples at different octreotide incubation periods, unleash EEF2K and KCNH5, which are upregulated in SSA incubated tissue samples, but downregulated in the octreotide incubated pituispheres.

On the other hand, when comparing PitNET tissue SSA therapy functional enrichment results (strength 0.63, FDR 0.05) with pituispheres combined octreotide incubation time point (strength 0.41, FDR 0.04) enrichment results (Additional file 1: Table S39–S42), only one pathway was found to be in both groups—‘Presynapse’, which specifically was detected at the octreotide 24 h incubation time point. Again, similarly to comparison of tissue SSA therapy and pituisphere cabergoline enrichment results, a larger number of DEGs was driving the enrichment in the tissue therapy group, supporting the assumption that a longer period of SSA therapy results in a larger number of genes with altered expression.

Pathway comparison of PitNET tissue SSA therapy enrichment results was also performed for the multiple time point results of octreotide incubated GH3 cells. Due to the low number of differentially expressed genes in the GH3 sample group, with even fewer identified as human orthologues none of them overlapped with PitNET tissue SSA therapy enriched pathways.

Discussion

In this study we investigated the transcriptomic landscape of GH-producing PitNETs to identify characteristic markers and follow alterations of the expression profiles induced by clinically used drugs in tumour tissue and available cell models. We have shown that GH-producing PitNETs have transcriptomic profiles with distinct perturbed growth hormone related pathways consistent to its functional status alongside changes in inner cell signalling, ion transport, cell adhesion and extracellular matrix characteristic pathways. Furthermore, we have provided additional evidence that pituispheres of primary cultures have higher similarity to the actual tumour tissue compared to MSCs which would limit MSC usability as a precise model for PitNET studies. Using the pituispheres we elucidated that treatment regimens (octreotide and cabergoline) affect specific cell proliferation gene expression (MKI67) and expression of members of core functionality pathways (RYR2, COL8A2, HLA-G, SLC2A1, ARFGAP1, TGFBR2). While using commercial GH3 cell line we observed that medication did not have transcriptomic effects similar to preoperative treatment in PitNET tissue or pituisphere model.

During the GH-producing PitNET transcriptome analysis we were able to identify 1595 DEGs that distinguished this group from non-functional PitNET, the expected result was that growth hormone related pathways were upregulated in GH-producing group, but alongside that we also detected other candidates among these only 41 have been mentioned in transcriptomic studies of PitNETs, therefore, there are distinct patterns that characterise GH-producing PitNETs that could be used for better understanding of tumour biology and development of management strategies [1]. We observed a large share of candidate gene expression implicated in ion flux regulation and inner cell signalling which is in concordance that these cells are metabolically active and producing large amounts of molecules for secretion.

Noteworthy, our pituisphere experiments showed that incubation in octreotide and cabergoline reduce expression levels of Ki-67 (MKI67) which is one of widely used markers in histology to assess tumour proliferation capacity therefore also aggressiveness [56]. This indicates that the pituisphere model is also affected by treatment in accordance with tumour biology in-vivo. Ki-67 has been previously downregulated in pre-operatively treated patients’ PitNETs [11]. On the other hand, we did not observe any significant downregulation of Ki-67 in medication treated GH3 cell line, for which further confirmation is needed but it indicates that pituispheres have similar response to drugs as PitNETs themselves.

In both modalities of pituisphere model (incubated in octreotide or cabergoline) we observed downregulation of RYR2, also slight upregulation of RYR1 and RYR3 was observed in GH-producing PitNET analysis (Additional file 1: Table S2). Altered expression of this factor has been found before in several transcriptomic studies of PitNET gene expression [51, 57, 58]. While another member of the ryanodine receptor family: RYR1 has been reported as carrier of somatic variants in at least three PitNET studies [21, 59, 60]. Ryanodine receptors (RYRs) are Ca²⁺ intracellular channels located in the endoplasmic reticulum (ER) being one of main triggers for signalling transduction coupled to calcium release from ER. These receptors have been implicated in the development of various human diseases with vascular or neuronal components—heart failure, arrhythmia, myopathies and neurodegenerative disorders [61, 62]. Calcium signalling has been previously linked to crosstalk with other signalling pathways of PitNET tumours therefore further studies of RYR family involvement in pathogenesis needs to be carried out [63, 64]. RYRs have been implicated in gonadotropin-releasing hormone signalling where in rat models the mRNA levels of RYRs have been found to regulated by gonadotropins, this opens novel interesting investigation window also for PitNET functional studies [65].

Several overlapping candidates related to neuronal system and ion channel regulation were observed in a group of cabergoline incubated pituispheres and pre-operatively treated patients’ PitNET (Table 6). As expected, we observed expressional changes of factor regulator of G protein signalling 4 (RGS4) that act by activating GTPase activating proteins (GAPs) that are crucial triggers in transmitting signals via G protein coupled receptors (GPCRs). This is not surprising as octreotide is targeting SSTRs and cabergoline D2R both belonging to GPCRs, treatment have been previously already reported to downregulate SSTR expression and our study helps to highlight specific factors involved in intracellular adaptation loop [66, 67]. Previously, only RGS16, also an RGS family member, has been mentioned in transcriptomic studies of PitNETs [68, 69]. Downstream signalling of SSTRs upon SSA stimulation have been widely studied in relation to potential SSA drug resistance mechanisms and primary roles of beta-arrestins and cytoskeleton protein filamin A have been investigated, however, the role of RGS candidates involvement in drug induced feed-back mechanism could also be investigated [70, 71]. We noticed slight upregulation of RGS7, 8, 16 and 17 in GH-producing PitNET transcriptomic results (Additional file 1: Table S2), but how this could be related to GH-producing functionality needs to be investigated.

Other overlapping factors found both in pituispheres and tissue analysis (Table 6) were eukaryotic elongation factor 2 kinase (EEF2K) and potassium voltage-gated channel subfamily H member 5 (KCNH5). According to Human Protein Atlas EEF2K is a protein kinase involved in downregulation of translation elongation, many studies have indicated the role of EEF2K in tumourigenesis and the role of EEF2K as anticancer drug target have been discussed, but so far, we have not found previous evidence of its involvement in PitNETs [72].

KCNH5 belongs to the family of ion channels that are involved in neurotransmitter and hormone release, that could be functionally linked to cell signal transduction pathways related to PitNET medication use. KCNH5 is predominantly expressed in various parts of the brain, but in smaller levels it has been detected also in adrenal and pituitary glands [73]. It has been demonstrated that in Cxcr2 knockout mice pituitary functions are altered and KCNH5 related pathways are downregulated, but how this could be related to PitNET development or response to therapy needs to be investigated [74].

Our data strongly supports evidence that extracellular matrix (ECM) is involved in molecular patterning of PitNET functionality and response to medical treatment. There is previous evidence that PitNETs can contain fibrotic tissue composed of collagens and alteration on collagen composition can have antitumoural effects [75,76,77]. It has also been demonstrated that tumours from patients pre-treated with SSA have softer tumour consistency it is not clear whether this could be related to the changes in ECM caused by medication [78, 79]. We have already previously observed collagens and other ECM factor modulation in PitNETs treated with SSA/DA [17]. In this study, we found significant downregulation of COL8A2 in octreotide incubated pituispheres, but in previous reports and our extended results here, COL8A2 was upregulated in patients with preoperative therapy. This indicates some consistency in tissue data but discrepancy with pituispheres, that could be explained by regulation of ECM also by inner biological environment of the tumour that is different from pituisphere propagation media (Table 7).

Table 7 Matching differentially expressed genes in detected significant enrichment terms between the PitNET tissue SSA therapy and combined cabergoline incubated PitNET pituisphere sample groups (DEG’s: L2FC >  ± 0.58, p-adjusted < 0.05; Enrichment terms: FDR < 0.05)

While we observed some overlap between our discovery cohort, [17] and the literature, the lack of co-occurring differentially expressed genes with the validation cohort from Neou et. al. is not surprising, given the fact that the compared subsample groups of SSA/DA therapy and non-therapy are of different sizes, resulting in differing statistical power, while such factors as sequencing technology, laboratory and sample level related batch effect could also play a major role [80].

We also wanted to compare drug effects on transcriptomic patterns of widely used rat cell model GH3 to pituispheres and tumour tissue of pre-operatively treated patients. To our surprise we did not observe similar patterns suggesting that at least on transcriptomic level these cell models differ considerably, and every experimental procedure needs to take into account this possibility when translating and extrapolating results of each model system. For our study we studied GH3 cell line that is widely used in PitNET research and was available in our cell culture facility [19, 20]. It would be interesting to specifically assess also other cell line models by transcriptomic analysis and compare the obtained results. For example, it would be specifically valuable to evaluate transcriptomic landscape of GH4C1 that compared to GH3 cells produce lower level of GH but higher of PRL and initially have been derived from the same GH3 clone [19]. This could help to assess transcriptomic changes in context of several cell line models and could help researchers in future to select more appropriate investigation model. We also propose that in future single cell transcriptome sequencing of spatial transcriptomics could be used to dissect complexity of various PitNET cell types and deduct their impact on functional experiments. Our study along with others highlights the added value of using NGS for tracing functionality of acquired primary cultures and these novel techniques could bring a novel dimension to PitNET model studies.

We consider that the most significant limitation of this study was the use of non-functional PitNETs as contrast for deduction of GH-producing PitNET transcriptome patterns. According to WHO Classification published in 2017 detailed histological characterization of PitNETs in subgroups according to cell lineage development transcription factors (SF-1, Tpit, Pit-1) is advised [81]. As our sample set dates back to 2010 our sample group included samples without proper assessment of these characteristics, nonetheless we specifically investigated patterns of GH-producing PitNETs compared to hormonally non-functional and non-secreting tumour groups and concentrated on therapy effects altering gene expression. We are well aware that our non-functional PitNET group is most likely heterogeneous that could affect the presented results, however, other molecular factors for example presence of genetic variants in USP8 and GNAS could introduce distinct expression patterns in affected tumours and wider scope on influential factors could affect groups selection alongside WHO PitNET Classification system [11].

Conclusions

We demonstrate that GH-producing PitNETs have distinct transcriptomics not only on growth hormone signalling pathways, but also in GH-producing PitNET functionality supportive signalling and regulatory pathways: ion transport, extracellular remodelling and calcium signalling. Our data indicate that pituispheres more closely represent transcriptomic profiles of tumour tissue while MSC have significantly altered gene expression. We show that pituispheres treated with octreotide and cabergoline are not directly comparable with pre-operatively treated patients’ tumour tissue, but both share similar patterns of transcriptomic alterations. Therapeutic effects on model cell line GH3 did not match the gene expression changes observed in tumour tissue nor pituispheres. This study highlights the importance of cell transcriptomic profiling for correct model system selection and data interpretation that could be achieved in future by incorporating NGS methods and detailed cell omics profiling in PitNET model research.