Whole Exome Sequencing Reveals Potential Single Nucleotide Polymorphisms and Copy Number Variations in 26 Sporadic Patients with Immature Teratoma.

Teratomas, a subtype of germ cell tumors, originate from 3 germ cell layers containing different types of tissue (such as hair, bone, and muscles). 1 , 2 Histopathologically, teratomas are divided into mature teratoma and immature teratoma. 3 Immature teratoma is a rare and malignant ovarian germ cell tumors related to a variable low potential of distant metastasis, which depends on the histological grade. 4 The lungs, liver and brain are the most distant sites of blood-borne metastasis of immature teratoma. 5 In addition, rare sites of distant spread include abdominal wall musculature and soft tissue of the thigh. 6 , 7 Clinically, immature teratoma is a curable disease. Adjuvant treatment with etoposide, bleomycin and cisplatin is the standard in stages greater than stage IA grade I. In addition, operative intervention combined with chemotherapy is required for advanced stage immature teratoma patients. 4 Interestingly, recurrent immature teratoma is known to undergo pathological retroconversion to contain mature elements following chemotherapeutic treatment, which is a unique characteristic of immature teratoma. 5 There are many causes of immature teratoma. It is suggested that copy-neutral loss of heterozygosity (resulting from meiotic errors) could be sufficient to generate immature teratomas from germ cells, which are characterized by extensive allelic imbalances, copy number alterations and a paucity of somatic mutations. 8 Currently, most studies on immature teratoma are retrospective clinical analyses, whereas studies focusing on its genetic and molecular mechanism remain limited. Therefore, investigations of single nucleotide polymorphisms (SNPs), copy number variations (CNVs) and potential drug targets may become key directions in future research. Meanwhile, whole exome sequencing (WES) has been increasingly applied in clinical genetic diagnosis and genomic research. Peripheral blood is a commonly used sample in WES and is easy to obtain. Tumor tissue from immature teratoma is clinically limited and difficult to obtain routinely. Furthermore, several studies on ovarian‑related tumors have also been conducted using peripheral blood samples. 9 , 10 In view of this, peripheral blood samples from 26 sporadic immature teratoma patients were collected for WES to identify new mutated genes, which may provide clues for the study of function and pathogenesis of immature teratoma.

After the screening, information of SNPs and indels of 26 sporadic patients with immature teratoma was obtained. Total variation type, type of base change and top 20 mutant genes is shown in Figure 1 . Nonsynonymous SNP was the major variation type ( Figure 1A ). C>T and G>A were the major base change types in SNP class ( Figure 1B ). It is noted that top 20 mutant genes were identified ( Figure 1C ), such as mucin 3A, cell surface associated (MUC3A), G protein regulated inducer of neurite outgrowth 2 (GPRIN2), mucin 2, oligomeric mucus/gel-forming (MUC2), Fc gamma binding protein (FCGBP), FSHD region gene 2 family member C (FRG2C). Interestingly, a total of 24 common mutated genes were identified in 26 sporadic immature teratoma patients. Among which, a total of 5 common mutated genes with common mutation sites were identified, including contacting associated protein family member 3 (CNTNAP3), Myb-related transcription factor, partner of profiling (MYPOP), FRG2C, GPRIN2 and MUC3A ( Table 2 ). Mutation type of FRG2C and MYPOP was frameshift_deletion. Mutation type of CNTNAP3, GPRIN2 and MUC3A was nonsynonymous_SNP. In addition, 2 common CNVs (deletion state) were screened out ( Table 3 ), which involving 6 genes, such as pre-mRNA processing factor 31 (RP11, also called PRPF31). Table 2 5 Common Mutated Genes with Common Mutation Sites in 26 Sporadic Immature Teratoma Patients Reference Gene Nucleotide Change AA Change CNTNAP3 Chr9- 39088492: T >A T1051S FRG2C Chr3- 75665673: AG >A R161Sfs*5 GPRIN2 Chr10- 46549538: A >G L400P MUC3A Chr7- 100958510: C>G T2244S MUC3A Chr7- 100957778: C >A T2000N MUC3A Chr7- 100965345: C > A A3149E MUC3A Chr7- 100957777: A >G T2000A MYPOP Chr19- 45890723: C >CG R367Pfs*83 Table 3 2 Common CNVs in 26 Sporadic Immature Teratoma Patients Chr Start_End Length Gene CNV_state 4 111,463,907-111,542,525 78618 PITX2:P|ENPEP:P|RP11-380D23.2:F Deletion 17 15,546,035-15,603,091 57056 RP11-640I15.1:F|TRIM16:P|ZNF286A:P|ZNF29P:F|RP11-385D13.1:P Deletion Figure 1 Total variation type ( A ), type of base change ( B ) and top 20 mutant genes ( C ) in 26 sporadic patients with immature teratoma. Three horizontal bar charts showing variant types, base changes and top 20 mutant genes. The image A showing the title “Variant Classification”. Horizontal axis label: (unclear), unit: (unclear), range: 0 to 15000 with ticks at 0, 5000, 10000, 15000. Vertical axis label: Variant Classification, unit: none. Bars and values: nonsynonymous single nucleotide variant about 17000 (extends past 15000 tick); unknown about 1000; nonframeshift deletion about 500; frameshift deletion about 400; nonframeshift insertion about 300; stopgain about 200; frameshift insertion about 100; startloss about 0; stoploss about 0. Pattern: nonsynonymous single nucleotide variant is the dominant category. The image B showing the title “Single nucleotide variant Class”. Horizontal axis label: (unclear), unit: (unclear), range: 0 to 4000 with ticks at 0, 1000, 2000, 3000, 4000. Vertical axis label: Single nucleotide variant Class, unit: none. Bars and values: cytosine to thymine 3813; guanine to adenine 3537; adenine to guanine 1749; thymine to cytosine 1729; guanine to cytosine 1106; cytosine to guanine 1031; cytosine to adenine 995; guanine to thymine 804; adenine to cytosine 546; thymine to guanine 545; thymine to adenine 476; adenine to thymine 425; cytosine guanine to thymine adenine 1; adenine cytosine to guanine thymine 1. Pattern: cytosine to thymine and guanine to adenine are the highest. The image C showing a horizontal stacked bar chart of top mutant genes with legend categories: frameshift insertion, stopgain, nonframeshift insertion, frameshift deletion, nonframeshift deletion, unknown, nonsynonymous single nucleotide variant. Horizontal axis label: (unclear), unit: (unclear), range: 0 to 1500 with ticks at 0, 500, 1000, 1500. Vertical axis label: gene names, unit: none. Genes listed: MUC3A, GPRIN2, MUC2, FCGBP, MUC19, PRSS2, PRAMEF18 times 3 and PRAMEF22, PRAMEF8, CNTNAP3, GOLGA6L7, NBPF8, PRAMEF19, FRG2C, ZNF280A, HERC2, OR10G2, OR4C46, NBPF26, LILRB1, LILRB3. Bar lengths: MUC3A about 1500; GPRIN2 about 1200; MUC2 about 300; FCGBP about 250; remaining genes each about 50 to 150. Relationship across A, B, C: A summarizes variant types, B breaks down single nucleotide variant base changes, C lists gene-level mutation counts by mutation type categories. 5 Common Mutated Genes with Common Mutation Sites in 26 Sporadic Immature Teratoma Patients 2 Common CNVs in 26 Sporadic Immature Teratoma Patients Total variation type ( A ), type of base change ( B ) and top 20 mutant genes ( C ) in 26 sporadic patients with immature teratoma. The physical/chemical properties of amino acid and conserved domain after mutations of CNTNAP3, FRG2C, GPRIN2, MUC3A and MYPOP were analyzed. Amino acid number, molecular weight and theoretical PI of CNTNAP3 are 1127, 124,115.80 and 6.67, respectively. In amino acid composition, Leu (9.80%) content is the highest. The total number of negatively and positively charged residues is 107 and 101, respectively. The extinction coefficient is 1.324/1.306. The instability index, fat index and the average hydrophilic are 37.61, 75.09 and −0.281, respectively. After mutation, the molecular weight (124101.78) decreased slightly, and the instability index changed to 38.39. The mutant amino acid is located in the Laminin_G_2 domain ( Figure 2A ). Amino acid number, molecular weight and theoretical PI of FRG2C are 282, 30,798.11 and 7.09, respectively. In amino acid composition, Ser (12.4%) content is the highest. The total number of negatively and positively charged residues is 35 and 35, respectively. The extinction coefficient is 0.571/0.551. The instability index, fat index and the average hydrophilic are 58.06, 54.40 and −0.886, respectively. After mutation, molecular weight (30729.00) and theoretical PI (6.75) decreased, and instability index became 57.19. The mutant amino acid was located in the FGR2 domain ( Figure 2B ). Amino acid number, molecular weight and theoretical PI of GPRIN2 are 458, 47,450.29 and 6.28, respectively. In amino acid composition, Ala (13.80%) content is the highest. The total number of negatively and positively charged residues is 43 and 38, respectively. The extinction coefficient is 1.090/1.075. The instability index, fat index and the average hydrophilic are 63.51, 70.24 and −0.320, respectively. After mutation, the molecular weight (47434.25) decreased slightly, and the instability index became 63.89. The mutant amino acid was located in the GRIN_C domain ( Figure 2C ). Amino acid number, molecular weight and theoretical PI of MUC3A are 3323, 345,126.77 and 5.14, respectively. In amino acid composition, Thr (29.40%) content is the highest. The total number of negatively and positively charged residues is 171 and 91, respectively. The extinction coefficient is 0.504/0.500. The instability index, fat index and the average hydrophilic are 43.76, 54.97 and −0.217, respectively. After mutation, the molecular weight (345154.58) increased slightly. The mutant amino acid was not in its conserved domain ( Figure 2D ). Amino acid number, molecular weight and theoretical PI of MYPOP are 399, 42,508.45 and 10.08, respectively. In amino acid composition, Pro (20.80%) content is the highest. The total number of negatively and positively charged residues is 39 and 53, respectively. The extinction coefficient is 0.755/0.752. The instability index, fat index and the average hydrophilic are 94.21, 66.17 and −0.677, respectively. After mutation, the molecular weight (345154.58) increased slightly. The mutant amino acid was not in its conserved domain ( Figure 2E ). Figure 2 Conserved domain analysis of mutated CNTNAP3 ( A ), FRG2C ( B ), GPRIN2 ( C ), MUC3A ( D ) and MYPOP ( E ) in 26 sporadic patients with immature teratoma. The red arrow in ( A – C ) represents the location of the mutant amino acid within the conserved domain. Screenshot collage of five conserved domain search result pages with graphical summaries and domain hit tables. A header area includes the text Protein Classification and Graphical summary. A horizontal sequence ruler is shown with a highlighted region and a downward arrow marker. A table titled List of domain hits shows columns Name, Accession, Description, Interval and Evalue. The image B showing a screenshot titled Conserved domains on a protein sequence. The page includes Graphical summary with a horizontal ruler and a long highlighted domain bar labeled FRG2 superfamily. A List of domain hits table is shown with the same columns. The image C showing a screenshot titled Conserved domains on a protein sequence. The Graphical summary highlights a domain labeled GRIN C and GRIN C superfamily. A List of domain hits table appears below. The image D showing a screenshot titled Conserved domains on a protein sequence. The Graphical summary shows a horizontal ruler with small marked regions near the right side. A List of domain hits table is shown. The image E showing a screenshot titled Conserved domains on a protein sequence. The Graphical summary highlights regions labeled Myb SANT like DNA binding domain and GLE superfamily. A List of domain hits table is shown with columns Name, Accession, Description, Interval and Evalue. Conserved domain analysis of mutated CNTNAP3 ( A ), FRG2C ( B ), GPRIN2 ( C ), MUC3A ( D ) and MYPOP ( E ) in 26 sporadic patients with immature teratoma. The red arrow in ( A – C ) represents the location of the mutant amino acid within the conserved domain. Based on GO analysis of 24 mutated genes, golgi lumen was the only one significantly enriched cytological component. A remarkably enriched biological process of mutated genes (involved MUC3A, MUC2 and MUC12) was identified, such as activation of the innate immune response. It is a pity that no enriched molecular function of mutated genes was found. According to KEGG analysis of 24 mutated genes, amoebiasis and gastric cancer were the only two significantly enriched signaling pathways of mutated genes. All significantly enriched biological processes, cytological components and signaling pathways are listed in Table 4 . In addition, PPI network was constructed based on 24 mutation genes ( Figure 3 ). The highest degree gene was MUC12. Moreover, MUC12 interacted with polypeptide N-acetylgalactosaminyltransferase 12 (GALNT12), MUC3A and FCGBP. Table 4 Functional Enrichment Analysis of 24 Mutated Genes in 26 Sporadic Immature Teratoma Patients ID Description P value FDR Gene ID Count BP GO:0016266 O-glycan processing 1.12E-06 3.01E-05 MUC3A/MUC2/MUC12 3 GO:0006493 Protein O-linked glycosylation 5.77E-06 4.54E-05 MUC3A/MUC2/MUC12 3 GO:0002223 Stimulatory C-type lectin receptor signaling pathway 7.39E-06 4.54E-05 MUC3A/MUC2/MUC12 3 GO:0002220 Innate immune response activating cell surface receptor signaling pathway 8.20E-06 4.54E-05 MUC3A/MUC2/MUC12 3 GO:0002758 Innate immune response-activating signal transduction 8.41E-06 4.54E-05 MUC3A/MUC2/MUC12 3 GO:0002218 Activation of innate immune response 1.62E-05 7.29E-05 MUC3A/MUC2/MUC12 3 GO:0045089 Positive regulation of innate immune response 5.51E-05 0.000212538 MUC3A/MUC2/MUC12 3 GO:0006486 Protein glycosylation 7.74E-05 0.000220905 MUC3A/MUC2/MUC12 3 GO:0043413 Macromolecule glycosylation 7.74E-05 0.000220905 MUC3A/MUC2/MUC12 3 GO:0002833 Positive regulation of response to biotic stimulus 8.50E-05 0.000220905 MUC3A/MUC2/MUC12 3 GO:0070085 Glycosylation 9.00E-05 0.000220905 MUC3A/MUC2/MUC12 3 GO:0045088 Regulation of innate immune response 0.000153637 0.000345683 MUC3A/MUC2/MUC12 3 GO:0009101 Glycoprotein biosynthetic process 0.000204504 0.000424738 MUC3A/MUC2/MUC12 3 GO:0031349 Positive regulation of defense response 0.000229828 0.000443239 MUC3A/MUC2/MUC12 3 GO:0002831 Regulation of response to biotic stimulus 0.000358931 0.000609945 MUC3A/MUC2/MUC12 3 GO:0009100 Glycoprotein metabolic process 0.000361449 0.000609945 MUC3A/MUC2/MUC12 3 GO:0002429 Immune response-activating cell surface receptor signaling pathway 0.000534341 0.000801511 MUC3A/MUC2/MUC12 3 GO:0002757 Immune response-activating signal transduction 0.000534341 0.000801511 MUC3A/MUC2/MUC12 3 GO:0030277 Maintenance of gastrointestinal epithelium 0.00813734 0.011563589 MUC2 1 GO:0010669 Epithelial structure maintenance 0.010714562 0.014464659 MUC2 1 GO:0022600 Digestive system process 0.03617272 0.046507782 MUC2 1 CC GO:0005796 Golgi lumen 2.22E-05 4.43E-05 MUC3A/MUC2/MUC12 3 KEGG hsa05146 Amoebiasis 0.012573964 0.01836785 MUC2 1 hsa05226 Gastric cancer 0.01836785 0.01836785 MUC2 1 Abbreviations : BP, biological processes; CC, cytological component; KEGG, Kyoto Encyclopedia of Genes and Genomes; FDR, false discovery rate. Figure 3 PPI network in 26 sporadic patients with immature teratoma. Light blue lines: Known interactions derived from curated databases.; Yellow lines: Associations supported by textmining evidence; Black lines: Associations derived from gene co-expression data; Light purple lines: Associations based on protein homology evidence. Nodes are connected if at least one line of evidence supports a functional or physical association between the corresponding proteins; the absence of an edge between two nodes indicates that no significant evidence for such an association was found in the analyzed database resources. Protein interaction network with nodes and lines indicating associations. Functional Enrichment Analysis of 24 Mutated Genes in 26 Sporadic Immature Teratoma Patients Abbreviations : BP, biological processes; CC, cytological component; KEGG, Kyoto Encyclopedia of Genes and Genomes; FDR, false discovery rate. PPI network in 26 sporadic patients with immature teratoma. Light blue lines: Known interactions derived from curated databases.; Yellow lines: Associations supported by textmining evidence; Black lines: Associations derived from gene co-expression data; Light purple lines: Associations based on protein homology evidence. Nodes are connected if at least one line of evidence supports a functional or physical association between the corresponding proteins; the absence of an edge between two nodes indicates that no significant evidence for such an association was found in the analyzed database resources. Based on the drugs-24 mutated genes interaction analysis, only compounds of CHEMBL35482, FLUOROURACIL, and DECITABINE interacted with MUC2. It is a pity that there was no interaction between drugs and the remaining 23 mutated genes. This suggested that MUC2 could be considered as a potential molecular target of CHEMBL35482, FLUOROURACIL, and DECITABINE in the therapy of immature teratoma. In order to explore the SNPs and CNVs of mutated genes in ovarian cancer, SNPs mutation data (involving 620 samples) and CNVs mutation data (involving 436 samples) of mutated genes in ovarian cancer were downloaded from TCGA dataset. In the SNPs data of ovarian cancer, 9 genes mutations were analyzed through “maftools” in the R package, including FCGBP, CNTNAP3, GPRIN2, MUC3A, GALNT12, MUC12, FRG2C, MYPOP and MUC2 ( Figure 4A ). The mutation frequency of FCGBP and CNTNAP3 was, respectively, 2% and 1%, which indicated rare mutation rate in ovarian cancer. In addition, the amplification (gain) and deletion (loss) frequency of FCGBP, CNTNAP3, GPRIN2, MUC3A, GALNT12, MUC12, FRG2C and MYPOP were calculated based on CNVs data of ovarian cancer ( Figure 4B ). FCGBP, CNTNAP3, GPRIN2, MUC3A and MUC12 had higher amplification frequency of CNV. FRG2C, GALNT12 and MYPOP had higher deletion frequency of CNV. Figure 4 Investigation of SNP ( A ) and CNV ( B ) of mutated genes in ovarian cancer. Different types of data visualizations such as an oncoprint and a stacked bar chart for gene variants. The image A showing an oncoprint with stacked bar charts summarizing single nucleotide polymorphism mutations by gene. The x-axis label and unit are unclear. The y-axis lists genes: FCGBP, CNTNAP3, GPRIN2, MUC3A, GALNT12, MUC12, FRG2C, MYPOP, MUC2. The top stacked bar chart has a y-axis range from 0 to 2 (unit unclear) and shows per-sample mutation counts. The right horizontal bar chart shows per-gene mutation frequency with labeled values: FCGBP 2 percent, CNTNAP3 1 percent, GPRIN2 0 percent, MUC3A 0 percent, GALNT12 0 percent, MUC12 0 percent, FRG2C 0 percent, MYPOP 0 percent, MUC2 0 percent. The bottom stacked bar chart shows per-sample mutation categories with legend entries: Missense Mutation, Multi Hit, Nonsense Mutation and single nucleotide variant classes labeled C to T, C to A, C to G, T to A, T to C, T to G (some legend text unclear). The image B showing a stacked bar chart of copy number variation frequency by gene. The x-axis label is unclear; x-axis categories are FRG2C, MUC3A, MUC12, CNTNAP3, GALNT12, GPRIN2, FCGBP, MYPOP. The y-axis label is CNV Frequency (100 percent) with a range from 0.00 to 0.35. Legend entries are GAIN and LOSS. Approximate total bar heights: FRG2C about 0.16, MUC3A about 0.15, MUC12 about 0.16, CNTNAP3 about 0.19, GALNT12 about 0.13, GPRIN2 about 0.17, FCGBP about 0.33, MYPOP about 0.30. Approximate segment heights (GAIN, LOSS): FRG2C about 0.07 and 0.09; MUC3A about 0.12 and 0.03; MUC12 about 0.13 and 0.03; CNTNAP3 about 0.13 and 0.06; GALNT12 about 0.00 and 0.13; GPRIN2 about 0.09 and 0.08; FCGBP about 0.21 and 0.12; MYPOP about 0.06 and 0.24. Relationship between A and B: A lists mutation frequencies for the same gene set, while B shows copy number variation gain and loss frequencies for genes overlapping that set. Investigation of SNP ( A ) and CNV ( B ) of mutated genes in ovarian cancer. To investigate additionally the mRNA expression of mutated genes in other reproductive-related tumors, GSE26712 dataset and GSE10615 dataset were used to detect mRNA expression of mutated genes in ovarian cancer ( Figure 5A ) and yolk sac tumor ( Figure 5B ), respectively. GALNT12 was significantly down-regulated in ovarian cancer compared with para-carcinoma tissue. MUC2 had an up-regulated tendency, and FCGBP and GPRIN2 had a down-regulated tendency in ovarian cancer without statistical difference. In addition, GALNT12, MUC2, FCGBP and GPRIN2 had an up-regulated tendency in yolk sac tumor compared with seminoma. Figure 5 MRNA expression of mutated genes in ovarian cancer ( A ) and yolk sac tumor ( B ). ****, p <0.0001. Illegible image: axis labels, units, tick values and box plot values cannot be read reliably. Illegible image: multiple box plots are visible but the x-axis label, y-axis label, units, tick marks, group labels and exact box plot statistics (median, quartiles, whiskers, outliers) are not readable, so exact plotted values cannot be extracted. Abbreviations : ns, no statistical significance; YST, yolk sac tumor. MRNA expression of mutated genes in ovarian cancer ( A ) and yolk sac tumor ( B ). ****, p <0.0001.

In the present study, a total of 26 sporadic patients with immature teratoma were enrolled. The peripheral blood samples of these individuals were collected for whole exome sequencing (WES) analysis. Detailed clinical information of 26 patients is shown in Table 1 . This study was approved by the ethics committee of the fourth hospital of hebei medical university. Written consent was obtained from enrolled individuals. Table 1 Clinical Information of 26 Sporadic Patients with Immature Teratoma Number Age (Years Old) Pelvic Mass Abdominal Pain Abdominal Distension Adnexal Mass Frequent Urination Excision Stage 1 30 Yes No Yes No No Yes III 2 46 No Yes No No Yes No II 3 44 No Yes No No Yes No II 4 21 No No No Yes No Yes IA 5 20 No No Yes Yes No No IA 6 77 No No Yes No No No IA 7 30 No No No No No Yes IA 8 73 No No No Yes No Yes IA 9 56 No No Yes No No No IIIC 10 13 Yes No No No No Yes IIB 11 28 No No No Yes No Yes IIIC 12 34 No Yes Yes Yes No No IA 13 28 No No No Yes No No IA 14 24 Yes No No No No No IA 15 35 No No No No No Yes IC1 16 28 No Yes Yes No No No IA 17 34 No No No No No Yes IA 18 15 Yes No No No No No IA 19 16 No No No Yes No No IC1 20 19 No No No No No No IA 21 23 Yes No No No No No IA 22 24 Yes No No No No No IA 23 39 No No No No Yes No IIIC 24 24 No No No Yes No No IA 25 25 Yes No No No No No IC 26 22 No No No No No No IA Clinical Information of 26 Sporadic Patients with Immature Teratoma Genomic DNA was extracted from peripheral blood samples. The DNA concentration was determined by microplate analyzer. The integrity and purity of DNA were determined by 1% of agarose gel. 1μg of genomic DNA was interrupted to 250–300bp by ultrasound and purified by AxyPrep Mag PCR Clean up Kit. The terminal repair reaction system was purified after a period of moderate temperature reaction. The reaction system with “A” at the end and joint connection was prepared. A certain amount of polymerase chain reaction (PCR) products were hybridized and hybridized with Agilent SureSelect Hybridization and Wash Kits. The eluted product was amplified and purified by the AxyPrep Mag PCR clean up Kit. Constructed libraries were tested for quality using Agilent 2100 Bioanalyzer and ABI StepOnePlus Real-Time PCR System. Qualified libraries were sequenced on Xten platform using 150PE sequencing. The original image data obtained by WES was transformed to sequence data by base calling, which is called raw data or raw reads. Raw reads were quality-controlled with fastQC and SeqPrep. The joint sequence in reads and the reads that have no inserted fragment (due to self-connection of the joint) was removed. The bases with low quality (mass value <20) at the end of the sequence (3’ end) were trimmed. Reads with an N ratio of more than 10% were removed. The adapter and sequences whose length is less than 20bp after quality trimming was discarded. The high-quality sequencing sequences (obtained after quality control) were compared with the designated reference genome (version: hg19) by BWA software ( https://github.com/lh3/bwa ). Picard software was used to remove repeated sequences. Through GATK software, the mismatch regions caused by indel were recalculated by BaseRecalibrator module. In addition, base quality values of reads in the BAM file were recalibrated by RealignerTargetCreator module. SNP is a phenomenon in which most nucleotide sequences of the same chromosome or the same site in different individuals are identical, but only one base is different. Indel is an insertion or deletion of multiple bases. Deletions or insertion at coding sites or splicing sites may alter protein translation. Strelka software was used for mutation detection of SNP and indel. ANNOVAR and tapes software were used to annotate mutation loci, including gene structure, genomic feature, the known mutation database, function of mutation-related gene and hazard prediction of non-synonymous mutation. Refseq and Gencode were used to annotate the gene structure of the mutated loci. The SIFT, PolyPhen, MutationAssessor, LRT and other methods were used to comprehensively evaluate the impact of non-synonymous mutations on the disease. The dbSNP database, SNP database of 1000 genome, cosmic database (involved knowing tumor somatic mutation) and ESP6500 mutation database were used for annotation. Candidate mutations were identified by eliminating mutations with a mutation frequency ≥0.01 in 1000 Genome or ExAC databases, synonymous mutations and mutations located in introns. In the hazard prediction, those mutations predicted to be benign/possibly benign were further filtered from ClinVar, InterVar, SIFT and Polyphen databases. In addition, in order to reduce the false positive rate of CNVs, fragments with length ≤300bp were filtered out. High-quality CNVs were filtered through annotation information (prediction quality score/default prediction score >40). To improve the accuracy of the results, only those genes in the CNVs region were retained. The physical/chemical properties of amino acid after gene mutation were predicted by using the online database protparam ( http://web.expasy.org/protparam/ ). Additionally, conserved domain analysis after gene mutation was performed through the Conserved Domain Search Service ( https://www.ncbi.nlm.nih.gov/cdd/?term ). In order to investigate the function of identified mutation genes, Gene Ontology (GO) and the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis were performed. False discovery rate (FDR) <0.05 was the threshold of significantly enriched GO and KEGG terms. Additionally, STRING ( https://string-db.org/ ) software was used to perform PPI analysis on the proteins encoded by mutated genes. Drug-Gene Interaction Database (DGIdb) is a database of drug–gene interactions that provides information about the association of genes with known or potential drugs. It has more than 14,000 drug–gene interactions involving 2,600 genes and 6,300 drugs that target these genes. It also includes 6,700 other genes. All these genes are likely to be considered potential drug targets. In this study, the interactions between mutated genes and drugs were searched from the database. In order to explore the SNPs and CNVs of mutated genes in ovarian cancer, SNPs mutation data (involving 620 samples) and CNVs mutation data (involving 436 samples) of mutated genes in ovarian cancer were downloaded from the Cancer Genome Atlas (TCGA) dataset ( https://gdc.xenahubs.net ). In the SNP data of ovarian cancer, the default parameters of the “maftools” in the R package were used to analyze the gene mutations. In addition, the amplification and deletion frequency of mutated genes were calculated based on CNVs data of ovarian cancer. To investigate additionally the mRNA expression of mutated genes in other reproductive-related tumors, GSE26712 and GSE10615 datasets were downloaded from the Gene Expression Omnibus database (GEO) database. GSE26712 dataset and GSE10615 dataset were used the mRNA expression of mutated genes in ovarian cancer and yolk sac tumor, respectively. For the GSE26712 dataset, the Wilcoxon test was used to analyze the statistical significance of differential expression. For the GSE10615 dataset, differential expression analysis was performed using the “limma” in the R package, and the results were visualized with the “ggplot2” in the R package. The mRNA expression of mutated genes was presented by box plot.

Previous multiregion exome sequencing study focused on somatic mutations and clonal origins in tumor tissues of immature teratoma. 8 However, germline variations from peripheral blood remain unclear. Our study was designed to fill this research gap. In 26 sporadic patients with immature teratoma, nonsynonymous SNP and C>T and G>A was the major variation type and base change types, respectively. A total of 24 common mutated genes were identified in 26 sporadic immature teratoma patients. Among which, 5 common mutated genes with common mutation sites were identified, including 2 frameshift mutant genes (MYPOP and FRG2C) and 3 nonsynonymous mutant genes (CNTNAP3, GPRIN2 and MUC3A). Moreover, molecular weight of these genes changed slightly after mutation. However, it is worth noting that this study is a preliminary exploratory investigation, and the results of this study still require further in-depth research and experimental verification in subsequent work to further confirm their reliability. It is found that MYPOP can regulate hematopoiesis and tumorigenesis. 11 In addition, a differentially methylated DNAm probe site of MYPOP is identified in endometrium. 12 FRG2C is more frequently mutated in estrogen receptor-negative breast cancer. 13 Additionally, mutation of FRG2C (c.A497G; p.E166G) is found in gestational choriocarcinoma. 14 Herein, frameshift mutations of MYPOP (R367Pfs*83) and FRG2C (R161Sfs*5) were found in immature teratoma. Interestingly, MYPOP and FRG2C had a higher deletion frequency of CNV in ovarian cancer. It is indicated that frameshift mutations in MYPOP (R367Pfs*83) and FRG2C (R161Sfs*5) may be associated with the tumorigenesis of immature teratoma. CNTNAP3 encodes the protein related to cell-recognition process within the nervous system. 15 CNTNAP3 is down-regulated in women with adenomyosis. 16 Chaiworapongsa T et al found that CNTNAP3 was up-regulated in late-onset preeclampsia patients 17 In peripheral blood samples of testicular germ cell tumors, loss of 9p13.1-p12 of CNTNAP3 was the only one novel CNV. 18 In this study, the nonsynonymous mutation of CNTNAP3 (T1051S) was identified in immature teratoma. In addition, CNTNAP3 had higher amplification frequency of CNV and had rare mutation frequency (1%) in ovarian cancer. It is suggested that nonsynonymous mutation of CNTNAP3 (T1051S) may be associated with germ cell development in immature teratoma. GPRIN2, encodes a G-protein regulated inducer of neurite overgrowth, is involved in formation and extension of neurite-like processes. 19 , 20 The mutation of GPRIN2 (rs112620425) is found in the breast cancer. 21 It is noted that CNV amplification of GPRIN2 is detected in fetuses with congenital malformations. In the present study, the nonsynonymous mutation of GPRIN2 (L400P) was identified in immature teratoma. GPRIN2 had a higher amplification frequency of CNV in ovarian cancer. In addition, mRNA expression levels of GPRIN2 were down-regulated and up-regulated in ovarian cancer and yolk sac tumor, respectively. It is assumed that GPRIN2 is involved in the embryonic neurodevelopment, the nonsynonymous mutation of GPRIN2 (L400P) may be associated with immature teratoma. Heinzelmann-Schwarz VA et al found that MUC3A was up-regulated in mucinous ovarian carcinoma 22 In this study, nonsynonymous mutations of MUC3A (T2244S, T2000N, A3149E and T2000A) were identified in immature teratoma. Furthermore, MUC3A had a higher amplification frequency of CNV in ovarian cancer. These nonsynonymous mutations of MUC3A may be associated with the process of immature teratoma. In the PPI network, the highest degree gene was MUC12. Furthermore, MUC12 interacted with FCGBP and GALNT12, except MUC3A. In breast cancer, MUC12 is primarily connected with the modulation and function of Janus kinase (JAKs). 23 FCGBP, located in the mucosa of endodermal-derived tissues, plays a crucial role in mucosal immunological defenses. 24 , 25 In ovarian cancer, FCGBP is significantly associated with overall survival. 26 GALNT12 is involved in embryonic development, organ development and organ morphology. 27 In healthy women, GALNT12 is associated with endometrial receptivity. 28 Herein, mutations of MUC12, FCGBP and GALNT12 were found in immature teratoma. Moreover, MUC12 and FCGBP had a higher amplification frequency of CNV, and GALNT12 had a higher deletion frequency of CNV in ovarian cancer. In addition, FCGBP was down-regulated and up-regulated in ovarian cancer and yolk sac tumor. These results indicated that MUC12, FCGBP, GALNT12 and MUC3A may act synergistically during the development of immature teratoma. Based on the drugs-24 mutated genes interaction analysis, only compounds of CHEMBL35482, FLUOROURACIL, and DECITABINE interacted with MUC2. Moreover, mRNA expression levels of MUC2 were up-regulated in both ovarian cancer and yolk sac tumor. MUC2 is expressed in villous goblet cells and immature crypt cells. 29–31 In ovarian cancer, MUC2 is a tumor promoter and is inversely associated with patient survival. 32 It is worth mentioning that MUC2 is expressed in ovarian teratoma. 33 In ovarian cancer, CHEMBL35482 is a potential drug interacted with cyclin-dependent kinase inhibitor 1B (CDKN1B). 34 It is found that the production of MUC2 can be blocked by CHEMBL35482 in pulmonary epithelial cells. 35 It is found that the adapted chemotherapy based on histology, such as FLUOROURACIL-based regimens for adenocarcinoma transformation, has recently been advocated to improve patient outcomes. 36 , 37 A patient with metastatic colonic-type adenocarcinoma arising in a teratoma has achieved response from FLUOROURACIL-based chemotherapy. 38 DECITABINE, a specific DNA methyltransferase inhibitor, can exert antitumor effects by re-expressing tumor suppressor genes. The early phase clinical trials in epithelial ovarian cancer have showed mixed results of the combination of DECITABINE and carboplatinum. 39–41 Good antitumor effect in combination DECITABINE with belinostat has been found in ovarian cancer. 42 It is suggested that MUC2 may be associated with the pathogenesis of immature teratoma and could be potentially modulated by CHEMBL35482, FLUOROURACIL, and DECITABINE in the treatment of immature teratoma. According to the GO analysis of 24 mutated genes, MUC3A, MUC2 and MUC12 were significantly enriched in activation of the innate immune response. The diminished innate immune response can help pluripotent stem cells maintain the rapid rate of cell proliferation, which is essential for the early embryo development. 43 In teratoma, targeting the immune response significantly inhibits the evolution of tumors progression. 44 Thus, it can be seen that the innate immune response may contribute to the pathological mechanism of immature teratoma by regulating the activity of MUC3A, MUC2 and MUC12. Moreover, previous studies have shown that MUC2 and MUC3A may regulate disease progression via the PI3K/Akt/mTOR and NF-kappaB signaling pathways, which are also involved in modulating macrophage-mediated innate immune responses. 45–48 MUC12 exerts an oncogenic role via the c-Jun/TGF-β signaling pathway, and TGF-β is involved in the regulation of innate immunity. 49 , 50 These studies suggest that the aberrations of MUC2, MUC3A and MUC12 identified in the study may participate in the development of immature teratoma by regulating the aforementioned signaling pathways and influencing innate immune responses. This provides a theoretical framework for elucidating the pathological mechanism of the disease, and also points out the potential direction for subsequent related research. In addition, 2 common CNVs (deletion state) were screened out in immature teratoma, which involving 6 genes, such as RP11. The expression of RP11 is decreased in human embryo competence and increased in ovarian germ cell tumor. 51 , 52 Maybe, CNVs deletion of RP11 is associated with the development of immature teratoma. In conclusion, this WES study based on peripheral blood was performed to explore germline variations in 26 patients with sporadic immature teratoma. We identified 24 common mutated genes. Among which, 5 common mutated genes with common mutation sites were identified, including 2 frameshift mutant genes (MYPOP and FRG2C) and 3 nonsynonymous mutant genes (CNTNAP3, GPRIN2 and MUC3A). In the PPI network, the highest degree gene was MUC12. MUC2 may be a potential target for CHEMBL3582, FLUOROURACIL and DECITABINE. Functional analysis showed that MUC3A, MUC2 and MUC12 were enriched in innate immune response. Two common deletion CNVs were also identified. These germline SNPs and CNVs may be associated with the development of immature teratoma and provide potential biomarkers and therapeutic targets. However, there are some limitations in this study. Firstly, the number of enrolled immature teratoma patients is small, larger pedigrees or sporadic cases are needed. Secondly, all the identified SNPs should be validated by Sanger sequencing. Finally, the investigation of the underlying molecular mechanism of SNPs and CNVs is further needed in animal model or cell experiments.