Identification of therapeutic targets in chordoma through comprehensive genomic and transcriptomic analyses.

Chordoma is a rare, orphan cancer arising from embryonal precursors of bone. Surgery and radiotherapy (RT) provide excellent local control, often at the price of significant morbidity because of the structures involved and the need for relatively high doses of RT; however, recurrence remains high. Although our understanding of the genetic changes that occur in chordoma is evolving rapidly, this knowledge has yet to translate into treatments. We performed comprehensive DNA (paired tumor/normal whole-exome and shallow whole-genome) and RNA (tumor whole-transcriptome) next-generation sequencing analyses of archival sacral and clivus chordoma specimens. Incorporation of transcriptomic data enabled the identification of gene overexpression and expressed DNA alterations, thus providing additional support for potential therapeutic targets. In three patients, we identified alterations that may be amenable to off-label FDA-approved treatments for other tumor types. These alterations include FGFR1 overexpression (ponatinib, pazopanib) and copy-number duplication of CDK4 (palbociclib) and ERBB3 (gefitinib). In a third patient, germline DNA demonstrated predicted pathogenic changes in CHEK2 and ATM, which may have predisposed the patient to developing chordoma at a young age and may also be associated with potential sensitivity to PARP inhibitors because of homologous recombination repair deficiency. Last, in the fourth patient, a missense mutation in IGF1R was identified, suggesting potential activity for investigational anti-IGF1R strategies. Our findings demonstrate that chordoma patients present with aberrations in overlapping pathways. These results provide support for targeting the IGF1R/FGFR/EGFR and CDK4/6 pathways as treatment strategies for chordoma patients. This study underscores the value of comprehensive genomic and transcriptomic analysis in the development of rational, individualized treatment plans for chordoma.

INTRODUCTION

Chordoma is a rare, locally aggressive, and invasive cancer arising from the notochord, which is a cartilage-like mesodermal rod. In vertebrates, the notochord is a precursor to the intervertebral disc. Chordoma is considered a malignant tumor of bone and accounts for 1%–4% of such tumors. The incidence is <1/1,000,000 per year and approximately 22 cases are diagnosed per year in the United States.

Chordoma may arise anywhere along the cranio-spinal axis. In adults, 50% of the cases involve the sacro-coccygeal region, 35% occur at the base of the skull, and 15% are found elsewhere in the vertebral column (Chugh et al. 2007). Gross total resection is typically limited by the cancer's proximity to critical neurovascular structures. Thus, adjuvant radiotherapy is used routinely. Patients have a median survival of 7.7 yr (Smoll et al. 2013); the survival rate is ∼70%–72% at 5 yr and falls to 40%–48% at 10 yr (Smoll et al. 2013; Chambers et al. 2014).

The brachyury (T, T-box transcription factor T) protein is essential for notochord differentiation. Expression of the T gene is used to differentiate chordoma from other chondroid tumors (Vujovic et al. 2006; Pillay et al. 2012). A germline single-nucleotide polymorphism (SNP) within the DNA-binding region of brachyury has been strongly associated with the risk of developing sporadic chordoma, whereas duplication of the T gene has been reported in the rare familial form of chordoma (Yang et al. 2009).

Previous analyses have led to the identification of pathological changes in chordomas. Such changes include transcriptomic changes (Bell et al. 2016), as well as point mutations and copy-number losses in genes such as ALK, CDKN2A, NRAS, PTEN, NBPF1, SETD2, and SMARCB1 (Hallor et al. 2008; Le et al. 2011; Choy et al. 2014; Fischer et al. 2015; Wang et al. 2016; Sa et al. 2017; Cote et al. 2018), whose loss characterizes poorly differentiated chordomas (Mobley et al. 2010). Of note, loss of INI1/SMARCB1 may also represent a discrete entity with a more aggressive phenotype more similar to rhabdoid tumors (Hasselblatt et al. 2016). As a rare and typically fatal condition with a complex etiology, and given the variety of molecular mechanisms that may be responsible for the development of this tumor, chordoma may be a candidate for individualized approaches to molecularly targeted treatment. Thus, in this study, we performed next-generation DNA and RNA sequencing of paired chordoma and germline samples to build on our understanding of the genomic and transcriptomic landscape of chordoma and to identify potentially actionable alterations that may be worth exploring further in clinical trials.

RESULTS

Across the four chordoma patients (Table 1), a mean target coverage of 165× was achieved for the tumors and 112× for the germline exomes. From the long-insert whole genomes (LIWGs), a mean clonal coverage of 43× was achieved for the tumors and 44× for the germline samples. A mean of 195 million paired reads was generated across the tumor RNAs using RNA-seq (RNA sequencing). The total number of somatic alterations identified in each patient's tumor (loss of function, missense, intronic, in-frame indel, splice site loss, splicing altered, untranslated region [UTR]) were 24, 32, 2, and 22, respectively for Patients 1–4. The brachyury (T) SNP rs2305089 was observed in all the tumors (Pillay et al. 2012) along with elevated expression of T in Patients 1, 2, and 4 with FPKMs of 155.5, 398.2, and 442.9, respectively. Table 2 lists the DNA or RNA alterations that were identified through sequencing and that are potentially targetable with an FDA-approved agent. For those genes associated with a drug–gene rule, any additional alterations that were detected in those genes across the cohort are also listed. Additional somatic alterations implicated in chordoma or other cancers (ALK, CDKN2A, NRAS, PTEN, SMARCB1, NBPF1, TP3) and potential pathogenic germline variants were also evaluated. Herein, we present clinical information, pathology results, and representative imaging data for each patient, as well as the results from sequencing and drug–gene analyses.

Table 1.

Summary of patients

Patient	Gender	Age at diagnosis	Ethnicity	Location	Tumor purity (%)a
1	M	70	Caucasian	Sacral	36
2	M	57	Caucasian	Sacral	74
3	M	29	Caucasian	Sacral	33
4	M	54	Caucasian	Petro-clival	44

Phenotypic information for each patient is shown.

aTumor purities were estimated in silico using Sequenza (Favero et al. 2015).

Table 2.

Selected aberrations and associated drug matches

Patient	Gene	Gene name	Variant typea	Chromosomal location	Alteration	Description	HGVS cDNA	Tumor DNA allele frequencyb	Constitutional DNA allele frequencyb	HGVS protein	ClinVar variation ID	ClinVar clinical significance	Predicted effect (MutationTaster/PolyPhen-2/FATHMM)c	Genotype (HM/HT)	Drug matches
1	FGFR1	Fibroblast growth factor receptor 1	Overexpression	Chr 8: 38,268,655	log2 fold = 3.3 (Qd = 0.03)	N/A	N/A	N/A	N/A	N/A	N/A	N/A	N/A	N/A	Ponatinib, pazopanib
	PIP4K2A	Phosphatidylinositol-5-phosphate 4-kinase type 2α	Substitution	Chr 10: 22,825,649	A>C	UTR	c.*481T>G	0.16 (5/31)	0.0 (0/19)	N/A	N/A	Not reported	N/A	HT	Palbociclib
	T	Brachyury	Substitution	Chr 6: 166,579,270	G>A	G177D; rs2305089	c.530G>A	0.48 (55/114)	0.54 (43/79)	p.Gly177Asp	N/A	Not reported	P/B/T	HT	N/A
2	CDK4	Cyclin-dependent kinase 4	Copy gain	Chr 12: 55,640,800	log2 fold = 1.5	N/A	N/A	N/A	N/A	N/A	441984	Pathogenic	N/A	N/A	Palbociclib, abemaciclib
	ERBB3	erb-b2 receptor tyrosine kinase 3	Copy gain	Chr 12: 55,640,800	log2 fold = 1.5	N/A	N/A	N/A	N/A	N/A	441984	Pathogenic	N/A	N/A	Afatinib, gefitinib, lapatinib
	GLI1	GLI family zinc finger 1	Copy gain	Chr 12: 55,640,800	log2 fold = 1.5	N/A	N/A	N/A	N/A	N/A	441984	Pathogenic	N/A	N/A	Arsenic trioxide
	GLI1	GLI family zinc finger 1	Substitution	Chr 12: 57,864,709	C>T	T729I	c.2186C>T	0.09 (31/359)	0.0 (0/101)	p.Thr729Ile	N/A	Not reported	N/B/T	HTe	N/A
	GLS2	Glutaminase 2	Copy gain	Chr 12: 55,640,800	log2 fold = 1.5	N/A	N/A	N/A	N/A	N/A	441984	Pathogenic	N/A	N/A	Temozolomide
	PIP4K2A	Phosphatidylinositol-5-phosphate 4-kinase type 2 alpha	Substitution	Chr 10: 22,898,627	A>T	V55D	c.164T>A	0.78 (25/32)	0.0 (0/54)	p.Val55Asp	N/A	Not reported	D/D/T	HM	Palbociclib
	T	Brachyury	Substitution	Chr 6: 166,579,270	G>A	G177D; rs2305089	c.530G>A	1.0 (109/109)	1.0 (75/76)	p.Gly177Asp	N/A	Not reported	P/B/T	HM	N/A
	TP53	Tumor protein p53	Copy loss	Chr 17: 5,003,600	log2 fold = −0.7	N/A	N/A	N/A	N/A	N/A	149158	Pathogenic	N/A	N/A	N/A
	IGF1R	Insulin-like growth factor 1 receptor	Substitution	Chr 15: 99,502,153	T>G	UTR	c.*1482T>G	0.28 (34/122)	0.00 (0/83)	N/A	N/A	Not reported	N/A	HT	Ganitumabf
3	ATM	ATM serine/threonine kinase	Germline substitution	Chr 11: 108,119,823	T>C	V410A; rs56128736	c.1229T>C	0.49 (19/39)	0.41 (23/56)	p.Val410Ala	127332	Conflicting (benign; likely benign; uncertain)	D/B/T	HT	Olaparib, niraparib
	CHEK2	Checkpoint kinase 2	Germline substitution	Chr 22: 29,091,178	G>T	D481Y; rs200050883	c.1441G>T	0.56 (79/141)	0.62 (66/107)	p.Asp481Tyr	128056	Conflicting (likely benign; uncertain)	D/D/T	HT	Olaparib, niraparib
	T	Brachyury	Substitution	Chr 6: 166,579,270	G>A	G177D; rs2305089	c.530G>A	1.0 (112/112)	1.0 (87/87)	p.Gly177Asp	N/A	Not reported	P/B/T	HM	N/A
4	IGF1R	Insulin-like growth factor 1 receptor	Substitution	Chr 15: 99,460,066	T>G	V721G	c.2162T>G	0.04 (9/206)	0.0 (0/142)	p.Val721Gly	N/A	Not reported	D/P/T	N/Ae	Ganitumabf
	PBRM1	Polybromo 1	Substitution	Chr 3: 52,668,696	A>G	F408S	c.1223T>C	0.24 (19/78)	0.0 (0/91)	p.Phe408Ser	N/A	Not reported	D/D/T	HT	N/A
	T	Brachyury	Substitution	Chr 6: 166,579,270	G>A	G177D; rs2305089	c.530G>A	0.59 (60/102)	0.37 (30/81)	p.Gly177Asp	N/A	Not reported	P/B/T	HT	N/A

aVariants listed are somatic unless indicated as germline.

bAllele frequencies are shown as the ratio of the # alternate reads/# total reads.

cPredicted effect: predicted functional impact based on MutationTaster, PolyPhen-2, and FATHMM. D, Damaging; P, possibly damaging; B, benign; N, neutral; T, tolerated.

dQ-value, corrected P-value using the Benjamini–Hochberg correction for multiple testing.

eSubclonal event identified using exome sequencing.

fInvestigational agent.

Sacral Chordomas

Patient 1

A 70-yr-old male patient presented with a 14 × 11 × 20-cm predominantly presacral mass and bony destruction of the distal sacrum and coccyx, as well as multiple amorphous calcifications. Surgery consisted of an S1–S3 laminectomy with en bloc resection of S3 and removal of the distal sacrum and coccyx (Fig. 1, upper row). Pathology reported the tumor to be chondromyxoid, with mild atypia and no mitoses or necrosis, as well as positive immunohistochemistry (IHC) for S100 and keratin AE1/AE3. Eleven months after diagnosis, a magnetic resonance imaging (MRI) showed recurrent local disease (Fig. 1, lower row), with pelvic and inguinal lymphadenopathy and pulmonary metastases. He began palliative treatment with intensity-modulated radiotherapy (IMRT) at that time with 15 × 2.5 = 37.5 Gy to the entire pelvis and an additional 5 × 2.5 = 12.5 Gy to the tumor, along with imatinib (400 mg daily). He developed intercurrent pneumonia and was transferred to hospice care.

Figure 1.

Patient 1 MRI images. All images are sagittal in midline. Panels in the left column (A,C) are T2-weighted, fast-spin echo, whereas those in the right column (B,D) are T1-weighted with contrast (T1 + Gd). Panels shown in the upper row (A,B) are postoperative, following initial gross total resection. Resection of the vertebrae caudal to S1 is indicated with the blue arrow. Panels in the lower row (C,D) were taken at time of recurrence. A large, heterogeneous mass extending rostrally almost to S1 and dorsally almost to the skin is present (blue arrows). The dimensions are almost equal to those of the tumor prior to initial resection.

Tumor DNA sequencing revealed no definitive therapeutic targets. However, a somatic point mutation in the UTR of PIP4K2A was identified and may thus impact transcriptional regulation of the gene. PIP4K2A has been reported to be essential for cell proliferation in acute myeloid leukemia, through AKT/mTOR activation and cell cycle regulation (Jude et al. 2015), and has also been found to be a target of the CDK4/6 inhibitor palbociclib (Sumi et al. 2015). Although PIP4K2A expression was not significantly up-regulated compared to a universal human control, expression of the gene was observed, with a FPKM of 17.8. Clinical trials are currently underway evaluating palbociclib as a single agent in advanced chordoma (NCT03110744), and in CDK4-amplified solid tumors (NCI Match, NCT02465060).

Overexpression of FGFR1 compared to a universal RNA control was also observed in Patient 1's tumor. The FGF pathway has been implicated in chordoma as a driver of brachyury overexpression (Nelson et al. 2012; Hu et al. 2014), which was observed in this patient's tumor. In preclinical models, chordoma cell lines express both FGFR2 and FGFR3 and are sensitive to FGFR inhibition (Hu et al. 2014). FDA-approved agents that target FGFRs include ponatinib, a pan-FGFR inhibitor, and pazopanib, a multityrosine-kinase inhibitor, which also has activity against VEGFR and PDGFR (Schutz et al. 2011; Gozgit et al. 2012). Several other FGFR inhibitors are also in clinical development.

With the availability of RNA-seq data, we further evaluated the expression of factors downstream from FGFR1 to determine if evidence for pathway activation is present. Genes downstream from FGFR demonstrated expression (FPKM > 1) with a mean FPKM of 42.7 (see Methods), but none demonstrated differential expression compared to a universal RNA control (corrected P < 0.05). Eight of 22 genes demonstrated higher levels of expression with FPKMs of >42.7 and ranging from 46.9 to 147.0. Although these data are not confirmatory of FGFR pathway activation, such data provide evidence that it is possible that pathway activation may have been present in this patient's tumor.

Patient 2

A 57-yr-old male presented with a 15 × 14 × 6-cm sacral mass. The two-stage surgery included a sacrectomy with S1 and S2 interspace embolization followed by closure of a myocutaneous flap. His postoperative MRI is shown in Figure 2. Pathology showed a tumor arranged in cords and lobules within a myxoid background and no mitoses or necrosis. A focus of lymphovascular invasion was identified. Four years after diagnosis, recurrence was noted along the borders of the resection cavity. He enrolled in a Phase 2 clinical trial (NCT02383498) of IMRT (70 Gy) with placebo versus a brachyury vaccine (GI-6301). Three years later, he remained free of further progression and continued to use the vaccine every 3 mo, according to the trial protocol.

Figure 2.

Patient 2 MRI images. Images for Patient 2 after his initial resection. (A) T2-weighted, with fat saturation via short-tau inversion recovery (STIR). (B) T1 + Gd. As with Patient 1 (Fig. 1A,B), resection of the vertebrae caudal to S1 is seen (blue arrow).

Sequencing revealed somatic DNA alterations clustered within regions of Chromosome 12 that are associated with targets for FDA-approved drugs, such as CDK4 and ERBB3. For instance, CDK4 amplification commonly occurs in liposarcomas, and the CDK4/6 inhibitor palbociclib showed activity in a Phase 2 trial for advanced liposarcoma (Dickson et al. 2013). In chordoma cell lines lacking p16, likely as a result of CDKN2A loss, CDK4/6 gene expression was observed to suggest the possibility of activation in this pathway (von Witzleben et al. 2015). Subsequent treatment of these lines with the CDK4/6 inhibitors palbociclib and abemaciclib inhibited cell growth to demonstrate the potential efficacy of this treatment strategy. In addition, overexpression of ERBB3 correlated with gefitinib sensitivity in a preclinical lung cancer model (von Witzleben et al. 2015) to suggest that ERBB inhibitors or monoclonal antibodies may have activity in this context. Although the expression levels of CDK4 and ERBB3 were not significantly altered compared to a universal RNA reference (P > 0.05), both genes demonstrated expression with FPKMs of 80.5 and 76.3, respectively. Furthermore, although variants associated with CDKN2A loss were not observed in the DNA, CDKN2A expression was low, with a FPKM of 3.5.

Sequencing also revealed alterations in GLI1 and GLS2, which encode candidate drug targets. For instance, arsenic trioxide binds GLI1 to inhibit transcription and cell growth in a cell culture model of Ewing sarcoma (Beauchamp et al. 2011). In Patient 2's tumor, a missense mutation in GLI1 was detected at a low DNA allele frequency (0.08), but there were no changes in the RNA expression. Furthermore, GLS2 overexpression is associated with decreased MGMT expression and increased sensitivity to temozolomide in a cell culture model of glioblastoma (Szeliga et al. 2012). In Patient 2's tumor, the GLS2 copy number was increased but there were no changes in gene expression. Unlike CDK4 and ERBB3, GLI1 and GLS2 demonstrated low expression, with FPKMs of 1.6 and 2.7, respectively.

Other potentially pathogenic somatic abnormalities were observed, although these do not suggest a specific treatment strategy. An expressed missense mutation in PIP4K2A was detected (RNA allele frequency [RAF] = 0.28), along with low-level loss of TP53 (log2 = −0.7). This specific PIP4K2A mutation has not been previously reported or characterized, and the functional consequence is unclear. However, as discussed above, PIP4K2A has been identified as a potential target of palbociclib, raising interest in evaluating the CDK4/6 inhibitor in this setting. A heterozygous IGF1R UTR point mutation was also observed. IGF1R has been implicated in chordoma, as ∼40% of chordomas demonstrate IGF1R/IR signaling activation and IGF1R is correlated with chordoma tumor volume (Sommer et al. 2010; Scheipl et al. 2012). Additionally, UTR variants have been proposed to regulate IGF1R expression. In this tumor, elevated expression of IGF1R was observed with an FPKM of 45.7, although the gene was not significantly overexpressed compared to a universal RNA reference (P > 0.05). The investigational IGF1R antibody, ganitumab, is currently being evaluated in clinical trials for other tumor types.

Patient 3

An incidental presacral mass was identified in a 29-yr-old male who had experienced chronic tailbone pain since age 18. The 4.5 × 8.5 × 4.5-cm tumor was removed en bloc (Fig. 3), and the chronic pain resolved. Pathology showed physaliferous cells arranged in cords, sheets, and nests, forming cyst-like masses separated by fibrocollagenous septa, within a myxoid background. Mild to moderate atypia was noted. IHC for S100, keratin AE1/AE3, and EMA confirmed the diagnosis. Focal invasion of the resection margins was noted. Proton beam radiotherapy was completed, 1.8 × 42 = 75.6 Gy, and 8 yr following resection, there were no signs of progression clinically or on MRI.

Figure 3.

Patient 3's MRI images. Postoperative images for Patient 3. (A) T2-weighted, fat saturated (STIR). (B) T1 + Gd. These images are largely normal; signs of resection of the presacral mass are present (blue arrows).

No somatic DNA or RNA alterations were found to be associated with the targets for FDA-approved drugs. However, evaluation of germline data led to the identification of rare missense mutations in CHEK2 (D481Y; population allele frequency [PAF] = 0.001) and ATM (V410A; PAF = 0.004). These genes encode proteins involved in DNA damage response and thus may have predisposed the patient to cancer. The same germline ATM mutation has been described in breast and ovarian cancer, as well as in acute lymphoblastic leukemia (Thorstenson et al. 2003; Liberzon et al. 2004; Broeks et al. 2008). The CHEK2 mutation we identified is located in the protein kinase domain. Although this mutation has not been reported, similar CHEK2 germline mutations have been reported in Li–Fraumeni syndrome (Bell et al. 1999) and breast cancer (Allinen et al. 2001; Domagala et al. 2012). Both mutations were found in the tumor, although loss of the other allele was not detected in the tumor for either gene. The overall expression of CHEK2 was low in the tumor (FPKM = 3.6), whereas ATM expression was elevated (FPKM = 52.4). Clinical trials are underway evaluating PARP inhibitors (olaparib, niraparib, talazoparib) in tumors with homologous recombination deficiency, including tumors with ATM or CHEK2 loss (e.g., TAPUR, NCT02693535; NCT02286687; NCT03375307). Inhibition of ATM is also being explored and such treatment has been shown to radiosensitize glioma cells (Golding et al. 2009). The patient's younger age supports the hypothesis that these mutations may be pathogenic (Chambers et al. 2014).

Chordoma of Clivus

Patient 4

A 54-yr-old male was diagnosed with a petro-clival mass and significant brainstem compression. Three-stage resection (Fig. 4, after Stage 1) was followed by adjuvant Cyber Knife radiosurgery to remove residual tumor. Pathology showed sheets of cells (physaliferous or epithelioid) within a chondromyxoid background. As with Patient 3, IHC for S100, keratin AE1/AE3, and EMA was confirmatory. The MIB-1 labeling index was 1%, and no mitoses were seen. The patient remains free of progression 18 mo after initial diagnosis.

Figure 4.

Patient 4's MRI images. Images recorded following the first stage of resection for Patient 4. Both panels shown are axial. (A) Tw-weighted (FLAIR, fluid-attenuated inversion recovery); (B) T1 + Gd; (C) sagittal T1 + Gd. Although the majority of the tumor has been resected (blue arrows), a significant component remains (orange arrow), arising from the clivus.

Exome sequencing revealed a subclonal missense IGF1R V721G mutation that was identified in 4% of reads. Expression of the mutation-containing transcript was also observed, with 11% of RNA reads supporting this event. A case study reported that an adult spinal chordoma patient who was treated with an investigational IGF1R/IR inhibitor linsitinib, in addition to the EGFR inhibitor erlotinib, remained stable on treatment for 5 yr (Aleksic et al. 2016). In a separate evaluation of dose escalation of linsitinib with erlotinib in solid tumors, a spinal chordoma patient also demonstrated a partial response to this combination (Macaulay et al. 2016). In Patient 4's tumor, both EGF and EGFR were expressed (FPKM = 50.9 and 70.0, respectively). Anti-IGF1R strategies continue to be explored in clinical trials.

A missense mutation in PBRM1 (F408S) was also expressed in the tumor (DNA allele frequency = 0.24, RAF = 0.40). PBRM1 encodes a component of complexes required for transcriptional activation via nuclear hormone receptors. Mutations in this gene are associated with renal cell carcinoma (Wang et al. 2017; de Velasco et al. 2018; Ricketts et al. 2018). In a targeted sequencing study, both deletion and missense mutation (S1315F) of PBRM1 have been reported in chordomas (Wang et al. 2016).

Pathway Analysis

Pathway analysis on differentially expressed genes (corrected P < 0.05) was additionally performed for each patient's tumor (Supplemental Table S1; Patient 1 = 482 genes, Patient 2 = 615 genes, Patient 3 = 561 genes, Patient 4 = 482 genes). Although analysis did not yield significant (P < 0.05) pathways that are specific to drug–gene matches for each patient, potentially interesting alterations were observed. In Patients 1 and 2, a potential immune component was revealed with the identification of T- and B-cell receptor pathways, and from a cohort-level analysis of all tumors, we observed additional commonly impacted processes across all tumors. These processes include ILK (integrin-linked kinase), integrin, and GTPase signaling, which are associated with the PI3K/Akt/mTor pathways; phospholipase C and calcium signaling pathways, which are associated with Fgf signaling; regulation of gene expression through the nuclear receptors FXR (farnesoid X receptor), TR (thyroid hormone receptor), LXR (liver X receptor), and RXR (retinoid X receptor); and what appear to be tumor microenvironment processes as defined by epithelial adherens junction signaling and tight junction signaling. Although ILK signaling has not been specifically reported in chordomas, obstructing β1 integrin function was suggested to be therapeutically relevant following radiation in chordomas (Harryman et al. 2017). Separately, the importance of the microenvironment of chordomas has been previously demonstrated through in vitro analyses, whereby hypoxia combined with exposure to CCN2, a connective tissue growth factor, up-regulated notochord progenitor markers (Patel et al. 2014). Although limited conclusions may be drawn from expression data alone from this small cohort, these findings shed additional light into the intricacies of chordomas.

DISCUSSION

Our work demonstrates that integrated and comprehensive DNA and RNA profiling is practical on archival chordoma samples and that such an approach aids in the identification of potential candidate therapeutic targets. Previous genetic studies of chordoma have reported somatic abnormalities in genes such as CDKN2A, PTEN, PIK3CA, and TP53 (Hallor et al. 2008; Choy et al. 2014; Fischer et al. 2015; Wang et al. 2016). Similar to previous studies, we did not observe highly recurrent somatic alterations; however, we identified alterations in genes that have been previously associated with chordoma, including PBRM1 (Wang et al. 2016) and IGF1R (Sommer et al. 2010; Scheipl et al. 2012). Incorporation of RNA-seq data was also beneficial, as it clarified which DNA alterations are expressed, and thus may serve as therapeutic targets. Based on comprehensive analyses, potential therapeutic strategies were identified across all patients. In the third patient, we identified germline base substitutions in CHEK2 and ATM, which are likely key driving events in this tumor. With the availability of RNA-seq data, we additionally performed pathway analyses and identified processes that may be associated with potentially targetable pathways, although further studies are needed to determine if and how such processes may impact the efficacy of different treatments.

Overall, the data from this study, along with the results of previously published chordoma studies, provide evidence of aberrations converging on overlapping pathways. Although chordoma appears to be a genetically heterogeneous tumor, activation of common pathways may represent key exclusive events in this disease. These include the Ras-Raf-Mek and mTor pathways that lie downstream from IGF1R, EGFR, and FGFR signaling (Patients 1, 2, and 4), as well as cell cycle progression pathways that lie downstream from TP53, which is activated by CHEK2 and ATM (Patient 3) and CDK4/6 (Patient 2). Inhibition of IGF1R, EGFR, FGFR, and CDK4/6 thus represents possible treatment themes for chordoma (Presneau et al. 2009; Sommer et al. 2010). Targeting of the mTor pathway has also been suggested as a potential therapeutic strategy for chordoma (Presneau et al. 2009).

Although we did not detect EGFR amplification or overexpression in the tumors analyzed here, targeting EGFR using erlotinib has shown promising results in two chordoma patients (Singhal et al. 2009; Houessinon et al. 2015), one of whom was refractory to imatinib. Imatinib inhibited tumor growth in a patient-derived chordoma xenograft (Siu et al. 2013), and in a separate Phase 2 study, imatinib inhibited tumor growth in previously progressive patients expressing PDGFB/PDGFRB (Stacchiotti et al. 2012; Hindi et al. 2015); both of these markers were expressed in the patients analyzed here but did not demonstrate significant overexpression in the tumors compared to a universal RNA control. Notably, the EGFR inhibitor afatinib is currently being evaluated in a Phase 2 trial for chordoma patients with EGFR-expressing tumors (NCT03083678) and also appears to be active in the central nervous system in non-small-cell lung carcinoma patients (Kelly et al. 2018). Last, for Patient 3, who appears to have germline disease, potential strategies for consideration include PARP inhibition and ATM inhibition combined with radiation, though ATM inhibitors are not yet FDA-approved.

Although our study sheds further light onto chordoma and potential treatment options, a number of limitations are apparent. The first is our small sample size, which is directly related to the rare nature of the malignancy and the availability of archival specimens. Second, given the diversity of mutations previously reported in chordoma and those reported here, additional analyses are needed to understand the role of identified mutations in altering functional pathways and how such perturbed pathways enable tumor development and progression in chordoma. Additional studies are also needed to evaluate efficacy and prioritization of candidate single-agent strategies and to explore potential use of combination therapies in this disease.

Even so, despite the rarity of chordoma, previous studies and our analysis here shed valuable light on our understanding of the disease. As the majority of previous chordoma sequencing studies have used targeted sequencing approaches that lack integrated transcriptomic analyses, we present analysis from a unique data set that integrates whole-genome, whole-exome, and whole-transcriptome data across patients. This analysis is thus the first to integrate both DNA and RNA sequencing data in sacral and clivus chordomas and also provides evidence of convergence of alterations in specific pathways. Notably, current clinical trials that are open to chordoma patients are focused on targeting brachyury (NCT02383498), PD-1/PD-L1 (NCT03173950, NCT02989636), SMARCB1/INI-1 (NCT02601950), and CDK4/6 (NCT03110744). Our analyses lend further credence to CDK4/6 inhibition and suggest that targeting IGF1R/EGFR/FGFR signaling may be a new candidate therapeutic strategy. Such convergence and continued molecular characterization of chordomas has the potential to improve our ability to find effective treatment strategies for patients.

METHODS

Clinical Details

De-identified clinical histories, pathology reports, and MRI results were reviewed for each patient. For patients who received care at another institution, some of the records were incomplete.

Sample Collection and Extraction

Chordomas, embedded in optimal cutting temperature compound, and matched whole blood were retrieved from the Biobank Core Facility of the Barrow Neurological Institute at St. Joseph's Hospital and Medical Center. The tissue was homogenized using the Bullet Blender (Next Advance, Inc.) and stainless-steel beads. The lysate was further homogenized using the QIAshredder (QIAGEN) and digested with proteinase K. For nucleic acid purification, we used the AllPrep DNA/RNA Mini Kit (QIAGEN), which allows for the simultaneous purification of genomic DNA and total RNA from the same biological sample. Germline DNA was extracted from whole blood using the Gentra Puregene Blood Kit, 3rd edition (QIAGEN), following the manufacturer's instructions, including RNase A treatment.

Next-Generation Sequencing

Paired tumor/normal exomes were constructed using SureSelectXT Human All Exon V5+UTR baits (Agilent). The KAPA Library Preparation Kit (Kapa Biosystems) was used for Patient 1 and the KAPA HyperPrep Kit (Kapa Biosystems) was used for Patients 2–4. Paired tumor/normal LIWGs were constructed using the KAPA HyperPrep Kit, as previously described (Liang et al. 2014, 2017). Tumor RNA was used to generate RNA libraries using the TruSeq RNA Access Library Prep Kit (Illumina). The libraries were paired-end sequenced using the TruSeq PE Cluster Kit v3 and the TruSeq SBS Kit v3 (Illumina) on the HiSeq 2500 (Illumina) for 2 × 109 bp reads with a 7-base index read.

Data Analysis

Binary base call files were converted to FASTQ files using bcl2fastq Conversion Software (Illumina) v1.8-2.17.1.14. The FASTQ files were aligned to the human reference genome (GRCh37) using the Burrows–Wheeler Aligner (Li and Durbin 2009). Somatic analyses were performed as previously described using an annotation pipeline built on GENCODE 3 (Ensembl; GRCh37.1) (Liang et al. 2017; Zerbino et al. 2018). To determine whether tumor single-nucleotide variants (SNVs) were present, we used the consensus of at least two of the three following tools: Seurat, a Bayesian framework (Christoforides et al. 2013) (threshold quality score of >30), MuTect, a Bayesian classifier (Cibulskis et al. 2013), and Strelka, a somatic small-variant tool (Saunders et al. 2012). Tumor cellularities, or purities, for each sample were estimated in silico using Sequenza (Favero et al. 2015). MutationTaster (Schwarz et al. 2010), PolyPhen-2 (Adzhubei et al. 2010), and FATHMM (Shihab et al. 2013) were used to predict the impact of identified base substitutions.

DNA copy number and breakpoint detection were assessed using LIWG data (Liang et al. 2014). For copy-number variation (CNV) detection, (1) read depths were determined every 100 bases across sequenced regions; (2) normalized log2 fold changes between the tumor and germline were calculated; and (3) a smoothing window was applied. The circular binary segmentation algorithm was used to partition chromosomes into sections of equal copy number to correct for noise in order to calculate log2 copy fold changes (Olshen et al. 2004).

RNA reads were aligned to GRCh37 using STAR (Dobin et al. 2013). Differential expression analysis was performed for each tumor sample against Universal Human Reference RNA (Agilent) using Cufflinks/Cuffdiff (Trapnell et al. 2012), which performs transcript assembly to identify differentially expressed genes, and DESeq2 (Anders and Huber 2010; Love et al. 2014), which uses the negative binomial distribution. False discovery rate adjustment of P-values for expression changes were performed to calculate a Q-value, or a corrected P-value, using the Benjamini–Hochberg correction for multiple testing. Expression levels are expressed as fragments per kilobase of transcript per million mapped reads (FPKMs). We defined a threshold FPKM value of >1 to indicate expression of a gene.

Pathway analysis was performed using differentially expressed genes (corrected P < 0.05) identified using both Cuffdiff and DESeq2. Analysis was performed using the Ingenuity Pathway Analysis (IPA) tool (version 01-12; QIAGEN) to identify potentially perturbed pathways in each tumor. Significant pathways were identified with P < 0.05. In a separate analysis for Patient 1, KEGG networks N00019 (FGF-FGFR-RAS-ERK signaling pathway) and N00037 (FGF-FGFR-PI3K signaling pathway) were used to identify factors downstream from FGFR; interrogated genes include GRB2, SOS1, SOS2, HRAS, KRAS, NRAS, ARAF, BRAF, RAF1, MAP2K1, MAP2K2, MAPK1, MAPK3, PIK3CA, PIK3CB, PIK3CD, AKT1, AKT2, AKT3, MTOR, RPS6KB1, and RPS6KB2.

Drug–Gene Matching

All somatic aberrations identified through LIWG, exome, and RNA sequencing data, including SNVs, small indels (insertion/deletions), CNVs, translocations, gene fusions, and RNA expression changes, were evaluated for each case. Potential therapeutic implications for RNA expression events are limited to genes that are differentially expressed compared to the Universal Human Reference RNA (corrected P < 0.05). For any abnormalities identified, we looked for potential therapies using a proprietary drug–gene database annotated with a pharmacopeia of FDA-approved agents. Genes and alterations associated with FDA-approved agents are reported here.

ADDITIONAL INFORMATION

Data Deposition and Access

The sequencing data generated from this study are available in the NCBI (National Center for Biotechnology Information) dbGaP (Database of Genotypes and Phenotypes) repository (https://www.ncbi.nlm.nih.gov/gap) under controlled access (accession number phs001643).

Ethics Statement

This work was approved by the IRB of St. Joseph's Hospital and Medical Center in Phoenix, Arizona. The clinical review was performed on de-identified records. Written informed consent was obtained from the patients to publish the work detailed here.

Acknowledgments

We would like to thank Jennifer Eschbacher, MD, and Catherine Seiler, PhD, of the Biobank Core Facility (St. Joseph's Hospital and Barrow Neurological Institute); Jessica Aldrich, Dr. Jonathan Keats, and Tyler Izatt for bioinformatics support (TGen): Dr. Michael Berens for support; and Cynthia Lechuga, Diana Davis, and Stephanie Buchholtz (TGen) for administrative support. Nancy Linford, PhD, provided editing assistance.

Author Contributions

W.S.L. and A.S.L. conceived the study. W.S.L. and C.D. wrote the manuscript. A.H., J.A., and L.C. performed nucleic acid extractions, prepared sequencing libraries, and performed sequencing. S.S. and D.E. performed bioinformatics analyses. S.B. performed drug–gene matching. S.B. and W.S.L. performed data interpretation. All authors reviewed the manuscript.

Funding

This work was supported by the Barrow Neurological Foundation and the Translational Genomics Research Institute. The funding bodies had no role in the study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Competing Interest Statement

The authors have declared no competing interest.

Referees

Robert Maki

Anonymous