Patient population involved
This retrospective study was conducted in agreement with the Reporting Recommendations for Tumor Marker Prognostic Studies (REMARK) criteria for reporting tumor biomarker prognostic studies [17]. A checklist for the criteria is provided (Additional file 1).
Set 1 (frozen tissue, University Medical Center Groningen (UMCG))
The discovery set consisted of prospectively collected chemo-naïve frozen tumor tissue of 18 patients with advanced-stage HGSOC operated by a gynecologic oncologist from the UMCG (Groningen, the Netherlands) in the period 1990–2008. All patients were staged according to the International Federation of Gynecology and Obstetrics (FIGO) guidelines. Standard treatment included debulking surgery followed by adjuvant chemotherapy consisting of platinum-based treatment regimens. After chemotherapy, patients were followed for up to 10 years with gradually increasing intervals. All the clinicopathological and follow-up data have been registered in an anonymous, password-protected database, in compliance with Dutch law. All patients gave informed consent. The responder group consisted of patients with advanced stage HGSOC, residual disease after primary surgery (> 2 cm), treated with adjuvant platinum-based chemotherapy and a PFS of more than 18 months. The non-responder group consisted of patients with advanced stage HGSOC, residual disease after primary surgery (> 2 cm), treated with adjuvant platinum-based chemotherapy and a PFS of less than 6 months. We have p53 and BRCA1/2 status information for 17 (8 responders and 9 non-responders) out of 18 discovery dataset patients. In this cohort, 16 were p53 mutated, except for one non-responder, and only two responders showed BRCA2 germline mutation. Detailed clinicopathological features are described in Additional file 2: Table S1.
Set 2 (mRNA dataset, UMCG)
This previously published gene expression data included 157 consecutive advanced-stage HGSOC patient samples from UMCG profiled using two-color oligonucleotide microarrays (35,000 Operon v3.0 probes), manufactured by The Netherlands Cancer Institute (Amsterdam, the Netherlands, https://www.nki.nl/topmenu/genomics-core-facility/, GSE 13876) as described by Crijns et al. [18]. For expression data integration, we used data from 11 patients (6 responders and 5 non-responders) who were also in the discovery set for MethylCap-Seq. Detailed clinicopathological features are described in Additional file 2: Table S1.
Set 3 (frozen tissue; UMCG + Innsbruck + Leuven)
The external validation cohort consisted of HGSOC patient tumors from 21 responders and 31 non-responders obtained from the UMCG, the Medical University of Innsbruck (Austria), and University Hospital Leuven (Belgium). All patients were selected based on the same inclusion criteria as the discovery set (Set 1). Detailed clinicopathological features are described in Additional file 2: Table S1.
Sets 4, 5, and 6 (Publicly available external cohort data)
For in silico validation of our findings, we used publicly available methylation and expression datasets from HGSOC patients. For methylation Set 4, Infinium 450K methylation array data of AOCS study group (http://www.aocstudy.org.) was extracted from the NCBI GEO portal using GEO accession no. GSE65820 [19] (http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE65820) and normalized using a beta-mixture quantile normalization as described previously [20]. The clinical data of patients was downloaded from the ICGC data portal (https://dcc.icgc.org/). Methylation probes for FZD10 (cg23054883), FAM83A (cg24833277), MYO18B (cg24035545), and MKX (cg14947429), which are related to same marker regions as identified with MethylCap-seq (shown in Fig. 2a for FZD10), have been used for univariate Mantel–Cox log-rank survival analysis to generate Fig. 4a, b and Additional file 3: Figure S4–S6. For methylation Set 5, Infinium 27K methylation array data of TCGA study group along with associated clinical information was extracted from Genomic Data Commons portal (https://gdc.cancer.gov/). Data was beta-mixture quantile normalized and the FZD10 methylation probe (cg23054883) was used for univariate Mantel–Cox log-rank survival analysis to generate Additional file 3: Figure S3.
Gene expression data from tumors of Set 6, the curated ovarian cancer database (Gyorffy et al. database, http://kmplot.com) [21], with carefully curated clinical annotations, was extracted in November 2015. We restricted our analysis to primary, high-grade (3), advanced stage (3 and 4), serous ovarian tumors, residual disease after surgery or suboptimal debulking, and platinum-containing therapy with available PFS and overall survival (OS). For candidate genes, we used the expression data of FZD10 (Probe ID: 219764_at), FAM83A (Probe ID: 239586_at), MYO18B (Probe ID: 1554579_a_at), and MKX (Probe ID: 239468_at) to perform univariate Mantel–Cox log-rank survival analysis and to generate Fig. 4c, d and Additional file 3: Figure S4–S6.
Detailed clinicopathological features are described in Additional file 2: Table S1.
DNA extraction and bisulfite modification
Histologic slides of patients were reviewed to confirm the diagnosis by an experienced gynecologic pathologist. Representative frozen blocks of each patient tumor were retrieved for DNA extraction. Frozen sections of 10 µm thickness were cut with periodic 4 µm sections prior to hematoxylin and eosin staining to assess the vital tumor cell percentage. For some samples, slides were macro-dissected to obtain more than 85% neoplastic cells. DNA was isolated using standard salt-chloroform extraction and isopropanol precipitation. Precipitated DNA was re-suspended in Tris-EDTA buffer (10 mM Tris; 1 mM EDTA, pH 8.0). Genomic DNA was amplified with multiplex PCR according to the BIOMED-2 protocol to check the DNA’s structural integrity [22]. DNA quantity was measured using Quant-iT™ PicoGreen® dsDNA Assay kit according to the manufacturer’s protocol (Invitrogen, Carlsbad, CA, USA). For DNA isolation from cell lines, the same standard method was followed. Bisulfite conversion was performed using EZ DNA methylationtm Kit (Zymo Research, Orange, CA, USA) according to the manufacturer’s protocol using 1 μg of DNA.
MethylCap-seq
MethylCap-seq was performed as described previously [23, 24]. Briefly, methylated DNA fragments were captured with methyl binding domains using the MethylCap kit according to manufacturer’s instructions (Diagenode, Liège, Belgium). The kit consists of the methyl binding domain of human MeCP2 as a C-terminal fusion with glutathione-S-transferase containing an N-terminal His6-tag. Before capturing, DNA samples (500 ng) were sheared to a size range of 300–1000 bps using a Bioruptor™ UCD-200 (Diagenode, Liège, Belgium) and fragments of approximately 300 bp were isolated. Captured DNA was paired-end-sequenced on the Illumina Genome Analyzer II platform according to the protocol (Illumina, San Diego, CA, USA). Results were mapped on the nucleotide sequence using Bowtie software [25], visualized using BioBix’ H2G2 browser (http://h2g2.ugent.be/biobix.html) and processed using the human reference genome (NCBI build 37). The paired-end fragments were unique and located within 400 bp of each other. The MethylCap-seq data have been deposited in the Gene Expression Omnibus under accession number GSE97128.
Bisulfite pyrosequencing
Based on the next-generation sequencing results of the discovery set (Set 1), all pyrosequencing primers were designed for the selected candidate differentially methylated regions (DMRs) of 45 genes using PyroMark Assay Design software (Qiagen, Hilden, Germany). Bisulfite pyrosequencing was performed as described previously [26]. Briefly, bisulfite treated DNA was amplified using PyroMark PCR kit (Qiagen). PCR reaction and cycling conditions were according to the kit manual. To generate the PCR product from bisulfite-converted DNA, we adopted the amplification protocol using a universal primer approach as described by Collela et al. [27]. The biotinylated PCR products were captured using 1 μL streptavidin-coated sepharose high-performance beads (GE Healthcare, Little Chalfont, UK). The immobilized products were washed with 70% alcohol, denatured with PyroMark denaturation solution (Qiagen) and washed with PyroMark wash buffer (Qiagen). The purified PCR product was then added to 25 μL PyroMark annealing buffer (Qiagen) containing 0.3 μM sequencing primers for specific genes (primer sequences are given in Additional file 4). Finally, pyrosequencing™ reactions were performed in a PyroMark Q24 MD System (Qiagen) according to the manufacturer’s instructions using the PyroGold Q24™ Reagent Kit (Qiagen). CpG site methylation quantification was performed using the methylation Software Pyro Q24 2.06 Version (Qiagen).
Cell line culturing
A panel of human ovarian cancer cell lines, A2780, C30, Cp70, SKOV3, OVCAR3, IGROV1, PEA1, PEA2, PEO14, and PEO23, was used for in vitro validation and functional analysis. The source, media, and culture conditions for the cell lines are shown in Additional file 2: Table S2. All cells were grown at 37 °C in a humidified atmosphere with 5% CO2 and were detached with 0.05% trypsin in phosphate-buffered saline (PBS, 0.14 mM NaCl, 2.7 mM KCl, 6.4 mM Na2HPO4, 1.5 mM KH2PO4, pH 7.4). Authenticity of all cell lines was verified by DNA short tandem repeat analysis (Baseclear, Leiden, The Netherlands) and mycoplasma testing was performed using in-house developed PCR-based method with specific primers (Invitrogen, NY) against various mycoplasma species. For global demethylation, cells at 40–50% confluency were treated with demethylating agent 5-aza-2′-deoxycytidine (DAC) at a final concentration of 1 μM for 72 h. Due to the low stability of DAC at 37 °C, the medium was replenished with DAC every 24 h. After 72 h, cells were trypsinized and processed for RNA and DNA isolation.
Total RNA isolation and quantitative reverse transcriptase PCR (qRT-PCR)
qRT-PCR was performed as described previously [26]. Total RNA was isolated from frozen tissue blocks and cell lines using the same procedure as described for DNA extraction. Total RNA was isolated using RNeasy mini kit (Qiagen) according to the manufacturer’s instructions. RNA was analyzed quantitatively using Nanodrop (Nanodrop Technologies, Rockland, DE), by using 1 μg of total RNA for cDNA synthesis by a RNase H+ reverse transcriptase using iScript cDNA synthesis kit (BioRad, Hercules, CA) according to the manufacturer’s instructions. qRT-PCR was performed in an ABI PRISM 7900HT Sequence Detector (Applied Biosystems, Foster City, CA) with the iTaq SYBR Green Supermix with Rox dye (Biorad). The reactions were analyzed with SDS software (Version 2.4, Applied Biosystems). The threshold cycles (Ct) were calculated and relative gene expression (∆Ct) was analyzed with GAPDH as the housekeeping gene (∆Ct = Ctgene – CtGAPDH) (primer sequences given in Additional file 4). The qRT-PCR primers used are available upon request. For the final analysis, data was imported to R to perform clustering and ggplot2 (http://ggplot2.org/) was used to create heat maps.
siRNA mediated silencing for in vitro experiments
Cells (1–3 × 105) were plated in a 6-well plate and grown overnight. FZD10 trisilencer-27 siRNAs (Origene Technologies, Rockville, MD) were used for transient knock-down using 20 nM of final concentration of siRNA (sequences given in Additional file 4). Scrambled and FZD10 targeted siRNAs were transfected using Oligofectamine (Invitrogen, NY) for 4 h with reduced growth factor serum-free opti-MEM media (Gibco, Life Technologies, CA). Subsequently, cell line-associated media (Additional file 2: Table S2) with 30% FCS were added to make a final FCS concentration of 10% for 48 h. Following 48 h after siRNA transfection, other functional assays (short- and long-term survival, migration, and apoptosis) were performed.
Short- and long-term survival assays
Short-term cellular viability was measured with the micro-culture tetrazolium assay (MTT) as described previously [28]. Briefly, in a 96-well culture plate, approximately 7500 SKOV3 cells, 15,000 OVCAR3 cells, 10,000 PEA2 cells, and 12,000 C-30 cells, either control or siRNA transfected, were seeded in 200 μL culture medium with or without cisplatin treatment. After 96 h, 20 μL of 3-(4,5-dimethythiazol-2-yl)-2,5-diphenyltetrazolium bromide (MTT, Sigma-Aldrich, St. Louis, MO, 5 mg/mL in PBS) was added and formazan production was measured colorimetrically using a Biorad iMark microplate reader at 520 nm wavelength.
For long-term assay, depending on the concentration of cisplatin, cells were seeded in 96-well plates at approximately 2000 cells per well for SKOV3 and 4000 cells per well for OVCAR3. After 8–10 h, indicated doses of cisplatin were added and allowed to grow for a set number of days. Finally, cells were fixed and stained in staining buffer (methanol (50%), acetic acid (20%) and 0.01% Coomassie brilliant blue), washed with water and dried, after which the plates were scanned. For quantification, 200 μL of 10% acetic acid was added to each well and left on a shaker for 30–60 min. Plates were read using a Biorad iMark microplate reader at 520 nm wavelength.
Wound healing assays
For the wound healing assays, cells were seeded in a 6-well plate at a density of 2 × 105 cells/well and grown overnight until confluency. A wound was created by manually scraping the cell monolayer with a 10 μL pipette tip and the medium was aspirated to remove the detached cells. Cells were then incubated with medium supplemented only with 10% FCS, and wound closure was observed within 24 h. Images were acquired with a Leica camera mounted on an inverted microscope and were processed using Image J software (NIH, Bethesda, MD; http://rsb.info.nih.gov/ij/). The distance cells migrated was determined by measuring the wound area at different time points followed by its correction from the wound area at time 0 h.
Western blot analysis
Various proteins in ovarian cancer cell lines were detected by the western blotting method as described previously [28]. Western blot membranes were probed overnight at 4 °C with primary antibodies (PARP, Cell Signaling Technology, Danvers, MA, #9532; Cleaved Caspase-3, Cell Signaling Technology, #9661). Afterwards, HRP-conjugated secondary antibodies (DAKO, Glostrup, Denmark) were used for detection using Lumi-LightPLUS Western Blotting Substrate (Roche Diagnostics, Hilden, Germany). Membranes were probed with β-actin antibody (mouse, A5441; Sigma-Aldrich, St. Louis, MO) to confirm equal loading.
Statistical analysis
MethylCap-seq
All methylation reads data were extracted using BioBix’ H2G2 browser (http://h2g2.ugent.be/biobix.html) for the broad promoter region (2000 bp upstream and 500 bp downstream of the transcription start site). Read counts were statistically compared between responder and non-responder groups using R/Bioconductor [29] package EdgeR [30], assuming that the data follow a negative binomial distribution, and ranked on P value.
Subsequently, integration of expression data was also performed using R-package LIMMA to find differentially expressed genes [31]. As an additional filter for further analysis, each candidate DMR had to be methylated (≥ 4 reads) in at least four samples of a specific response group. Given the fact that putatively relevant loci were selected based on both differential methylation and expression, and that several rounds of subsequent independent biological validation were performed, a relatively permissive error rate control cut-off (P = 0.05) was used for expression as well as validation.
Bisulfite pyrosequencing
Methylation percentage results were analyzed using statistical software IBM SPSS 21 (SPSS Inc., Chicago, IL) and a non-parametric statistical test (Mann–Whitney U test) was performed to find differences between responder and non-responder groups. P values of less than 0.05 were assumed to be statistically significant for all tests. To present data as a heatmap, all methylation percentage data were imported to Genesis software (Graz University of technology, genome.tugraz.at/genesis) for clustering and heatmap visualization.
In silico validation of candidate markers
For prognostic validation of candidate gene methylation, methylation data of the AOCS and TCGA study groups was extracted and normalized as mentioned in ‘patient population involved’ Set 4 and Set 5, respectively. Low and high methylation cut-offs were based on median beta value. This resulted in 89 patients for PFS analysis (proxy for sensitivity to platinum containing chemotherapy) and 91 patients for OS analysis in AOCS data (Set 4). For the TCGA cohort (Set 5), we used 91 patients for PFS analysis and 105 patients for OS analysis. To handle the missing data, we used the listwise deletion methodology.
For marker expression, data (Set 6) was derived for analysis using KM plotter [21] in November 2015, in which we selected only advanced stage (3 and 4) HGSOC cancer patients with suboptimal debulking surgery, all of whom had received platinum therapy. This resulted in 200 patients for PFS and 208 patients for OS analysis using univariate Mantel–Cox log-rank survival analysis with FZD10 probe (Probe ID: 219764_at), and 100 patients for PFS and 102 patients for OS analysis with FAM83A (Probe ID: 239586_at), MYO18B (Probe ID: 1554579_a_at), and MKX (Probe ID: 239468_at). With an expression range of probes for different genes, an auto cut-off value for PFS and OS analysis was used, based on computation of upper and lower quartiles with the default portal settings [21].
To review the gene expression of FZD10 in other cancer types, we used the TCGA data from the TCGA FIREHOSE pipeline (http://gdac.broadinstitute.org/) [32]. To predict FZD10 expression across 41 tumor types, we used their functional genomic mRNA (FGmRNA) profiles as described earlier [33, 34]. In this methodology, non-genetic transcriptional components were used as covariates to correct microarray expression data and the residual expression signal (i.e., FGmRNA profile) was found to capture the downstream consequences of genomic alterations on gene expression levels [33]. We quantified the percentage of samples across 41 tumor types with a significantly increased FGmRNA signal (i.e., a proxy for underlying gene amplification). For each of the 19,746 tumor samples, FZD10 was marked as significantly amplified when the FGmRNA signal was above the 97.5th percentile threshold as defined in the non-cancer samples [33].
In vitro experiments
Statistical significance was calculated by two-sided Student’s t test between two groups, unless otherwise mentioned in the figure legends. P values of less than 0.05 were defined as statistically significant for all tests.