All human studies performed in the current study were approved by the Medical Ethics Committee of Kunming Children’s Hospital before collection (approval number: 20190822002), and the guardians of the children signed informed consent forms.

We retrieved the transcriptome data of MRTK from the TARGET database (https://ocg.cancer.gov/programs/target).

A total of 200 genes in the gene set of HALLMARK_MTORC1_SIGNALING were retrieved from the Molecular Signatures Database (MSigDB).

Clinical tissue specimen: The tissue specimens of six children with MRTK cancer and adjacent tissues in our hospital were selected. The relevant studies have been approved by the relevant ethics committee, and informed consent was signed by their families.

A total of 50 MRTK samples and six normal samples were applied to conduct GSVA, of which h.all.v7.2.symbols.gmt was considered to be the preset pathway. Then, we used the “limma” package to analyze the differences of GSVA scores between the MRTK and control groups. The screening conditions were |t-value| > 2 and p-value < 0.05. In the GSVA, the screened pathways with t-value > 0 were considered to be activated in the MRTK group, while screened pathways with t-value < 0 were considered to be activated in the normal group.

The genes were extracted from the identified HALLMARK_MTORC1_SIGNALING pathway. Thereafter, we used the “limma” R package to perform the differential analysis for identifying the differentially expressed genes (DEGs) between the MRTK samples and normal samples. The screening conditions of DEGs were |log₂FC| > 1 and p-value < 0.05. We combined the DEGs with 200 genes of the HALLMARK_MTORC1_SIGNALING gene set, obtaining the differentially expressed mTORC1 pathway-related genes (DE-mTRGs).

A total of 49 MRTK samples containing complete clinical information were applied to construct an mTRG-related gene signature. The univariate Cox regression analysis and Kaplan–Meier (K-M) survival curves were performed by employing the survival R package to identify the DE-mTRGs related to the overall survival (OS) of MRTK (p < 0.05). DE-mTRGs with p < 0.05 were input into a multivariate Cox regression analysis for the construction of an mTRG-related gene signature. The risk score of each MRTK sample is obtained as the following formula: Risk score = h0(t)*exp (β1X1+β2X2+…+βnXn). In this formula, β is the regression coefficient, and h0(t) is the benchmark risk function. A total of 49 MRTK samples were classified into high-risk and low-risk groups with the boundary of the median risk score. The K-M survival curve was generated via the “survdiff” functions of the R package to compare the OS of two risk groups. Thereafter, a time-independent receiver operating characteristic (ROC) curve was conducted to demonstrate the effectiveness of this mTRG-related gene signature for OS in MRTK patients.

The relationship between risk score and clinicopathological features was evaluated using the Chi-square test. Then, we executed univariate and multivariate Cox regression analyses to evaluate whether these clinicopathological features were independent predictors for RKT prognosis. A nomogram containing independent prognostic factors was generated using the “rms” R package, and the corresponding calibration plot was further established to evaluate the efficiency of the nomogram.

Gene set enrichment analysis (GSEA) of prognostic mTRGs was analyzed by the “clusterProfiler” R package (version 3.18.0) to explore the potential signaling pathways that were involved. The “limma” package was executed to screen the DEGs in high- and low-risk groups with the screening conditions being |log₂FC| > 1 and p-value < 0.05. In addition, functional enrichment analyses including Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses were further conducted using Metascape software based on the DEGs.

The tumor stem cell index of MRTK samples was calculated using the OCLR algorithm. The correlations of tumor stem cell index and the prognostic factors were analyzed by the Spearman correlation analysis. Additionally, we searched the target drugs of these prognostic factors from The Comparative Toxicogenomics Database (CTD). The interactions of target drugs and prognostic genes were visualized by Cytoscape software.

This study has been approved by the Ethics Committee of Kunming Children’s Hospital, and all patients and their parents signed informed consent before joining the study. A total of six cases of fresh tumor and paratumor tissues were collected from children with MRTK admitted to Kunming Children’s Hospital from July 2014 to June 2020. There were 2 male and 4 female cases aging from 4 months to 4 years and were not treated with radiotherapy before surgery. The tumor tissues were excised as soon as possible after isolation in 0.5 cm³ volume, and relatively normal kidney tissues were taken from the distal cut edge of pathologically confirmed paratumor-infiltrated tissues. The samples were snap-frozen with liquid nitrogen and transferred to −80 °C for storage prior to use.

Total RNA was extracted from MRTK tumor and paratumoral tissues using a TRIzol (T9424, Sigma, Germany) reagent following the manufacturer’s instructions. The RNA concentration and purity were determined using a NanoDrop 2000 spectrophotometer (Thermo Fisher Scientific). Then, mRNA was reversely transcribed into cDNA using a PrimeScript RT reagent kit with gDNA Eraser (RR047A, Takara, Japan). qRT-PCR was performed using a GoTaq qPCR and RT-qPCR kit (A6001, Promega, United States) using an ABI ViiATM7 real-time PCR system (Applied Biosystems, United States). The PCR program was started with an initial step of 95°C for 3 min, followed by 40 cycles of 15 s denaturation at 95°C and 30 s extension at 60°C. The mRNA expression levels were calculated by the 2^–ΔΔt method [△△ Ct = (Ct target gene-Ct internal reference gene) experimental group-(Ct target gene-Ct internal reference gene) control group] and normalized to GAPDH. The experiments were performed in triplicate for at least three samples of each group. The primers were synthesized by Sangon Biotech (Shanghai) and are listed in Table 1.

The protein was extracted from tumor and paratumor tissues using RIPA buffer (Beyotime, Shanghai, China). BCA protein quantification kit (Shanshan Jinqiao Company, Beijing, China) was used for the quantification of protein samples. About 30 μg protein of each sample was loaded on SDS-PAGE gels. After transferring to the polyvinylidene fluoride membrane (Millipore, United States), the blots were probed with the appropriate primary antibody (anti-P4HA1, 1:1000, ProteinTech; anti-MLLT11, 1:1000, ProteinTech; anti-AURKA, 1:1000, Wanleibio; anti-GOT1, 1:1000, Wanleibio; anti-β-actin, 1:1000, Zen-Bioscience). After incubation with the primary antibodies, the membrane was incubated with the corresponding secondary antibodies and monitored with an enhanced chemiluminescence (ECL, Affinity, United States) substrate kit (Amersham, Biosciences, United States). The protein bands were photographed using an Image Lab system (Bio-Rad Laboratories, Inc., United States). After developing, the image was saved, and the Image Lab (Bio-Rad Laboratories, Inc., United States) analyzed the gray value.

All bioinformatics analyses involved in this research were conducted using the R software (version 3.6.3). The relationship between risk score and clinicopathological features was demonstrated using the Chi-square test. Moreover, the independent prognostic factors for MRTK identified the clinicopathological features by univariate and multivariate Cox regression analyses. GraphPad8.0 statistical software was used for data analysis of mRNA and protein expression. The t-test was used to compare the difference between the tumor and paratumor groups. p < 0.05 was considered statistically significant unless specified.

Through GSVA analysis, we screened 9 activated pathways in the MRTK group and 15 activated pathways in the normal group. Interestingly, we observed that the mTORC1 signaling pathway was activated in the MRTK samples (Figure 1A; Supplementary Table S1). Subsequently, we extracted the genes from the mTORC1 signaling pathway to screen the DEGs between the MRTK and normal samples. As a consequence, a total of 4787 DEGs were screened in the MRTK samples compared with normal samples, of which 1294 were upregulated and 3493 were downregulated (Figure 1B; Supplementary Table S2). After intersection with the 200 genes of the HALLMARK_MTORC1_SIGNALING gene set obtained from the MsigDB, we identified 70 DE-mTRGs in MRTK (Figure 1C). The expressions of 70 DE-mTRGs between the MRTK and normal samples were displayed in a heatmap. Among the 70 DE-mTRGs, 22 genes were downregulated and 48 genes were upregulated in patients with MRTK (Figure 1D).

To assess whether these 70 DE-mTRGs were related to survival, we first performed the K-M survival curves and univariate Cox regression analysis. As a result, four DE-mTRGs were demonstrated to be related to the OS of MRTK (Figure 2A, Supplementary Table S3, p < 0.05). A total of four DE-mTRGs including P4HA1, MLLT11, AURKA, and GOT1 were obtained from the multivariate Cox regression analysis, which were used to develop an mTRG-related gene signature (Figure 2B). Moreover, the survival analyses of each gene uncovered that patients with high expression of these genes were all related to a poor prognosis (all with p < 0.05, Supplementary Figure S1). The risk score of each MRTK patient was obtained as the formula mentioned in Materials and Methods. The MRTK patients were stratified into a high-risk group (n = 25) and a low-risk group (n = 24) by the boundary of the median risk score. The risk curves and patient survival scatter plots showed that the number of patient deaths increased with the risk score (Figure 2C). The K-M curves demonstrated that the risk score significantly differentiated clinical outcomes in MRTK patients (p < 0.0001), with a high-risk score implying a poor prognosis (Figure 2D). The area under curve (AUC) was 0.811 and 0.878 at 1 and 3 years, respectively, indicating the validity of the 4-gene-based prognostic signature (Figure 2E). Furthermore, the heatmap showed that all 4 model genes possessed relatively high expression levels in the high-risk group (Figure 2F).

Next, a Chi-square test was implemented to assess the distribution of clinical characteristics of MRTK patients in the high- and low-risk groups. The results showed that the age distribution of patients in the high- and low-risk groups was significantly different (p = 0.032), with 72.0% of patients in the high-risk group being ≥365 days old and 62.5% of patients in the low-risk group being less than 365 days old (Table 2). Furthermore, the Wilcoxon rank sum test indicated that the level of the risk score was significantly higher in MTRK patient’s ≥365 days than in MTRK patients’ <365 days (p < 0.01); however, the risk score was not statistically significant in the gender and stage subgroups (Figure 3A). Subsequently, we performed an independent prognostic analysis by Cox regression to assess whether the risk score could affect the clinical outcome of OS in MRTK patients independent of the clinicopathological characteristics (age, sex, and stage). The univariate Cox regression analysis showed that the risk score and stage were significantly associated with OS in patients with MRTK (p < 0.05; Table 3). Finally, multivariate Cox regression analysis pointed out that the risk score and stage were independent prognostic factors for MRTK patients (Table 4). Thereafter, we constructed a nomogram based on two independent prognostic factors that predicted patients’ 1-year and 3-year OS (Figure 3B), which had a C-index of 0.767. The calibration curve indicated that the predictive value of the nomogram model for OS in MRTK patients was similar to the actual observed value (Figure 3C), implying that the combined model of the two independent prognostic factors may be more clinically applicable.

We conducted GSEA to further explore the involved signaling pathways of the four prognostic biomarkers. The results indicated that AURKA was mainly involved in cell cycle, DNA replication, non-alcoholic fatty liver disease, oxidative phosphorylation, proteasome, ribosome, and ribosome biogenesis in eukaryotes, RNA transport, spliceosome, and thermogenesis (Figure 4A; Supplementary Table S4). Similarly, GOT1 was significantly associated with various signaling pathways such as ECM–receptor interaction, carbon metabolism, human papillomavirus infection, AGE–RAGE signaling pathway in diabetic complications, focal adhesion, insulin resistance, lysosome, MAPK signaling pathway, peroxisome, and tight junction (Figure 4B; Supplementary Table S5). Interestingly, the expression of MLLT11 was correlated with several diseases, including amyotrophic lateral sclerosis, coronavirus disease-COVID-19, Huntington’s disease, Parkinson’s disease, and prion disease (Figure 4C; Supplementary Table S6). In addition, P4H41 was also mainly involved in bile secretion, chemical carcinogenesis, ECM–receptor interaction, lysosome, proteasome, RNA degradation, and RNA polymerase (Figure 4D; Supplementary Table S7).

To elucidate the molecular mechanism that was involved in the risk score, we conducted GO and KEGG functional enrichment analyses on DEGs of the two risk groups (Figures 5A,B). The results showed that the DEGs mainly focused on several important biological processes including central nervous system neuron differentiation and developmental process involved in reproduction, forebrain development, forebrain neuron differentiation, positive regulation of cell development, reproductive structure development, and reproductive system development (Figure 5C). KEGG analysis indicated that these DEGs were mainly associated with basal cell carcinoma, pathways in cancer, proteoglycans in cancer, melanogenesis, PI3K–AKT signaling pathway, gastric cancer, signaling pathways, mTOR signaling pathway, regulating pluripotency of stem cells, hepatocellular carcinoma, human papillomavirus infection, Hippo signaling pathway, cGMP-PKG signaling pathway, breast cancer, Cushing syndrome, and Wnt signaling pathway (Figure 5D).

It has been reported that mTOR exerts important functions in cancer stem cells through the specific functions related to stemness (Matsubara et al., 2013). We calculated the tumor stem cell index (mRNAsi) of MRTK samples based on the OCLR algorithm (Zhang et al., 2020) (Figure 6A). The correlation analysis showed that the expressions of AURKA and GOT1 were positively correlated with the mRNAsi in MRTK (Figure 6B). Then, we searched the corresponding target drugs of AURKA and GOT1 from the CTD. A network containing 262 nodes and 349 edges showed the interactions between the target drugs and prognostic factors that were positively correlated with mRNAsi (Figure 6C). These target drugs predicted by AURKA and GOT1 may play an essential role in the therapy of MRTK.

RT-qPCR results showed that the abundance of P4HA1, MLLT11, AURKA, and GOT1 were significantly elevated in tumor tissues (Figure 7A). Consistent with the RT-qPCR results, Western blot also showed that the protein levels of P4HA1, MLLT11, AURKA, and GOT1 were markedly higher in the tumor tissues than those in the paratumor tissues (Figure 7B). These results were in accordance with the RNA-sequencing results, indicating the reliability of our analysis.