NELL2-PAX7 Transcriptional Cascade Suggests Activation Mechanism for RAD52-Dependent Alternative Lengthening of Telomeres During Malignant Transformation of Malignant Peripheral Nerve Sheath Tumors: Elongation of Telomeres and Poor Survival

TL;DR

This study identifies a gene cascade that helps cancerous nerve tumors maintain telomeres, leading to worse patient outcomes.

Contribution

The discovery of a NELL2-PAX7-RAD52 transcriptional cascade activating telomere elongation in MPNSTs is novel.

Findings

NELL2 activates PAX7, which in turn activates RAD52 to elongate telomeres in MPNSTs.

RAD52-dependent telomere elongation is associated with poor survival outcomes in MPNST patients.

Increased expression of PAX7 and RAD52 correlates with worse metastasis-free and overall survival.

Abstract

Background: In eukaryotes with a double-stranded linear DNA genome, the loss of terminal DNA during replication is inevitable due to an end-replication problem; here, telomeres serve as a buffer against DNA loss. Thus, the activation of the telomere maintenance mechanism (TMM) is a prerequisite for malignant transformation. Methods: We compared neurofibroma (NF, benign) and malignant peripheral nerve sheath tumors (MPNSTs) occurring in the same patient with type 1 neurofibromatosis, where each NF–MPNST pair shared the same genetic background and differentiation lineage; this minimizes the genetic bias and contrasts only those changes that are related to malignant transformation. A total of 20 NF–MPNST pairs from 20 NF1 patients were analyzed. Whole-transcriptome sequencing (WTS) was conducted to profile the transcriptional relationship, and whole-genome sequencing (WGS) was performed to measure the telomere length. Results: We identified 22 differentially expressed genes (DEGs) during the malignant transformation of MPNSTs. Among them, NELL2 activated PAX7, which sequentially activated RAD52, the recombinase of RAD52-dependent alternative lengthening of telomeres (ALT). RAD52 elongated MPNSTs–telomeres (p = 0.017). Otherwise, neither NELL2 nor PAX7 affected telomere length (p = 0.647 and p = 0.354, respectively). RAD52 increased MPNSTs–telomeres length, independently of NELL2 and PAX7 in multiple analyses (p = 0.021). The group with increased telomere length during the malignant transformation showed inferior overall survival (OS) (HR = 3.809, p = 0.038) to the group without increased telomere length. Accordingly, the group with increased PAX7 showed inferior OS (HR = 4.896, p = 0.046) and metastasis-free survival (MFS) (HR = 9.129, p = 0.007) in comparison to the group without increased PAX7; the group with increased RAD52 showed inferior MFS (HR = 8.669, p = 0.011) in comparison to the group without increased RAD52. Conclusions: We suggest that the NELL2-PAX7 transcriptional cascade activates RAD52-dependent ALT to increase telomere length during the malignant transformation of MPNSTs, resulting in a poor prognosis.

Figures

A). We focused on the paired box (PAX) genes because they are essential for neural development and we hypothesized that differentially expressed genes (DEGs) activate PAX genes and turn on their unorthodox function, driving oncogenesis during the malignant transformation of MPNSTs. As a result, we suggest that the *NELL2*-*PAX7* transcriptional cascade increases telomere length during the malignant transformation of MPNSTs by activating RAD52-dependent ALT, which results in poor survival.

Full text

1. Introduction

Telomeres are essential structures on the ends of the chromosomes of eukaryotes with double-stranded linear DNA genomes, and they assist in overcoming end-replication problems [1]. Due to the 5′ to 3′ unidirectionality of DNA polymerization, an Okazaki fragment with an RNA primer is required on lagging strands [2] and DNA loss corresponding to the terminal RNA primer on the lagging strand is inevitable [1]. The cumulative effects of DNA loss during DNA replication are responsible for replicative senescence. Telomeres are located at both ends of the chromosomes and serve as a buffer against DNA loss through DNA replication and they maintain genomic stability.

The Fédération Nationale des Centres de Lutte Contre Le Cancer (FNCLCC) grading system consists of three categories: tumor differentiation, mitotic count, and tumor necrosis [3,4]. The mitotic count is the underlying determinant for tumor dedifferentiation and necrosis, and unlimited mitosis is the most fundamental cellular feature of malignant tumors. The activation of the telomere maintenance mechanism (TMM) is essential for overcoming the Hayflick limit and achieving immortalization. The first identified TMM was a type of reverse transcriptase, later called telomerase [5,6]. Human telomerase is composed of TERT, TERC, DKC1, and TEP1. Shelterin, which includes TRF1, TRF2, TPP1, POT1, TIN2, and RAP1, supports telomerase recruitment to telomeres and the protection of telomeres by forming secondary telomeric structures, such as T- and D-loops [7]. The H/ACA snoRNP complexes associated with DKC1, NHP2, NIP10, GAR1, and NAF1, together with hTERT, constitute the catalytic core of telomerase [8]. Approximately 85 to 90% of human cancers exhibit increased telomerase activity, whereas the remaining 10 to 15% lack telomerase activity and adopt the alternative lengthening of telomeres (ALT) [9,10,11]. In addition to the ALT, which is based on homologous recombination (HR), break-induced replication (BIR) [12,13,14] and mitotic DNA synthesis (MiDAS) [13] have also been reported to be associated mechanisms by which RAD52-dependent ALT overcomes deficient telomerase activity.

Type 1 neurofibromatosis (NF1) is an autosomal dominant disorder including a tumor suppressor gene, neurofibromin 1, located on 17q11.2, which encodes a GTPase-activating protein involved in the RAS/MAPK pathway [15]. The average global prevalence of NF1 is approximately 1 in 3000 individuals, and individuals harboring NF1 have a life expectancy of 8 to 21 years less than the general population [16]. The most fatal clinical manifestation of NF1 is a malignant peripheral nerve sheath tumor (MPNST), which involves the malignant transformation of preexisting neurofibromas. The cumulative lifetime risk for a malignant transformation in NF1 patients is estimated to be approximately 8–13%. MPNSTs are rare sarcomas that account for approximately 5–10% of all soft-tissue sarcomas. Notably, both telomerase activity and ALT have been observed in MPNSTs, instead of the ALT replacing telomerase activity [17,18], which indicates that MPNSTs are optimal tumors in comprehensive analyses of the TMM.

The ALT is a major TMM of malignant tumors originating from the mesenchyme and neuroepithelium. It has been reported that there are two major independent ALT pathways with distinct mechanisms and times of action: the RAD51- and RAD52-dependent ALT pathways. RAD51-dependent ALT promotes semiconservative homologous recombination (HR) with high fidelity [19], whereas RAD52-dependent ALT promotes conservative HR as well as breaks-induced replication-related DNA synthesis [20,21]. This conservative DNA synthesis mechanism is observed only in the late phase of the cell cycle, G2/M [12,13]. These findings suggest that RAD51- and RAD52-dependent ALT are regulated independently of each other. Since the molecular mechanisms of the ALT share a DNA repair process and DNA damage signals, DNA damage signaling has been characterized as an initiator of the ALT. Indeed, the association between the DNA damage response and RAD51-dependent ALT has been well documented [22,23]. However, the activation signals for RAD52-dependent ALT have not been elucidated and remain an outstanding question in telomere biology.

New strategies were employed in this study as follows. First, we designed a genetic analysis model conceptually similar to a “subtractive cDNA library”. We focused on the MPNSTs associated with NF1 and compared the MPNSTs and NFs from the same NF1 patient (Figure 1A). Since the MPNSTs associated with NF1 arose from preexisting NFs, each NF–MPNST pair in this model shared the same genetic background and differentiation lineage, minimizing the bias due to genetic variation and differentiation. Furthermore, the comparison of each NF–MPNST pair contrasted only the malignant transformation-related changes, structurally subtracting those changes from the normal cell to the benign tumor (Figure 1B). Second, we designed a combined statistical analysis method in which we analyzed the interaction effects on the telomere length as well as the transcriptional relationships for all the gene pairs, distinguishing the TMM-related transcriptional relationships and filtering out those that were unrelated under a strict false-positive control with FDR = 0.1 (Figure 1C). A total of 20 NF–MPNST pairs from 20 NF1 patients were analyzed to evaluate the activation of the TMM during an MPNST malignant transformation. Whole-transcriptome sequencing (WTS) was performed to profile the genome-wide gene expression and whole-genome sequencing (WGS) was performed to measure the telomere length.

The dysregulation of the major transcription factors (MTFs) involved in development may mediate the aberrant gene expression associated with oncogenesis in cancer cells [24]. Although they are essential for differentiation during development (orthodox function), the unexpected activation of MTFs in cancer cells has recently been reported to induce the expression of oncogenes in various cancers (unorthodox function) (Figure 3A). We focused on the paired box (PAX) genes because they are essential for neural development and we hypothesized that differentially expressed genes (DEGs) activate PAX genes and turn on their unorthodox function, driving oncogenesis during the malignant transformation of MPNSTs. As a result, we suggest that the NELL2-PAX7 transcriptional cascade increases telomere length during the malignant transformation of MPNSTs by activating RAD52-dependent ALT, which results in poor survival.

2. Materials and Methods

2.1. Patients and Samples

We analyzed 20 patients with NF1 who underwent surgery for both NFs and MPNSTs from November 2000 to July 2017. The clinical characteristics of the patients and tumors (NFs and MPNSTs) are described in Supplementary Material S1. Only 4 of the NF–MPNST pairs were derived from fresh-frozen tissues, and the remaining 16 NF–MPNST pairs were derived from formalin-fixed paraffin-embedded (FFPE) tissues. The FFPE samples were prepared as slides for histopathological qualification.

2.2. Whole-Genome Sequencing (WGS)

The integrity of the genomic DNA was checked by agarose gel electrophoresis and gDNA was quantified using Quant-IT PicoGreen (Invitrogen, Waltham, MA, USA). The sequencing libraries were prepared according to the instructions for the TruSeq DNA Nano Library Prep Kit. Briefly, genomic DNA (100 ng) was fragmented using adaptive-focused acoustic technology (Covaris) and end-repaired to create 5-phosphorylated blunt-ended dsDNA molecules. Following end-repair, the size of the DNA was selected via a bead-based method. DNA fragments were further processed by the addition of a single ‘A’ base, and ligation with TruSeq indexing adapters. The purified libraries were quantified using qPCR according to the qPCR Quantification Protocol Guide (KAPA Library Quantification Kit for Illumina sequencing platforms) and quantified using an Agilent Technologies 2200 TapeStation (Agilent Technologies, Santa Clara, CA USA). Then, paired-end (2 × 150 bp) sequencing was performed by Macrogen (Seoul, Korea) using the NovaSeq platform (Illumina). WGS was used to measure the telomere length of 38 samples (19 each for NFs and MPNSTs).

2.3. Whole-Transcriptome Sequencing (WTS)

The total RNA concentration was calculated by Quant-IT RiboGreen (Invitrogen, Waltham, MA, USA). To determine the DV200 (% of RNA fragments > 200 bp) value, samples were run on a TapeStation RNA ScreenTape platform (Agilent). Total RNA (100 ng) was subjected to sequencing library construction using a TruSeq RNA Access Library Prep Kit (Illumina, San Diego, CA, USA) according to the manufacturer’s protocols. Briefly, total RNA was first fragmented into small pieces using divalent cations at an elevated temperature. The cleaved RNA fragments were copied into first-strand cDNA using SuperScript II reverse transcriptase (Invitrogen, #18064014) and random primers, followed by second-strand cDNA synthesis using DNA Polymerase I, RNase H, and dUTP. These cDNA fragments then underwent an end-repair process, involving the addition of a single ‘A’ base, and the ligation of adapters. These products were then purified and enriched via PCR to create a cDNA library. All libraries were normalized, and six libraries were pooled into a single hybridization/capture reaction. The pooled libraries were incubated with a cocktail of biotinylated oligos corresponding to coding regions of the genome. The targeted library molecules were captured via hybridized biotinylated oligo probes using streptavidin-conjugated beads. After two rounds of hybridization/capture reactions, the enriched library molecules were subjected to a second round of PCR amplification. The captured libraries were quantified using a KAPA Library Quantification Kit for the Illumina sequencing platforms according to the qPCR Quantification Protocol Guide (KAPA BIOSYSTEMS, Wilmington, MA, USA, #KK4854), and qualified using the TapeStation D1000 ScreenTape platform (Agilent Technologies, # 5067–5582). The indexed libraries were then submitted to sequencing on the Illumina NovaSeq platform (Illumina) and paired-end (2 × 100 bp) sequencing was performed by Macrogen. WTS was used to profile the transcriptomes of 40 samples (20 each for NFs and MPNSTs). All TPM values (NFs and MPNSTs) of the 131 genes evaluated in this study are shown in Supplementary Material S2.

2.4. DEG Selection

The raw RNA-seq data (fastq) of all the samples were mapped using STAR (v2.7.3) and quantified by individual genes with a normalized count value, transcripts per million (TPM), using RSEM (v.1.2.31). Noncoding genes, including pseudogenes and long noncoding RNA, were excluded from the analysis and only coding genes were used (20,315 coding genes). The expression value was log2 transformed (log_2_(TPM + 1)), where the entire expression profile was proportionally transformed, and the 0 conversions were corrected by adding 1 to the substitution. By comparing 20 NF–MPNSTs, the DEGs between the NFs and MPNSTs were selected on the basis of an average fold-change greater than two and a t-test p-value < 0.001. To achieve a fairly high confidence level in DEG selection, we adopted a p-value < 0.001 to select the most reliable and appropriate number of genes. Using the above filtering criteria, 22 significant DEGs were selected. Statistical testing for each gene was achieved using paired t-tests and a p-value < 0.05 was considered statistically significant.

2.5. Measuring the Length of Telomeres

Telomeric reads were calculated using TelomereHunter [25]. The program extracts telomeric reads with nonconsecutive repeats, and the extracted telomeric reads are categorized according to their alignment coordinates and mapping quality. The reads are classified into different telomeric regions. The telomere content is given as the fraction of intratelomeric reads per million reads. The program also accounts for GC bias and GC is used to correct the telomeric content by dividing the intratelomeric reads by the GC number of reads, with a GC content between 48% and 52%. In this study, GC-corrected telomeric reads were used.

2.6. Statistical Analysis

Statistical analyses were performed using SPSS (SPSS, Inc., Chicago, IL, USA). Transcriptional relationships were analyzed with LR. For each gene pair, transcriptional activation of the B gene by the A gene was defined as follows: the ΔA transcripts during NF–MPNST transformation were correlated with both the ΔB transcripts and the MPNST-B transcripts, where the significance of the LR results between the ΔA transcripts and the ΔB transcripts was confirmed using the Benjamini–Hochberg procedure. The interaction effects of genes A and B on the telomere length were also analyzed with a GLM (Figure 1C). Correlations with the histologic grade were evaluated using logistic regression, and MFS and OS were evaluated using Cox regression.

3. Results

3.1. Differentially Expressed Genes During MPNST Malignant Transformation

In this study, 20 NF–MPNST pairs from 20 NF1 patients were analyzed. The patient demographics and tumor characteristics are described in Table 1. The mean latency to MPNST presentation was 36.55 years, and 30.0% of the patients were in the metastatic stage at MPNST presentation. The OS rate was 35.0%, and the metastasis rate was 71.4%. In terms of the tumors, 50% of the MPNSTs were high-grade, with an FNCLCC grade of 3. There were no significant differences in tumor sizes (p = 0.430). There was a significant difference in the telomere length between the NFs (999.96 ± 423.69 intratel_reads × 1,000,000/total_reads_with_tel_gc) and MPNSTs (727.32 ± 490.75 intratel_reads × 1,000,000/total_reads_with_tel_gc) in the paired t-tests (p = 0.043) (Table 1). Although there were statistical differences between the NFs and MPNSTs telomere lengths, the differences did not show a clear distinction, suggesting that MPNTs achieve immortalization through adequate TMM activation.

Twenty-two differentially expressed genes (DEGs) were identified from the WTS during the MPNST malignant transformation by selecting the genes with an average fold-change of greater than twofold. The Z-score heatmap (Figure 2A) and volcano plot (Figure 2B) are illustrated. All the NF and MPNST groups do not have the same profile, and the level of gene expression appears to be different in the heatmap. This is caused by the fact that the cellular heterogeneity of tissues may vary slightly from sample to sample, which may lead to slight differences in the gene expression profiles. Eleven upregulated genes were mainly classified as related to (1) cell signaling—NELL2 and DLGAP5; (2) cytokinesis—ASPM and BUB1; (3) forkhead-box transcription factors—FOXG1 and FOXM1; and (4) transcription factors:—SIX1 and SOX11, whereas eleven downregulated genes were mainly classified as related to (1) cell adhesion—CHL1, CDH19, CLDN1, SORBS1, and DMD; and (2) synapse-related—PRIMA1, GRIK3, and SLITKR2. The details of the 22 DEGs are schematically shown in Figure 2C.

Among the upregulated DEGs, BUB1 promotes telomere replication during the S phase in HeLa cells [26]. Although the NELL2-Robo3 complex is involved in the activation of axonal guidance [27], NELL2 is also involved in various cancers, for instance, the NELL2/cdc42-BAF complex in Ewing’s sarcoma cells [28], fat mass and obesity-associated protein (FTO)/E2F1/NELL2 in non-small-cell lung cancer (NSCLC) [29], and NELL2/N-cadherin in embryonic carcinoma cells [30]. FOXM1 is highly expressed and a marker of poor prognosis in several cancers [31], including MPNST [32]. SOX11 is also highly expressed in nervous system neoplasms [33]. Among the downregulated DEGs, PRIMA1 has been reported to be a tumor suppressor by restoring mutant p53 [34] and inducing tumor cell death [35]. DMD is also a tumor suppressor in sarcomas, hematologic malignancies, and nervous system tumors [36]. SORB1 suppression promotes lung adenocarcinoma [37].

3.2. NELL2 Activates PAX7 During Oncogenesis of MPNST

The paired box (PAX) genes encode transcription factors with highly conserved N-terminal DNA-binding domains, known as paired domains, which have been reported to be important in neural development and oncogenesis. In addition to their orthodox functions in normal neural development, several studies have described their unorthodox oncogenic functions in various cancers (Figure 3A). Mutations in the PAX genes are associated with congenital human diseases related to eye development and deafness, including Aniridia and Peter’s anomaly (PAX6) and Waardenburg syndrome (PAX3), suggesting that the PAX genes play a central role in the development of the nervous system [38]. Moreover, the PAX genes have also been reported to be involved in oncogenesis [39]. The nine PAX genes are subclassified into four groups according to their homeodomain and octapeptide domain: Group I (PAX1 and PAX9), Group II (PAX2, PAX5, and PAX8), Group III (PAX3 and PAX7), and Group IV (PAX4 and PAX6). Group III contributes to sarcomas (Ewing’s sarcoma and rhabdomyosarcoma) and neural crest-derived tumors [40,41] (Figure 3B).

We analyzed both the transcriptional relationship and the interaction effects on the telomere length for 198 DEG-PAX gene pairs, comprising 22 DEGs and 9 PAX genes. The transcriptional activation of a PAX gene by a DEG was defined as the significant positive correlation of the ΔDEG transcript with both the ΔPAX transcript and the MPNST*-PAX* transcript according to a linear regression (LR). In addition, to screen for only TMM-related transcriptional activation, the interaction effects of each gene pair on the telomere length (MPNST–telomere) were also analyzed using a generalized linear model (GLM). All 594 p-values from the LRs for ΔPAX and MPNST-PAX and the GLM were verified to control the false discovery rate (FDR) using the Benjamini-Hochberg (B-H) procedure, with a false discovery rate (FDR) of 0.1. Conclusively, the transcriptional activation of the TMM activation pathway was defined as a significant result for both the LR and GLM, following verification by the B-H procedure (Figure 1C). NELL2 transcriptionally activates PAX7 and SORBS1 transcriptionally activates PAX5 (Supplementary Material S3). The transcriptional activation of PAX7 by the 22 DEGs is summarized in Table 2.

3.3. PAX7 Activates RAD52-Dependent ALT During Oncogenesis of MPNST

An in-depth search of the literature revealed 100 genes that play important roles in the TMM. These genes were classified into four main categories based on their function in the TMM: (1) telomere machinery, (2) telomere topology, (3) regulatory genes involved in the DNA damage response and cell cycle/checkpoint, and (4) TMM-effector genes. The abbreviations and details for the 100 genes are listed in Supplementary Material S1. To evaluate the downstream regions of PAX7 and PAX5, we analyzed both the transcriptional relationship and interaction effects on the telomere length of 200 pairs comprising these 2 genes and 100 TMM-related genes. In addition, we evaluated the process of TMM activation caused by DNA damage signals. The DNA damage markers H2AFX, PARP1, P53, and RB1 were also analyzed in the same way as above for 396 pairs comprising these 4 genes and the remaining 96 TMM-related genes. Finally, all 1788 p-values from the two LRs and GLM for the 596 pairs were verified to control the false discovery rate (FDR) by the Benjamini-Hochberg (B-H) procedure (FDR = 0.1). PAX7 activates RAD52 and SLX4IP [24], whereas PAX5 has no effect on TMM activation. In the DNA damage signaling pathway, H2AFX activates RAD51 and P53 activates FEN1, but those results did not pass verification by the B-H procedure. (Supplementary Material S4).

3.4. RAD52-Dependent ALT Elongates Telomeres During Oncogenesis of MPNST, Which Results in Poor Survival

We evaluated the biological and clinical manifestations of RAD52-dependent ALT induced by the NELL2-PAX7 transcriptional cascade. We first evaluated their effects on telomere length during the malignant transformation. The correlations between ΔNELL2, ΔPAX7, ΔRAD52, ΔH2AFX, and ΔRAD51 and the telomere length are below. The Δtelomere length is the telomere length change during a malignant transformation from a benign (NF) to malignant tumor (MPNST) and the MPNST–telomere length represents the telomere length in an MPNST. They were evaluated using an LR analysis. RAD52 was the only gene that correlated with the MPNST–telomere length in the simple analysis (B = 114.894, p = 0.017), and RAD52 was independent of NELL2 and PAX7 in the multiple analysis (B = 161.989, 0.021). Although RAD51-dependent ALT is also a major component of TMM, RAD51 was not associated with the telomere length (B = −23.229, p = 0.633 for Δtelomere length and B = −34.669, p = 0.603 for MPNST–telomere length). This finding strongly suggests that NELL2-PAX7-induced RAD52-dependent ALT plays a dominant role in telomere maintenance and overcomes replicative senescence to achieve the malignant transformation of MPNST (Table 3).

To evaluate the association between telomere length and prognosis, we divided all the factors into two groups according to their changes during the malignant transformation: increased (Δtranscript > 0) and not increased (Δtranscript ≤ 0) groups. The overall survival (OS) and metastasis-free survival (MFS) were compared between the two groups using a Cox regression. The group with increased telomere length showed an inferior OS to the group without increased telomere length (Figure 4A right panel, HR = 3.809 (1.076 to 13.479), p = 0.038), suggesting that achieving telomere maintenance leads to a poor prognosis. Additionally, the group with an increased PAX7 transcript showed an inferior MFS (Figure 4C left panel, HR = 9129, p = 0.007) and OS (Figure 4C right panel, HR = 4.896, p = 0.046) in comparison with the group without an increase in the PAX7 transcript. The group with an increased RAD52 transcript showed an inferior MFS to the group without an increased RAD52 transcript (Figure 4A left panel, HR = 8.669 (1.654 to 45.432), p = 0.011). However, NELL2 was not associated with survival (Figure 4).

4. Discussion

We conducted this study with smart strategies. First, we designed a new analysis model optimized for the study of malignant transformation: we compared NFs and MPNSTs in a single NF1 patient. Since MPNSTs develop from preexisting NFs, they must have the same genetic background, which minimizes the genetic bias and the comparisons of the NF–MPNST pairs indicate only the changes that develop during the malignant transformation. Second, we designed a combined statistical analysis method: we analyzed the interaction effects on the telomere length and the transcriptional relationship to distinguish only the TMM-related transcriptional relationships and filter out those that are unrelated. Using these strategies, we suggested that the NELL2-PAX7 cascade activates RAD52-dependent ALT during the malignant transformation of MPNSTs. This answers an outstanding question in telomere biology: what signal activates RAD52-dependent ALT? In addition, we also determined the biological and clinical manifestations of RAD52-dependent ALT. RAD52-dependent ALT activated by the NELL2-PAX7 transcriptional cascade increases telomere length, leading to a poor prognosis for the MPNSTs.

Transcriptional activation is the most universal and fundamental process in gene activation and is involved in numerous signaling pathways and developmental processes. However, although studies of the TERT promoter [42] and transcriptional regulation of the TMM have been attempted [43,44], the transcriptional activation of the TMM is still unknown. No notable TMM activation signals have been identified, except for those related to DNA damage. Under these challenging conditions, we have suggested a new NELL2-PAX7 transcriptional cascade that activates RAD52-dependent ALT. This discovery could provide a new systematic understanding of the overall architecture of the TMM activation process, providing insights into the fundamental pathways for TMM activation via the transcriptional cascades. Although PAX7 is known to be essential for neural development (orthodox functions), it has been reported to play unorthodox functions in oncogenesis. At the molecular level, there are similarities between its orthodox and unorthodox functions. A major function of PAX7 during normal development is to maintain cell proliferation with self-renewal, which is comparable to overcoming the replicative senescence during a malignant transformation. Additionally, resistance to apoptosis and cell migration are required for both normal development and oncogenesis (metastasis). Therefore, the unorthodox functions of PAX7 may be considered a reinvocation of their orthodox functions during oncogenesis, albeit with different biological consequences.

Limitations of This Study

Despite these results, there are some limitations that should be considered. First, functional validation through a molecular biology experiment was not performed. Since this study was based on the changes in mRNA transcripts rather than protein–protein interactions, we were unable to analyze the genes for epigenetic changes, such as phosphorylation and SUMOylation (small ubiquitin-related modifier). In this study, H2AFX showed no association with ATM or ATR, which encode kinases and convert H2AX to its active form, γH2AX, which is a DNA damage marker. ATRX/DAXX, which belongs to the SWI/SNF family, suppresses the ALT by depositing H3F3A onto telomeres and resolving G-quadruplex (G4) loops [45]. ATRX encodes an ATP-dependent helicase that undergoes cell cycle-dependent phosphorylation and DAXX encodes a potent transcriptional repressor that binds to SUMOylated transcription factors. These genes showed no association with the NELL2-PAX7 cascade in this study. In addition, the SMC5/SMC6 complex and NSMCE2, which are representative examples of SUMOylation regulation in the TMM [46], also showed no significant interactions in this study. Therefore, in this study, the significance of DNA damage signaling (MDC1-ATM/ATR-CHEK1 pathway), the telomere topology (ATRX/DAXX), and telomere recruitment to the ALT-associated promyelocytic leukemia nuclear body (APB) (SMC complex and NSMCE2) may have been partly obscured and underestimated due to methodological reasons. The findings in this study need to be further validated at the protein–protein interaction level by wet-laboratory research.

Second, this study was conducted entirely through statistical analyses and consisted of a relatively small number of cases, which may be a clear weakness. However, NF1-associated MPNST is very rare, with an incidence of approximately 0.0005% in the general population. Therefore, it is difficult to secure a large cohort. To maximize the reliability of the statistical results of this study, we focused on false-positive control. To minimize false positives in the statistical results, extremely high p-value controls (FDR = 0.1) were applied to all the statistical results in this manuscript.