ΔNp63 drives epithelial differentiation in glioma

Dear Editor, Glioma is known for its morphologic and molecular diversity. Clinically documented glioma with epithelial differentiation is rare and linked to poor prognosis with unknown mechanisms.1,2 We identified an epigenetically modulated and stress-related glioma epithelial transcriptional program showing correlation with aggressiveness, recurrence and prognosis in glioma patients, and highlighted ΔNp63 as a likely driver of this cellular plasticity. The study flow diagram is shown in Figure S1. We applied the GENIE3 algorithm3 to transcriptomes of 212 heterogeneous Gene Expression Omnibus brain samples with genes prefiltered by correlation (Pearson’s R > 0.8) for regulatory network reconstruction. Regulons, namely regulatory relationships, were established concerning 37 core transcription factors (TFs) (top 1% most confident regulations and ≥10 targets). Based on connections among TFs (Table S1), we aggregated the regulons to eight clusters denoted C1 to C8. Nearly all these regulon clusters were univariately associated with prognosis in TCGA glioma patients (Figure S2A). However, only one regulon cluster, the C6 containing 8 TFs and 76 targets, showed independent prognostic significance after adjustment for established glioma prognosticators (Figure S2B). Markers as well as regulators for epithelial differentiation were enriched in C6 genes. We further deconvoluted C6 into two major subsets of genes (Figure S3), with the larger (n= 37) as C6a, whichmainly consisted of end targets including epithelial/keratinocyte markers (e.g., cytokeratins), and the smaller (n = 18) as C6b, which mainly consisted of upstream regulators including TP63 (key regulator of basal/squamous subtype carcinoma) and MACC1 (coactivator of the receptor tyrosine kinase gene MET also involved in basal/squamous subtype carcinoma)4 (Table S2). Higher C6b expression was found in glioblastoma (GBM) versus low grade glioma and in recurrent versus primary tumors (Figure S2C). High C6b was an independent predictor for unfavorable prognosis in glioma (Figure S2D). We developed a network-based approach to prioritize key C6 regulators using scRNA-seq data (n = 658). MET


ΔNp63 drives epithelial differentiation in glioma
Dear Editor, Glioma is known for its morphologic and molecular diversity. Clinically documented glioma with epithelial differentiation is rare and linked to poor prognosis with unknown mechanisms. 1,2 We identified an epigenetically modulated and stress-related glioma epithelial transcriptional program showing correlation with aggressiveness, recurrence and prognosis in glioma patients, and highlighted ΔNp63 as a likely driver of this cellular plasticity. The study flow diagram is shown in Figure S1.
We applied the GENIE3 algorithm 3 to transcriptomes of 212 heterogeneous Gene Expression Omnibus brain samples with genes prefiltered by correlation (Pearson's R > 0.8) for regulatory network reconstruction. Regulons, namely regulatory relationships, were established concerning 37 core transcription factors (TFs) (top 1% most confident regulations and ≥10 targets). Based on connections among TFs (Table S1), we aggregated the regulons to eight clusters denoted C1 to C8. Nearly all these regulon clusters were univariately associated with prognosis in TCGA glioma patients ( Figure S2A). However, only one regulon cluster, the C6 containing 8 TFs and 76 targets, showed independent prognostic significance after adjustment for established glioma prognosticators ( Figure S2B). Markers as well as regulators for epithelial differentiation were enriched in C6 genes. We further deconvoluted C6 into two major subsets of genes ( Figure  S3), with the larger (n = 37) as C6a, which mainly consisted of end targets including epithelial/keratinocyte markers (e.g., cytokeratins), and the smaller (n = 18) as C6b, which mainly consisted of upstream regulators including TP63 (key regulator of basal/squamous subtype carcinoma) and MACC1 (coactivator of the receptor tyrosine kinase gene MET also involved in basal/squamous subtype carcinoma) 4 (Table S2). Higher C6b expression was found in glioblastoma (GBM) versus low grade glioma and in recurrent versus primary tumors ( Figure S2C). High C6b was an independent predictor for unfavorable prognosis in glioma ( Figure S2D).
We developed a network-based approach to prioritize key C6 regulators using scRNA-seq data (n = 658). MET and its canonical ligand HGF, as well as a classical regulator for glioma stemness and metastasis, CD44, were considered in addition to C6b genes. 5 Single cells were clustered to three subpopulations ( Figure 1A, left). The highest C6b expression and rate of cells expressing CD44, TP63, MET, and MACC1 were found in one single-cell cluster ( Figure 1A, middle and right; Table S3). Coexpression dependency network analysis demonstrated TP63, MET, and MACC1 as the top three regulators ( Figure 1B, left; Table S4). Correlation dynamics analysis revealed higher connectivity among these three genes in samples presenting higher expression, suggesting positive regulatory circuit ( Figure 1B, right).
Genetic alterations of TP63, MET, and MACC1 were rare in glioma ( Figure S4). Thus, an epigenetically mediated activation was suspected. We examined the CpG methylation of TP63, MET, and MACC1 in TCGA samples from top and bottom quarter regarding C6b expression. We found for all these genes a significantly lower CpG methylation level in the top quarter ( Figure 1C, left). However, in GSC models, there was no difference in CpG methylation between the three models showing highest overall expression (GSC17, GSC18, and GSC20) and the one expressing none of the three genes (GSC36), suggesting a role of other epigenetic mechanisms ( Figure 1C, right). We found multiple H3K27ac modification peaks detected around the TP63, MET, and MACC1 genomic regions in GSC17, GSC18, and GSC20, but not in GSC36 ( Figure 1D). These observations suggest coordinated epigenetics for switching on the expression of the key regulators.
Differential expression and functional impact of p63 isoforms have been reported. 6 We found the ΔNp63 was the primary isoform expressed in glioma bulk tumors (ΔNp63, n = 589; TAp63, n = 46; Figure S5). Significant positive correlation was found between expression of epithelial markers and ΔNp63 but not TAp63 (Figure 2A), and higher expression of these markers were presented in tumors expressing ΔNp63 compared to those expressing TAp63 at matched level ( Figure 2B). No co-occurrence of these two isoform types was revealed at single cell level, and cells expressing ΔNp63 showed significantly  . F, Association between TP63 genetic dependency (CRISPR/CAS9-mediated knockout effect on cell viability), TP63 expression, and TP63 isoform constitution in glioma cell lines with adequate TP63 expression (log2(TPM+1) > 1.5; n = 10). The two cell lines presenting high TP63 expression with ΔNp63 as dominant showed significant decrease in cell viability following TP63 knockout, while almost no effect in others decreased expression and altered exon usage of genes involved in cell cycle and RNA splicing compared to cells expressing TAp63, indicating a stressed phenotype as according to a previous report that stress responses induced by ɣ-irradiation, hypoxia, and chemotherapy are all associated with downregulation and splicing alteration of genes involved in cell cycling as well as pre-mRNA splicing machinery in cancer cells ( Figure 2C-E). 7 We also found in bulk tumors a significant positive correlation between expression of C6 genes involved in response to stress and ΔNp63 but not TAp63 ( Figure S6). Thus, selection pressure by hypoxia and chemoradiation may via stress response cascade lead to resistance linked with lineage plasticity such as epithelial differentiation. In glioma cell lines with adequate TP63 expression (n = 10), we observed a positive correlation between expression of epithelial markers and ΔNp63 but not TAp63 transcripts, coherent with findings in bulk tumors ( Figure S7). CRISPR/Cas9-mediated functional screening revealed that only the two cell lines presenting high TP63 expression with ΔNp63 as dominant showed significant decrease in cell viability following TP63 knockout, while almost no effect in others ( Figure 2F; Table S5). These data altogether demonstrated stress-related ΔNp63's isoform-specific role as a master regulator driving epithelial differentiation in glioma.
Tumor heterogeneity and phenotypic plasticity are important aspects of cancer hall marks. 8 Glioma is well known for its morphologic and molecular diversity. 9 Our findings suggest an epithelial transcriptional program that exists in glioma and show correlation with disease progression, recurrence, and patients' survival. Our data highlighted the ΔN isoforms of TP63 as a likely master regulator driving this cellular plasticity. Potential roles of epigenetic remodeling and cell stress response was also suggested, consistent with similar findings in previous studies. 10 Glioma cells with epithelial trans-differentiation potential may contribute to the repertoire of tumor progression, treatment resistance, and metastasis, under natural or therapeutic selection pressure. The novel insights brought by this study form a new perspective and basis for further mechanistic and translational studies aiming to improve the understanding and management of glioma.
One limitation of this work is the retrospective nature of data used. Our multilevel multi-omics strategy and rigorous statistical design could mitigate its potential effects. Besides, mechanistic exploration is limited at current stage. Nevertheless, the preliminary functional observations reported open a new path to fundamental glioma research.
In summary, we for the first time reported an epigenetically modulated and stress-related glioma epithelial transcriptional program showing clinical significance, with ΔNp63 as a likely driving master regulator that interacts with MET signaling. Our findings deepen the understanding of glioma plasticity and provide novel insights for its management and research.

A C K N O W L E D G E M E N T
This work was supported by the Improvement Project for Theranostic ability on Difficulty Miscellaneous disease (Tumor) of Zhongnan Hospital of Wuhan University (ZLYNXM202008), and National Natural Science Foundation of China (No. 81672114). This work was also funded by Medical talented youth development project in Health Commission of Hubei Province. XYM was supported by a fellowship from ITMO Cancer AVIESAN within the framework of Cancer Plan.

C O N F L I C T O F I N T E R E S T
The authors declare that they have no conflict of interest.

A U T H O R C O N T R I B U T I O N S
Jia-Hao Lu, Shuo Li, Wen-Jie Zhang, and Qi Zhang were involved in study concept and design, acquisition of data, analysis and interpretation of data, and drafting of the manuscript. Chun-Ming Liu, Shu-Yang Jiang, and Qiao-Hua Xiong were involved in critical revision of the manuscript for important intellectual content and material support. Xiang-Yu Meng and Fu-Bing Wang were involved in study concept and design, drafting of the manuscript, critical revision of the manuscript for important intellectual content, statistical analysis, obtaining funding, material support, and study supervision. All authors read and approved the final manuscript.

D ATA AVA I L A B I L I T Y S TAT E M E N T
All the data obtained and/or analyzed associated with the current study were available from the corresponding authors upon reasonable request.