Spatial transcriptomics reveal PI3K-AKT and metabolic alterations in aggressive, treatment-resistant lactotroph pituitary neuroendocrine tumors

Abstract

Clinically aggressive lactotroph pituitary neuroendocrine tumors (PitNET) are invasive tumors with an unusually rapid growth rate despite maximally tolerated doses of dopamine agonist (DA). We aimed to unravel the molecular heterogeneity of lactotroph PitNET and to identify biomarkers of aggressiveness and resistance to pharmacological treatment. A total of 13 patients harboring DA-resistant lactotroph PitNET were included in this study. Visium Spatial Transcriptomics (ST), whole transcriptome sequencing (WTS), and whole exome sequencing (WES) were performed in tumors from 4 of these patients; WTS and WES was carried out in 5; tumors from two patients underwent ST and WES and tumors from two other patients underwent only ST. Tumors were classified as null or partial responders according to their response to DA treatment. The eight PitNET analyzed by ST exhibited significant intratumoral heterogeneity, with clones showing alterations in PI3K/AKT and lipid metabolism pathways, particularly inositol phosphate, glycerophospholipid, and sphingolipid metabolism. The cell-cell communication analysis showed FGF-FGFR ligand receptor interaction whilst the transcription factors RXRA and CREM showed participation in both groups. A trajectory exploration was performed by including all PitNET together in a single analysis to determine whether there was a tendency or molecular pathway showing a differentiation pattern that would guide the transition from a partially responsive PitNET to a completely unresponsive one. We did not observe any such pattern. All of these findings were corroborated in the cohort of DA-resistant PitNETs in which only bulk WTS and WES were performed. The bulk WTS corroborated lipid metabolism and PI3K-AKT pathway alteration in PitNET, whereas the WES showed only SF3β1 and TP53 variants in one tumor each. Our work suggests that the PI3K/AKT pathway may constitute a molecular target at which to aim therapeutic strategies designed to treat aggressive and DA-resistant lactotroph PitNET.

Introduction

Pituitary neuroendocrine tumors (PitNET) are epithelial neoplasms arising from adenohypophyseal cells and constitute the second most frequent intracranial tumor [18, 19, 22]. Lactotroph PitNET (PRL-PitNET) are the most common type of functioning neoplasm of the pituitary gland and account for 50% of all pituitary tumors in both women and men [24]. These highly heterogeneous lesions can behave invasively in 40–45% of cases, with up to 15% being clinically aggressive [20]. Aggressive lactotroph PitNET are invasive tumors with an unusually rapid growth rate despite maximally tolerated doses of dopamine agonist (DA) and thus, frequently require adjuvant therapy, as they are at risk for persistence and recurrence [21, 24].

DA such as cabergoline suppress PRL production, reduce cell proliferation and constitute the pharmacological treatment of choice for lactotroph PitNET [15]. In case of persistence or recurrence, the therapeutic alternatives include surgical resection, radiation therapy, and in very aggressive lesions or when metastases are present, chemotherapy with the alkylating agent temozolomide, with or without capecitabine [24]. Predicting the response to any of these treatment modalities remains a challenge [26]. In this context, the clinical, pathological, and molecular definitions of malignant and aggressive lactotroph PitNET remain to be clearly defined, as primary lactotroph carcinomas are similar to aggressive adenomas except for the presence of craniospinal or distant metastases [21].

To ascertain the molecular heterogeneity of lactotroph PitNET and to identify biomarkers of aggressiveness and resistance to pharmacological treatment, we carried out spatial transcriptomics as well as whole exome and whole transcriptome sequencing of tumoral tissues from patients with aggressive lactotroph PitNET resistant to DA who underwent transsphenoidal surgical resection.

Materials and methods

Patients

A total of 13 patients harboring DA-resistant lactotroph PitNET were included in this study. Visium Spatial Transcriptomics (ST), whole transcriptome sequencing (WTS), and whole exome sequencing (WES) were performed in tumors from 4 of these patients; WTS and WES was carried out in 5; tumors from two patients underwent ST and WES and tumors from two other patients underwent only ST. One of the patients presented to the ear, nose and throat department with a mass in the sphenoid sinus, which was biopsied and turned out to be a PRL-secreting lesion documented by MRI, so she was treated with cabergoline. The remaining 12 patients had been subjected to transsphenoidal surgery because of lack of response to DA treatment.

Response to cabergoline treatment was defined as follows: complete response, normal PRL levels and > 50% reduction of tumor mass; partial or stable response, > 50% reduction in PRL levels and any reduction in tumor mass; and null response, less than 50% reduction in PRL levels and no tumor mass reduction. All tumors included in the study were sporadic and were collected from patients diagnosed, treated and followed at the Endocrinology Service and the Neurosurgical department of Hospital de Especialidades, Centro Médico Nacional Siglo XXI of the Instituto Mexicano del Seguro Social. All participating patients were recruited with signed informed consent and ethical approval from the Comisión Nacional de Ética e Investigación Científica of the Instituto Mexicano del Seguro Social in accordance with the Helsinki declaration.

Immunophenotyping of PitNET: immunohistochemistry (IHC) for hormones and transcription factors (TF)

Immunohistochemistry was performed as previously described [35]. Briefly, paraffin-embedded, formalin-fixed tissue blocks were obtained and 3 μm sections were stained with hematoxylin-eosin and reviewed by a neuropathologist. Tumors were represented with a 2-fold redundancy. Sections were cut and placed onto coated slides. Immunostaining was performed by means of the HiDef detection HRP polymer system (Cell Marque, CA, USA), using specific antibodies against each pituitary hormone (TSH, GH, PRL, FSH, LH and ACTH) and the lineage specific transcription factors (TF) TBX19 (T-PIT), POU1F1 (PIT-1) and NR5A1 (SF1), as previously described. Interpretation of IHC for pituitary hormones and TF was carried out by two independent observers.

Visium Spatial transcriptomics

Pathologist reviewed the hematoxylin and eosin-stained (H&E) slides and selected representative areas from each tumor. Visium Spatial Gene Expression for FFPE (10X Genomics) instructions were followed as provided by manufacturers. The FFPE tissue blocks were cut 5 μm-thick as recommended. Slices were placed onto Visium Spatial Gene Expression Slide, and the tissue was deparaffinized at 60ºC for 2 h and xylene baths followed by decreasing concentrations of ethanol and finally molecular grade water. Once the tissue was deparaffinized and re-hydrated it was stained with H&E. The tissue images were digitalized with the Aperio CS2 by Leica. Decrosslinking was performed with 0.1 N HCl as recommended. Libraries were constructed with the Visium Spatial Gene Expression Reagent Kit as follows, FFPE probes were hybridized overnight and subsequently were ligated, the RNA was digested, and the probes released for extension and eluted for amplification and index ligation and finally library cleanup. The libraries were pair-end sequenced in a NextSeq2000.

Visium bioinformatic analysis

The bioinformatic analysis was conducted using Ubuntu 22.04.5 LTS, based on Linux. Statistical analyses and figure constructions were carried out using R version 4.4.1, unless otherwise specified as in the following sections.

Data preprocessing was performed using Space Ranger 2.1.1 with GRCh38 as the reference genome. For further analysis, the Seurat package was used, whereby data was normalized with the “SCTransform” function using default parameters on the “Spatial” assay, followed by dimensionality reduction and clustering. Variable features were calculated using the “FindAllMarkers” function and “FindSpatiallyVariableFeatures” for spatial variogram variable features.

Cell-cell communication analysis was conducted using CellChat v. 2.1.2 with the full CellChatDB.human database. Communication probability/strength was calculated using the “computeCommunProb” function with default parameters, overexpressed genes and filtering of overexpressed interactions, and network centrality was also calculated by means of “netAnalysis_computeCentrality.”

Enrichment analysis was performed using the “clusterProfiler” package, v. 4.12.6. The “enrichKEGG” function was used for enrichment with default parameters, and results from “FindAllMarkers” were filtered by absolute avg_log2FC > 0.5 and an adjusted p value of < 0.01. Term similarity and figures were generated using the same package. Figures over tissue images were constructed by adding an assay to individual Seurat objects.

The “monocle3” package, v. 1.3.7, was used for trajectory analysis on the aggregate, utilizing its functions “as.cell_data_set,” “cluster_cells,” “learn_graph,” and “plot_cells,” as default parameters.

Gene regulatory network analysis using “Scenic” required transforming the Seurat object to loom using the “build_loom” function from “SCopeLoomR” v. 0.13.0. Loom files were then processed with the “grn,” “ctx,” and “aucell” functions from “pyscenic” v. 0.12.1, using the following references: “hg38__refseq_r80__500bp_up_and_100bp_down_tss.mc9nr.genes_vs_motifs.rankings.feather” and “hg38__refseq-r80__10kb_up_and_down_tss.mc9nr.genes_vs_motifs.rankings.feather.” Python v. 3.12.14 and the “binarize” function from “pyscenic” were used to construct the binary matrix.

Analysis of metabolic activity was carried out by means of the “sc.metabolism.Seurat” function of the “scMetabolism” package v. 0.2.1, using the “AUCell” method and KEGG metabolism type. The results were added to the original Seurat objects as new assays for figure construction.

Velocity analysis (gene transcriptional patterns) was performed using “velocyto.R” v. 0.6 and “SeuratWrappers” v. 0.3.5, with the “RunVelocity” function. Figures were constructed with “show.velocity.on.embedding.cor” using standard parameters.

For aggregate analysis, the Seurat-proposed method was used with the “merge” function on raw assays. Normalization, dimensionality reduction, clustering, and visualization were performed using the same methods as for individual assays.

DNA purification

Pituitary tissue was lysed in proteinase K solution. After lysis, 300 µL of 5 M ammonium acetate was added to precipitate proteins and cellular components. The aqueous phase was transferred to a fresh tube and 600 µL of isopropanol were added and the mixture incubated overnight at -20 °C. The mixture was then centrifuged at 14 000 rpm for 30 min. The resulting DNA pellet was washed with 1 mL 75% ethanol and centrifuged at 10 000 rpm for 5 min; the pellet was air-dried, and the DNA resuspended in nuclease-free water.

RNA purification

Total RNA was extracted from tumor tissues using the miRNAeasy Mini Kit (Qiagen Inc, CA, USA) according to manufacturer’s instructions. Tissue samples were disrupted and homogenized in 700 µL Qiazol Lysis Reagent. They were then incubated at room temperature for 5 min. Next, 200 µL of chloroform were added, and samples were incubated at room temperature for 3 min. The mixture was centrifuged at 12 500 rpm for 15 min at 4 °C. The aqueous phase was transferred to a fresh tube and mixed with an equal volume of 70% ethanol. Samples were then transferred to an RNAeasy Column in a 2 mL tube, and centrifuged at 10 000 rpm for 15 s. After centrifugation, 700 µl of RW1 buffer was added and the mixture was centrifuged at 10 000 rpm for 15 s. Flow-through was discarded and 500 µL of RPE buffer were added to the membrane and then centrifuged at 10 000 rpm for 15 s and a second time for 1 min. The column was transferred to a new collection tube adding 30 µL of RNAse free water and centrifuged for 1 min at 10 000 rpm. RNA was quantified using a Nanodrop-ND-1000 spectrophotometer (Thermo Scientific, DE, USA); RNA integrity was evaluated by a Bioanalyzer 2100.

Whole RNA sequencing

RNA integrity of each sample was assessed using a R1 RNA Cartridge for QSep 400 (BiOptic, New Taipei City, Taiwan), RNA concentration was measured with the Qubit RNA HS Assay Kit (Invitrogen, Carlsbad, CA, United States) and purity was analyzed with a NanoDrop 1000 (Thermo Fisher Scientific, Wilmington, DE). Transcriptome libraries were prepared with the TruSeq Stranded Total RNA Library Prep with Ribo-Zero Gold (Illumina, San Diego CA, United States). Fragmentation times were adjusted based on RIN. Transcriptome libraries were quantified with Qubit dsDNA HS Assay Kit (Invitrogen, Carlsbad, CA, United States), library size was analyzed in S2 Standard DNA Cartridge for QSep 400 (BiOptic, New Taipei City, Taiwan), and sequencing was performed in a NovaSeq 6000 sequencer (Illumina, San Diego CA, United States) in a 150 bp pair-end configuration.

Whole exome sequencing (WES)

Genomic DNA (gDNA) was shipped to the Genomics Core Lab of the Instituto Tecnológico y de Estudios Superiores de Monterrey for exome sequencing. gDNA was quantified using Qubit dsDNA BR Assay Kit (Invitrogen, Carlsbad, CA, USA). Quality was determined spectrophotometrically using a Nanodrop One spectrophotometer (Thermo Fisher Scientific, Waltham MA, USA). WES libraries were prepared using Illumina DNA Prep with Exome 1.0 Enrichment (Illumina, San Diego CA, United States). All libraries were quantified with the Qubit dsDNA BR Assay Kit (Invitrogen, Carlsbad, CA, USA), libraries size was analyzed in S2 Standard DNA Cartridge for Sep 400 (BiOptic, New Taipei City, Taiwan), and sequencing was performed in a NovaSeq 6000 sequencer (Illumina, San Diego CA, United States) in a 150 bp pair-end configuration.

FastQC and preprocessing

Quality assessment of spatial transcriptomics, RNAseq and exome libraries was carried out using FastQC (Babraham Bioinformatics) to determine the quality of the sequencing. All raw sequences passed the initial quality filter. Adapters were removed and a quality and length filter were performed with Trimmomatic 0.40.

Differential expression

Preprocessed reads were aligned using STAR against the human reference sequence (hg38). Once the BAM files were obtained, HT-Seq package was used to estimate gene expression whereby a count table is obtained that can be used to perform differential expression analysis. Statistical analysis of differential gene expression (DGE) among the groups was performed using the DESeq2 R package, version 2.13 (R Foundation for Statistical Computing, Vienna, Austria). The false discovery rate was set at (FDR) < 0.01 and a threshold normalized absolute log 2-fold change > 1.0.

WES computational analysis

Preprocessed sequences were aligned to the human reference sequence (hg38) using the Illumina-Dragen Enrichment pipeline (llumina, San Diego CA, United States). This pipeline was set to produce copy number variants (--enable-cnv true). The BAM files resulting from the enrichment were removed from PCR duplicates using Picard Tools (http://broadinstitute.github.io/picard). Each BAM file was used to obtain the somatic variants using the GATK pipeline, and variants were annotated using ANNOVAR based on the following databases: Clinvar, gnomAD, refGene, cytoBand, exac03, avsnp147, dbnsfp30a. The somatic variants were then transformed to MAF using Funkotator from GATK. Additionally, converted annotated variant files were analyzed with Maftools package from R programming language to visualize the landscape of critical mutations.

The mutation analysis was carried out using maftools, from this package the merge_mafs tools were used to combine samples. The mutations were filtered using subsetMaf with the parameter “Variant_Classification = = Missense_Mutations”. Graphs of mutations in genes of interest were constructed using lollipopPlot to observe mutations in general and lollipopPlot2 to compare mutations by study subgroup.

Results

Clinical and demographic characteristics of the patients

Table 1 summarizes the demographic, clinical and biochemical characteristics of the patients included in the study. Median age was 25 years (IQR 22–50), 9 were women and 4 men. Median PRL at diagnosis was 538 ng/mL (IQR 193–944), median maximum tumor diameter at diagnosis was 25 mm (IQR 16–39), median PRL after DA treatment was 698 ng/mL (IQR 419–1972), and median maximum cabergoline dose was 3 mg/week (IQR 3–4). The Ki67 proliferative index varied widely among tumors, from 0 to 20% and did not correlate with tumor size or invasiveness. According to our definition, only three of the 13 patients were categorized as being partially responsive to DA treatment, the remaining 10 were considered resistant. Spatial transcriptomics was carried out in a total of 8 patients, whereas bulk WTS and WES were performed in 9 and 11 patients respectively. A total of four patients underwent ST, WTS and WES; two underwent ST and WES, five WTS and WES and two ST only.

Table 1 Demographic, clinical and hormonal features of patients

Spatial transcriptomics shows heterogenous transcriptome within lactotroph PitNET

A total of 8 unique lactotroph PitNET from 8 different patients were evaluated by means of Visium (10X Genomics) spatial transcriptomics. From spatial transcriptomics we sequenced a total of 11,054 spots with 107,473 mean reads per spot. Spot deconvolution was carried out by both, gene expression and by manual cell identification of immunohistochemistry slides evaluated by two expert neuropathologists. In spots where the two predictions juxtaposed the corresponding cell type was tagged.

All tumors analyzed were composed by tumor cells (PRL, POU1F1), endothelial cells (VWF, CD34), immune response cells, particularly macrophages (PTPRC, CD68) and stem cells (SOX2, LGR4 and RBPMS). Alongside, we carried out Ki67 staining to evaluate proliferative index. The three partial responders showed 0%, 1% and 10% Ki67 staining positive cells. Whereas the null responders 0%, 2%, 10%, 10% and 20% Ki67 staining positive cells, showing high heterogeneity among all resistant tumors.

All lactotroph PitNET included in the ST study showed a high intratumoral cellular heterogeneity. Each individual tumor was found to have 3 to 10 cell clusters comprising the tumor mass. One of the three patients categorized as having a partial response to cabergoline harbored a tumor with only three cell clusters, whereas one of the five null responders had tumors consisting of up to 10 cell clusters (Fig. 1).

Fig. 1

Panels a-h) Different cell clusters that comprise tumor mass and the cell deconvolution of each tumor analyzed by spatial transcriptomics. They are organized by their response to pharmacological therapy as completely unresponsive and partially responsive. Tumor complexity increases as response to therapy worsens. Tumor cluster range from 3 to 10 cell clusters. The tumor that was treated primarily with cabergoline showed the lowest number of clusters (a), the partially responsive tumors showed five and eight cell clusters and finally the completely unresponsive tumors showed five, six, seven, eight and ten cell clusters. Tumors are mainly comprised by tumor cells, stem cells, macrophages and endothelial cells, as noted by gene expression deconvolution, H&E analysis and immunohistochemistry

Most of the spots clustered together regardless of the tissue of origin or whether the tumor came from a partially or totally resistant patient (Fig. 2a-d). The pathway enrichment analysis showed that both, lactotroph PitNET from null and partially responsive patients presented altered gene expression in pathways such as PI3K − AKT signaling (PIK3C3 p = 1.599e-07 and AKT3 p = 1.382e-170), calcium signaling (GRM5 p = 4.075e-40 and ADCY1 p = 1.140e-11), AGE − RAGE signaling in diabetic complications (NOX4 p = 2.453e-43 and PRKCZ p = 1.722e-31) and MAPK signaling (IGF1R p = 0.000 and MAPK1 p = 1.278e-17) (Fig. 2e-f, Supplementary Fig. 1) among others. Partial responders showed a slightly different enrichment profile with proteoglycans in cancer and miRNAs in cancer pathways (Fig. 2e-f). Tumors from the three patients who were partially responsive to cabergoline treatment showed expression of CCN1 (p = 0.000), FFAR4 (p = 0.000) and IGFBP5 (p = 0.000). Tumors from the five nonresponsive patients showed expression of FNDC9 (p = 4.641e-222), SCGN (p = 1.196e-52), and DLK1 (p = 0.000) (Supplementary Fig. 2a-b).

Fig. 2

Panels a) and b) portray the UMAP and heatmap of the spots comprising all analyzed tumors clustered by sample, showing shared gene expression profiles as most samples group close together, whereas c) and d) depict the UMAP and heatmap clustering by therapy response, again showing shared gene expression profiles as most samples group close together. Panels e-f) show the altered pathways in null and partial responders, respectively

Since several of the altered pathways are related to metabolic traits, we performed a more detailed metabolic analysis. Interestingly, inositol phosphate metabolism was overrepresented in the null response tumors, both as a group and in the individual lesions. Null response lactotroph PitNET showed overrepresentation of glycerophospholipid, sphingolipid and pyruvate metabolism. Other altered metabolic traits in tumors from null responders included oxidative phosphorylation and the metabolism of nucleic acid precursors such as pyrimidines and purines, which positively correlated with tumor size. The tumors from the three patients with partial response to cabergoline showed overrepresentation of fatty acid degradation, ether lipid and glycerolipid metabolism, as well as pentose phosphate pathway, Krebs cycle and drug metabolism. All metabolic alterations observed as a group were corroborated in the individual lesions showing the presence of metabolically heterogeneous clones comprising the PitNET tumor mass (Fig. 3a-i).

Fig. 3

Panel a) Bubble plot of the altered metabolic pathways. Tumors from totally unresponsive patients showed alterations in inositol phosphate, glycerolipid, and pyruvate metabolism. Tumors from partially responsive patients displayed representation of fatty acid elongation, fatty acid degradation and ether lipid metabolism. The metabolism of nucleic acid precursors such as pyrimidine metabolism and pentose phosphate pathways are represented in both lesions. Panels b-i) show the spots with lipid metabolism from the null and partially responsive lesions

We then explored the potential ligand-receptor interaction of each group. Lactotroph PitNET from unresponsive patients showed overrepresentation of FGF9-FGFR1 and DLK1-NOCTH3, whereas in lesions from partially responding patients HSPG2-DAG1 and COL6A2-CD44 were among the predominant ligand-receptor interactions found. Interestingly, MIF-(CD74 + CD44) and APP-CD74 were present in all tumors (Fig. 4a-h). We explored gene regulatory networks by identifying transcription factors involved in the altered pathways. RXRA, RXRG and ASCL1, are represented in the null and partial responders lactotroph PitNET. CREM, ETV1 and HMX3 were also present in both, but with a slightly higher expression in lesions from partially responding patients (Supplementary Fig. 3).

Fig. 4

Panels a-h) show the cell-cell interaction networks among clusters from all analyzed PitNET, showing COL6A1-CD44, FGF9-FGFR1, HSPG2-DAG1, MIF-CD74, APP-CD74 and PTN-NCL

A trajectory analysis exploration was performed by interrogating all lactotrph PitNET together in a single evaluation to determine whether there was a molecular pathway showing a differentiation pattern that would guide the transition from the PitNET with a partial response to the PitNET with null response. We did not observe any such pattern. Velocyto analysis showed neither a differentiation pattern nor the presence of transitions in the lesions (Supplementary Fig. 4a), which agrees with the detailed morphological analysis carried out by our neuropathologists. To corroborate this result, we performed a second analysis using Monocle3 algorithm which yielded the same results (Supplementary Fig. 4b-c).

In the H&E-stained slides, the PitNET is composed of cells arranged in a solid architectural pattern. These cells exhibit a soft cytology, polygonal morphology, and moderate cytoplasm ranging from eosinophilic to pale eosinophilic. The nuclei are ovoid with granular “salt and pepper” chromatin and inconspicuous nucleoli. The cells are monotonous and display mild cytological atypia, with abundant small-caliber blood vessels with irregular shapes. None of the PitNET showed marked atypia, mitotic figures, necrosis, or signs of anaplasia (Supplementary Fig. 4d). This lack of morphologic and potentially molecular differentiation was independent of Ki67 expression (Supplementary Fig. 4e-f).

Whole transcriptome sequencing analysis corroborates lipid metabolism alterations found by Spatial transcroptomics

WTS was carried out in tumors from 8 patients, 4 of whom had also been subjected to ST. WTS results confirmed most of the ST findings. Fatty acid degradation, sphingolipid metabolism and phospholipase D signaling pathway are among the metabolism-related molecular events altered in these lesions, which all but one had been categorized as null response to cabergoline treatment (Fig. 5a-b, Supplementary Fig. 5a). To characterize metabolically these lactotroph PitNET we performed Metaflux analysis through transcriptomic means, and again, tumors showed a high metabolic heterogeneity, pointing to potential clone heterogeneity. Concordant with ST findings, most metabolic abnormalities were related to lipid metabolism the PI3K/AKT signaling pathway (Fig. 5c, Supplementary Fig. 5a).

Fig. 5

Panels a) and b) show the PI3K-AKT and fatty acid metabolism gene expression in bulk RNAseq from therapy resistant tumor showing expression in both null and partial responders. Panel c) Metaflux analysis results showing mostly alteration in lipid metabolism

The most frequent molecular finding in patients with lactotroph PitNET resistant to cabergoline therapy is the reduced expression of dopamine D2 receptors (D2R) and low FLNA expression. In our study, the expression of genes encoding D2R and FLNA was found to be down regulated in both, the ST and bulk WTS experiments (Supplementary Fig. 2c-d).

Whole exome profiling

WES was carried out in 11 tumors, two from partially resistant and 9 from totally unresponsive patients. Few relevant genetic variants were identified in these lesions, the most prominent being in SF3β1 (c.C1873T: p.R625C). One PitNET showed a TP53 pathogenic variant (c.C817T: p.R273C) (Supplementary Fig. 5c), one tumor showed a variant of uncertain significance in DRD2 (c.A979G: p.K327E), and another one in PRLR (c.A298G: p.I100V). Mutational signatures were retrieved using COSMIC SBS database. The most represented signatures were SBS1 and SBS5, which can be considered as background of endogenous cellular damage of DNA; these two signatures are ubiquitously present across many cell types, both normal as well as tumoral. Interestingly, three other signatures, SBS6, SBS31 and SBS87, were represented in the therapy resistant lactotroph PitNET. The SBS6 is related to defective DNA mismatch repair, whereas SBS31 and SBS87 are related to prior platinum- and thiopurine-based chemotherapy (Supplementary Fig. 5b).

These SBS mutational signature could indicate that PitNET tumorigenesis could not be driven by genetic variants. We did not find any copy number alteration in our cohort. The calculated tumor mutation burden (TMB) was not significantly different in totally unresponsive patients compared to those who responded partially to the DA, and there were no correlations between TMB and clinical characteristics of the patients. Nevertheless, further research on DNA mutations is needed.

Discussion

Here we present the heterogeneous and complex biology and architectural cell structure of aggressive lactotroph PitNET resistant to pharmacological therapy. Leveraging spatial transcriptomics, we have described the spatial heterogeneity of lactotroph PitNET through the integration of their histological features and their molecular characteristics. We confirmed our findings by performing WTS and WES in two thirds of our tumor samples subjected to ST. Although somewhat expectedly, we found DA-resistant lactotroph PitNET to be made of transcriptionally heterogeneous clones with altered lipid metabolism and PI3K-AKT signaling alongside with other metabolic abnormalities that could potentially drive tumor aggressiveness and pharmacological resistance. Heterogeneity provides the fuel for resistance; therefore, an accurate assessment of tumor heterogeneity is essential for the development of effective therapies [2].

During the course of the disease, tumors generally become more heterogeneous. As a result of this intratumoral heterogeneity (ITH), the tumor mass might include a diverse collection of cells harboring distinct molecular signatures with varying degrees of sensitivity to pharmacological treatments. Also, as a results of ITH there is a non-uniform distribution of key molecularly distinct tumor-cell subpopulations across regions of the tumor (spatial heterogeneity). Many solid tumors adopt a branched pattern of evolution, whereas a linear pattern is invoked in phylogenetic depictions of certain hematological malignancies [2]. The findings of multiple studies demonstrate that higher levels of ITH predispose patients to inferior responses to anti-tumor therapies, including targeted agents [2]. The major contributors to intratumoral heterogeneity are (i) genetic variation, (ii) stochastic processes, (iii) microenvironment and (iv) cell and tissue plasticity [30]. Each of these factors impacts on drug sensitivity [30].

All normal cells share nearly identical wild-type genomes, which encode highly diverse phenotypic manifestations. Phenotypic diversity in normal cells reflects developmental processes triggered by responses to microenvironmental cues. Normal tissues are organized in functional and structural units with near-equal access to vasculature, which provides oxygen and nutrients while removing waste products [17]. Moreover, blood and lymphatic vasculature in tumors are disorganized with significant functional, spatial, and temporal heterogeneity [17]. In addition to diversification of contextual cues, structural and microenvironmental disorganization in carcinogenesis also creates relatively consistent new microenvironments [17]. Pre-existing ITH increases the likelihood that at least some tumor cells will survive therapy-induced elimination, while ongoing diversification of tumor cell phenotypes during treatment enables tumor cells to adapt to therapy-imposed selective pressures, leading to de novo resistance and eventual relapse [17]. ITH has been shown to be associated with poor outcome and decreased response to cancer treatment by multiple groups in multiple human cancer types implying a universal role in therapeutic resistance [17]. As observed in our data, the partially responsive lesions showed less ITH when compared to the completely unresponsive tumors. It remains to be elucidated whether ITH exists from the beginning, or it is induced by clone selection through pharmacological exposure. At least in some cases, phenotypic transition toward drug tolerance might be induced by treatment, although the ability to undergo phenotypic transition is likely to be limited to only certain subpopulations of tumor cells [17].

Plasticity between different signaling pathways, specifically the adaptive activation of bypass signaling, is a potential manifestation of temporal and spatial heterogeneity under therapeutic selective pressure [2]. Drug tolerance is often thought as a distinct, well-defined phenotypic state that might be defined by similar molecular underpinnings across multiple cancer types, such as dependency on certain metabolic pathways [17].

Partially responsive, and to a lesser extent totally unresponsive PitNET showed expression of CCN1 which is an extracellular matrix protein that does not contribute to matrix formation but rather participates in cell-matrix interactions and regulates survival, angiogenesis, adhesion, and chemotaxis in a cell-dependent manner [23]. Unresponsive tumors expressed GREB1, which is a glycosyltransferase that participates in the hexosamine biosynthetic pathway and temporally modified ER⍺ acting as is coactivator and implicated in hormone-dependent tumors such as breast, ovarian and prostate cancer [32]. Genes encoding both of these proteins are closely related to the PI3K/AKT pathway [7, 12]. Other molecular markers such as CCNB1has been previously correlated to recurrence and progression of lactotroph PitNET [26].

Cell clones from the unresponsive PitNET showed abnormalities in lipid metabolism and PI3K/AKT pathway. PI3K is a family of intracellular signal transduction enzymes capable of phosphorylating a hydroxyl group in the inositol ring of phosphatidylinositol, whereas AKT is a serine-threonine kinase [3]. It has been described that the PRL gene promoter could be regulated by the PI3K/AKT/CREB [8]. Furthermore, the PI3K/AKT/CREB pathway not only could influence mRNA expression but also it may have a role in PRL secretion [27]. Accordingly, the clones with alterations in the PI3K/ATK pathway could lead to an up-regulation of PRL gene expression. Not only PI3K/AKT pathway could lead to augmented PRL levels, but also it could drive apoptosis inhibition [4], cell migration, invasion and epithelial-mesenchymal transition in lactotroph PitNET [4, 6]. Previous studies have shown alterations in the PI3K/AKT pathway in PitNET, including lactotroph lineage [20, 29]. Evidence also suggests that PI3K/AKT pathway is one of the mechanisms by which cabergoline act in the responsive lactotroph PitNET inducing cell death by activating ERK, but also regulating autophagic events [28, 34]. The PI3K/AKT pathway contributes to occurrence and progression of tumor, but also to therapeutic resistance in several tumor types such as breast, ovary, lung, and hepatic cancer as well as in melanoma, activating survival signaling and limiting antitumor drug cytotoxic effects as well as favoring tumor stem cell characteristics [16]. PI3K/AKT pathway could mediate these effects by interacting with several other proteins such as MDM2, TSC2, GSK3, FOXO and mTOR and through these proteins modulate growth, survival, pharmacological therapy resistance and metabolism [9]. Through their interaction with GSK3, ACLY, SREBP and mTOR it can influence glucose and lipid synthesis and metabolism as well as the pentose phosphate pathway which are all required for ATP production [10]. We and others have previously described altered metabolic events, particularly those related to lipid metabolism, in aggressive and recurrent PitNET [35, 37]. In the present study, we found the same alterations in lipid metabolism in our DA-resistant lactotroph PitNET. Lipid metabolism can influence therapy resistance by reducing membrane fluidity based on the predominance of saturated fatty acyl chains in membrane lipids, which disrupts the drug uptake via passive diffusion and/or endocytosis [11]. It can induce formation of detergent-resistant membrane domains, which can activate membrane bound ATP-binding cassette (ABC) multidrug efflux transporters such as ATP-dependent translocase and finally, it can render cells less sensitive to toxic lipid peroxidation which can trigger apoptosis and ferroptosis, which occurs in response to the oxidative stress due to pharmacological therapy [11].

Among the activation cues of PI3K/AKT there is the FGF-FGFR signaling [36], which is also related to therapy resistance [33]. According to our ligand-receptor analysis, DA-resistant lactotroph PitNET display abnormalities in FGF-FGFR signaling which can potentially contribute to pharmacological resistance.

Finally, the RXR⍺ transcription factor gene regulatory network, which is known to interact with the PI3K/AKT signaling pathway, was represented in both the null and the partial responsive PitNET [31]. We also found an over-expression of CREM regulons, which is a cAMP response element modulator [25]; cAMP is related to synthesis and secretion of prolactin [5]. The interaction of PI3K/AKT pathway with lipid metabolism [13] could be targeted as a therapeutic strategy [14].

The potential use of molecules targeting PI3K/AKT pathway was previously described, NVP-BKM120 a pan-PI3K inhibitor and NVP-BEZ235 and PI3K/mTOR inhibitor which show apoptosis induction as well as cytostatic effects [1], also several other molecules targeting PI3K-AKT pathway have been used as an alternative treatment such as XL765 which also induced apoptosis are reported [3].

Conclusions

The understanding of the spatial structure of resistant PitNET may help the development of new targeted therapies. The intratumoral heterogeneity is more readily identifiable in pharmacological resistant PitNET than in partially responsive lesions with lipid metabolism alterations alongside with PI3K/AKT pathway alterations which may contribute to pharmacological resistance and therefore represent an attractive therapy target for these tumors.