DNA Extraction and Whole-Exome Sequencing

GenRead™ DNA FFPE Kit (Qiagen; Germany) was used to extract genomic DNA from formalin-fixed paraffin-embedded tissues and paired normal blood DNA. Then the genomic DNA was fragmented and hybridized following instructions on manufacturer’s protocol. These two methods were combined to verify the quality of isolated genomic DNA: 1) measuring DNA concentration using Qubit^® DNA Assay Kit in Qubit^® 2.0 Flurometer (Invitrogen, USA); 2) monitoring DNA degradation and contamination using 1% agarose gels. The Agilent SureSelect Clinical Research Exome Version 2 Probe set kit was applied to perform exome capture. Then, the Illumina HiSeq X platform was used to sequence the captured DNA as instructed in the manufacturer’s protocol, with a paired-end run of 2 × 151 bp in Nantong ZhongKe Medical Laboratory. HCS (HiSeq Control Software v2.2) and RTA (Real Time Analysis. v1.18) were used to analyze initial alignment and quality control. Illumina package bcl2fastq (v1.8.4) was applied to convert the binary base call (BCL) to a FASTQ format. Related sequencing data in this article have been uploaded to the GenBank Data Libraries with accession number as PRJNA684940.

Bioinformatics Analysis

Short reads (Raw data) in FASTQ format were obtained from fluorescence images on Hiseq platform by base calling. Fastp was used for quality control by discarding low quality paired reads, including reads that contain adapter contamination or more than 10% of uncertain bases in either one read or over 50% low-quality bases in either one read or less than 20 bases after trimming (18). Mapping results were obtained using Burrows-Wheeler Aligner (BWA) and stored in BAM format (19, 20). PCR duplicates of mapped reads were removed with Sambamba (21).

Germline mutations, including single-nucleotide polymorphisms (SNPs) and small insertions and deletions (InDels) of each exome, were detected using the HaplotypeCaller tool of Genome Analysis Toolkit (GATK) software following the protocol (22). To avoid false positives, these criteria were applied to filter.vcf outputs: SNP: “QD< 2.0 || FS > 60.0 || MQ< 40.0 || MQRankSum< -12.5 || ReadPosRankSum< -8.0”, InDel: “QD< 2.0 || FS > 200.0 || ReadPosRankSum< -20.0”. Annotation of filtered SNPs was conducted using ANNOVAR (23) and snpEFF (24). Only mutations that met the following criteria were considered for further analysis (1): the mutation locates in the exon or its flank region ( ± 2bp) (2); its allele frequency is less than 0.01 in GnomAD, 1,000 genomes and in house database (25, 26) (3); it is not a synonymous mutation (4); more than two softwares of SIFT, Polyphen2HDIV, Polyphen2HVAR, MutationTaster2, predict it to be harmful (27–30) (5). Its host gene is predicted to be associated with PC by phenolyzer software (31).

Mutect2 was used to detect somatic mutations, including single-nucleotide variants (SNVs), small insertions and deletions (InDels) (32). Filter Mutect Calls, a tool of GATK, was employed to assess the qualities of all potential variants and only the site with tag “PASS” in the “FILTER” filed were remained. With Vcf2maf (https://github.com/mskcc/vcf2maf) and VEP, vcf files were converted maf files and annotated somatic mutations (33). Cnvkit was used to identify the copy number variants (CNVs) (34).

The relative weights of mutational signatures of in each sample were calculated according to all somatic SNVs. The R package “deconstructSigs,” built on the basis of mutational signatures in COSMIC (http://cancer.sanger.ac.uk/cosmic/signatures), was applied to quantify the contribution of each signature for each tumor (35). The microsatellite instability (MSI) score was calculated withmsisensor2 (36).

The ABSOLUTE algorithm (37) was employed to assess the tumor purity and ploidy of all samples with the default parameters, which took segmented copy-number and MAF profiles as input. The dNdScv R package (38), which is a group of maximum-likelihood dN/dS methods, was used to quantify selection in cancer and somatic evolution in this study.

Cell Culture and Transfection

BxPC-3 cells were cultured in 1640 medium supplemented with 10% fetal bovine serum (FBS) (Gibco, Thermo Fisher Scientific, Inc, Waltham, MA, USA) and maintained in a humidified 5% CO2 atmosphere at 37˚C. WT (wild type) or Mut (a specific mutation) plasmids of SMAD3 and SMAD4 were purchased from Genechem (Shanghai, China). Briefly, the mutation plasmids were constructed as the mutation information found in the sample (exon2: c.220C>T for SMAD3 gene and exon9: c.1138A>T for SMAD4 gene) by chemical synthesis. To reduce the interference of endogenous expression of SMAD3 or SMAD4, lentiviruses containing SMAD3 or SMAD4 shRNA, purchased from Genechem, were first utilized to knockdown their expression before WT and Mut plasmids transfection in BxPC-3 cells. The siRNAs targeting P21 or FOXO3a were also purchased from Genechem (Shanghai, China).

Cell Viability Assay and LDH Release Assay

BxPC-3 cells transfected with WT or Mut plasmid were seeded in 96-well plates. We then utilized Cell Counting Kit-8 (CCK-8) from MedChemExpress (Madison, WI, USA) and LDH assay kit (Promega, Madison, WI) to examine cell viability and cytotoxicity detection as the manufacturer’s instructions, respectively. The results were normalized against the Mock group from three independent experiments done in triplicate.

Colony Formation and Clonogenic Survival

Cell colony formation detection was performed in a 6-well plate. Cells in different groups were seeded at optimal density. Fresh medium was changed every three days. After two weeks, we fixed the cells with methanol and then stained them with crystal violet. The pictures were photographed by a digital scanner. The surviving fraction (SF) was represented by calculating a ratio of the number of colonies formed and the product of the number of cells plated and the plating efficiency.

Western Blot and Real-Time PCR

Protein extraction, protein determination, Western blot, and real-time PCR were performed as described previously (39).

Sequencing Data Summary

Tumor tissues and paired blood samples of six individuals were collected, respectively to conduct the whole exome sequencing. After removing the low-quality, too short and adaptor containing reads, about 48.5 G bases per library were obtained, and the sequence depth ranged from 330X to 533X. Otherwise, the bases and depth of paired samples were 15.4 G and 149X, averagely. Then, all reads were mapped to the UCSC Homo sapiens genome (version 19), with the total mapping rate more than 99.5% in all samples. Clean data have the Q30 value more than 90% in our research and fits the need of following identification of somatic mutations (details in Supplemental Material S1).

Somatic Mutations in PDAC

A total of 2,380 somatic variations (median, 403.5; range, 74–735; SD, 241.3) were identified in the six pancreatic carcinomas samples, including 1,989 single nucleotide polymorphisms and 391 small insertions and deletions (Figure 1A). As expected, most variations are located in the exonic region (35.03%) and intronic area (44.33%) of the genomes (Figure 1B). We also found that cytosine (C) to thymine (T) transitions were the predominant nucleotide changes (Figure 1C). The somatic mutation profiles were analyzed to address the possible mechanism in the occurring and progressing of PDAC (Figure 1D). Signature 1, age, has the high weight in all six individuals, and signature 3 enriched in insensitive group members (COSMIC Signatures version 2).

Open in a separate window

Figure 1

The somatic mutation profiles of six primary pancreatic tumor samples. (A) Total number of SNVs and InDels detected in each sample. (B) The statistics of somatic mutations in the genome. (C) The distribution of mutation type for each sample. (D) The mutation signature distribution based on COSMIC signature database.

The TMB, MSI, and Clonality Are Not Related to Radiosensitivity in PDAC

It is known to all that tumor mutational burden and microsatellite instability were closely related with the response to immunotherapy in various types of cancers. To determine whether there are also correlations between tumor mutational burden and radiosensitivity in PDAC patients, we examined the number of all mutations in the target region as tumor mutation burden in this tumor type. As shown in Table 2, our results showed that there was no significant correlation between overall survival time (OS) and MSI (p=0.9194). Otherwise, the TMB was negatively correlated with the OS with a coefficient of -0.66 (p=0.175). To determine whether radiosensitivity was correlated with clonality, the purity and copy number alterations were then analyzed in this study. Our results indicated that there was no obvious association between OS and purity or ploidy (r = 0.029424 and r = 0.2, respectively) (Table 2). These results showed that the overall survival is not significantly correlated with the number of missense mutations, the instability of microsatellite, and clonality in PDAC patients treated with radiotherapy. The possibly influencing factors, age and tumor stage, were also analyzed (Table 2). However, the results showed that all these factors are not the core factors causing radioresistance in PDAC.

Differences Between Sensitive and Insensitive Individuals in PDAC

Somatic mutations were filtered according to the criterions mentioned in the method to expose the candidate deleterious sites. 459 in 2301 variants remained, which is located in 441 genes, such as TP53, KRAS, PTEN, and so on (Figure 2A, detail in Supplemental Material S2). We next counted the number of variants for each gene in sensitive and insensitive members. Among them, 434 genes were only mutated in group N or group S, which were defined to be radiosensitivity-related genes. To address the possibly functions that these radiosensitivity-related genes took part in, Gene Ontology and KEGG enrichment analysis were conducted (Figure 2B). The most significant GO terms contained response to extracellular stimulus, mitotic cell cycle checkpoint, DNA catabolic process, and so on (detail in Supplemental Material S3). Pathway terms, like cell cycle, TGF-beta signaling pathway, and FoxO signaling pathway, were significantly enriched in our study, which might be closely related with radiosensitivity in PDAC (detail in Supplemental Material S4). As shown in Figures 2C, E, no significant differences on total number of genes and variation were observed between group N and group S. There are eight common genes in both groups, as well as 258 and 175 unique genes belonged to the group N and group S, respectively (Figure 2D). Moreover, two common mutations were shared in both groups, but 270 and 187 unique mutations belonged to the group N and group S, respectively (Figure 2F). Germline mutations were also analyzed in our study. All the mutations were filtered by the procedures mentioned in the method, and 754 sites in 639 genes remained (Details in Supplemental Material S5).

Open in a separate window

Figure 2

Specific mutated genes between group N and group S. (A) The oncoprint plot of all somatic mutations. (B) The bubble chart of kegg terms enriched for all diff-mutated genes. (C) No significant difference on total number of genes was observed between group N group and group S. (D) There are 8 common genes in both groups, as well as 258 and 175 unique genes belonged to the group N and group S, respectively. (E) No significant difference on total number of variations was observed between group N and group S. (F) 2 common mutations were shared in both groups, but 270 and 187 unique mutations belonged to the group N and group S, respectively.

Pathway-Based Protein Interaction Network of Radiosensitivity-Related Genes

Based on previous studies, 30 RI-related pathways were selected, including cell cycle, homologous recombination, nucleotide excision repair, and so on. Combined with the protein-protein interactions, a pathway-based protein-protein interaction (PPI) network was constructed with the RI-related genes (Table 3). Seventeen genes were involved. SMAD3, SMAD4, ARID1B, CRKL, and CASP10 played a center role in the network (Figure 3A). Moreover, the results of qPCR showed that the mRNA expression of these five genes was not affected by ionizing radiation in BxPC-3 cells (Figure 3B). It is worth noting that both SMAD3 and SMAD4 belong to the SMAD proteins family, which is a group of intracellular signal transducer proteins. The synonymous mutations of somatic mutations are listed in Supplemental Material S6.

Open in a separate window

Figure 3

(A) Protein-protein interaction network of RI-related genes. (B) The mRNA expression of SMAD3, SMAD4, ARID1B, CRKL, and CASP10 was not affected by ionizing radiation (4Gy) in PDAC cells.

Mutations of SMAD3 and SMAD4 Affect the Radiosensitivity in PDAC

The ectopic wild-type (WT; without the variation) and mutant (Mut; with the variation) SMAD3 or SMAD4 cells were constructed in BxPC-3 cells (a PDAC cell line) to determine whether the mutation of SMAD3 or SMAD4 affected the radiosensitivity. As shown in Figures 4A, B, results of qPCR showed that there were no obvious differences on the expression of SMAD3 or SMAD4 between the WT and Mut cells, respectively. We then examined cell survival and death after ionizing radiation (IR) by utilizing CCK-8 and LDH release assay. The results showed that transfection of WT SMAD3 repressed the decreases of cell viability and the increase of LDH release induced by ionizing radiation (IR, 4Gy), whereas mutation of SMAD3 facilitated the effects of IR on cell survival and death (Figures 4C, D). Moreover, transfection of SMAD3 mutation significantly promoted IR-induced decreases of colony formation and surviving fraction in PDAC cells compared with the WT group (Figures 4E, F ). On the contrary, the effects of IR on cell survival, cell death, colony formation, and surviving fraction were inhibited by the mutation of SMAD4 compared with the WT group (Figures 4G–J). These results indicate that SMAD3 mutation or SMAD4 mutation leads to the different response to radiotherapy in PDAC.

Open in a separate window

Figure 4

The radiosensitivity is affected by the mutation of SMAD3 or SMAD4 in PDAC. (A, B) Results of qPCR showed that there were no obvious differences on the expression of SMAD3 (A) or SMAD4 (B) between the WT and Mut cells, respectively. (C, D) The regulatory effects of IR on cell viability (C) and LDH release (D) were facilitated by the SMAD3 mutation. (E, F) The decreases of colony formation (E) and surviving fraction (F) induced by IR were enhanced by the mutation of SMAD3 compared with the WT group. (G, H) SMAD4 mutation mitigated the effects of IR on cell survival (G) and LDH release (H). (I, J) The decrease of colony formation (I) and surviving fraction (J) induced by IR were antagonized by SMAD4 mutation compared with the WT group. Cells infected with the empty vector were defined as the mock group. *P < 0.05.