In this study, three individuals with grade 3 immature ovarian teratomas, as defined by the WHO (17), were enrolled ( Table 1 ). Additionally, a validation cohort consisting of 7 patients across grades 1, 2 and 3 was included ( Supplementary Table 4 ). This study was approved by the Institutional Review Board of West China Second University Hospital (IRB No. 2020112), and informed consent was obtained from patients (18).

Immature ovarian teratomas were dissected and sheared on ice after being washed with Hanks’ balanced salt solution (HBSS). The tissues were then digested for 1 hour at 37°C with collagenase I (2 mg/ml) (Gibco 1710-0017), collagenase IV (1 mg/ml) (Gibco 1710-4019), and 0.25% pancreatic enzymes (Gibco 25200-056). To remove cell debris and large clumps, the digested tissues were filtered through a 40-mm strainer. At 4°C, the suspensions were pelleted at 500g for 5 minutes. After lysis with 10 RBC Lysis Buffer (Thermo Fisher Scientific, 00-4300-54), the pellets were resuspended in HBSS with 0.04% bovine serum albumin (BSA) to determine cell viability and concentration (Counting Star, Aber Instruments Ltd.). The remaining cells were pelleted at 500g for 5 minutes at 4°C before being stored at 80°C. Finally, single-cell suspensions were diluted to about 5000 cells per milliliter. All tissues were obtained from West China Biobanks at Sichuan University’s West China Hospital’s Department of Clinical Research Management.

According to the directions provided by the manufacturer in the Chromium Single Cell 3’ Reagents Kits v2 User Guide, the Chromium Single Cell 3’ Library & Gel Bead Kit v2 (PN-120237), Chromium Single Cell 3’ Chip Kit v2 (PN-120236), and Chromium i7 Multiplex Kit (PN-120262) were utilized. Phosphate-buffered saline (PBS) + 0.04% BSA was used to wash the single-cell suspension twice. With the use of the TC20 Automated Cell Counter, cell quantity and concentration were verified. A 10x Genomics Chromium Controller machine was used to generate gel beads in emulsion (GEM) from the cells right away. A 10x Genomics Chromium Single Cell 3’ reagent kit (V2 chemistry) was used to prepare barcoded complementary DNAs (cDNAs), which were then recovered, purified, and amplified to provide enough for library construction. With the use of an Agilent Bioanalyzer 2100, library concentration and quality were evaluated. For PE150 sequencing, libraries were performed on Illumina’s NovaSeq 6000 platform.

Analysis on single cell transcriptome was performed with the pipeline as we described previously (19, 20). Briefly, the raw data was processed using Cell Ranger (v3.0) (21). Following Cell Ranger processing, a raw unique molecular identifier (UMI) count matrix was generated, which was then converted into a Seurat object by the R package Seurat (22). All Cell Ranger commands were executed with the default settings. Cells with UMI numbers below 500, gene numbers below 300 or greater than 6,000, or mitochondrial-derived UMI counts of more than 15% were considered low-quality and were removed. Furthermore, we used Scrublet (v0.2) to calculate doublet scores and predicted potential doublets for each cell, with expected doublet rate = 0.06 (23). The predicted doublets were then removed. After quality control, 6274, 8435, and 7444 cells were left for each patient for subsequent analysis.

The Seurat v4 data integration pipeline was used to integrate the data (22). To summarize, we first used total-counts normalization to normalize the filtered counts matrix before log-transforming the result. We identified highly variable genes in each sample and chose the top 3000 genes that appeared to be overdispersed across the most teratomas. We then used sctransform (24) to normalize, scale, and correct for mitochondrial and ribosomal read percentages. Nonbiological batch effects were removed by performing a canonical correction analysis on individual samples using the Seurat functions FindIntegrationAnchors and IntegrateData.

Following data integration, the RunPCA function was used to perform PCA, with the first 40 principal components being used. The RunUMAP function was used to generate UMAP. The Seurat functions FindNeighbors and FindClusters were used to find clusters with a resolution of 2.0. We then manually examined the marker genes for each cluster and annotated the cell type based on canonical marker expression ( Supplementary Table 2 ).

Signature genes in human fetal Cell Atlas (25) and CTD teratomas (15) were defined as the top 50 DEGs with the lowest p values for each cell type. The signature scores for each cell type were then calculated using Seurat’s AddModuleScore function. The signature scores that resulted were interpreted as correlation assessments.

In this study, we compared the heterogeneity of CDT and PDT cell types using ROGUE (26), an entropy-based universal metric for assessing the purity of single cell populations. The ROGUE index has been scaled from zero to one. One represents a completely pure subtype with no significant genes, and zero represents the population’s most heterogeneous state.

For PCA reduction projection, expression matrix from PDT and CDT was extracted from SCT assay in Seurat objects. To reduce batch effects, we only select shared genes across PDT and CDT, and their expression matrix were normalized as Z-scores. PCA embedding were calculated on combined expression matrix by R prcomp function. In CDTs, the PCA similarity distance was calculated as the Euclidean distance between each CDT sample to the lower left CDT sample using the Pythagorean theorem.

Monocle 2 was used to perform pseudotime analysis to determine the dramatic developmental trajectory and translational relationships of neural cell types (27). 1500 significantly changed features were identified by the FindVariableFeatures function in Seurat with the vst method and used for cell ordering. The default parameters were used to calculate CytoTRACE scores, and the results were used to infer the origin of neural cells (28).

BEAM analysis (q-value < 1e-50) was used to identify cell differentiation-related genes, which were then visualized using the plot genes branched heatmap function.

SCENIC analysis was carried out as previously described (29). The pySCENIC package (version 0.10.3) was used, which is a lightning-fast Python implementation of the SCENIC pipeline. The Wilcoxon rank sum test was used to identify the differentially activated TFs of each subcluster in comparison to all other cells of the same cell type in PDTs.

The M1/M2 phenotype of each macrophage cell was defined as the mean expression of gene signatures (30). M1 macrophage signature includes IL23, TNF, CXCL9, CXCL10, CXCL11, CD86, IL1A, IL1B, IL6, CCL5, IRF5, IRF1, CD40, IDO1, KYNU, and CCR7. the M2 macrophage signature includes IL4R, CCL4, CCL13, CCL20, CCL17, CCL18, CCL22, CCL24, LYVE1, VEGFA, VEGFB, VEGFC, VEGFD, EGF, CTSA, CTSB, CTSC, CTSD, TGFB1, TGFB2, TGFB3, MMP14, MMP19, MMP9, CLEC7A, WNT7B, FASL, TNFSF12, TNFSF8, CD276, VTCN1, MSR1, FN1 and IRF4.

Phagocytosis genes contain MRC1, CD163, MERTK, and C1QB (30). Antitumor cytokines include TNF, IFNB1, IFNA2, CCL3, TNFSF10, and IL2 (31). Protumor cytokines contain IL4, IL5, IL13, IL10, GATA3 and CCR4 (31). Angiogenesis genes include VEGFA, VEGFB, VEGFC, PDGFC, CXCL8, CXCR2, FLT1, PGF, CXCL5, KDR, ANGPT1, ANGPT2, TEK, VWF, and CDH5 (30). The Effector T cells traffic phenotype was defined as the mean expression of CXCL9, CXCL10, CXCL11, CX3CL1, CCL3, CCL4, CX3CR1, CCL5, and CXCR3 (31). T cell cytotoxic signature contains GZMA, GZMB, GZMK, IFNG, NKG7, PRF1, CST7, CCL4 and TYROBP (32, 33). T exhaustion signature contains PDCD1, LAG3, TIGIT, HAVCR2 and CTLA4. The corresponding mouse signature were defined by their human homologous genes by function convert_mouse_to_human_symbols in R package nichenetr.

The quantification of ligand-receptor pairs among different cell types was used to assess cell-cell communication. CellphoneDB or CellChat (34, 35) software was fed gene expression matrices and metadata with cell type annotations. CellPhoneDB’s default database and parameters were used. PDT02 and PDT01/03 were analyzed separately for CellChat, and the netVisual_diffInteraction function was used to visualize the differential interaction strength between the two groups by using cycling neural progenitors as target cells.

The velocyto CLI (36, 37) was used to count the spliced and unspliced UMIs for each gene in each cell, and scvelo (37) was used to perform the subsequent analyses. Specifically, the function scv.pp.moments was used to compute the moments for each cell with the default parameters of normalized spliced/unspliced counts. These moments aided RNA velocity estimation in the function scv.tl.velocity, with the mode set to ‘stochastic’. Using the function scvelo.tl.velocity graph, a velocity graph representing transition probabilities was constructed based on estimated velocities. The velocity graph was then used to embed RNA velocities in the UMAP as streamlines using the scv.pl.velocity embedding stream function.

All involved surgical tumor tissues from 7 additional PDT patients were processed into paraffin blocks. Multiplex mIF staining was performed as our previous described (19), using the Opal 7-color kit (Akoya Bioscience, NEL801001KT). Tissues were sliced into 4μm sections and heat in the retrieval solution of citric acid or EDTA to recover antigens. All the relative antibodies, CD163 (ab182422, Abcam, 1:500, Opal 620), CD68 (ab269587, Abcam, 1:1000, Opal 570), SOX2 (11064-1-AP, Proteintech, 1:100, Opal 690), Doublecortin (ab18723, Abcam, 1:1000, Opal 520), were evaluated via immunohistochemistry. Briefly, the sections were dewaxed with xylene and ethanol was used for rehydration. Microwave treatment was performed for antigen retrieval with buffer (pH 9.0 or 6.0). Next, all sections were blocked using an antibody diluent/block (72424205; Akoya Bioscience), Slides were incubated with a primary antibody followed by secondary reagents and tyramide signal amplification reagents at room temperature for 10 min (Opal 480, Opal520, Opal 620, and Opal 690, Akoya Bioscience, 1:100). MWT antigen retrieval was performed until all markers were stained. Finally, DAPI-based nuclear staining was performed at room temperature for 5 min. Then the sections were mounted using anti-Fade fluorescence mounting medium (ab104135, Abcam), Samples were scanned by Vectra Polaris Automated Quantitative Pathology Imaging System and analyses were performed with QuPath software.

To explore the cellular composition of teratoma, we generated single-cell gene expression matrices of clinically annotated immature ovarian teratomas from 3 patients ( Table 1 and Supplementary Figure 1 ). After quality control, we acquired transcriptomes of 22,153 cells, which can be divided into 28 clusters based on cluster-specific markers obtained by gene differential expression analysis, including immune cells, epithelial cells, fibroblasts, and neuron cells ( Figure 1A and Supplementary Table 1 ). These clusters were refined and manually annotated using canonical cell type markers, such as meningeal fibroblast (ZIC1, CXCL12), schwann cells (MPZ, PLP1), chondrogenic fibroblast (COL2A1, SOX9), and immature neuron (DCX, MAP2) ( Figures 1A, B and Supplementary Table 2 ). Also, the expression of canonical marker genes was visualized for each cell type to assess the robustness of cell type annotations ( Supplementary Figure 2A ). To further validate the cell type annotations, we assessed the relationship between the expression signature of each teratoma cell type and that of the human fetal cell atlas (25), revealing that each teratoma cell type generally correlates with at least one human fetal cell type ( Supplementary Figure 2B ). For instance, immature neurons from the teratoma correlate with fetal human intestine neurons, airway epi/ciliated cells correlate with human fetal lung ciliated epithelial cells, whereas pericytes present at intervals along the walls of capillaries and contain contractile fibers as vascular smooth muscle thus correlate accordingly (38). Such correlation confirmed the similarity between different components of teratomas and that of fetal cells.

Moreover, given injection of hPSC into immunodeficient mice can generate cell-derived immature teratomas (CDTs), which has been profiled through scRNA-seq ( Supplementary Table 3 ) (15), we thus compared the expression signatures of their correlated cell types with PDTs ( Supplementary Figure 2C ). Not surprisingly, most of the cell types in PDTs are strongly correlated with the expected cell types in CDTs, while some discrepancies were observed possibly induced by the differences in developmental stage and xenograft-specific expression ( Figure 1C ).

To further investigate the similarities and differences between PDTs and CDTs, we first compared the proportions of 31 cell types. Around half of the cell types (15 out of 31) were shared by both teratoma types (e.g., radial glia, immature neuron, and schwann cells), whereas eight were distinct in PDTs (e.g., embryonic stem cells, endothelial cells and pericytes) and eight were distinct in CDTs (e.g., retinal neurons, malenoblasts, and cardiac/skeletal muscle) ( Figure 2A ). Considering the germ layer, most of the cells from both PDTs and CDTs are derived from mesoderm and ectoderm but less from endoderm ( Figure 2B ), suggests that endoderm is less prevalent during early development. Despite distinct cell types observed in PDTs and CDTs, there is no statistically significant difference between these two teratomas types in terms of germ layer components ( Figure 2C ). We next compared the heterogeneity between PDTs and CDTs in their shared 15 cell types and found that all cell types but not immune cell from PDTs have significant lower purity than that from CDTs ( Figure 2D ), suggesting the heterogeneity and complexity of PDTs. Moreover, cell-cell interaction network of the shared cell types was also estimated and compared between PDTs and CDTs, revealing higher overall interactions among PDT-derived cell types ( Figure 2E ). Three interaction modules (e.g., high, medium, and low) were defined based on the cell-cell interaction counts in both PDTs and CDTs ( Figure 2E ). Interactions between all five fibroblast cell types and other cell types are found in high/medium modules, particularly meningeal fibroblast and cardiomyocyte/myofibroblast are in high modules in both teratoma types, whereas interactions between all four epithelial cell types and other cell types are found in low/medium module, particularly retinal epithelial are in low modules in both teratoma types ( Figure 2E ). For neuroectodermal cells, immature neuron exhibited the least cell-cell interactions with all the rest cell types, whereas cycling neural progenitor and radial glia cells increase the intensity of their interactions from low in CDTs to medium in PDTs ( Figure 2E ), suggesting more activated neural development in malignant teratomas.

Given immature and embryonic tissue from all three germ layers can be found in immature teratomas, and the amount of primitive neuroectodermal tissue is highly predictive of prognosis (3, 39), we first focused on the neuroectodermal cell types. In this study, PDT02 had a lower proportion of ectodermal cells than the other two patients ( Figure 3A and Supplementary Figure 3A ), consistent with her event-free survival status and low AFP level ( Table 1 ). Three major neuroectodermal cell types were characterized by specific markers ( Supplementary Table 2 ), including radial glia (SOX2, NES), cycling neural progenitor (HMGB2, MKI67), and immature neuron (STMN2, DCX) ( Figure 3B ), which is perfectly matched to the cell types with the same marker genes in CDTs ( Supplementary Figure 3B ). As validation by multiplexed immunofluorescence (mIF) staining, DCX⁺ and SOX2⁺ cells were either present in DCX-/SOX2-specific region or expressed in the mixed region in mutually exclusive manner, and exhibited different morphological features ( Figure 3C and Supplementary Figure 3C ). According to the guideline of gynecologic pathology and WHO-defined histopathology, immature teratomas are characterized by immature neuroepithelium (rosettes and tubules), and neuroectodermal elements can be highlighted by neural markers, including NSE, GFAP, and SOX2, among them, SOX2 is more specific for immature neural tissue and diagnosis of immature teratoma (4). Therefore, enriched presentation of these genes indicated the dominate role of radial glia to determine immature teratoma.

Next, we performed PCA analysis to estimate the transcriptional similarity of neuroectodermal cell types and found that all three PDTs were clustered together based on the transcriptome profile of both radial glia and cycling neural progenitor, whereas immature neuron from PDT02 shared less similarity with that from PDT01/PDT03 ( Figure 3D and Supplementary Figure 3D ), suggesting the inter-heterogeneity among PDTs can be featured by immature neuron cells. Moreover, significant association of transcriptome disparity in H1-based CDTs (i.e., distance between ter5 and other CDTs in PCA) with the cell proportions was observed in immature neuron but not in radial glia and cycling neural progenitor clusters ( Figures 3D, E and Supplementary Figure 3E ). Next, we screened the differentially expressed genes (DEGs) of immature neuron cells between PDT02 and PDT01/PDT03, and found that a series of the top up-regulated genes (e.g., STMN2, STMN1, CD24, SOX11, SOX4, and DCX) were canonical markers for immature neurons and contribute to neuronal differentiation and migration during neurogenesis ( Figure 3F ) (40–44). For instance, SOX4 and SOX11, which encode two neurogenesis-related transcription factors, are functionally essential for hippocampal neurogenesis (42). Consistently, some of these genes exhibited similar positive correlations between their expression levels and immature neuron proportions in CDTs ( Figure 3G ), with SOX4 and DCX reached statistical significance in terms of linear and/or categorical approaches ( Figure 3H ). To further validate the findings, we performed mIF staining of DCX in validation cohort (n = 7) ( Supplementary Table 4 ), and identified the positive correlations between the overall expression levels of DCX and teratoma grade ( Figure 3I and Supplementary Figure 3F ). However, due to the limited number of patients enrolled in the validation cohort, the correlation did not reach statistical significance ( Supplementary Figure 3F ). On the other hand, many top down-regulated genes are related to gliogenesis and glial differentiation (e.g., OLIG1, OLIG2, APOD, and S100B) (45) ( Figures 3F, G ), while expression of S100B and SCRG1 exhibited negative correlation with immature neuron proportions in CDTs ( Supplementary Figure 3G ). In contrast, no trend has been observed for the correlation between markers gene expressions and the cell proportions in radial glia and cycling neural progenitor ( Supplementary Figures 3H, I ).

Given that neurogenesis process in teratomas could be critical to patients’ prognosis, we conducted Monocle-based cell trajectory analysis to investigate the evolutionary transitions of three neuroectodermal cell types described above. Two distinct differentiated paths were revealed along with the evolutionary trajectory from cycling neural progenitors to radial glia/immature neuron, enhancing the definition of three states ( Figure 4A ), including State1 (cycling neural progenitors dominant, defined as root), State2 (radial glia dominant) and State3 (immature neuron dominant) ( Figure 4B ). Similar trends were observed in CDTs ( Supplementary Figures 4A, B ). Neuroectodermal cells from PDT01/PDT02 are more enriched in State1/State2 compared with PDT03 ( Figure 4C ), which may contribute to the metastatic ability of PDT03 according to our clinical follow-up ( Table 1 ). Finally, tracing the gene fluctuation along biforked trajectories, genes that differ as a function of pseudotime were identified as incrementally up-regulated or down-regulated from root to State2 or State3, and can be grouped into three clusters (Cluster1, Cluster2 and Cluster3). Genes in Cluster1 were highly expressed by immature neurons and associated with nervous system development and neuron actin remodeling ( Figure 4D ). Interestingly, many of them are up-regulated genes identified in comparison between PDT01/PDT03 and PDT02 ( Figures 3F , and 4D ). Genes in Cluster2 were up-regulated in radial glia cells, consisting of many down-regulated genes described above and enriched in several neuron development related pathways, including glial cell differentiation ( Figures 3F , and 4D ). Finally, genes in Cluster3 were up-regulated in cycling neural progenitors and were enriched in chemokine and cytokine signaling regulation ( Figure 4D ).

We also conducted RNA velocity-based trajectory inference to delineate cellular differentiation of neuroectodermal cells. Consistently, cycling neural progenitors were identified as the origin of the global differentiation processes during neuroectoderm diversification to functional neural subtypes ( Figure 4E ). Moreover, samples from PDT02 exhibited suppressed neural cells differentiation from cycling neural progenitor to immature neuron/radial glia ( Figure 4E ), which may contribute to the lower abundance of immature neuron/radial glia in PDT02 than that in PDT01/PDT03.

To elucidate transcriptional states across different cell types, SCENIC was performed to identify potential regulatory elements ( Supplementary Figure 5 and Supplementary Table 5 ). For neuroectodermal cell types, BHLHE22, ONECT22, and NHLH1 are the specific regulons in immature neurons ( Figure 5A ), consistent with their reported role in cell-type determination, proliferation, and differentiation within the developing nervous system (46, 47). Interestingly, only the top transcription factors enriched in immature neuron are present at a significant lower level in PDT02 than PDT01/PDT03 ( Figures 5B, C ), suggesting the relative inactivation status of neuroectodermal cell in PDT02. Consistently, the top transcription factors are specifically over-presented in immature neuron of H1-derived CDTs ( Figure 5D ), but fail to correlate with the immature neuron proportions.

Next, we investigated the cell-cell crosstalk among different cell types through CellphoneDB and found more prolific overall interactions in PDT01/PDT03 than those in PDT02 ( Figure 5E and Supplementary Figure 6A ). For neuroectodermal cell types, radial glia has much fewer interactions with fibroblasts/epithelial cell types in PDT02 than those in PDT01/PDT03, whereas the differences were less pronounced in cycling neural progenitor and immature neuron ( Figure 5E ). To focus on the potential ligands that drive the transcriptome difference in three PDTs, we performed Cellchat to estimate the number of ligand-receptor interactions between neuroectodermal cells and other cell types. Consistent with the observations in CellphoneDB, PDT02 showed decreased number of interactions with primarily fibroblasts and epithelial cells in all three neuroectodermal cell types ( Figure 5F and Supplementary Figure 6B ), which might contribute to the better prognosis in PDT02.

Large differences were observed in terms of the cell type and proportion of immune cells between PDTs and CDTs. For instance, myeloid cells were dominantly identified in CDTs, compared with multiple immune cell types (e.g., myeloid, T cells, NK cells, B cells, and neutrophils) identified in PDTs ( Figures 1A , 6A , and Supplementary Figure 2C ), while the immune cells exhibited distinct copy number variation (CNV) patterns compared to other cells across all PDT samples ( Supplementary Figure 7 ), suggesting that most of the immune cells in PDTs were infiltrated rather than teratoma-derived. Additionally, immune cells were the only cell type that exhibited no difference between CDT and PDT in terms of cell purity ( Figure 2D ), further suggesting the altered origin of the immune cells in PDT. Because CDTs were generated in mice, we examined the mouse-derived cells from the scRNA-seq data of CDTs, and identified a large number of immune cells (mostly macrophage and dendritic cells because the mice are immunodeficient) ( Figure 6B ). Given myeloid is the most abundant cell types in both PDTs and CDTs, as well as the reported role of macrophages in aiding the initiation and progression of teratomas (48), we compared the expression profile of myeloid cells between PDT02 and PDT01/03, the DEGs were enriched in multiple immune related pathways, including major histocompatibility complex class II (MHC II) signaling ( Figure 6C ), compared to non-immune related pathways enriched in other cell types ( Supplementary Figure 8A ). Expression of multiple MHC-II components encoding genes were upregulated in myeloid cells from PDT02 compared to PDT01/PDT03, and the correlation of their expression with immature neuron proportion is in line with that of infiltrated myeloid cells of CDTs ( Figures 6D, E ), but not myeloid cells derived from CDTs ( Supplementary Figure 8B ), thus further confirmed the potential role of infiltrating immune cells in teratomas progression. Furthermore, we determined the gene signature based myeloid functional scores in each PDT, we found that PDT02, which has a lower proportion of immature neurons, showed an enrichment of M1 polarization and antitumor cytokine scores ( Figure 6F ). Conversely, PDT01/PDT03, with a higher proportion of immature neurons, showed an enrichment of M2 polarization, phagocytosis, protumor cytokine, and angiogenesis scores ( Figure 6F ). This was further confirmed by mIF staining in the validation cohort, which showed a positive association between the percentage of M2 macrophages and teratoma grading, whereas a negative correlation between the percentage of M1 macrophages as well as M1/M2 ratio and teratoma grading ( Figures 6G-I and Supplementary Figure 9A ). Interestingly, the infiltrated macrophages were more abundant in DCX-specific and mixed region than in SOX2-specific region, indicating different tumor immune microenvironment among neuroectodermal subtypes ( Figure 6J and Supplementary Figure 9B ).

On the other hand, we investigated potential role of other potential infiltrated immune cells in PDTs, which is undetectable in CDTs due to immunodeficiency of the injected mice, such as T and NKT cells. PDT02 has higher proportions of T cells and NK cells than PDT01/PDT03 ( Figure 3A ). Although no pronounced difference of interactions between infiltrated immune cells and other cell types was identified among three PDTs ( Figures 5E, F ), PDT02 exhibited significant increased signature-estimated cytotoxicity and mobility but not exhaustion features of NK and T cell compared to PDT01/PDT03 ( Supplementary Figures 10A-C ), thus may contribute to the better prognosis of PDT02.