Zihao Wang1,2†
Mengqi Chang1,2†
Yanruo Zhang3,4†
Gang Zhou1,2†
Penghao Liu1,2
Jizhong Lou3,4
Yuekun Wang1,2Yuan Zhang1,2
Xiaopeng Guo1,2
Yaning Wang1,2
Xinjie Bao1,2Wei Lian1,2
Yu Wang1,2
Renzhi Wang1,2
Wenbin Ma1,2
Bing Xing1,2*Jun Gao1,2*- 1Department of Neurosurgery, Peking Union Medical College Hospital, Chinese Academy of Medical Sciences and Peking Union Medical College, Beijing, China
- 2State Key Laboratory of Complex Severe and Rare Diseases, Peking Union Medical College Hospital, Chinese Academy of Medical Science and Peking Union Medical College, Beijing, China
- 3Key Laboratory of RNA Biology, CAS Center for Excellence in Biomacromolecules, Institute of Biophysics, Chinese Academy of Sciences, Beijing, China
- 4University of Chinese Academy of Sciences, Beijing, China
Purpose: Pituitary adenomas (PAs) are the second most common intracranial neoplasms. Total surgical resection was extremely important for curing PAs, whereas tumor stiffness has gradually become the most critical factor affecting the resection rate in PAs. We aimed to investigate the molecular mechanisms of tumor stiffening and explore novel medications to reduce stiffness for improving surgical remission rates in PA patients.
Methods: RNA sequencing, whole-genome bisulfite sequencing, and whole exome sequencing were applied to identify transcriptomic, epigenomic, and genomic underpinnings among 11 soft and 11 stiff PA samples surgically resected from patients at Peking Union Medical College Hospital (PUMCH). GH3 cell line and xenograft PA model was used to demonstrate therapeutic effect of sunitinib, and atomic force microscopy (AFM) was used to detect the stiffness of tumors.
Results: Tumor microenvironment analyses and immunofluorescence staining indicated endothelial cells (ECs) and cancer-associated fibroblasts (CAFs) were more abundant in stiff PAs. Weighted gene coexpression network analysis identified the most critical stiffness-related gene (SRG) module, which was highly correlated with stiff phenotype, ECs and CAFs. Functional annotations suggested SRGs might regulate PA stiffness by regulating the development, differentiation, and apoptosis of ECs and CAFs and related molecular pathways. Aberrant DNA methylation and m6A RNA modifications were investigated to play crucial roles in regulating PA stiffness. Somatic mutation analysis revealed increased intratumoral heterogeneity and decreased response to immunotherapy in stiff tumors. Connectivity Map analysis of SRGs and pRRophetic algorithm based on drug sensitivity data of cancer cell lines finally determine sunitinib as a promising agent targeting stiff tumors. Sunitinib inhibited PA growth in vitro and in vivo, and also reduced tumor stiffness in xenograft PA models detected by AFM.
Conclusion: This is the first study investigating the underlying mechanisms contributing to the stiffening of PAs, and providing novel insights into medication therapy for stiff PAs.
Introduction
Pituitary adenomas (PAs) are the second most common type of intracranial tumors, accounting for approximately 15% of primary central nervous system tumors (Ostrom et al., 2014). The mass effect of PA and secondary hypopituitarism and the multisystem complications caused by excessive secretion of hormones seriously reduce the quality of life and increase the mortality of PA patients (Molitch, 2017). During the past few decades, transsphenoidal surgery has been the first-line therapy for PAs, and only a few patients gain limited benefits from radiotherapy and medical treatments (Tabaee et al., 2009; Molitch, 2017; Almutairi et al., 2018). However, gross total resection can only be achieved in 66–78% of PA patients (Tabaee et al., 2009; Almutairi et al., 2018). Despite concurrent treatments, many patients still suffer repeated recurrence of tumors, with a 10-years recurrence rate as high as 7–12% (Reddy et al., 2011; Salomon et al., 2018). Hence, pursuing total resection of tumors during surgery and exploring new targeted drugs have become hopeful directions for reducing recurrence and curing PA patients. Due to the complicated anatomic structure in the sellar region and the limited operative field of view, it is very difficult for neurosurgeons to completely remove pituitary tumors with stiff texture, large size, and cavernous sinus invasion. In particular, tumor stiffness has become the most critical factor that affects the surgical resection rate despite the rapid progression of endoscopic surgery systems (Zhao et al., 2010; Sughrue et al., 2011). Soft tumors, even with larger size and cavernous sinus invasion, can be easily curetted through suctioning and usually have a better surgical outcome. Hence, achieving a better understanding of the underlying mechanisms of tumor stiffness and exploring novel medications to transform stiff tumors into soft ones are important for improving remission in PA patients.
Multiple studies have found that changes in tissue mechanical properties can both precede and drive disease treatment (Kai et al., 2016; Northcott et al., 2018), with tumor stiffness correlating with prognosis in several tumor types, including colorectal cancer, breast cancer and pancreatic ductal adenocarcinoma (Colpaert et al., 2001; Paszek et al., 2005; Kai et al., 2016; Laklai et al., 2016; Northcott et al., 2018; Anlaş and Nelson, 2020; Shen et al., 2020). At the tissue level, stiffness is governed by the cell cytoskeleton (Grady et al., 2016) and the extracellular matrix (ECM) (Humphrey et al., 2014). Fibrillar collagens are the most abundant ECM scaffolding proteins and contribute significantly to tissue stiffness (Mouw et al., 2014). Aberrant ECM remodeling with collagen I (COL-I) enrichment have been identified as major causes of tissue stiffening during cancer progression (Levental et al., 2009; Pickup et al., 2014). As the major source of ECM, CAFs further modify the tumor mechanical environment by expressing lysyl oxidase (LOX), an amine oxidase that initiates the process of covalent intramolecular and intermolecular crosslinking of collagen (Kagan and Li, 2003; Levental et al., 2009). In experimental models, inhibiting matrix stiffening via LOX inhibition ameliorates tumor growth and improves therapy (Levental et al., 2009). Thus, CAFs are regarded as a promising therapeutic target for limiting cancer progression (Pickup et al., 2014). Despite these findings, the molecular mechanisms that regulate the mechanical properties and contribute to stiffening in PAs still await elucidation.
Additionally, PA stiffness can be predicted through magnetic resonance imaging (MRI). Several studies testing the ability of MRI strategies utilizing machine learning algorithms to predict the tissue consistency of PAs have been performed (Hughes et al., 2016; Fan et al., 2019; Yao et al., 2020). Such strategies can help to determine individualized therapeutic schemes for PA patients. For those patients predicted to have stiff tumors, preferential use of medications that reduce tumor stiffness preoperatively might significantly improve the surgical resection rate.
In this study, by performing comprehensive multi-omics analyses of transcriptomic, genomic, DNA methylation, and m6A RNA methylation data from soft and stiff tumors, we aimed to explore the specific alterations associated with the unique biology of stiff tumors and the underlying mechanisms contributing to the stiffening of PAs. Furthermore, the responses to targeted therapy and immunotherapy were also explored to demonstrate the potential therapeutic value of treating stiff PA tumors. Sunitinib, the most critical predicted drug, was applied to treat PAs in vitro and in vivo in order to explore whether sunitinib could inhibit tumor growth and reduce stiffness.
Materials and Methods
Human Samples
We prospectively enrolled 30 patients with nonfunctioning PAs who underwent transsphenoidal adenectomy at Peking Union Medical College Hospital (PUMCH) between April 2018 and June 2018. All patients received pituitary endocrinological and neuroradiological evaluations before and after surgery and were diagnosed with nonfunctioning PAs confirmed by postoperative histopathology. The inclusion criteria were as follows: 1) intact and available perioperative endocrine examinations and sellar MRI; and 2) no history of pituitary surgery, radiotherapy or medical treatment (e.g., bromocriptine or cabergoline) preoperatively. During the surgical procedure, the consistency of the pituitary tumor was classified as stiff or soft by at least two of four experienced neurosurgeons together (Bahuleyan et al., 2006). In contrast to that of other solid tumors, the stiffness of PAs can be easily distinguished in most instances. Soft pituitary tumors are pasty, loose, and amenable to being suctioned out piecemeal with an aspirator, whereas stiff pituitary tumors grow to be almost firm spherical ellipsoids and cannot be removed easily with suction despite vigorous movement of the suction tips (Bahuleyan et al., 2006). Postoperatively, fresh tumor specimens were immediately frozen in liquid nitrogen and then stored at −80°C. After rigorous screening according to the inclusion criteria, 22 samples, including 11 soft and 11 stiff pituitary tumors, were used for integrated transcriptomic, genomic, and epigenomic profiling analyses. This study was approved by the Institutional Ethics Committee and Institutional Review Board of PUMCH (No. S-K431), and informed consent was obtained from all participants for the publication of any potentially identifiable images or data included in the study.
DNA and RNA Extraction
Pituitary tumor samples (10–50 mg) were powdered under liquid nitrogen. DNA was extracted and purified following proteinase K digestion using a TIANamp Genomic DNA Kit (TIANGEN Biotech, Catalog No. DP304-03). DNA concentrations were determined using a Qubit DNA HS Assay Kit (Thermo Fisher, Catalog No. Q32854). The DNA quality was assessed with Agilent 2,200 TapeStation Genomic DNA Analysis (Agilent Technologies).
Total RNA was extracted using TRIzol reagent (Thermo Fisher, Catalog No. 15596018). The RNA quality was checked by Agilent 2,200 TapeStation RNA Analysis (Agilent Technologies) to determine the RNA integrity number (RIN). The RNA was considered acceptable for cDNA library construction when the RIN was >7.0.
RNA Sequencing
cDNA libraries were constructed for each RNA sample using the TruSeq Stranded mRNA Library Prep Kit (Illumina, Inc.) according to the manufacturer’s protocols. The quality of the cDNA libraries was assessed with the Agilent 2,200 system, and the libraries were sequenced with the HiSeq X Ten system (Illumina, Inc.) with a 150 bp paired-end run. Raw reads were filtered with FAST-QC. Before read mapping, clean reads were obtained from the raw reads by removing the adaptor sequences and low-quality reads. The clean reads were then mapped to the human genome hg19 sequence (GRCh37) using HISAT2 (Kim et al., 2015). HTseq was used to generate gene counts, and the RPKM method was used to determine gene expression (Anders et al., 2015).
Whole-Genome Bisulfite Sequencing
Genomic DNA was bisulfite converted with the EZ DNA Methylation-Gold Kit (Zymo Research) and then processed with the TruSeq DNA Methylation Kit (Illumina, Inc.) according to the manufacturer’s instructions for WGBS library construction. The tagged WGBS libraries were used for 150 bp paired-end sequencing in a single lane of the HiSeq X Ten system (Illumina, Inc.) with 10-15% phi-X for base balance. Before read mapping, clean reads were obtained from the raw reads by removing the adapters and sequences at the 5′ and 3′ ends with methylation bias by using Trim Galore. The bisulfite mapping of methylation sites, including alignment to human genome hg19 (GRCh37), removal of duplicates, and extraction of clean reads to a CpG count matrix, was performed by utilizing Bismark and then indexed by Bowtie2 (Krueger and Andrews, 2011; Langmead and Salzberg, 2012).
Exome Sequencing and Variant Calling
To generate standard exome capture libraries, we used the Agilent SureSelectXT Reagent Kit (Agilent, Catalog No. G9611B) protocol for the Illumina HiSeq Paired-end Sequencing Library, SureSelectXT Human All Exon v6 system (Agilent, Catalog No. 5190-8864), followed by sequencing of libraries using paired-end mode (2 × 75 bp) on the HiSeq X Ten (Illumina, Inc.). Raw reads were filtered and evaluated with FAST-QC. Before read mapping, clean reads were obtained from the raw reads by removing the adaptor sequences, reads with >5% ambiguous bases (noted as Ns) and low-quality reads containing more than 20% bases with qualities of <20. The clean reads were then aligned to the human genome hg19 sequence (GRCh37) using BWA-MEM with the default settings (Houtgast et al., 2018). Variant calling was performed by using the GATK3 (version 3.6) standard pipeline with default parameters based on the mapping bam file, which was sorted, indexed and deduplicated with SAMtools and recalibrated with GATK tools (McKenna et al., 2010; Li, 2011). Mutation sites were annotated by the VEP analysis pipeline and filtered by the following criteria: allele frequency <0.01 in common frequency databases, such as gnomAD; IMAPCT over moderate; mutation frequency >0.1; sequencing depth >10; intolerant or damaging impact on protein structure and function predicted by the SIFT and PolyPhen-2 databases (McLaren et al., 2016; Amadori et al., 2020).
Analysis of the Tumor Microenvironment and Immunogenomic Patterns of PAs
To clarify the cellular components of the PA TME, the xCell algorithm was employed to accurately quantify the abundances of 64 cell types within the admixtures of tumor samples by using the xCell package in R (Aran et al., 2017). The enrichment levels of 7 subgroups of epithelial cells, nine subgroups of hematopoietic progenitors, 21 subgroups of lymphoid cells, 13 subgroups of myeloid cells, and 14 subgroups of stromal cells were estimated via single-sample gene set enrichment analysis (ssGSEA) based on the gene expression profiles of PAs (Hänzelmann et al., 2013). The Estimation of Stromal and Immune Cells in Malignant Tumors using Expression Data (ESTIMATE) algorithm was utilized to assess factors of the overall TME, including the abundances of intratumoral stromal and immune cells and tumor purity, based on the transcriptomic profiles of PA samples (Yoshihara et al., 2013). The immune and stromal scores were used to represent the abundances of immune and stromal cells, the ESTIMATE score were used to represent general nontumoral components, and tumor purity was used to reflect general tumoral proportions. CIBERSORT, a deconvolution algorithm based on linear support vector regression, was also utilized to quantify the abundances of 22 subtypes of tumor-infiltrating immune cells (TIICs) in PAs (Newman et al., 2015). Furthermore, 31 immune signatures (gene sets) introduced by He et al. (2018) were utilized to reflect the overall immune activity of solid tumors in terms of the types, related immune functions and molecular pathways of TIICs, and these factors were quantified by the ssGSEA algorithm (Supplementary Table S1). Then, based on the enrichment scores of those immunogenomic signatures, unsupervised hierarchical clustering was performed to categorize PA patients into different clusters according to immune subtype. The high immunity group included “immune hot” tumors, which had the highest enrichment scores; the low immunity group included “immune cold” tumors, which had the lowest enrichment scores; and the medium immunity group included “immune altered” tumors, which had the potential to transform into hot or cold tumors (Wang et al., 2021). Additionally, to elucidate the functional heterogeneity of cancer cells at single-cell resolution, the gene signature profiles of 14 functional states of cancer cells obtained from CancerSEA were evaluated by ssGSEA using transcriptomics data from PA samples (Supplementary Table S1) (Yuan et al., 2019).
Immunofluorescence Analysis
Immunofluorescence (IF) staining was performed on formalin-fixed, paraffin-embedded (FFPE) sections of tumor tissues. Briefly, FFPE tissues were cut into 4 μm sections, followed by deparaffinization and rehydration using xylene and ethanol. Next, the slides were incubated in EDTA antigen retrieval buffer at subboiling temperature and then in blocking solution (BSA; G5001, Servicebio) for 30 min at room temperature. The slides were incubated overnight with primary antibodies, including anti-CD31 (ab28364, Abcam), anti-von Willebrand factor (VWF) (bs-10048R, Bioss), anti-smooth muscle actin-α (αSMA) (A2547, Merck), and anti-S100A4 (bs-3759R, Bioss) antibodies, followed by incubation with fluorochrome-conjugated secondary antibodies for 50 min at room temperature. DAPI (G1012, Servicebio) was used to stain cell nuclei. Pictures were taken with an Ortho fluorescence microscope (ECLIPSE C1, NIKON). All staining was quantified using NIH ImageJ 1.51s analysis software with the same threshold for each stain.
Gene Set Variation Analysis
GSVA was utilized to evaluate the 50 most significantly enriched hallmark pathways obtained from the Molecular Signatures database (MSigDB) (Subramanian et al., 2005) in soft and stiff tumors by using the GSVA package in R (Hänzelmann et al., 2013). Differential analysis of the enrichment scores of molecular pathways between two groups was performed with the limma package in R (Smyth et al., 2005). The hallmark pathways with |t value| > 4, indicating a false discovery rate (FDR) < 0.05, were considered the most differentially enriched molecular pathways between the two groups (Lambrechts et al., 2018).
Weighted Gene Coexpression Network Analysis
The differentially expressed genes (DEGs) between soft and stiff tumors were identified by using the edgeR package in R based on the raw count data (Robinson et al., 2010). The criteria for selecting DEGs were Benjamini–Hochberg corrected FDR <0.05 and |fold change (FC)| > 2. Next, the coexpression network for the DEGs was constructed by the WGCNA package in R based on the RPKM data (Langfelder and Horvath, 2008). First, sample clustering was performed to detect outliers. The pickSoftThreshold method was used to select an appropriate soft threshold (power) to achieve a scale-free topology fit index >0.85 and maintain optimal mean connectivity. Afterwards, the adjacency matrix was transformed into a topological overlap matrix (TOM) to define gene coexpression similarity, and gene hierarchical clustering for TOM-based dissimilarity was performed to obtain the hierarchical clustering dendrogram. The Dynamic Tree Cut package was used to identify the modules with a minimum gene size of 50, and then the similarity cut-off was set to 0.75 to merge the modules after calculating the dissimilarity of module eigengenes (MEs), representing the overall expression profiles of each module. The adjacency of the MEs of all modules was determined by Pearson correlation analysis.
Identification of Clinically Significant Modules
Pearson correlation analyses between MEs and clinical traits were performed to determine the key gene signatures associated with the stiffness of tumors. Furthermore, gene significance (GS) was calculated as the absolute value of the correlation between each gene within MEs and each trait, and module membership (MM) was calculated as the correlation of the gene expression profile and each ME. A high correlation between GS and MM suggested a highly significant association between the modules and clinical traits (Zhang and Horvath, 2005).
Construction of the Protein-Protein Interaction Network and Functional Enrichment Analysis of Key MEs
The genes within the key module associated with stiffness were mapped with the STRING database (http://string-db.org) to evaluate their functional associations, and a combined score >0.4 was defined as significant (Mering et al., 2003). The PPI network, representing the topology of the interactions between SRGs, was constructed and visualized by Gephi software. To further explore the biological properties and molecular mechanisms of the SRGs, gene ontology (GO) enrichment analyses were performed with the ClueGO plug-in of Cytoscape (Bindea et al., 2009). An FDR <0.05 was considered significant.
Identification of Differentially Methylated Regions and DMR-Associated Genes
DMRs between soft and stiff tumors were identified by using the DSS package in R (Park and Wu, 2016). First, the DMLtest function was used to perform statistical tests for differentially methylated loci (DML), including estimating mean methylation levels for all CpG sites and dispersions at each CpG site and performing the Wald test. Based on the DML results, the callDMR function was applied to identify DMRs consisting of at least four statistically significant CpG sites. A |delta| > 0.1 and p < 0.01 were considered the thresholds for determining DMLs and DMRs. Genes overlapping with DMRs in the whole genome were considered DMR-associated genes.
Evaluation of m6A Modification
To explore the potential roles of RNA modifications in regulating the stiffness of PAs, three types of m6A regulators, including 2 demethylases (erasers), 7 methyltransferases (writers), and 11 RNA binding proteins (readers), were compared between soft and stiff tumors (Li et al., 2019). The post-methylation regulation of stiffness-related mRNAs was evaluated by Pearson analysis of the correlation between the expression of readers and SRGs. The overall effects of m6A in regulating the stiffness of PAs were comprehensively assessed as reported in the literature (He et al., 2019).
Somatic Mutation Analysis
The different genomic variations between soft and stiff tumors were explored by somatic mutation analysis. The mutation types and frequencies of the top mutated genes were visualized by a waterfall plot using the maftools R package (Mayakonda et al., 2018). The differentially mutated genes (DMGs) between soft and stiff tumors were detected by the mafComapre function, which used Fisher’s test to compare all genes between the two groups. An FDR <0.05 was considered statistically significant. As reported by previous studies, mutation of transcription factors (TFs) plays a critical role in the development and progression of multiple cancers. The differentially mutated TFs were visualized by a lollipop plot. The target genes of TFs were obtained from the Gene Transcription Regulation database (GTRD, http://gtrd.biouml.org/) based on the ChIP-seq data. Tumor mutation burden (TMB) was defined as the total number of nonsynonymous mutations in the coding region per megabase (Budczies et al., 2018). Tumor heterogeneity was inferred by clustering variant allele frequencies (VAFs) using the inferHeterogeneity function, which clustered variants to infer clonality. According to the VAF clustering, clones were defined as either mutations that occurred in most cancer cells or that only existed in a small number of cells (subclones). The mutant-allele tumor heterogeneity (MATH) score is a novel quantitative measure of intratumoral genetic heterogeneity that calculates the width of the VAF distribution (Mroz et al., 2015). Tumors with high heterogeneity and a high subclonal fraction tend to experience immune evasion and resist immunotherapy (Guan and Shastri, 2018).
Prediction of Targeted Drugs for PAs and Therapeutic Response
The Connectivity Map (CMap) database (https://clue.io/) was employed to explore potential compounds targeting the molecular pathways and genes associated with the stiffness of PAs (Subramanian et al., 2017). The SRGs of stiff tumors were considered potential targets of compounds used to query the CMap database. The enrichment scores of compounds were calculated, and compounds with enrichment scores < -95 and p < 0.05 were considered potential therapeutic drugs for stiff PAs. The most enriched mode of action (MoA) and corresponding drugs were selected for further analysis. Drug sensitivity data of cancer cell lines obtained from the Cancer Cell Line Encyclopedia (CCLE) project were obtained from the Genomics of Drug Sensitivity in Cancer database (GDSC, https://www.cancerrxgene.org/) (Geeleher et al., 2014a). These data provide the half maximal inhibitory concentration (IC50) as the measure of drug sensitivity, and lower IC50 values indicate increased sensitivity to compounds. After integrating the gene expression profiles of cell lines (as the training set) and PA samples (as the test set), the IC50 values of drugs in PA patients were estimated by ridge regression analysis via the pRRophetic R package, and the prediction accuracy was assessed by 10-fold cross-validation (Geeleher et al., 2014b).
Prediction of the Response of PAs to Immunotherapy
The Tumor Immune Dysfunction and Exclusion (TIDE, http://tide.dfci.harvard.edu/) algorithm quantifies T cell dysfunction signatures by testing how the expression of each gene in solid tumors affects the infiltration of cytotoxic T lymphocytes (CTLs) to influence the immunotherapy response (Jiang et al., 2018). TIDE scores, which are suggestive of the clinical response to immunotherapy, were calculated based on the gene expression profiles of PA samples. A TIDE score <0 was considered to indicate sensitivity to immunotherapy, and >0 was considered to indicate resistance to immunotherapy. In addition, an unsupervised subclass mapping (https://cloud.genepattern.org/gp/) method was utilized to predict the clinical response to anti-PD1 and anti-CTLA4 therapy in soft and stiff tumors (Hoshida et al., 2007). An FDR <0.05 was considered to indicate a significant response to immune checkpoint inhibitors (ICIs).
Cell Culture
The GH3 cell line were purchased from National Infrastructure of Cell Line Resource (Beijing, China) and cultured in DMEM/F12 medium supplemented with 15% horse serum, 2.5% fetal bovine serum, and an antibiotic-antimycotic solution. Cell cultures were maintained in an incubator at 37°C and 5% CO2.
Cell Viability Assay
For Cell Counting Kit-8 (CCK-8) cell viability assays, cells were seeded in 96-well plates at a density of 5,000 cells per well and then incubated with dimethyl sulfoxide (DMSO; control) or a series of 2-fold-diluted concentrations of sunitinib (cat. no. S7781, Selleck) for 2 days. Six parallel wells were used for each concentration of the drugs. After incubation, CCK-8 assays (cat. no. CCK-8, Dojindo) were performed to measure cell viability according to the manufacturer’s instructions.
Xenograft Experiment
A total of 12 Wistar Furth rats were used for generation of the GH3 xenograft model. Four-week-old female rats were subcutaneously injected with GH3 cells (5 × 105) into the left lumbar area. Tumor size was measured every 3 days, and tumor volume was calculated as width2 × length × 0.5. Twelve rats were randomly divided into 2 groups, including drug treatment group and control group. sunitinib was intragastrically administered at a dose of 40 mg/kg once daily when the tumor volumes approached approximately 100 mm3. All rats were sacrificed upon completion of the 12-days experiment, and the tumors were excised, weighed, and maintained in PBS at 4°C for stiffness analysis. All animal experiments and euthanasia were approved by Animal Care and Use Committee of PUMCH.
Atomic Force Microscopy-Based Young’s Modulus Measurement of Tumor Samples
Resected xenograft PA samples were cut to proper size by a scalpel and immobilized on 6 cm-cell culture dishes (Thermofisher, cat. no. 15462) with two slices of 1mm× 5 cm-parafilm (Bemis, PM-996). Then the tissues were immersed in 1× PBS, and applied to the AFM force-measuring setting. A homemade AFM from Institute of Biophysics, Chinese Academy of Sciences was used in this study. Using constant force mode, more than 2000 force-to-distance traces (force curves) were recorded and more than 300 traces were selected for each group. The slope K (pN/nm) of every trace was calculated and then transformed to Young’s modulus (kilopascal, kPa) by using JPK Data Processing software. The final Young’s modulus of each tumor was calculated taking into account all traces recorded for a single tissue sample. The tumor stiffness between sunitinib treatment and control group were compared by mean Young’s modulus of each group, which was calculated by the Gaussian fitting via the amplitude version of Gaussian peak function (GaussAmp) in Origin software.
Statistical Analyses
Independent Student’s t test was utilized for continuous variables and the χ2 test was utilized for categorical variables when making comparisons between two groups. The Mann-Whitney U test was used to compare categorical variables and nonnormally distributed variables between two groups. The Kruskal–Wallis test was used to compare multiple groups. Correlation analysis was performed by the Pearson correlation test, and a p value <0.05 and |correlation coefficient| > 0.3 were considered to indicate significant correlation. The statistical analyses in this study were performed with R 3.6.1 software. A two-tailed p value <0.05 was considered to indicate statistical significance.
Results
Clinicopathological Features of Soft and Stiff PAs
The demographics and clinicopathological features of 22 PA patients are summarized in Table 1. In general, almost all the clinical variables did not differ significantly between soft and stiff tumors. Gross total resections were achieved in all soft PAs (100%), whereas 18.2% of stiff tumors were only subtotally resected via transsphenoidal surgery due to firm texture. After a long-term follow-up (2.35 ± 0.06 years), recurrence occurred in three patients (all from the stiff tumor group), with an average recurrence time of 1.30 years. The recurrence in two of these patients was due to the active growth of residual tumors.
TABLE 1. Demographics and clinicopathological features of 22 patients diagnosed with pituitary adenomas.
ECs and CAFs Were More Abundant in Stiff PAs
The overall workflow of this integrated transcriptomic, genomic, and epigenomic profiling analysis is displayed in Figure 1. First, the xCell algorithm was employed to quantify the cellular components of PA samples, including epithelial, hematopoietic progenitor, lymphoid, myeloid, and stromal cell clusters (Figure 2A). Then, the ESTIMATE algorithm was employed to assess the overall TME of PAs. Compared with the respective scores in soft tumors, the stromal and ESTIMATE scores were significantly higher and the tumor purity score was significantly lower in stiff PAs (all p < 0.05), and stiff PAs demonstrated high abundances of general stromal cells and low tumor purity (Figure 2B). However, the immune score did not differ significantly between soft and stiff tumors, and the CIBERSORT algorithm also demonstrated no significant difference in the infiltration abundances of the majority of TIICs between the two groups (Figure 2C). Next, the enrichment levels of 31 immune signatures, representing the overall immune activity of PAs, were quantified by ssGSEA, and the 22 PA patients were classified into three immune subtypes by unsupervised hierarchical clustering (Figure 2D). The distributions of soft and stiff tumors in the high-, medium-, and low-immunity groups did not differ significantly (p = 0.659).
FIGURE 1. The overall workflow of the integrated transcriptomic, genomic, and epigenomic profiling analyses in this study.
FIGURE 2. Estimation of tumor immune microenvironment patterns associated with the stiffness of PAs. (A) Heatmap illustrating the cellular components of PA samples, including epithelial, hematopoietic progenitor, lymphoid, myeloid, and stromal cell clusters, quantified by the xCell algorithm. (B) Comparisons of stromal score, immune score, ESTIMATE score and tumor purity between soft and stiff tumors. (C) Comparisons of the abundances of 22 immune cells between soft and stiff tumors. (D) Heatmap illustrating the immune subtypes of PA patients, which were categorized based on the overall immune activity of tumors. (E) Comparisons of ECs, lymphatic ECs, microvascular ECs, and CAFs between soft and stiff tumors. (F) Correlation analysis of ECs, lymphatic ECs, microvascular ECs, and CAFs. (G) Correlation analysis of ECs, CAFs and 14 functional processes of cancer cells. (H) Left panels: Immunofluorescence staining of CAF markers (αSMA and S100A4) and EC markers (CD31 and VWF) (green) in soft and stiff tumors. Cell nuclei were counterstained with DAPI (blue). Scale bar: 50 μm. Right panels: Quantification of the fluorescence intensity of four cell markers between soft and stiff PAs. a. u., arbitrary unit. (I) The differential hallmark pathways associated with stiff and soft tumors. * means p < 0.05, ** means p < 0.01, and *** means p < 0.001.
Among all the cellular components of PA samples, general ECs, lymphatic ECs, microvascular ECs, and CAFs showed significantly higher abundances in stiff tumors than in soft tumors (Figure 2E), and they were highly positively correlated (Figure 2F). Furthermore, to elucidate the functional associations between ECs, CAFs and 14 functional processes of cancer cells, Pearson correlation analysis was used, and the results demonstrated that general, lymphatic, and microvascular ECs were negatively correlated with the cell cycle, DNA damage/repair, hypoxia, and stemness, whereas CAFs were positively correlated with angiogenesis, apoptosis, epithelial-mesenchymal transition (EMT), hypoxia, inflammation, invasion, metastasis, and proliferation (Figure 2G). To verify the difference in stromal cells between soft and stiff tumors, we performed immunofluorescence staining of CAF markers (αSMA and S100A4) and EC markers (CD31 and VWF) (Figure 2H). Significantly higher expression levels of αSMA, S100A4, CD31 and VWF were observed in stiff tumors than in soft tumors, suggesting that there are higher abundances of CAFs and ECs in stiff PAs.
GSVA was further performed to explore the hallmark pathways and underlying molecular mechanisms associated with the stiffness of PA samples. A total of 23 differentially enriched molecular pathways were identified, including 18 pathways positively correlated with stiff tumors and five pathways positively correlated with soft tumors (Figure 2I). EMT was the most significant hallmark pathway related to stiff PAs, suggesting a potential molecular mechanism underlying the stiff phenotype.
All these findings suggest that stromal cells, especially ECs and CAFs, and EMT might play critical roles in contributing and promoting the stiffening of PAs, whereas immune cells might not be associated with the stiffness of PAs.
Identification of the Top Stiffness-Related Gene Module of PAs by WGCNA
A total of 1,288 DEGs between soft and stiff PAs were identified, including 759 upregulated and 529 downregulated genes in stiff tumors, and the results are displayed in the volcano plot (Figure 3A). To determine the phenotypic relevance of the DEGs, WGCNA was performed to identify the most critical gene module related to the stiffness of tumors. By using the pickSoftThreshold method, 30 was selected as the soft-thresholding power needed to achieve a scale-free topology fit index >0.85 and maintain optimal mean connectivity (Supplementary Figure S1). A hierarchical clustering dendrogram was obtained, and 12 gene modules were ultimately generated by employing the Dynamic Tree Cut package (Figure 3B). The gray module, including all the genes that could not be enrolled into any other modules, was excluded from the subsequent analysis. The network heatmap demonstrated that MEs were highly correlated within each module, suggesting that highly coexpressed eigengenes in the same module (Supplementary Table S2) may possess similar biological significance and function together (Figure 3C). In addition, to explore the coexpression similarity of all modules, the modules were mainly divided into two clusters in the hierarchical clustering dendrogram and the eigengene adjacency heatmap according to their correlations with each other (Figure 3D). Furthermore, analyses of the correlations between MEs and clinical traits were performed, and the turquoise module, consisting of 131 genes, was positively correlated with stiffness (R = 0.84, p = 1 × 10−05), general ECs (R = 0.90, p = 2 × 10−07), lymphatic ECs (R = 0.94, p = 4 × 10−10), microvascular ECs (R = 0.94, p = 2 × 10−09), and CAFs (R = 0.70, p = 0.003) (Figure 3E). There were high correlations of MM for genes in the turquoise module and GS with stiffness, ECs and CAFs, which also indicated the critical roles of the turquoise module in promoting the stiffness of PAs (Supplementary Figure S2). Hence, the 131 MEs in the turquoise module were defined as SRGs.
FIGURE 3. Identification of the top stiffness-related gene module of PAs by WGCNA. (A) Volcano plot displaying the DEGs between stiff and soft tumors. Red dots represent the upregulated genes, and green dots represent the downregulated genes in the stiff PAs. (B) Cluster dendrogram of the DEGs. Upper panel: Each branch represents one single gene. Lower panel: Each color represents one coexpression module. (C) Heatmap illustrating the interactions of coexpressed genes. The brightness of yellow in the heatmap represents the degree of connectivity of different modules, with a lighter color indicating greater overlap. The colors of the horizontal and vertical axes represent different modules, and the branches represent different genes. (D) Upper panel: Hierarchical clustering of module genes. Lower panel: Heatmap of the adjacencies in the module gene network. (E) Heatmap displaying the correlations between module eigengenes and the clinical traits of patients with PA. The turquoise module was the most critical module and was positively correlated with stiffness, general ECs, lymphatic ECs, microvascular ECs, and CAFs. (F) PPI network including 131 SRGs from the turquoise module. A change in dot color from lighter to darker indicates an increasing degree of gene connectivity. Red circles represent TFs, and purple circles represent TF cofactors. (G) Functional enrichment analysis of the SRGs. Different dot colors indicate different GO clusters.
To better illustrate the functional associations and topology of the interactions between SRGs, we constructed a PPI network (Figure 3F). Three TFs (SOX18, HES1, and ZNF358) and 4 TF cofactors (CCDC85B, ENG, PYCARD, and TRIP6) were identified in the network. Functional annotations of the 131 SRGs indicated that they were mainly enriched in 14 GO clusters, representing 30 GO terms (Figure 3G). Vascular smooth muscle cell development, including 11 terms (36.7%), was the most important functional cluster, followed by lymph vessel development (13.3%), endothelial cell apoptotic process (6.7%), regulation of VEGFR signaling pathway (6.7%), etc. These findings suggest that the SRGs might regulate the stiffness of pituitary tumors by regulating the development, differentiation, and apoptosis of ECs and CAFs and related molecular pathways.
Associations Between DMRs and the Stiffness of PAs
To investigate epigenomic associations with the stiffness of PAs, genome-wide methylation profiles were utilized to determine DMRs between soft and stiff tumors. A total of 188 significant CpG sites and 38 DMRs were identified, including 25 regions with higher methylation in soft tumors and 13 regions with higher methylation in stiff tumors (Figure 4A and Supplementary Table S3). The hypermethylated DMRs in stiff tumors were more enriched within exons (p = 0.037) than were those in soft tumors, whereas the percentage of hypermethylated DMRs in the promoter regions did not differ significantly between the two groups (Figure 4B). DMR-associated genes were identified for all the DMRs, and only two of them (C5orf66-AS1 and CAVIN3) belonged to the SRGs. The methylation levels of the promoter regions of C5orf66-AS1 and CAVIN3 were significantly lower in stiff tumors than in soft tumors, and in contrast, their gene expression was higher in the stiff tumors (Figures 4C,D). The strong negative correlations between promoter methylation and mRNA expression levels of C5orf66-AS1 and CAVIN3 suggest that aberrant DNA methylation of the SRGs might be a crucial process contributing to the stiffening of PAs.
FIGURE 4. Associations of DNA methylation and m6A RNA methylation with the stiffness of PAs. (A) Heatmap illustrating the 38 DMRs between soft and stiff tumors. (B) Comparisons of hypermethylated DMRs between soft and stiff tumors. (C,D) Left panel: Comparisons of the methylation levels of CpG sites within the DMRs of promoter regions between soft and stiff PAs. Middle panel: Comparisons of the methylation levels of promoter regions between soft and stiff PAs. Right panel: Comparisons of the gene expression levels between soft and stiff PAs. (E) Heatmap comparing 20 m6A regulators between soft and stiff tumors. (F) Analysis of the correlation between reader and SRG expression. Blue dots represent SRGs, and purple rhombi represent readers. Red lines represent negative correlations, and green lines represent positive correlations. (G) The overall role of m6A in regulating the stiffness of PAs is believed to be two-sided, including promotion and inhibition of stiffness.
Dynamic Bilateral Regulation of Stiffness by m6A
To explore the potential roles of RNA modifications in regulating the stiffness of PAs, we assessed the associations between m6A regulators and SRGs. Twenty m6A regulators were compared between soft and stiff tumors, and we found that FTO and RBM15 were significantly downregulated in stiff tumors, suggesting low m6A levels in these PAs (Figure 4E). Analysis of the correlation between reader and SRG expression demonstrated that the post-methylation regulatory effects on target genes were mainly negative, but some readers showed positive effects (IGF2BP1 and IGF2BP3) (Figure 4F). m6A has been previously reported to have a dual role in regulating the stiffness of PAs, mostly promoting stiffening and rarely inhibiting stiffening of pituitary tumors (Figure 4G).
High Intratumoral Heterogeneity in Stiff PAs
Somatic mutation analysis was performed to investigate distinct genomic alterations between soft and stiff tumors that have been reported to regulate the TME and immunity. A total of 12 DMGs were identified between the two groups (Figure 5A), and the mutation frequency of AR, the only TF, was significantly higher in the soft tumor group than in the stiff tumor group (63.6 vs. 18.2%, p = 0.035). A lollipop plot was generated and illustrated that the mutation sites and types of AR were distinct in the soft and stiff tumors (Figure 5B). TF mutation has been reported to play a crucial role in pathogenesis; hence, we performed correlation analysis between AR and its target genes. SELENBP1, one of the SRGs and one of the target genes of AR, was positively correlated with stiff tumors (R = 0.65, p = 0.029) (Figure 5C). In addition, as shown in Figure 5D, patients with stiff tumors had lower TMB (p = 0.11), higher MATH score (p = 0.023), lower clonal fraction (p = 0.019), and higher subclonal fraction (p = 0.019) than patients with soft tumors, indicating that there is higher intratumoral heterogeneity in stiff PAs (Figure 5E). All these findings suggest underlying differences in immune evasion and immunotherapy response between soft and stiff tumors. The general TME patterns of soft and stiff tumors are displayed in Figure 6.
FIGURE 5. Comparisons of somatic variations between soft and stiff tumors. (A) Waterfall plots showing the mutation types and frequencies of 12 DMGs between soft and stiff tumors. (B) Lollipop plot illustrating the mutation sites and types of AR in soft and stiff PAs. aa: amino acid. Androgen_recep: androgen receptor (6-449 aa); zf-C4: zinc-finger, C4 type (two domains) (558-626 aa); Hormone_recep: ligand-binding domain of nuclear hormone receptor (692-877 aa). (C) Analysis of the correlation between AR and its target gene (SELENBP1) in soft and stiff PAs. An AR motif is displayed in the upper panel. (D) Comparisons of TMB, MATH score, clonal fraction, and subclonal fraction between soft and stiff tumors. (E) Schematic diagram displaying the intratumoral heterogeneity of soft and stiff PAs.
FIGURE 6. Schematic diagram illustrating the general TME patterns of soft and stiff tumors.
Identification of Potential Drugs for Treating Stiff PAs
CMap analysis of the SRGs and related molecular pathways was performed to explore potential compounds that could be used to treat stiff PAs. MoA analysis revealed 35 molecular pathways targeted by 36 compounds in the stiff tumors (Figure 7A and Supplementary Table S4). Regarding the most enriched and critical MoAs, there were 14 compounds sharing the same MoA as VEGFR inhibitors, 7 compounds sharing the same MoA as PDGFR inhibitors, and three compounds sharing the same MoA as FGFR inhibitors in the stiff tumors. Hence, the VEGF, PDGF and FGF signaling pathways might serve as potential therapeutic targets for stiff PAs. As shown in Figure 7B, the enrichment levels of the VEGF, PDGF and FGF signaling pathways were significantly higher in the stiff tumors than in the soft tumors, indicating that there is activation of these molecular pathways in stiff PAs and that they have potential roles in promoting the stiffening of PAs. Regarding the structurally related factors and receptors of the VEGF family, only VEGFR3 (FLT4) was observed to be significantly higher in the stiff PAs (Figure 7C). In addition, the therapeutic responses of stiff PAs to axitinib, pazopanib, sorafenib, and sunitinib were evaluated by using the pRRophetic algorithm based on GDSC data. By integrating the gene expression profiles of cell lines and PA samples, we estimated the IC50 values of the four drugs in each PA patient using ridge regression analysis. The estimated IC50 value of sunitinib was significantly lower in stiff PAs than in soft PAs (p = 0.003), indicating that PA patients with stiff tumors tended to be more sensitive to sunitinib therapy (Figure 7D).
FIGURE 7. Prediction of PA targeted therapy response. (A) CMap-based exploration of candidate drugs and molecular pathways that can be used for treatment of stiff PAs based on the SRGs. The MoA analysis revealed 35 molecular pathways targeted by 36 compounds. (B) Comparisons of the ssGSEA enrichment scores of the VEGF, PDGF, FGF and RTK signaling pathways between soft and stiff tumors. (C) Comparisons of the gene expression of structurally related factors and receptors of the VEGF family between soft and stiff tumors. (D) Comparisons of the IC50 values of axitinib, pazopanib, sorafenib, and sunitinib as treatment between soft and stiff PA samples, as estimated by the pRRophetic algorithm based on GDSC data. (E) Comparisons of the gene expression of ICMs between soft and stiff tumors. (F) Comparisons of the proportions of nonresponders and responders to immunotherapy between the soft and stiff tumor groups. NR, nonresponder; R, responder. (G) Comparisons of the TIDE score between soft and stiff tumors. A TIDE score <0 indicates sensitivity to immunotherapy, and a TIDE score >0 indicates resistance to immunotherapy. (H) Subclass mapping analysis for predicting the likelihood of response to ICI treatments in soft and stiff PAs.
Stiff PAs Were More Resistant to Immunotherapy
Regarding immune checkpoint molecules (ICMs), the expression levels of PD-1 and TIGIT were significantly higher in the soft PAs than in the stiff PAs, whereas other ICMs did not differ between the two groups (Figure 7E). Then, the TIDE algorithm, which quantifies T cell dysfunction signatures, was applied to predict the likelihood of immunotherapy response in PA patients. The proportion of responders to immunotherapy in the patients with soft PAs was two times greater than that in the stiff tumors (54.5 vs. 18.2%, p = 0.183) (Figure 7F). The TIDE scores, calculated based on gene expression profiles, were significantly higher in the stiff tumors than in the soft tumors (p = 0.039), suggesting that the patients with soft PAs were more sensitive to immunotherapy than were those with stiff tumors (Figure 7G). In addition, an unsupervised subclass mapping method was utilized to predict the clinical response to ICIs, including PD1 and CTLA4 inhibitors, of soft and stiff tumors. Patients with soft tumors were more likely to respond to anti-PD-1 therapy than those with stiff tumors (FDR = 0.011), whereas neither group was more sensitive to anti-CTLA4 therapy (Figure 7H). Generally, the higher TMB, lower intratumoral genetic heterogeneity, lesser subclonality, and higher ICM expression of soft PAs might explain why they are more sensitive to immunotherapy, especially anti-PD-1 treatment, than stiff PAs.
Sunitinib Inhibited PA Growth in Vitro and in Vivo, and Reduced Tumor Stiffness in Xenograft PA Models Detected by AFM
To further investigate the sensitivity of PAs to sunitinib treatment, we tested the cell viability of GH3 cell line treated by sunitinib. After 2 days of treatment, sunitinib exhibited promising anti-proliferative effects in GH3 cells, with a half-maximal inhibitory concentration value of 41.81 µM (Figure 8A). Subsequently, to evaluate the impact of the sunitinib on GH3 cells in vivo, we generated a xenograft PA model by transplanting GH3 cells subcutaneously into the flanks of Wistar Furth rats. Once the tumor volume approached approximately 100 mm3, treatment was started with intragastric administration of 40 mg/kg sunitinib or vehicle control once daily for 12 days. All rats were sacrificed after completion of the 12-days experiment. Sunibtinib treatment significantly inhibited tumor growth with regard to tumor volume (59% inhibition, p < 0.05) and tumor weight (37% inhibition, p < 0.05) in comparison with the control regimen (Figures 8B–D). Additionally, sunibtinib treatment showed minimal effects on rats’ body weights, demonstrating its safety (Figure 8E). Then, AFM was further applied to examine the mechanical properties of the resected tumor samples from two groups. The overall working schematic of the AFM setup used for mechanical property measurement, especially Young’s modulus changes, was shown in Figures 8F,G. The distributions of Young’s modulus of all traces recorded for the tumors in two groups were displayed by the statistic histograms (Figures 8H,I). Mean Young’s modulus (xc ± SE), calculated by the Gaussian fitting, was 0.85 ± 0.34 kPa for the sunitinib treatment group and 0.90 ± 0.03 kPa for the control group. Kolmogorov-Smirnov test demonstrated that the Young’s modulus of collected traces was lower in the sunitinib treatment group than that in the control group (p < 0.0001; Figure 8J). All these findings indicated that sunitinib can inhibit tumorigenesis and reduce the stiffness of pituitary tumors, which suggested sunitinib could serve as a potential candidate drug for stiff PAs.
FIGURE 8. Sunitinib inhibited PA growth in vitro and in vivo, and reduced tumor stiffness in xenograft PA models examined by AFM. (A) Cell viability of GH3 cells treated with sunitinib or DMSO for 2 days. Three independent experiments were repeated for each result. *** means p < 0.001. Bar represents mean ± SD. (B) Tumor volume of rats treated with sunitinib or vehicle control once daily, measured for 12 days. (C) Tumor weight of rats treated with sunitinib or vehicle control. ** means p < 0. 01. (D) Macroscopic image of the resected tumors. Scale bar: 1 cm. (E) Body weight of rats measured for 12 days “NS” means no significant difference between groups. (F) Working schematic of the AFM setup used for mechanical property measurement of tumor samples. (G) Slope changes deliver Young’s modulus changes of different interaction surfaces. Black curve represents force vs. distance curve recorded for the chamber surface (reference), and the red one represents force vs. distance curve for the tumor tissue. Statistic histograms show the distributions of Young’s modulus of all traces recorded for tumor tissues in the sunitinib treatment group (H) and control group (I). Mean Young’s modulus (xc ± SE), calculated by the Gaussian fitting, was 0.85 ± 0.34 kPa for the sunitinib treatment group and 0.90 ± 0.03 kPa for the control group. E (kPa) means the Young’s modulus of the traces in each sample. (J) Comparisons of the Young’s modulus of all collected traces between two groups. **** means p value <0.0001 via using the Kolmogorov-Smirnov test.
Discussion
As the first-line and even the only therapy for PAs, total resection of tumors in the transsphenoidal surgery was extremely important for curing PAs. Despite rapid progression of endoscopic surgery systems, tumor stiffness has become the most critical factor that affects the surgical resection rate in invasive tumors. Particularly, total resection of stiff PAs invading the cavernous sinus is the most challenging for neurosurgeons (Bao et al., 2016; Kim et al., 2018). Consistent with this, in this study, stiff tumors had a significantly lower gross total resection rates than soft tumors, and the recurrence rate of stiff tumors was higher. These results suggest that the basis of this investigation is well founded and that the samples used in the study are robust. However, the molecular mechanisms that regulate the mechanical properties and contribute to stiffening in PAs were still unknown.
Tissue stiffness can be increased by collagen content and fiber organization, which lead to angiogenesis (Bayer et al., 2019; Shen et al., 2020). In our study, the xCell algorithm revealed that ECs and CAFs were more abundant in stiff tumors than in soft tumors. EMT might play critical roles in the stiffening of PAs. Furthermore, immunofluorescence staining of CAF markers (αSMA and S100A4) and EC markers (CD31 and VWF) verified the differences in stromal cells. Therefore, CAFs and ECs play important roles in tumor stiffness. In addition, the SRGs identified by the transcriptome analysis were enriched in EC and CAF regulation pathways. These results imply that these SRGs probably lead to a high abundance of ECs and CAFs, which might lead to stiffening of tumors.
DNA methylation, in which alterations of CpG dinucleotides block the transcriptional mechanism and silence gene expression, is the most frequently studied epigenetic phenomenon (Pease et al., 2013). By studying the correlation between DNA methylation data with SRG expression data, we found that the DNA methylation level was negatively correlated with the RNA expression of SRGs, including C5orf66-AS1 and CAVIN3, which indicates that SRGs are probably regulated by DNA methylation. C5orf66-AS1 is a long noncoding RNA (lncRNA) that suppresses the development and invasion of pituitary null cell adenomas (Yu et al., 2017). Two groups also reported aberrant methylation-mediated downregulation of the lncRNA C5orf66-AS1. One study was in gastric cardia adenocarcinoma (Guo et al., 2018), and the other study was in head and neck squamous cell carcinoma (Sailer et al., 2018). The relationships between DNA methylation, SRGs, and ECs and CAFs are probably responsible for the regulation of PA stiffening.
RNA methylation is an important regulatory factor in different physiological and pathological processes. including development (Yoon et al., 2017; Li et al., 2018; Ma et al., 2018; Wang et al., 2018), neurogenesis (Li et al., 2017), innate immunity (Zheng et al., 2017), and tumorigenesis (Zhang et al., 2016; Zhang et al., 2017a; Cui et al., 2017; Ma et al., 2017; Chen et al., 2018; Chen et al., 2019; Han et al., 2019; Lan et al., 2019; Niu et al., 2019; Jin et al., 2020). m6A has been shown to play a role in fate determination of cells (Batista et al., 2014; Aguilo et al., 2015; Chen et al., 2015; Zhang et al., 2017b). In this study, we found that m6A regulators (FTO and RBM15) were significantly downregulated in stiff tumors, suggesting aberrant m6A levels in stiff PAs. Moreover, the expression levels of m6A “readers” and SRGs were mainly decreased, while some readers showed increased expression. m6A has been shown to regulate RNA expression levels by influencing RNA degradation (Wang et al., 2014; Ke et al., 2017). Therefore, m6A might play a role in the posttranscriptional regulation of SRGs and further regulate EC and CAF production. Consistent with this, m6A was found to be critical for the development of cardiac fibrosis (Li et al., 2021). In another study, the m6A-mediated MALAT1/miR-145/FAK pathway was found to be involved in renal fibrosis (Liu et al., 2020). Three ECM components (COL6A1, LAMA5, and FN1) are target genes of the m6A reader IGF2BP3 and can be regulated in an m6A-dependent manner (Gu et al., 2020). Together, these mechanisms provide strong evidence that m6A regulation might play an important role in ECM component production, which is important for tissue stiffening.
Consistent with the increased levels of ECM components seen in stiff tumors, tumor heterogeneity is also increased in stiff tumors. Compared with soft tumors, stiff tumors showed higher MATH scores, indicating that there is higher intratumoral heterogeneity in stiff PAs. Increased heterogeneity can lead to tumor evolution (McGranahan and Swanton, 2017; Pelham et al., 2020), drug resistance (Lim and Ma, 2019), and immune evasion (Vinay et al., 2015). In addition, immune markers were decreased in stiff tumors, and stiff PAs were likely to show a decreased response to immunotherapy. Taken together, these data suggest that stiff tumors have higher heterogeneity, greater subclonality, and fewer immunotherapy targets than soft tumors, which make them harder to treat with targeted therapy. Therefore, transforming stiff tumors into soft tumors can make tumor treatment easier.
We used CMap to obtain compounds that can inhibit SRG expressions. MoA analysis revealed compounds targeting the most enriched and critical molecular pathways. Sunitinib was ultimately selected as the drug to which stiff PAs would be most sensitive using the pRRophetic algorithm based on GDSC database. In the GH3 cell cultures, sunitinib significantly decrease the cell viability. In addition, since the pituitary gland sites outside the blood-brain barrier, a rat xenograft flank model was applied to explore the efficacy of sunitinib on PA in vivo. Sunitinib significantly inhibited GH tumorigenesis but did not show toxicity to rats when orally administered for 12 days. In addition, AFM-based Young’s modulus measurement also demonstrated the stiffness of PA samples were significantly reduced after sunitinib treatment, which suggested sunitinib could serve as a potential candidate drug for stiff PAs. As reported in the literature, sunitinib was developed to inhibit growth factor receptor tyrosine kinases in ECs and pericytes, which are implicated in angiogenesis (Tran et al., 2016). Sunitinib is a standard-of-care first-line therapy for patients with advanced renal cell carcinoma (Motzer et al., 2019). Sunitinib was also reported to block survival in the rat pituitary cell line GH4C (Chenlo et al., 2019). In summary, sunitinib is a potential therapeutic agent for stiff PAs. Prospective clinical studies investigating sunitinib will be necessary to assess its efficacy.
In conclusion, we assessed the molecular landscape during PA stiffening and discovered unique features of cell components and gene regulation in stiff tumors. Our results indicate that the epigenome and epitranscriptome are essential in the regulation of tumor stiffness-related RNA expression. In vitro and in vivo also provided evidence that sunitinib has the potential to reverse PA stiffening. Therefore, characterizing the mechanism by which PAs become stiff will offer a new theoretical basis for developing novel therapies for individualized treatment.
Data Availability Statement
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://ngdc.cncb.ac.cn/gsa-human/, HRA000954, HRA000955, and HRA000956.
Ethics Statement
The studies involving human participants were reviewed and approved by the Institutional Ethics Committee and Institutional Review Board of PUMCH. The patients/participants provided their written informed consent to participate in this study. The animal study was reviewed and approved by the Institutional Ethics Committee and Institutional Review Board of PUMCH. Written informed consent was obtained from the individual(s) for the publication of any potentially identifiable images or data included in this article.
Author Contributions
JG, and BX designed the study, ZW, MC, and GZ analyzed and interpreted the data, MC performed the immunofluorescence analysis, YRZ and JL performed the AFM analysis, ZW and MC wrote and edited manuscript, GZ, PL, YKW, YZ, XG, YW, XB, and WL collected the samples, YW, WM, and RW assisted to design the study, and all authors read and approved the manuscript.
Funding
This study was supported by the Chinese Academy of Medical Sciences Innovation Fund for Medical Sciences (grant number: 2017-I2M-1-001), the China Postdoctoral Science Foundation (grant number: BX2019045), the National Natural Science Foundation of China (82103302 to C. M.), and the Graduate Innovation Fund of the Chinese Academy of Medical Sciences and Peking Union Medical College (grant number: 2019-1002-73).
Conflict of Interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Publisher’s Note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Acknowledgments
Wang is grateful for the invaluable support received from his parents and from BX over the years.
Supplementary Material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fcell.2022.820562/full#supplementary-material
Supplementary Figure S1 | Clustering of PA samples and determination of soft-thresholding power. (A) The hierarchical clustering dendrogram of PA samples. (B) Analysis of the scale-free fit index for various soft-thresholding powers. The red line indicates a correlation coefficient of 0.85. (C) Analysis of the mean connectivity for various soft-thresholding powers.
Supplementary Figure S2 | Scatter plots of Module Membership (x-axis) in the turquoise module vs. Gene Significance for stiffness (A), general ECs (B), lymphatic ECs (C), microvascular ECs (D), and CAFs (E) (y-axis).
References
Aguilo, F., Zhang, F., Sancho, A., Fidalgo, M., Di Cecilia, S., Vashisht, A., et al. (2015). Coordination of M 6 A mRNA Methylation and Gene Transcription by ZFP217 Regulates Pluripotency and Reprogramming. Cell Stem Cell 17, 689–704. doi:10.1016/j.stem.2015.09.005
Almutairi, R. D., Muskens, I. S., Cote, D. J., Dijkman, M. D., Kavouridis, V. K., Crocker, E., et al. (2018). Gross Total Resection of Pituitary Adenomas after Endoscopic vs. Microscopic Transsphenoidal Surgery: a Meta-Analysis. Acta Neurochir 160, 1005–1021. doi:10.1007/s00701-017-3438-z
Amadori, E., Scala, M., Cereda, G. S., Vari, M. S., Marchese, F., Di Pisa, V., et al. (2020). Targeted Re-sequencing for Early Diagnosis of Genetic Causes of Childhood Epilepsy: the Italian Experience from the 'beyond Epilepsy' Project. Ital. J. Pediatr. 46, 92. doi:10.1186/s13052-020-00860-1
Anders, S., Pyl, P. T., and Huber, W. (2015). HTSeq--a python Framework to Work with High-Throughput Sequencing Data. Bioinformatics 31, 166–169. doi:10.1093/bioinformatics/btu638
Anlaş, A. A., and Nelson, C. M. (2020). Soft Microenvironments Induce Chemoresistance by Increasing Autophagy Downstream of Integrin-Linked Kinase. Cancer Res. 80, 4103–4113. doi:10.1158/0008-5472.CAN-19-4021
Aran, D., Hu, Z., and Butte, A. J. (2017). xCell: Digitally Portraying the Tissue Cellular Heterogeneity Landscape. Genome Biol. 18, 220. doi:10.1186/s13059-017-1349-1
Bahuleyan, B., Raghuram, L., Rajshekhar, V., and Chacko, A. G. (2006). To Assess the Ability of MRI to Predict Consistency of Pituitary Macroadenomas. Br. J. Neurosurg. 20, 324–326. doi:10.1080/02688690601000717
Bao, X., Deng, K., Liu, X., Feng, M., Chen, C. C., Lian, W., et al. (2016). Extended Transsphenoidal Approach for Pituitary Adenomas Invading the Cavernous Sinus Using Multiple Complementary Techniques. Pituitary 19, 1–10. doi:10.1007/s11102-015-0675-0
Batista, P. J., Molinie, B., Wang, J., Qu, K., Zhang, J., Li, L., et al. (2014). m6A RNA Modification Controls Cell Fate Transition in Mammalian Embryonic Stem CellsA RNA Modification Controls Cell Fate Transition in Mammalian Embryonic Stem Cells. Cell Stem Cell 15, 707–719. doi:10.1016/j.stem.2014.09.019
Bayer, S. V., Grither, W. R., Brenot, A., Hwang, P. Y., Barcus, C. E., Ernst, M., et al. (2019). DDR2 Controls Breast Tumor Stiffness and Metastasis by Regulating Integrin Mediated Mechanotransduction in CAFs. Elife 8, e45508. doi:10.7554/eLife.45508
Bindea, G., Mlecnik, B., Hackl, H., Charoentong, P., Tosolini, M., Kirilovsky, A., et al. (2009). ClueGO: a Cytoscape Plug-In to Decipher Functionally Grouped Gene Ontology and Pathway Annotation Networks. Bioinformatics 25, 1091–1093. doi:10.1093/bioinformatics/btp101
Budczies, J., Seidel, A., Christopoulos, P., Endris, V., Kloor, M., Győrffy, B., et al. (2018). Integrated Analysis of the Immunological and Genetic Status in and across Cancer Types: Impact of Mutational Signatures beyond Tumor Mutational burden. Oncoimmunology 7, e1526613. doi:10.1080/2162402x.2018.1526613
Chen, M., Wei, L., Law, C.-T., Tsang, F. H.-C., Shen, J., Cheng, C. L.-H., et al. (2018). RNA N6-Methyladenosine Methyltransferase-like 3 Promotes Liver Cancer Progression through YTHDF2-dependent Posttranscriptional Silencing of SOCS2. Hepatology 67, 2254–2270. doi:10.1002/hep.29683
Chen, T., Hao, Y.-J., Zhang, Y., Li, M.-M., Wang, M., Han, W., et al. (2015). m6A RNA Methylation Is Regulated by MicroRNAs and Promotes Reprogramming to PluripotencyA RNA Methylation Is Regulated by microRNAs and Promotes Reprogramming to Pluripotency. Cell Stem Cell 16, 289–301. doi:10.1016/j.stem.2015.01.016
Chen, Y., Peng, C., Chen, J., Chen, D., Yang, B., He, B., et al. (2019). WTAP Facilitates Progression of Hepatocellular Carcinoma via m6A-HuR-dependent Epigenetic Silencing of ETS1. Mol. Cancer 18, 127. doi:10.1186/s12943-019-1053-8
Chenlo, M., Rodriguez-Gomez, I. A., Serramito, R., Garcia-Rendueles, A. R., Villar-Taibo, R., Fernandez-Rodriguez, E., et al. (2019). Unmasking a New Prognostic Marker and Therapeutic Target from the GDNF-RET/PIT1/p14ARF/p53 Pathway in Acromegaly. EBioMedicine 43, 537–552. doi:10.1016/j.ebiom.2019.04.007
Colpaert, C., Vermeulen, P., Van Marck, E., and Dirix, L. (2001). The Presence of a Fibrotic Focus Is an Independent Predictor of Early Metastasis in Lymph Node-Negative Breast Cancer Patients. Am. J. Surg. Pathol. 25, 1557–1558. doi:10.1097/00000478-200112000-00016
Cui, Q., Shi, H., Ye, P., Li, L., Qu, Q., Sun, G., et al. (2017). m 6 A RNA Methylation Regulates the Self-Renewal and Tumorigenesis of Glioblastoma Stem CellsA RNA Methylation Regulates the Self-Renewal and Tumorigenesis of Glioblastoma Stem Cells. Cel Rep. 18, 2622–2634. doi:10.1016/j.celrep.2017.02.059
Fan, Y., Hua, M., Mou, A., Wu, M., Liu, X., Bao, X., et al. (2019). Preoperative Noninvasive Radiomics Approach Predicts Tumor Consistency in Patients with Acromegaly: Development and Multicenter Prospective Validation. Front. Endocrinol. 10, 403. doi:10.3389/fendo.2019.00403
Geeleher, P., Cox, N., and Huang, R. S. (2014). pRRophetic: an R Package for Prediction of Clinical Chemotherapeutic Response from Tumor Gene Expression Levels. PLoS One 9, e107468. doi:10.1371/journal.pone.0107468
Geeleher, P., Cox, N. J., and Huang, R. (2014). Clinical Drug Response Can Be Predicted Using Baseline Gene Expression Levels and In Vitro Drug Sensitivity in Cell Lines. Genome Biol. 15, R47. doi:10.1186/gb-2014-15-3-r47
Grady, M. E., Composto, R. J., and Eckmann, D. M. (2016). Cell Elasticity with Altered Cytoskeletal Architectures across Multiple Cell Types. J. Mech. Behav. Biomed. Mater. 61, 197–207. doi:10.1016/j.jmbbm.2016.01.022
Gu, Y., Niu, S., Wang, Y., Duan, L., Pan, Y., Tong, Z., et al. (2020). DMDRMR-mediated Regulation of m6A-Modified CDK4 by m6A Reader IGF2BP3 Drives ccRCC Progression. Cancer Res. 81, 923–934. doi:10.1158/0008-5472.can-20-1619
Guo, W., Lv, P., Liu, S., Xu, F., Guo, Y., Shen, S., et al. (2018). Aberrant Methylation-Mediated Downregulation of Long Noncoding RNA C5orf66-AS1 Promotes the Development of Gastric Cardia Adenocarcinoma. Mol. Carcinogenesis 57, 854–865. doi:10.1002/mc.22806
Han, J., Wang, J.-z., Yang, X., Yu, H., Zhou, R., Lu, H.-C., et al. (2019). METTL3 Promote Tumor Proliferation of Bladder Cancer by Accelerating Pri-miR221/222 Maturation in m6A-dependent Manner. Mol. Cancer 18, 110. doi:10.1186/s12943-019-1036-9
Hänzelmann, S., Castelo, R., and Guinney, J. (2013). GSVA: Gene Set Variation Analysis for Microarray and RNA-Seq Data. BMC Bioinformatics 14, 7.
He, L., Li, H., Wu, A., Peng, Y., Shu, G., and Yin, G. (2019). Functions of N6-Methyladenosine and its Role in Cancer. Mol. Cancer 18, 176. doi:10.1186/s12943-019-1109-9
He, Y., Jiang, Z., Chen, C., and Wang, X. (2018). Classification of Triple-Negative Breast Cancers Based on Immunogenomic Profiling. J. Exp. Clin. Cancer Res. 37, 327. doi:10.1186/s13046-018-1002-1
Hoshida, Y., Brunet, J.-P., Tamayo, P., Golub, T. R., and Mesirov, J. P. (2007). Subclass Mapping: Identifying Common Subtypes in Independent Disease Data Sets. PLoS One 2, e1195. doi:10.1371/journal.pone.0001195
Houtgast, E. J., Sima, V.-M., Bertels, K., and Al-Ars, Z. (2018). Hardware Acceleration of BWA-MEM Genomic Short Read Mapping for Longer Read Lengths. Comput. Biol. Chem. 75, 54–64. doi:10.1016/j.compbiolchem.2018.03.024
Hughes, J. D., Fattahi, N., Van Gompel, J., Arani, A., Ehman, R., and Huston, J. (2016). Magnetic Resonance Elastography Detects Tumoral Consistency in Pituitary Macroadenomas. Pituitary 19, 286–292. doi:10.1007/s11102-016-0706-5
Humphrey, J. D., Dufresne, E. R., and Schwartz, M. A. (2014). Mechanotransduction and Extracellular Matrix Homeostasis. Nat. Rev. Mol. Cel Biol 15, 802–812. doi:10.1038/nrm3896
Jiang, P., Gu, S., Pan, D., Fu, J., Sahu, A., Hu, X., et al. (2018). Signatures of T Cell Dysfunction and Exclusion Predict Cancer Immunotherapy Response. Nat. Med. 24, 1550–1558. doi:10.1038/s41591-018-0136-1
Jin, D., Guo, J., Wu, Y., Yang, L., Wang, X., Du, J., et al. (2020). m6A Demethylase ALKBH5 Inhibits Tumor Growth and Metastasis by Reducing YTHDFs-Mediated YAP Expression and Inhibiting miR-107/lats2-Mediated YAP Activity in NSCLCA Demethylase ALKBH5 Inhibits Tumor Growth and Metastasis by Reducing YTHDFs-Mediated YAP Expression and Inhibiting miR-107/lats2-Mediated YAP Activity in NSCLC. Mol. Cancer 19, 40. doi:10.1186/s12943-020-01161-1
Kagan, H. M., and Li, W. (2003). Lysyl Oxidase: Properties, Specificity, and Biological Roles inside and outside of the Cell. J. Cel. Biochem. 88, 660–672. doi:10.1002/jcb.10413
Kai, F., Laklai, H., and Weaver, V. M. (2016). Force Matters: Biomechanical Regulation of Cell Invasion and Migration in Disease. Trends Cel Biol. 26, 486–497. doi:10.1016/j.tcb.2016.03.007
Ke, S., Pandya-Jones, A., Saito, Y., Fak, J. J., Vågbø, C. B., Geula, S., et al. (2017). m6A mRNA Modifications Are Deposited in Nascent Pre-mRNA and Are Not Required for Splicing but Do Specify Cytoplasmic Turnover. Genes Dev. 31, 990–1006. doi:10.1101/gad.301036.117
Kim, D., Langmead, B., and Salzberg, S. L. (2015). HISAT: a Fast Spliced Aligner with Low Memory Requirements. Nat. Methods 12, 357–360. doi:10.1038/nmeth.3317
Kim, E. H., Oh, M. C., Chang, J. H., Moon, J. H., Ku, C. R., Chang, W.-S., et al. (2018). Postoperative Gamma Knife Radiosurgery for Cavernous Sinus-Invading Growth Hormone-Secreting Pituitary Adenomas. World Neurosurg. 110, e534–e545. doi:10.1016/j.wneu.2017.11.043
Krueger, F., and Andrews, S. R. (2011). Bismark: a Flexible Aligner and Methylation Caller for Bisulfite-Seq Applications. Bioinformatics 27, 1571–1572. doi:10.1093/bioinformatics/btr167
Laklai, H., Miroshnikova, Y. A., Pickup, M. W., Collisson, E. A., Kim, G. E., Barrett, A. S., et al. (2016). Genotype Tunes Pancreatic Ductal Adenocarcinoma Tissue Tension to Induce Matricellular Fibrosis and Tumor Progression. Nat. Med. 22, 497–505. doi:10.1038/nm.4082
Lambrechts, D., Wauters, E., Boeckx, B., Aibar, S., Nittner, D., Burton, O., et al. (2018). Phenotype Molding of Stromal Cells in the Lung Tumor Microenvironment. Nat. Med. 24, 1277–1289. doi:10.1038/s41591-018-0096-5
Lan, T., Li, H., Zhang, D., Xu, L., Liu, H., Hao, X., et al. (2019). KIAA1429 Contributes to Liver Cancer Progression through N6-methyladenosine-dependent post-transcriptional Modification of GATA3. Mol. Cancer 18, 186. doi:10.1186/s12943-019-1106-z
Langfelder, P., and Horvath, S. (2008). WGCNA: an R Package for Weighted Correlation Network Analysis. BMC Bioinformatics 9, 559. doi:10.1186/1471-2105-9-559
Langmead, B., and Salzberg, S. L. (2012). Fast Gapped-Read Alignment with Bowtie 2. Nat. Methods 9, 357–359. doi:10.1038/nmeth.1923
Levental, K. R., Yu, H., Kass, L., Lakins, J. N., Egeblad, M., Erler, J. T., et al. (2009). Matrix Crosslinking Forces Tumor Progression by Enhancing Integrin Signaling. Cell 139, 891–906. doi:10.1016/j.cell.2009.10.027
Li, H. (2011). A Statistical Framework for SNP Calling, Mutation Discovery, Association Mapping and Population Genetical Parameter Estimation from Sequencing Data. Bioinformatics 27, 2987–2993. doi:10.1093/bioinformatics/btr509
Li, L., Zang, L., Zhang, F., Chen, J., Shen, H., Shu, L., et al. (2017). Fat Mass and Obesity-Associated (FTO) Protein Regulates Adult Neurogenesis. Hum. Mol. Genet. 26, 2398–2411. doi:10.1093/hmg/ddx128
Li, M., Zhao, X., Wang, W., Shi, H., Pan, Q., Lu, Z., et al. (2018). Ythdf2-mediated m6A mRNA Clearance Modulates Neural Development in Mice. Genome Biol. 19, 69. doi:10.1186/s13059-018-1436-y
Li, T., Zhuang, Y., Yang, W., Xie, Y., Shang, W., Su, S., et al. (2021). Silencing of METTL3 Attenuates Cardiac Fibrosis Induced by Myocardial Infarction via Inhibiting the Activation of Cardiac Fibroblasts. FASEB J. 35, e21162. doi:10.1096/fj.201903169r
Li, Y., Xiao, J., Bai, J., Tian, Y., Qu, Y., Chen, X., et al. (2019). Molecular Characterization and Clinical Relevance of m6A Regulators across 33 Cancer Types. Mol. Cancer 18, 137. doi:10.1186/s12943-019-1066-3
Lim, Z.-F., and Ma, P. C. (2019). Emerging Insights of Tumor Heterogeneity and Drug Resistance Mechanisms in Lung Cancer Targeted Therapy. J. Hematol. Oncol. 12, 134. doi:10.1186/s13045-019-0818-2
Liu, P., Zhang, B., Chen, Z., He, Y., Du, Y., Liu, Y., et al. (2020). m6A-induced lncRNA MALAT1 Aggravates Renal Fibrogenesis in Obstructive Nephropathy through the miR-145/FAK pathwayA-Induced lncRNA MALAT1 Aggravates Renal Fibrogenesis in Obstructive Nephropathy through the miR-145/FAK Pathway. Aging 12, 5280–5299. doi:10.18632/aging.102950
Ma, C., Chang, M., Lv, H., Zhang, Z.-W., Zhang, W., He, X., et al. (2018). RNA m6A Methylation Participates in Regulation of Postnatal Development of the Mouse Cerebellum. Genome Biol. 19, 68. doi:10.1186/s13059-018-1435-z
Ma, J. z., Yang, F., Zhou, C. c., Liu, F., Yuan, J. h., Wang, F., et al. (2017). METTL14 Suppresses the Metastatic Potential of Hepatocellular Carcinoma by Modulating N 6 ‐methyladenosine‐dependent Primary MicroRNA Processing. Hepatology 65, 529–543. doi:10.1002/hep.28885
Mayakonda, A., Lin, D.-C., Assenov, Y., Plass, C., and Koeffler, H. P. (2018). Maftools: Efficient and Comprehensive Analysis of Somatic Variants in Cancer. Genome Res. 28, 1747–1756. doi:10.1101/gr.239244.118
McGranahan, N., and Swanton, C. (2017). Clonal Heterogeneity and Tumor Evolution: Past, Present, and the Future. Cell 168, 613–628. doi:10.1016/j.cell.2017.01.018
McKenna, A., Hanna, M., Banks, E., Sivachenko, A., Cibulskis, K., Kernytsky, A., et al. (2010). The Genome Analysis Toolkit: a MapReduce Framework for Analyzing Next-Generation DNA Sequencing Data. Genome Res. 20, 1297–1303. doi:10.1101/gr.107524.110
McLaren, W., Gil, L., Hunt, S. E., Riat, H. S., Ritchie, G. R. S., Thormann, A., et al. (2016). The Ensembl Variant Effect Predictor. Genome Biol. 17, 122. doi:10.1186/s13059-016-0974-4
Mering, C. v., Huynen, M., Jaeggi, D., Schmidt, S., Bork, P., and Snel, B. (2003). STRING: a Database of Predicted Functional Associations between Proteins. Nucleic Acids Res. 31, 258–261. doi:10.1093/nar/gkg034
Molitch, M. E. (2017). Diagnosis and Treatment of Pituitary Adenomas. JAMA 317, 516–524. doi:10.1001/jama.2016.19699
Motzer, R. J., Penkov, K., Haanen, J., Rini, B., Albiges, L., Campbell, M. T., et al. (2019). Avelumab Plus Axitinib versus Sunitinib for Advanced Renal-Cell Carcinoma. N. Engl. J. Med. 380, 1103–1115. doi:10.1056/nejmoa1816047
Mouw, J. K., Ou, G., and Weaver, V. M. (2014). Extracellular Matrix Assembly: a Multiscale Deconstruction. Nat. Rev. Mol. Cel Biol 15, 771–785. doi:10.1038/nrm3902
Mroz, E. A., Tward, A. M., Hammon, R. J., Ren, Y., and Rocco, J. W. (2015). Intra-tumor Genetic Heterogeneity and Mortality in Head and Neck Cancer: Analysis of Data from the Cancer Genome Atlas. Plos Med. 12, e1001786. doi:10.1371/journal.pmed.1001786
Newman, A. M., Liu, C. L., Green, M. R., Gentles, A. J., Feng, W., Xu, Y., et al. (2015). Robust Enumeration of Cell Subsets from Tissue Expression Profiles. Nat. Methods 12, 453–457. doi:10.1038/nmeth.3337
Niu, Y., Lin, Z., Wan, A., Chen, H., Liang, H., Sun, L., et al. (2019). RNA N6-Methyladenosine Demethylase FTO Promotes Breast Tumor Progression through Inhibiting BNIP3. Mol. Cancer 18, 46. doi:10.1186/s12943-019-1004-4
Northcott, J. M., Dean, I. S., Mouw, J. K., and Weaver, V. M. (2018). Feeling Stress: the Mechanics of Cancer Progression and Aggression. Front. Cel Dev. Biol. 6, 17. doi:10.3389/fcell.2018.00017
Ostrom, Q. T., Gittleman, H., Liao, P., Rouse, C., Chen, Y., Dowling, J., et al. (2014). CBTRUS Statistical Report: Primary Brain and central Nervous System Tumors Diagnosed in the United States in 2007-2011. Neuro Oncol. 16 Suppl 4 (Suppl. 4), iv1–63. doi:10.1093/neuonc/nou223
Park, Y., and Wu, H. (2016). Differential Methylation Analysis for BS-Seq Data under General Experimental Design. Bioinformatics 32, 1446–1453. doi:10.1093/bioinformatics/btw026
Paszek, M. J., Zahir, N., Johnson, K. R., Lakins, J. N., Rozenberg, G. I., Gefen, A., et al. (2005). Tensional Homeostasis and the Malignant Phenotype. Cancer Cell 8, 241–254. doi:10.1016/j.ccr.2005.08.010
Pease, M., Ling, C., Mack, W. J., Wang, K., and Zada, G. (2013). The Role of Epigenetic Modification in Tumorigenesis and Progression of Pituitary Adenomas: a Systematic Review of the Literature. PLoS One 8, e82619. doi:10.1371/journal.pone.0082619
Pelham, C. J., Nagane, M., and Madan, E. (2020). Cell Competition in Tumor Evolution and Heterogeneity: Merging Past and Present. Semin. Cancer Biol. 63, 11–18. doi:10.1016/j.semcancer.2019.07.008
Pickup, M. W., Mouw, J. K., and Weaver, V. M. (2014). The Extracellular Matrix Modulates the Hallmarks of Cancer. EMBO Rep. 15, 1243–1253. doi:10.15252/embr.201439246
Reddy, R., Cudlip, S., Byrne, J. V., Karavitaki, N., and Wass, J. A. H. (2011). Can We Ever Stop Imaging in Surgically Treated and Radiotherapy-Naive Patients with Non-functioning Pituitary Adenoma? Eur. J. Endocrinol. 165, 739–744. doi:10.1530/eje-11-0566
Robinson, M. D., McCarthy, D. J., and Smyth, G. K. (2010). edgeR: a Bioconductor Package for Differential Expression Analysis of Digital Gene Expression Data. Bioinformatics 26, 139–140. doi:10.1093/bioinformatics/btp616
Sailer, V., Charpentier, A., Dietrich, J., Vogt, T. J., Franzen, A., Bootz, F., et al. (2018). Intragenic DNA Methylation of PITX1 and the Adjacent Long Non-coding RNA C5orf66-AS1 Are Prognostic Biomarkers in Patients with Head and Neck Squamous Cell Carcinomas. PLoS One 13, e0192742. doi:10.1371/journal.pone.0192742
Salomon, M. P., Wang, X., Marzese, D. M., Hsu, S. C., Nelson, N., Zhang, X., et al. (2018). The Epigenomic Landscape of Pituitary Adenomas Reveals Specific Alterations and Differentiates Among Acromegaly, cushing's Disease and Endocrine-Inactive Subtypes. Clin. Cancer Res. 24, 4126–4136. doi:10.1158/1078-0432.ccr-17-2206
Shen, Y., Wang, X., Lu, J., Salfenmoser, M., Wirsik, N. M., Schleussner, N., et al. (2020). Reduction of Liver Metastasis Stiffness Improves Response to Bevacizumab in Metastatic Colorectal Cancer. Cancer Cell 37, 800–817. e7. doi:10.1016/j.ccell.2020.05.005
Smyth, G. K., Michaud, J., and Scott, H. S. (2005). Use of Within-Array Replicate Spots for Assessing Differential Expression in Microarray Experiments. Bioinformatics 21, 2067–2075. doi:10.1093/bioinformatics/bti270
Subramanian, A., Narayan, R., Corsello, S. M., Peck, D. D., Natoli, T. E., Lu, X., et al. (2017). A Next Generation Connectivity Map: L1000 Platform and the First 1,000,000 Profiles. Cell 171, 1437–1452. doi:10.1016/j.cell.2017.10.049
Subramanian, A., Tamayo, P., Mootha, V. K., Mukherjee, S., Ebert, B. L., Gillette, M. A., et al. (2005). Gene Set Enrichment Analysis: a Knowledge-Based Approach for Interpreting Genome-wide Expression Profiles. Proc. Natl. Acad. Sci. 102, 15545–15550. doi:10.1073/pnas.0506580102
Sughrue, M. E., Chang, E. F., Gabriel, R. A., Aghi, M. K., and Blevins, L. S. (2011). Excess Mortality for Patients with Residual Disease Following Resection of Pituitary Adenomas. Pituitary 14, 276–283. doi:10.1007/s11102-011-0308-1
Tabaee, A., Anand, V. K., Barrón, Y., Hiltzik, D. H., Brown, S. M., Kacker, A., et al. (2009). Endoscopic Pituitary Surgery: a Systematic Review and Meta-Analysis. Jns 111, 545–554. doi:10.3171/2007.12.17635
Tran, T. A., Leong, H. S., Pavia-Jimenez, A., Fedyshyn, S., Yang, J., Kucejova, B., et al. (2016). Fibroblast Growth Factor Receptor-dependent and -independent Paracrine Signaling by Sunitinib-Resistant Renal Cell Carcinoma. Mol. Cel Biol 36, 1836–1855. doi:10.1128/mcb.00189-16
Vinay, D. S., Ryan, E. P., Pawelec, G., Talib, W. H., Stagg, J., Elkord, E., et al. (2015). Immune Evasion in Cancer: Mechanistic Basis and Therapeutic Strategies. Semin. Cancer Biol. 35 (Suppl. l), S185–S198. doi:10.1016/j.semcancer.2015.03.004
Wang, C.-X., Cui, G.-S., Liu, X., Xu, K., Wang, M., Zhang, X.-X., et al. (2018). METTL3-mediated m6A Modification Is Required for Cerebellar Development. Plos Biol. 16, e2004880. doi:10.1371/journal.pbio.2004880
Wang, X., Lu, Z., Gomez, A., Hon, G. C., Yue, Y., Han, D., et al. (2014). N6-methyladenosine-dependent Regulation of Messenger RNA Stability. Nature 505, 117–120. doi:10.1038/nature12730
Wang, Z., Guo, X., Gao, L., Wang, Y., Guo, Y., Xing, B., et al. (2021). Classification of Pediatric Gliomas Based on Immunological Profiling: Implications for Immunotherapy Strategies. Mol. Ther. - Oncolytics 20, 34–47. doi:10.1016/j.omto.2020.12.012
Yao, A., Rutland, J. W., Verma, G., Banihashemi, A., Padormo, F., Tsankova, N. M., et al. (2020). Pituitary Adenoma Consistency: Direct Correlation of Ultrahigh Field 7T MRI with Histopathological Analysis. Eur. J. Radiol. 126, 108931. doi:10.1016/j.ejrad.2020.108931
Yoon, K.-J., Ringeling, F. R., Vissers, C., Jacob, F., Pokrass, M., Jimenez-Cyrus, D., et al. (2017). Temporal Control of Mammalian Cortical Neurogenesis by m6A Methylation. Cell 171, 877–889. doi:10.1016/j.cell.2017.09.003
Yoshihara, K., Shahmoradgoli, M., Martínez, E., Vegesna, R., Kim, H., Torres-Garcia, W., et al. (2013). Inferring Tumour Purity and Stromal and Immune Cell Admixture from Expression Data. Nat. Commun. 4, 2612. doi:10.1038/ncomms3612
Yu, G., Li, C., Xie, W., Wang, Z., Gao, H., Cao, L., et al. (2017). Long Non-coding RNA C5orf66-AS1 Is Downregulated in Pituitary Null Cell Adenomas and Is Associated with Their Invasiveness. Oncol. Rep. 38, 1140–1148. doi:10.3892/or.2017.5739
Yuan, H., Yan, M., Zhang, G., Liu, W., Deng, C., Liao, G., et al. (2019). CancerSEA: a Cancer Single-Cell State Atlas. Nucleic Acids Res. 47, D900–D908. doi:10.1093/nar/gky939
Zhang, B., and Horvath, S. (2005). A General Framework for Weighted Gene Co-expression Network Analysis. Stat. Appl. Genet. Mol. Biol. 4, Article17. doi:10.2202/1544-6115.1128
Zhang, C., Chen, Y., Sun, B., Wang, L., Yang, Y., Ma, D., et al. (2017). m6A Modulates Haematopoietic Stem and Progenitor Cell Specification. Nature 549, 273–276. doi:10.1038/nature23883
Zhang, C., Samanta, D., Lu, H., Bullen, J. W., Zhang, H., Chen, I., et al. (2016). Hypoxia Induces the Breast Cancer Stem Cell Phenotype by HIF-dependent and ALKBH5-Mediated m6A-Demethylation of NANOG mRNA. Proc. Natl. Acad. Sci. USA 113, E2047–E2056. doi:10.1073/pnas.1602883113
Zhang, S., Zhao, B. S., Zhou, A., Lin, K., Zheng, S., Lu, Z., et al. (2017). m 6 A Demethylase ALKBH5 Maintains Tumorigenicity of Glioblastoma Stem-like Cells by Sustaining FOXM1 Expression and Cell Proliferation ProgramA Demethylase ALKBH5 Maintains Tumorigenicity of Glioblastoma Stem-like Cells by Sustaining FOXM1 Expression and Cell Proliferation Program. Cancer Cell 31, 591–606. doi:10.1016/j.ccell.2017.02.013
Zhao, B., Wei, Y.-K., Li, G.-L., Li, Y.-N., Yao, Y., Kang, J., et al. (2010). Extended Transsphenoidal Approach for Pituitary Adenomas Invading the Anterior Cranial Base, Cavernous Sinus, and Clivus: a Single-center Experience with 126 Consecutive Cases. Jns 112, 108–117. doi:10.3171/2009.3.jns0929
Keywords: pituitary adenoma, stiffness-related gene, DNA methylation, m6A, sunitinib
Citation: Wang Z, Chang M, Zhang Y, Zhou G, Liu P, Lou J, Wang Y, Zhang Y, Guo X, Wang Y, Bao X, Lian W, Wang Y, Wang R, Ma W, Xing B and Gao J (2022) Multi-Omics Investigations Revealed Underlying Molecular Mechanisms Associated With Tumor Stiffness and Identified Sunitinib as a Potential Therapy for Reducing Stiffness in Pituitary Adenomas. Front. Cell Dev. Biol. 10:820562. doi: 10.3389/fcell.2022.820562
Received: 23 November 2021; Accepted: 01 March 2022;
Published: 15 March 2022.
Edited by:
Mantang Qiu, Peking University People’s Hospital, ChinaCopyright © 2022 Wang, Chang, Zhang, Zhou, Liu, Lou, Wang, Zhang, Guo, Wang, Bao, Lian, Wang, Wang, Ma, Xing and Gao. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Jun Gao, gaojunpumc@163.com; Bing Xing, xingbingemail@aliyun.com
†These authors have contributed equally to this work