Single Cell Metabolic Landscape of Pituitary Neuroendocrine Tumor Subgroups and Lineages

Pituitary neuroendocrine tumors (PitNETs) are common intracranial tumors comprising numerous subtypes whose metabolic profiles have yet to be fully examined. The present in silico study analyzed single-cell expression profiles from 2311 PitNET cells from various lineages and subtypes to elucidate differences in metabolic activities. Gonadotroph tumors exhibited high activities with histidine metabolism, whose activity is low in lactotroph tumors. Somatotroph tumors enriched for sulfur and tyrosine metabolism, while lactotroph tumors were enriched metabolism of nitrogen, ascorbate, and aldarate. PIT-1 lineage tumors exhibited high sulfur and thiamine metabolism. These results set precedence for further translational studies for subgroup/lineage specific targeted therapies.


Introduction
Pituitary neuroendocrine tumors (PitNETs, known previously as pituitary adenomas) are predominantly benign tumors of the pituitary gland, which may manifest with infertility, hyperprolactinemia, or Cushing disease, among other possible clinical manifestations. The clinical representation may vary due to functioning and non-functioning tumors [1][2][3]. The 2012-2016 CBTRUS (Central Brain Tumor Registry of the United States) Statistical Report, which provides data on primary brain and other CNS tumors, indicated that pituitary tumors represent 16.8% of all CNS tumors (n = 68020), of which 0.2% were malignant [4]. Current treatment guidelines have wavered in recent years owing to the many PitNET lineages and subtypes. For instance, non-functioning PitNETs -marked by lower Ki-67 index value and absence of hormone hypersecretory syndrome, compared to functioning PitNETs [5] -may benefit from active surveillance or, in some cases, evidence aggressive growth and resistance to drug therapies [6,7]. In TSH-secreting PitNETs, albeit comparatively rare, radiation therapy is complicated by resultant hypopituitarism, whereas medical therapy with somatostatin analogues and dopamine receptor agonists is highly effective [6]. When standard therapies fail to treat aggressive PitNET subtypes, administration of temozolomide (TMZ) is indicated, despite progressive disease in 26% of recipient PitNET patients (median 9 cycles, n = 116, p = 0.01) [8]. Suboptimal response to TMZ also underscores the need to explore alternative therapies for aggressive PitNETs [9,10].
Current World Health Organization (WHO) classification of PitNETs includes distinction based on the transcription factor involved and hormone expression in the tumor cells, while previous standards delineated by immuno-histochemical type, granulation subtype, and proliferative markers. For example, gonadotrophs are sparsely granulated PitNETs recognized by GATA3 immunostaining, whereas somatotrophs can be densely granulated and report positive immunostaining for low weight cytokeratins. PitNETs can also be classified according to the cell type-specific hormone genes and lineage-specific transcription factors. These include SF-1 (NR5A1) lineage comprised of gonadotrophs, T-PIT (TBX19) lineage comprising corticotrophs, PIT-1 (POU1F1) lineage comprising somatotrophs, lactotrophs, and thyrotrophs, along with null-cell and plurihormonal pituitary adenomas (PPA) [11][12][13][14]. A predominant fraction of pituitary adenomas is sporadic and generally consists of monoclonal origin [15]. Although certain recurrent somatic driver mutations have been found (GNAS in 40-60% of somatotroph tumors and USP8 in 40-60% of corticotroph tumors), most tumors contain no identifiable somatic mutations. Additionally, mechanistic links between copy number variations (CNVs) and pathogenesis have yet to be established [16][17][18][19].
The pathobiology of a large fraction of PitNETs remains poorly understood in part to the heterogeneity of tumor subtypes. To explore other avenues for novel therapies, it is critical to examine the pathobiology of these tumors, namely their metabolism. It is known that tumor metabolism greatly contributes to disease heterogeneity. Such heterogeneity can be seen in pituitary adenomas and cells from tumors of different subgroups may have unique metabolic demands shaped by their environments, ligand-receptor interactions, and transcriptional regulation factors [11][12][13][14]. At a cellular level, metabolism is dictated by metabolite fluxes and concentrations. Therefore, it is complicated to quantify these measurements in vivo.
However, the advent of single-cell gene expression profiling allows for an indirect method of examining cellular metabolism, especially within single cells across whole tumors [20]. Analyzing the single-cell transcriptome in lieu of the metabolome has been shown to accurately represent metabolite dynamics and has provided novel insights into cellular metabolism [21,22].
As the cellular metabolic landscape may contribute to pathogenesis and has yet to be specifically examined across different pituitary adenomas, the present study uses previously published single-cell sequencing data encompassing a multitude of pituitary adenoma subtypes and lineages to characterize inter-tumoral metabolic heterogeneity.

Expression Data
Previously published single-cell tagged reverse transcription (STRT) sequencing data derived from surgically resected pituitary adenomas from 21 patients was retrieved from the Genome Sequence Archive project PRJCA002946 [23]. For the present analysis, gene expression data from the entire population of malignant cells, as annotated in the original study, were used (n = 2311). To be consistent with the 2017 WHO classification used in the original study from which the data was described, 862 cells were from somatotroph tumors (n = 4), 327 cells derived from silent corticotroph tumors (n = 4), 289 cells from PIT-1-positive, plurihormonal pituitary adenomas (PIT l + PPA) (n = 3), 94 cells came from a pituitary adenoma with unusual immunohistochemical combinations (PAwUIC) (n = 1), 627 cells were from gonadotroph tumors (n = 8), and 112 cells were derived from a lactotroph tumor (n = 1) [23].
This article does not contain any studies with human participants performed by the author. Our research complies with the guidelines for human studies and was conducted ethically in accordance with the World Medical Association Declaration of Helsinki (also includes research on identifiable human material and data).

Preprocessing and Normalization
Gene expression values were gene length normalized to TPM (transcripts per million) and log 2 transformed. For genes associated with multiple transcripts, the longest transcript length was used. Metabolic gene lists were downloaded from the KEGG database (http:// www.kegg.jp/) for genes involved in metabolic pathways. Clustering visualization was performed in R version 4.0.4 using the 'Rtsne' package with perplexity set to 30 and using 20 dimensions retained in initial principal component analysis (PCA) step [24,25]. Deconvolution normalization using the scran package computed tumor subgroup-specific size factors [26]. Read counts were divided by size factors corresponding to tumor subgroups and then transformed back to TPM. Only genes with dropout rate <0.75 were used as reference genes for normalization to avoid noise from lowly expressed genes. Code is available at https://github.com/Batchu-Sai/PitNET_metabolism.

Pathway Activity and Metabolic Heterogeneity Analysis
The pathway score analysis from Xiao et al. was used with the inputted deconvolution normalized values [22]. The pathway activity score is defined as the relative gene expression values averaged over all genes in a specific pathway and all tumor subgroups or lineage cells of this type [22]. Statistical significance of pathway activities in specific subgroups was calculated by random permutation test, where subgroup and lineage labels were randomly shuffled 1000 times to simulate null distribution, and then compared pathway activity scores to original scores. To assess intra-subgroup metabolic heterogeneity, PCA was applied to the log 2 transformed TPM values without imputation. For each metabolic gene, a PCA score was derived from the sum of absolute values of loadings in the top principal components that account for at least 80% of the variance across cells in each subgroup. The genes were sorted according to these descending PCA scores, and gene set enrichment analysis (GSEA) was applied to these lists to identify pathways enriched with most variable genes. GSEA was applied with javaGSEA (http://software.broadinstitute.org/gsea/downloads.jsp).

Results
First analyzing the global structure through dimensionality reduction, pituitary adenoma subgroups were seen to cluster distinctly according to their subgroup type based on expression levels of 1566 metabolic genes used ( Figure 1A). Correspondingly, coloring the cells by their lineage shows a sharp clustering pattern ( Figure 1B).
To identify the overall landscape of metabolic pathways amongst each subgroup, a pathway activity score was calculated for each subgroup, which compares the relative gene expression value averaged over all genes in the specific pathway set and all cells of the specific subgroup [22]. From the 80 metabolic pathways analyzed, 77 pathways were more highly or lowly expressed in at least one subgroup, suggesting that most metabolic pathway activities are subgroup specific (Figure 2A). Distinct patterns of metabolic activities were also seen. Silent corticotroph tumors exhibited high taurine and hypotaurine metabolism along with degradation of ketone bodies. The pituitary adenoma with unusual immunohistochemical combinations, PAwUIC, reported uniquely high glycosaminoglycan biosynthesis and phosphonate and phosphinate metabolism. The tumor cells resected from PAwUIC also evidenced low neomycin, kanamycin, and gentamicin biosynthesis. Moreover, the collection of silent corticotroph tumor cells revealed upregulation (activity score >1) across 55% of the analyzed metabolic pathways, compared to the lactotroph (11%), PIT-1positive plurihormonal pituitary adenoma (11%), and somatotroph (6%) tumor cells.
The pathway analysis was performed again on the cells aggregated by lineage ( Figure  2B). From the 80 metabolic pathways analyzed, 67 pathways were more highly or lowly expressed in at least one of the three lineages, suggesting that certain metabolic pathway activities are also lineage specific but also that more are shared between lineages than between subgroups. T-PIT-positive tumor cells reported uniquely low expression of metabolic genes associated with histidine and thiamine metabolism. The T-PIT-positive tumor cells also exhibited distinctly high levels of nitrogen metabolism and biosynthesis of unsaturated fatty acids. Activity levels for these pathways reflected opposite directionality in the NR5A1 lineage: SF-1-positive tumor cells reported high histidine and thiamine metabolism, low nitrogen metabolism and biosynthesis of unsaturated fatty acids. PIT-1positive tumors showed unique enrichment for sulfur metabolism and folate biosynthesis, among several others. The collection of PIT-1-positive tumors revealed downregulation (activity score <1) across 50% of the analyzed metabolic pathways, compared to T-PITpositive tumors (10%) and SF-1-positive tumors (35%).
The described gene analysis supports the heterogeneity of the PitNET metabolic landscape, potentially targeting high-profile gene sets (i.e., HAL and HDC in histidine metabolism) for immunogenomic therapy of aggressive PitNETs. Notably, metabolic gene analysis identified high levels of histidine metabolism in gonadotroph tumors, and correspondingly low histidine metabolism in lactotroph tumors (Figure 2A). Somatotroph PitNETs enriched sulfur and tyrosine metabolism, while lactotroph tumors enriched pathways associated with metabolism of nitrogen, ascorbate, and aldarate. Metabolic profiling of PIT-1 lineage tumors discovered high, differential expression of genes whose protein products related to sulfur, propanoate, thiamine, and tyrosine metabolism, among others ( Figure 2B).
GSEA was done to determine which metabolic pathways were enriched in genes contributing to most variability within each subgroup and lineage. Analysis revealed that oxidative phosphorylation was the only significant metabolic pathway across all tumor subgroups and lineages, indicating its involved genes' unique contributions to PitNET intratumoral heterogeneity (Table S1, see https://www.genome.jp/pathway/map00190 for gene-level granularity).

Discussion
PitNETs are mostly slow-growing tumors of the anterior pituitary gland constituent of multiple subtypes. While specific PitNETs may be managed with dopamine and somatostatin receptor drugs, the long-term efficacy of these drugs for the treatment of PitNETs can be limited: 10% of prolactinomas resist cabergoline treatment and 30% of acromegaly cases resist octreotide treatment, for instance [27,28]. At the time of writing, aggressive PitNETs which do not lend themselves well to conventional therapies (i.e., surgery, chemotherapy, and radiotherapy) are treated with TMZ or TMZ concomitantly with radiotherapy [9]. Radiological response to TMZ therapies has shown promise [29], but a subset of PitNETs can develop TMZ resistance via MGMT-mediated repair of DNA adducts and tumors with high MGMT immunostaining are not expected to respond to TMZ [30]. The adduct produced by TMZ can also be recognized by mismatch repair which can correct it or potentially generate lethal DNA breaks. A working mismatch repair system is required for TMZ susceptibility; therefore, MSH6 mutations may confer resistance [30]. As such, the development of advanced molecular therapies which target PitNET subgroupand lineage-specific metabolic pathways is warranted, necessitating an understanding of the PitNET metabolic repertoire.
Taken in combination with previous studies showing that methotrexate sensitivity is affected by histidine catabolism [31,32], our finding that the gonadotroph PitNET subtype demonstrated increased levels of histidine metabolism relative to the lactotroph subtype may be of interest for future targeted therapy research (Figure 2A). Recent studies have shown that the expression of the LAT1 transporter involved in the metabolism of histidine and other amino acids may be regulated by pro-carcinogenic transcription factors [33], offering a possible mechanistic explanation for the dysregulated histidine metabolism found in our analysis. However, consideration should be given to the unequal number of cells available for the current analysis with 627 cells from gonadotroph tumors and 112 cells from one lactotroph tumor.
Observing metabolic pathway activities across the three, assessed PitNET lineages revealed additional distinctions ( Figure 2B). For instance, SF-1 positive tumor cells were enriched for thiamine metabolism and oxidative phosphorylation. Previous research indicated that thiamine supplementation modulates tumor proliferation, suggesting a role for intracellular thiamine in PitNET oncogenesis and aggressiveness [34,35]. The PIT-1 family of PitNETs evidenced high thiamine metabolism coupled with subtype-specific enrichment for sulfur, propanoate, tyrosine, and sugar metabolism, each of which whose dysregulation has been linked to tumor malignancy [36][37][38][39]. Previous transcriptome analyses of PIT-1 lineage PitNETs also identified differentially expressed genes associated with fatty acid and nitrogen metabolism [39,40], in addition to enriched pentose phosphate pathway activity and fatty acid degradation [41,42]. However, while Figure 2B reports distinct gene activity related to nitrogen metabolism (depletion), it does not evidence regulation of fatty acid metabolism, nor gene-level enrichment of the pentose phosphate pathway among assessed POU1F1 tumor cells. Continuing, pathway enrichment analyses of ACTH-secreting tumors discovered differential expression connected to drug, tryptophan and pyrimidine metabolism [40], corroborating metabolic analysis of TBX19 lineage PitNETs.
The present study is limited by the relatively small sample size for each tumor subtype. Additionally, functional corticotroph subtypes were not included in the analysis. It is also worth considering that the characterization of single cell population metabolism is currently restrained by the quantity of single cells that may be profiled at one time due to the innate noise associated with single cell gene expression profiling and cell type diversity. Although the described findings may therefore not be fully representative of PitNET subtype heterogeneity, the findings discussed here are valuable in order to start understanding the metabolic microenvironment of PitNETs at the single-cell level for targeted therapy applications.

Conclusion
The above-described findings offer novel insight into the metabolic landscape of pituitary neuroendocrine tumors and further validate single-cell sequencing for profiling tumor cells. Here we employed a benchmarked approach to a collection of 2311 PitNET tumor cells, which represent six subtypes and three lineages. Striking subtype-specific results include gonadotroph PitNET enrichment for histidine metabolism, which is correspondingly depleted in lactotroph PitNETs; lineage-specific results also include evidence of unique enrichment for sulfur metabolism in PIT-1 lineage tumors and unique depletion of histidine metabolism in TBX-19 lineage tumors. The authors hope the study findings will encourage future, metabolic exploration of treatment-resistant tumors and inform the adoption of alternative therapies for aggressive PitNET subtypes.