Thirteen TGCT tissue samples (11 from SEM and 2 from non-SEM patients with predominantly malignant embryonal carcinoma) and 4 para-carcinoma tissues were obtained from the Department of Urology, the Affiliated Cancer Hospital of Xiangya School of Medicine of Central South University ( Table 1 ). Fresh tissues were collected and snap frozen in liquid nitrogen. The study was approved by the Ethics Committee, and informed consent was obtained from each patient.

Total RNA was extracted from tissue samples using Trizol reagent (Thermo, USA), and 1µL per sample was used to determine the concentration using a NanoDrop 1000 spectrophotometer. After removing rRNAs (Ribo-Zero Gold rRNA Removal Kit, MRZG12324, Illumina), the remaining mRNAs and ncRNAs were purified using Agencourt RNA Clean XP Beads. The enriched mRNA and ncRNA fragments were cut into shorter sequences using a fragmentation buffer and reverse-transcribed using the Transcriptor First-strand cDNA synthesis Kit (Roche, USA). The cDNA fragments were purified to repair the ends add poly (A) and ligate the Illumina sequencing joints (QiaQuick PCR extraction kit). The second-strand cDNA was digested with uracil-N-glucanase, and amplified by PCR. The Sensitivity DNA assay Kit (Agilent Technologies) was used to analyze the quality of the library. Finally, ABI StepOnePlus Real-Time PCR System (Life Technologies) was used to quantify and pool the sequences, and the DNA sequenced on the Illumina HiSeqTM 4000 platform (Gene Denovo, Guangzhou, China).

The differentially expressed mRNAs between the two groups were screened using the edgeR software with FDR (false discovery rate) < 0.01 and log2FC (fold change) >1 as the thresholds for significant differences. The DEmRNAs were functionally annotated by gene ontology (GO) using the R software ClusterProfile package, and KEGG pathways using the Metascape software (http://metascape.org/gp/index.html#/main/step1) (15).

The DElncRNAs and DEmRNAs between the non-tumor and SEM, non-tumor and non-SEM, and non-tumor and non-SEM+SEM groups were identified using the edgeR software as described above. The common DElncRNAs and DEmRNAs among the three comparison groups were screened, and those with Spearman correlation coefficient r2>0.9999 and p<0.001 were selected for network construction. The target miRNAs of the lncRNA-mRNA interaction pairs were predicted using the R software miRNAtap package, TargetScan (16) and Miranda (17). The miRNAs obtained from all three programs were used to construct the lncRNA-miRNA-mRNA network, which was visualized using Cytoscape_v3.8.0.

The DEmRNAs were validated using the GEPIA2 database (http://gepia2.cancer-pku.cn/#index) (18). The RFS of the different groups (total cases 134) were plotted using the Kaplan-Meier Plotter (http://kmplot.com/analysis/index.php?p=background), and compared in terms of P values and hazard ratios (HR) with 95% confidence intervals (CI) by the Log-rank test and univariate Cox proportional hazard regression (19).

RNA-sequencing data (level 3) of 134 tumors and corresponding clinical information were obtained from The Cancer Genome Atlas (TCGA) dataset (https://www.cancer.gov/). The predictive accuracy and risk score of genes were compared by TimeRoc analysis. The minimum absolute shrinkage and selection operator (Lasso) regression algorithm was used for feature selection with 10-fold cross-validation using the R software package GLMnet. All statistical analyses were performed using R software (R Foundation for Statistical Computing, 2020) version V4.0.3. P <0.05 was considered statistically significant.

The genomic gene set analysis in TGCT was used the Gene Set Cancer Analysis (GSCA) (http://bioinfo.life.hust.edu.cn/GSCA/#/), which is an integrated database for genomic and immunogenomic gene set cancer analysis (20). Immunogenomic analysis was performed by ImmuCellAI (Immune Cell Abundance Identifier) algorithm with 24 immunes cells (http://bioinfo.life.hust.edu.cn/ImmuCellAI#!/) (21). P <0.05 was considered statistically significant.

A total of 11,469 mRNAs were differentially expressed between the non-tumor and non-SEM samples, of which 3,894 were up-regulated and 7,575 were down-regulated in the latter ( Figures 1A, D ). Furthermore, there were 21,506 DEmRNAs in SEM relative to the non-tumor samples, of which 5974 were up-regulated and 15,632 were down-regulated ( Figures 1B, D ). Finally, the SEM and non-SEM tumors exhibited 20,197 DEmRNAs, including 5,468 up-regulated and 14,729 down-regulated mRNAs compared to the non-tumor samples ( Figures 1C, D ). The expression levels of the randomly selected ADAM19, DNMT3L, BRCC3, ZMIZ1, BRAP, DCAF5, FBF1 and OPTN were verified using the GEPIA2 database ( Figure 1E ).

To determine the contribution of DEmRNAs to TGCT, we analyzed the top 100 mRNAs with high expression levels in the non-tumor (vs SEM + non-SEM), SEM (vs non-tumor + non-SEM), and non-SEM (vs non-tumor + SEM) samples ( Figure 2A ). The highly expressed genes in the non-tumor samples were enriched in the GO terms such as “structural constituent of ribosome”, “structural molecule activity” and “nucleoside-triphosphatase activity” ( Figure 2B ), and the “Tight junction”, “Gap junction” and “Apoptosis” pathways ( Figure 2E ). In the non-SEM samples, the significantly enriched GO terms were “molecular function”, “MHC protein complex binding” and “translation factor activity, RNA binding” ( Figure 2C ), whereas the main KEGG pathways were enriched in immune-related signaling pathways, such as “B cell receptor signaling pathway”, “Antigen processing and presentation” and “intestinal immune network for IgA production” ( Figure 2F ). The highly expressed genes in SEM showed significant enrichment in “MHC protein complex binding” “structural constituent of ribosome” and “molecular function” ( Figure 2D ), and the immune-related signaling pathways, such as “Antigen processing and presentation”, “intestinal immune network for IgA production” and “Th17 cell differentiation” ( Figure 2G ). Furthermore, the top 100 highly expressed mRNAs in all groups were associated with multiple transcription factors, including SMARCA5 and SMARCC1, which were overexpressed in all three groups ( Figure 2H ).

As shown in Figure 3A , there were many down-regulated mRNAs in SEM and non-SEM samples compared to the controls ( Figure 3A ). Specifically, 3047 mRNAs were up-regulated and 7,883 mRNAs were downregulated in the SEM and total TGCT (SEM + non-SEM) samples relative to the non-tumor controls ( Figures 3B, C ). The non-SEM and combined TGCT samples showed upregulation of 196 mRNAs and downregulation of 88 mRNAs compared to the non-tumor group. Four upregulated and 18 downregulated mRNAs were common to both SEM and non-SEM groups compared to the non-tumor control ( Figures 3B, C ). A pairwise comparison of all groups revealed that 2,091 upregulated mRNAs and 6,270 downregulated mRNAs were common. In addition, 1,603 and 1,191 mRNAs were respectively upregulated and downregulated only in the non-SEM group relative to the non-tumor group, whereas 732 upregulated and 1,461 downregulated mRNAs were specific to the non-tumor-vs-SEM comparison, and 134 upregulated and 488 downregulated mRNAs were only detected when comparing non-tumor samples to all TGCT samples ( Figures 3B, C ). GO analysis revealed that the DEmRNAs common to all three groups were enriched in “cell motility”, “movement of cell or subcellular component” and “cell cycle” ( Figure 3D ), and the top 20 pathways were “cell adhesion molecules”, “adherence junction” and “Tight junction” ( Figure 3E ). Furthermore, the signaling pathways significantly associated with the top 100 common DEmRNAs included “apoptosis signaling pathway”, “positive regulation of cellular component movement”, “leukocyte activation involved in immune” and “cytokine signaling in immune system”, and were interconnected ( Figure 3F ). The protein-protein interaction (PPI) network of the top 100 common DEmRNAs revealed 3,098 connections ( Figure 3G ).

A total of 2,924 DElncRNAs and 4,842 DEmRNAs were common to the non-tumor vs SEM, non-tumor vs non-SEM+SEM, and non-tumor vs non-SEM comparisons ( Figures 4A, B ). The correlated mRNA-lncRNA pairs were screened with R2 > 0.9999 and p < 0.001 as the criteria. Nine lncRNAs were identified, including SPATA42, AL450226.1, TCONS_00053325, DPYD-AS2, ELW67898.1, DIAPH1-AS1, Z97205.2, ISPD-AS1 and TCONS_00022110. The nine mRNAs were FAM71E2, HEPACAM2, RGSL1, DPYD-AS2, ZPBP2, DYNLRB2, PCYT2, BANF2, and GRK4. The regulatory miRNAs of the candidate mRNAs and lncRNAs were predicted by the R software miRNAtap package, TargetScan (16), and Miranda (17). We identified 8 lncRNA–miRNA–mRNA regulatory networks that included 9 lncRNAs nodes, 9 mRNAs nodes, and 82 miRNAs nodes with 166 possible interactions ( Figure 4C ). The DEmRNAs common to non-tumor vs non-SEM and non-tumor vs SEM pairs were associated with “movement of cell or subcellular component”, “reproduction” and “cell motility” both in the groups of ( Figures 4D, E ), whereas “male gamete generation”, “cell differentiation” and “reproduction progress” were the functions enriched in the DEmRNAs between non-tumor vs non-SEM + SEM ( Figure 4F ). The top 20 enriched pathways among the DEmRNAs between non-tumor vs non-SEM groups were “focal adhesion” and “cell adhesion molecules” ( Figure 4G ), that for non-tumor vs SEM were “biosynthesis of secondary metabolism”, “steroid hormone biosynthesis” and “cell adhesion molecules” ( Figure 4H ), and “cell adhesion molecules” and “metabolic pathway” were enriched for non-tumor vs non-SEM + SEM ( Figure 4I ).

To explore the clinical value of the nine mRNAs (see above) in TGCTs, we obtained the raw RNA-Seq data (Level 3) of 134 TGCTs and the corresponding clinical information from TCGA dataset (https://www.cancer.gov/) ( Supplementary Table 1 ). A polygenic risk score (PRS) was then calculated based on the correlation between the expression levels of the nine mRNAs and prognosis using the LASSO cox regression analysis ( Figures 5A, B ). Accordingly, the patients were divided into the high-risk and low-risk groups ( Figure 5C , top). Overexpression of GRK4, PCYT2 and RGSL1 was associated with increased mortality rates ( Figure 5C , middle) and a higher risk score ( Figure 5C , bottom). The prognostic index was calculated as follows:

KM survival analysis also confirmed that the high expression level of this gene set portended poor prognosis (log-rank p=0.0404, HR=2.077, Figure 5D ). Furthermore, the area under the time-dependent ROC curve of the gene set was greater than 0.5, indicating good predictive power ( Figure 5D ). To further establish the prognostic value of GRK4, PCYT2 and RGSL1, we analyzed the correlation between their expression levels with RFS using Kaplan-Meier Plotter. PCTY2 (P=0.049, Figure 5F ) expression level showed a significant negative correlation with RFS, while that of GRK4 (P=0.12, Figure 5E ) and RGSL1 (P=0.064) were negatively correlated with RFS without statistical significance ( Figures 5F, G ).

The potential therapeutic relevance of GRK4, PCYT2 and RGSL1 was evaluated by GSCA (http://bioinfo.life.hust.edu.cn/GSCA/#/) (20) and Spearman correlation analysis. As shown in Figure 6A , the expression levels of this gene set in TGCT showed a significant negative correlation with DNA damage, androgen receptor (AR), estrogen receptor (ER) and other pathways, and positive correlation with RASMAPK, TSCmTOR and other pathways ( Figure 6A ). Furthermore, the gene set was also negatively correlated with the infiltration of immune cells such as Cytotoxic, TFH, B cell, CD4_T and CD8_T, and positively with Monocyte, Th2 and Macrophage ( Figure 6B , Supplementary Tables 2, 3 ).

Finally, we analyzed the correlation between drug sensitivity of cancer cell lines in GDSC (Genomics of Drug Sensitivity in Cancer), and the GRK4, PCYT2 and RGSL1 expression levels (22). GRK4 was correlated with the target drug 17-AAG and increased sensitivity to RDEA119, PD-0325901, CX5461, Trametinib and Selumetinib. PCYT2 showed a significant positive correlation with sensitivity to targeted drugs such as CX-5461, PI-103 (PI3K and mTOR inhibitors) and MPS-1-IN-1 ( Figure 6C ).