Introduction

Breast Cancer (BC) is the most common malignancy with an incidence of 24.2% and mortality of 15% among all cancers in women. The mortality rates are higher in transitioning countries including Africa and South Asia, indicating the requirement for effective population based screening programs to detect BC at the earliest1. Mammography and ultrasonography are commonly used for detection of breast abnormalities but the sensitivity of these imaging modalities is moderate and results in false negatives (4–34%) particularly in women (under 40 years) with dense breast tissue2,3. Additionally, in low-resource setting access to Mammogram can be limited and hence a non-invasive test could help identify those who would need further investigations. Thus, there is an unmet need to develop a strategy that circumvents the shortcomings of the existing methods.

In the past decade, technological advances have facilitated the mining of blood-based BC biomarkers at genomic, transcriptomic, proteomic and metabolomic scale4,5,6,7. Carcinoembryonic antigen (CEA)8,9, human epidermal growth factor receptor-2 (HER2)10,11, the oncogenic protein RS/DJ-112, tissue polypeptide antigen (TPA)13 tissue polypeptide specific antigen (TPS)14 and CA 15–315,16 have been recommended as biomarkers for metastatic BC, as they were deficient in early BC detection. Amongst the potential multianalyte biomarker panels for early detection of BC17,18,a panel of 5 complex autoantigens [B11 (LGALS3), B18 (PHB2), B119 (MUC1), B130 (GK2), and CA15-3], had a high sensitivity of 87% and specificity of 76%19.

Parallelly, the analysis of global and gene-specific DNA methylation profiles in circulating cell-free DNA (cfDNA) has led to the discovery of epigenetic markers for BC detection20. Liu et al., carried out targeted sequencing of 9223 frequently methylated CpG sites in plasma cfDNA samples of cancer patients. Methylation scores could accurately detect colorectal cancer (96.3%), followed by breast cancer (91.7%), melanoma (81.8%) and NSCLC (61.1%)21. The PanSeer assay can simultaneously evaluate a set of differentially methylated genes to detect five cancers (stomach, oesophagus, colorectum, lung or liver). It displayed an overall sensitivity of 95% in detecting the five cancers 0–4 years prior to clinical diagnosis22. In BC, the circulating levels of differentially methylated genes such as BRCA123, ATM24, CDKN2A25, ERα26, APC27, RASSF1A and RARβ228,29 have also been extensively evaluated in single or combination. The factors limiting the translation of these markers include lack of desired sensitivity and specificity, limited sample size, lack of prospective independent validation, high cost and absence of significant clinical correlation30. Thus, the present study aims to develop a combinatorial biomarker assay that serves as a minimally invasive, cost-effective adjunct for BC diagnosis. To achieve this, a transcriptomic profile using microarray was obtained from 41 primary breast tumours (T), 18 paired normal (PN) (adjacent normal tissue to the tumour confirmed to be morphologically normal with no tumour cell present in the sample) and 6 apparent normal (AN) (patients with benign breast disease had breast tissue sampled well away from the lesion) tissues procured from South Asian women. The differentially expressed genes were validated using Quantitative Real-Time Polymerase Chain Reaction (qRT-PCR). Fifteen candidate markers were chosen, and their protein levels were first assessed in tissues [T, PN, AN] and then in plasma of patients with breast cancer (stage 0–4), benign breast abnormalities and in healthy volunteers. The protein marker panel lacked the desired diagnostic power, hence 3 candidate genes which were found to be downregulated in our BC microarray study and confirmed to be methylated in tumours but not in normal breast tissue or in peripheral blood mononuclear cells (PBMC) were chosen for epigenetic marker evaluation. Detection of these methylated genes in cfDNA was performed in the same set of subjects recruited for the protein marker study. The circulating methylation biomarkers achieved 100% sensitivity and up to 90% specificity in discrimination of BC from benign cases and controls.

Results

Study population

The clinicopathological details of the samples used for the microarray are given in Supplementary Table S1 and that used in antibody array are furnished in Supplementary Tables S2 and S3 respectively.

Identification of potential biomarkers using BC gene expression profiling

Transcriptomic profiling was performed on 41 T,18 PN and 6 AN tissue to obtain a BC gene expression signature (Supplementary Table S1). The raw data analysis using the Class comparison module in the BRB-Array Tools software v 3.7.0 (criteria used were P-value = 0.001 and twofold or higher difference) identified 57 genes overexpressed in T compared to either AN or PN, 146 genes with elevated expression in PN compared to AN or T and 173 genes with higher expression in AN compared to PN or T [indicative of downregulated genes in T] (Fig. 1a). Validation of 67 differentially expressed genes (short-listed from microarray analyses) using TaqMan® Low density arrays (TLDA) resulted in 56 genes being identified as differentially expressed, resulting in more than 80% concordance with the microarray data. Fifteen genes were highly overexpressed in T compared to AN and PN; 9 genes were overexpressed in PN compared to AN and T; 32 genes were grossly downregulated in tumours compared to PN or AN samples (Fig. 1b, c). Further, we confirmed the differentially expressed genes from our series using an independent GEO Microarray Dataset [GSE 22820] which included 176 primary breast tumour samples and 10 normal breast tissues from reduction mammoplasties.

Figure 1
figure 1

Identification of candidate biomarkers for breast cancer by microarray analyses. (a) Unsupervised hierarchical clustering of primary breast tumour, paired normal (PN) and apparently normal (nor) tissues. Each column represents an individual sample and each row a gene. Upregulated genes are indicated in red and downregulated genes in green. (b) List of genes upregulated in tumour and paired normal tissues relative to apparently normal (> 2.5-fold differential expression, P < 0.0001) in gene expression profiling. (c) List of genes downregulated in tumour and paired normal tissues relative to apparently normal (> 2.5-fold differential expression, P < 0.0001).

Based on the differential expression, 15 markers [CFD, WIF-1, aFGF, bFGF, DKK-3, sFRP-3/FRZB1, LEP, IL-17B, IGF-I, IP-10/CXCL10, LOX-1/OLR1, MIG/CXCL9, MIP-1d/CCL15, OPN/SPP1, SDC1] which were either secreted or likely to be secreted were chosen for further validation.

Validation of protein biomarkers in BC tissues and plasma

To determine the protein levels of 15 markers, Quantibody array (RayBiotech, USA) based multiplexed sandwich ELISA was first done in tissues (46 T, 23 PN and 6 AN). Among the 15 markers the levels of CFD (P = 0.0117) and LEP (P < 0.0001) were significantly downregulated in T compared to AN and PN. In contrast, 7 out of 15 proteins including, IL17B (P = 0.0001), IP10 (P = 0.0135), LOX-1 (P < 0.0001), MIG (P < 0.0001), MIP-1d (P < 0.0001), SDC1 (P < 0.0001) and OPN (P < 0.0001) were significantly upregulated in tumour tissues compared to normal (AN and PN) (Table 1).

Table 1 The expression levels of the 15 candidate biomarkers in breast cancer tissues and normal, estimated by quantibody array.

Further, to quantify the levels of these markers in plasma, Quantibody array analysis was carried out on the plasma from patients with invasive (n = 202), non-invasive (n = 16), benign (n = 37) and healthy control subjects (n = 203). (Supplementary Fig. S1). Plasma levels of CFD (P = 0.002), LEP (P = 0.0629) and DKK3 (P = 0.0116) were found to be downregulated in patients with Ductal carcinoma in-situ (DCIS) and invasive BC when compared to controls (healthy controls and patients with benign breast conditions). The markers significantly upregulated in BC relative to controls were FGF1 (P = 0.001), FGF2 (P = 0.001), sFRP3 (P = 0.007), IGF-I (P = 0.001), IL-17B (P = 0.001), IP10 (P = 0.001), LOX-1 (P = 0.0026), MIG (P = 0.0001), MIP1d (P = 0.0001), OPN (P = 0.0237), SYN-1 (P-0.0001), WIF1(P = 0.037) Table 2. Of the 15 markers, the levels of FGF2, sFRP3 and IL17B were highest in the invasive cases relative to DCIS and they gradually declined in benign and control groups. In contrast, LEP levels were lowest in invasive group and the levels gradually increased in DCIS, benign and controls (Supplementary Fig. S2).

Table 2 The median (range) concentration of each protein in plasma of cases and controls along with their area under the curve (AUC) value of single markers.

Diagnostic Evaluation of the protein markers in plasma of BC patients

To analyse the diagnostic properties of the protein markers, the plasma samples were split into training (N = 371) and test (141) sets (Supplementary Table S4). The individual diagnostic performance of the 15 markers was determined using the receiver operating characteristic (ROC) analysis. The specificity was fixed at 70% and sensitivities were generated. The markers which had an AUC < 0.5 were CFD (35.7 ng/ml, AUC 0.494, 31%), LEP (45.7 ng/ml, AUC 0.455, 29%) and DKK3 (124.6 ng/ml, AUC 0.418, 30%). Markers with AUC > 0.5 were FGF2 (0.0068 ng/ml, AUC 0.628, 47%), IL17B (0.67 ng/ml, AUC 0.655, 52%) and SDC1 (4.565 ng/ml, AUC 0.647, 49%). To improve the sensitivity, multiparameter diagnostic evaluation was performed. The Median + 1 MAD (Median Absolute Deviation) value of each marker across cases and controls set as cut-off to generate sensitivity and specificity. Leave one out biomarker analysis was performed to choose combination with maximum sensitivity. (Supplementary Table S5). Overall, the combination of 6 proteins CFD, LEP, DKK3, IL17B, SDC1, FGF2 showed the best diagnostic potential with sensitivity and specificity of 65% and 79% in training set (AUC 0.661), and 70% and 85% in test set (AUC 0.678). respectively (Table 3). Although the sensitivity of multiple proteins improved the sensitivity of the model, it did not satisfy our objective, hence we evaluated additional markers.

Table 3 The diagnostic performance of different multi-marker models in detecting breast cancer.

Identification of Circulating DNA methylation markers

The top downregulated genes [WIF1, DACT2, ID4, TP63, SOX10] in tumours in our transcriptomic profile were regulators of Wnt signalling, cell differentiation and morphogenesis and acted as inhibitors of tumour growth in BC31,32,33,34. Epigenetic silencing of various Wnt signalling antagonists promote abnormal cell proliferation in BC, since they have a tumour suppressive role35. Hence, we chose WIF1 and DACT2 since they showed the greatest downregulation in T relative to AN (Fig. 1c). Additionally, SOSTDC1 found to be downregulated in tumours in our gene expression study, had been confirmed to be downregulated in 98.2% (56/57) of breast tumour tissues and cell lines due to promoter hypermethylation in our earlier study, was also included36.

Epigenetic silencing of WIF1, DACT2 and SOSTDC1 in cell lines and primary breast carcinomas

To validate the downregulation of WIF1, DACT2 and SOSTDC1 in various subtypes of BC, qRT-PCR was performed in breast-cancer derived cell lines: ZR751, MCF7, T47D (ER/PR+ HER2-), SKBR3 (ER/PR- HER+), MDAMB231, MDAMB468 (ER/PR- HER2-) and a non-tumorigenic breast cell line HBL10037. The mRNA levels of all three genes were found to be decreased in BC cell lines with greater downregulation in triple negative MDAMB231, MDAMB468 cells and HER2 overexpressing SKBR3 cells when compared to others (Fig. 2a). We then did Methylation Specific PCR (MSP) analysis and found complete methylation of WIF1 and SOSTDC1, and partial methylation of DACT2 in all BC cell lines (Fig. 2b).

Figure 2
figure 2

Downregulation of SOSTDC1, DACT2, WIF1 in breast cancer cell lines is associated with methylation. (a) Expression levels of SOSTDC1, DACT2, WIF1 in breast cancer cell lines relative to non-tumorigenic HBL100, analysed by qRT-PCR. The CT values were normalized against the internal control GAPDH. Each sample was run in triplicates and mean ± SD of 2^∆∆CT values were plotted in log scale, (b) MSP results of SOSTDC1, DACT2, WIF1 in breast cancer cell lines, (c) Representative figure of bisulphite gene sequencing of the putative promoter region of DACT2 and WIF1 in MDAMB231 cells using gene specific BSP primers. The DACT2 and WIF1 amplicons were 257 bp and 493 bp respectively. The black symbols () in chromatogram represent the methylated CpG sites.

Bisulfite Sequencing PCR (BSP) analysis of the putative promoter regions of WIF1 and DACT2 uncovered methylated CpG sites, concordant with the MSP analysis (Fig. 2c). The BSP sequencing of SOSTDC1 was not repeated since, its promoter region was already characterized with the promoter region having a 54 bp CpG island upstream of 5’UTR encompassing 4 CpG sites38. Further, the methylation status of these genes was validated in 10 primary breast carcinoma tissues (having > 70% tumour cells) and 10 paired normal tissues. In concordance with the cell lines, SOSTDC1 and WIF1 were completely methylated but, DACT2 was hemi-methylated in tumours. All three markers were unmethylated in the paired normal tissues and lymphocytes confirming that the methylated alleles of SOSTDC1, DACT2 and WIF1 originate from the breast tumour (Fig. 3a–c).

Figure 3
figure 3

Differential methylation of candidate gene markers in breast cancer. MSP analysis of (a) DACT2, (b) SOSTDC1, and (c) WIF1 in primary breast tumour tissues, M/U represents the ratio of methylated allele and unmethylated allele, L-lymphocytes, T-breast tumour tissues, N-paired normal tissues (d) MSP analysis of DACT2, SOSTDC1 and WIF1 in plasma samples from patients with invasive BC (C = Cancer), non-invasive BC (D = DCIS), benign breast disease (B = benign) and healthy volunteers (O = Controls). The amplicons were resolved in 2% agarose gel. POS- positive control (100% in-vitro methylated and unmethylated DNA standards), NTC-no template control, M-methylated allele, U-unmethylated allele. (e) Densitometric analysis results of the amplicons obtained from methylation sensitive-PCR. Relative intensity of each band was calculated with 100 bp DNA ladder as reference. M/U represents Ratio of relative intensity of Methylated allele/Unmethylated allele respectively. The median with range intensity was plotted and Kruskal Wallis statistic was used to test the differential methylation intensity. The full-length gels of above images are provided in Supplementary Fig. S4.

WIF1, DACT2 and SOSTDC1 as potential circulating biomarkers in BC

For the validation of SOSTDC1, DACT2 and WIF1 in plasma samples, cfDNA was isolated from patients (invasive BC = 202, DCIS = 16, benign = 37) and controls (n = 203) recruited for the case control study. Increased concentration of cfDNA was detected in plasma from patients with invasive BC and DCIS when compared to plasma from patients with benign breast diseases and healthy controls (P < 0.0001). cfDNA levels were associated with stage of the disease and were significantly elevated (P < 0.0001) in stage II (Median [M]: 437 ng/ml, Range [R]:236–889), stage III (M: 763 ng/ml, R: 310–2262) and stage IV (M: 1687 ng/ml, R: 421–2447) patients when compared to stage I (M:297 ng/ml, R:279–320) and DCIS (M:221 ng/ml R:178–513) patients. Although the difference in cfDNA levels were not statistically significant between healthy volunteers (M:174 ng/ml, R:66–317), patients with benign breast disease (M:181 ng/ml, R:98–295), DCIS (M:221 ng/ml R:178–513) and early BC (M:297 ng/ml, R:279-320 ng/ml) an increasing trend in the cfDNA concentration was observed (Supplementary Fig. S3). This lack of significance could also be due to the smaller number of DCIS and stage 1 patients included in the study.

MSP detection of SOSTDC1, DACT2 and WIF1 was carried out with bisulphite modified cfDNA samples isolated from the cases and controls. The amplicons were visualised in 2% agarose gel. (Fig. 3d, Supplementary Fig. 4). SOSTDC1 was methylated in 95.5% (193/202) of invasive BC patients and 75% (12/16) of DCIS patients, while it was unmethylated in 83.7% (31/37) benign patients and (202/203) in 99.5% of healthy controls. DACT2 was methylated in 91.08% (184/202) of invasive BC and 75% (12/16) of DCIS patients, and it was unmethylated in 89% (33/37) benign patients and (193/203) in 95.07% of healthy controls. The methylated allele of WIF1 was detected in 99% (200/202) of invasive BC and 75% of DCIS patients and unmethylated allele was detected in 83.7% (31/37) of benign patients and 98.5% (200/203) of controls. Densitometric analysis was used to obtain the ratio of relative intensity of the methylated and unmethylated alleles (M/U) in each sample. The methylation intensity of SOSTDC1, DACT2 and WIF1 was significantly elevated (P < 0.0001) in the invasive BC and DCIS, when compared to benign and controls. The benign cases were poorly methylated or unmethylated but a slight increase in methylation intensity was observed in benign lesions with proliferative potential (data not shown). SOSTDC1 and DACT2 did show prominent differential methylation in DCIS and invasive BC, but methylated WIF1 was significantly increased (P = 0.031) in DCIS tumours (Fig. 3e, Supplementary Table S6).

Association of SOSTDC1, DACT2 AND WIF1 methylation with clinicopathological variables

The association of SOSTDC1, DACT2 AND WIF1 methylation with the clinicopathological parameters of 218 patients (202 invasive BC, 16 DCIS) in the case control study was studied. Elevated methylation intensity of SOSTDC1, DACT2 and WIF1 was observed in patients who had advanced tumour stage, positive nodes and local or distant metastasis. In addition, BC patients with grade III had high levels of methylated SOSTDC1 and DACT2. However, there was no correlation identified between the methylation patterns of the genes with age, menstrual status and receptor status, with the exception of methylated DACT2 which was high in patients with HER2- BC (Table 4).

Table 4 Association of SOSTDC1, DACT2 and WIF1 methylation with clinicopathological parameters of 202 invasive BC patients and 16 DCIS patients.

Diagnostic evaluation of SOSTDC1, DACT2 AND WIF1 as a multi-marker panel for BC

To analyse the diagnostic performance of the predicted diagnostic model, invasive BC, DCIS patients were clubbed as cases and the benign and healthy controls were combined as controls. The samples were randomly split into training and test sets. Univariate ROC analysis of the markers in the training set (151 cases, 166 controls) revealed an AUC of 0.993 (95% CI 0.9866–0.9997) for SOSTDC1, AUC of 0.994 (95% CI 0.9898–0.9987) for DACT2, AUC of 0.9966 (95% CI 0.9930–1.000) for WIF1. Multivariate ROC analysis of the 3 gene panel in the training set displayed an AUC of 1.000 (95% CI 1.000–1.000). The markers were validated in an independent set of samples (67 cases, 73 controls) for testing the accuracy of the markers. In the test set, univariate ROC analysis showed an AUC of 0.994 (95% CI 0.9873–1.002) for SOSTDC1, AUC of 0.983 (95% CI 0.9687–0.9976) for DACT2 and AUC of 0.995 (95% CI 0.9902–1.001) for WIF1. The combined AUC of 3 genes was 1.000 (95% CI 1.000–1.000) (Supplementary Fig. S5). The overall performance of the markers in single and in combination were concordant in the training and test sets. Hence, individually the markers indicated excellent discriminatory power and when combined, the 3 markers were able to detect BC with 100% efficiency.

To confirm the accuracy of the markers, an arbitrary cut-off of ‘1’ was assigned and leave one out analysis was performed using the training set. In cases, if M/U (ratio of methylated allele and unmethylated allele) value is greater than 1, a score of ‘1’ was assigned and if M/U value is lesser than 1, a score of ‘0’ was assigned. In controls, if M/U value is greater than 1, a score of ‘-1’ was assigned and if M/U value is lesser than 1, a score of ‘0’ was assigned. The diagnostic parameters were calculated using MEDCALC® tool. In the training set, individually all three markers showed > 90% sensitivity and specificity, with an accuracy of 94% for SOSTDC1 (91.17- 96.60), 94% for DACT2 (90.42–96.10) and 97% for WIF1 (93.88–98.26). The combination of DACT2 with SOSTDC1 displayed a sensitivity and specificity of 100% (97.59–100.00) and 91% (85.53–94.85); DACT2 with WIF1 showed a sensitivity of 100% (97.59–100) and specificity of 91% (86.25–95.31). The 3 gene signature showed a sensitivity and specificity of 100% (97.59–100.00) and 88% (82.70–92.97). In the test set, SOSTDC1 and WIF1 could discriminate BC cases with approximately 95% sensitivity and specificity, but DACT2 exhibited only 85% sensitivity. All the markers individually displayed > 90% specificity in the test set, in concordance with the training set data. In the multivariate analysis, a combination of SOSTDC1, DACT2 and WIF1 yielded 100% (94.64–100.00) sensitivity, 89% (79.80–95.22) specificity and accuracy 94.33% (89.13–97.52) (Table 5).

Table 5 Discriminant analysis of single gene and multi-gene biomarkers in breast cancer.

Thus, combining multiple markers improved the overall sensitivity of the predicted diagnostic model.

Discussion

In the last few years several attempts have been made to develop a non-invasive assay for the early diagnosis of BC. In recent times, the use of multiple markers for BC detection is gaining interest because, the combination of functionally important genes and proteins involved in breast tumorigenesis, have exhibited better diagnostic accuracy than single candidate markers30. A 3-protein diagnostic model named Mastocheck was able to predict BC with a sensitivity, specificity of 71.58%, 85.25% and AUC 0.832339,40. A tumour suppressor methylation panel (APC, BIN1, BMP6, BRCA1, CST6, ESR-b, GSTP1, P16, P21 and TIMP3) showed 91.7% sensitivity in distinguishing BC patients from normal individuals41. A 6 epimarker panel (SFN, P16, hMLH1, HOXD13, PCDHGB7, RASSF1a) discriminated BC and benign cases from controls with 78% sensitivity42 and combinatorial methylation panel comprising GSTP1, RARβ2, RASSF1 and APC detected early-stage tumours with a sensitivity of 33%. In addition, SNiPER gene panel showed a sensitivity of 63% for DCIS and 51% early invasive tumour detection at 80% specificity43. Apart from these, attempts at targeted sequencing of ubiquitously methylated CpG sites have made it possible to detect and classify multiple cancer types from plasma cfDNA21,22.

In the present study, we have aimed at developing a minimally invasive, multi-marker BC diagnostic model with improved sensitivity and specificity in South Asian women. Transcriptomic profiling followed by qRT PCR validation identified 56 candidate markers for biomarker evaluation. Fifteen of the 56 genes were then analysed at the protein level using a multiplexed sandwich ELISA platform (Quantibody array), first in tissues and later in plasma. Eight markers (FGF2, IL17B, IP10, MIG, MIP1d, LOX1, OPN and SDC1) found to be upregulated in primary tumours were present in higher concentration in plasma of breast cancer patients. Similarly, 3 markers (CFD, LEP and DKK3) were downregulated both in tumour lysates and plasma from patients with cancer. Hence, the expression of 11 out 15 markers showed similar trend in tumour tissues and in plasma. IGF-1 induces Ras-Raf-MAPK and PI3K-AKT signalling components to promote tumour progression and elevated levels of serum IGF-1 was associated with poor prognosis44,45,46. Consistent with the existing literature, IGF-1 mRNA and protein levels were found to be higher in BC in our study population. The mitogenic growth factors FGF1 and FGF2, known to possess potent angiogenic properties47,48, were increased in patients with malignant breast tumours and decreased in patients with benign breast lesions indicating their possible role in malignant cell transformation. IL17B, is produced by induction of memory T lymphocytes and plays an important role in inflammatory responses by binding to the membrane receptor IL17RB. IL17B-IL17RB interaction triggers the activation of NF-kB signalling cascade leading to production of anti-apoptotic Bcl-249. IL17B levels were significantly higher in our DCIS and invasive mammary carcinoma patients and diminished in benign and control groups. The protein SDC1 was also significantly increased in our BC cases, on par with previous findings which have demonstrated that high serum levels of SDC1 were associated with aggressive phenotype, poor prognosis and decreased response to chemotherapy50,51,52. In majority of studies, increased levels of CFD and LEP have been linked to obesity associated BC progression via enhanced TGFβ signalling and MMP modulation53,54,55. In contrast, our results indicate decreased levels of the adipokines, CFD and LEP in BC tissues as well as plasma when compared to controls but it is similar to a transcriptome profiling study in Arabian women that has reported the downregulation of leptin and other downstream leptin metabolism genes in BC56.

The dysregulation of WNT signalling cascade in tumorigenesis has been well documented and WNT antagonists sFRP3, DKK3 and WIF1 are frequently downregulated in BC contributing to constitutive activation of oncogenic growth factors35. Our data depicts decreased plasma levels of DKK3 in cases which correlate with the previous gene expression data but the median protein levels of WIF1 were higher in BC tissues and plasma, contrary to our microarray results. Inter-individual variability of WIF1 was observed to be high in our data, and WIF1 levels were below the detection limit in 30–40% of the plasma samples, hence the median and normal range was difficult to establish making it less reliable. Univariate ROC analysis of the markers showed limited sensitivity and specificity, so markers which had AUC above or below 0.5 were chosen for multivariate analysis. Combination of 6 proteins including CFD, IL17B, SDC1, DKK3, FGF2, LEP showed a sensitivity of 65%, 70% and specificity of 79% and 85% in training and test sets respectively. The multi-protein panel possessed higher discriminant performance than single markers but did not achieve the desired sensitivity and specificity.

Next, we assessed the methylation status of SOSTDC1, DACT2 and WIF1 in cfDNA samples isolated from reserved aliquots of plasma from cases and controls recruited. All three genes are known negative modulators of the canonical Wnt signalling cascade. The deficiency of SOSTDC1 correlated with greater tumour size and treatment with recombinant SOSTDC1 effectively blocked WNT signalling components which contribute to cell proliferation, indicating its antagonistic role against Wnt pathway57. In addition, our previous study reported that SOSTDC1 is downregulated in 98.1% of breast tumour tissues and it coincides with DNA methylation of CpG sites in promoter region of SOSTDC136. Epigenetic inactivation of WIF1 contributes to constitutive activation of WNT signalling pathway in breast tumorigenesis58. The differential methylation of WIF1 was also able to predict the clinical efficacy of neoadjuvant chemotherapy (docetaxel, pirarubicin and cyclophosphamide) in sera of locally advanced BC patients59. DACT2 acts as a tumour suppressor by inhibiting canonical WNT signalling and is frequently silenced by promoter hypermethylation in BC. Overexpression of DACT2 inhibited the expression of β-catenin target genes associated with tumour growth and metastasis60,61. Although these 3 genes were known to be frequently methylated in breast tumour tissues, their methylation status in circulating DNA remained unknown. Therefore, SOSTDC1, DACT2 and WIF1 seemed candidate markers for the development of a potential methylation-based diagnostic model.

Our data showed reduced mRNA levels and methylation of the 3 genes in breast cancer cell lines and tumour tissues consistent with previous studies which have reported that promoter methylation mediated silencing of WIF1 and DACT2 was observed in 63.3% (95/150)62 and 73% (107/147)33 of primary BC tissues.

BSP analysis was then carried out to confirm the methylation status in a methylation-independent manner. Consequently, the putative promoter methylation of the 3 markers was analysed in cfDNA isolated from plasma of patients diagnosed with either non-invasive or invasive BC, benign breast diseases and healthy individuals. The markers of interest were found to be methylated in more than 90% of invasive and 75% of pre-invasive BC cases. They were either unmethylated or weakly methylated in more than 80% of the benign cases and 95% of healthy controls. The clinicopathological correlation revealed that rise in methylation intensity of the said genes were significantly associated with advanced tumour stage, high grade, nodal status, and metastases. These results signify the putative tumour suppressive role of the markers and possible role in progression and severity of BC.

The receptor status of breast tumours did not show any influence on the methylation levels except for DACT2 which displayed significant variation among tumours which were HER2 positive and negative. The HER2 receptor status of 50% of the patients had a score of 2 + and FISH confirmation was not done for all due to financial issues (Supplementary Table S7). Diagnostic evaluation of the single markers displayed a high sensitivity range of 92%-98% and specificity of 95% in the training set. In the test set, the sensitivity and specificity of the single markers ranged from 85 to 95% and 92–99% respectively. Combination of 2 or 3 biomarkers generated 100% sensitivity and 89–91% specificity in both the training and test datasets. Hence, the predicted diagnostic model was robust and reliable since it could discriminate BC (non-invasive and invasive) from plasma obtained from patients with benign breast disease and healthy controls with a better sensitivity than the previously proposed models.

Despite, the promising results of various biomarker panels reported previously, the current study has some strengths. Although multiple markers show better diagnostic potential than single markers, it is imperative to establish an assay comprising of minimal number of markers with high diagnostic ability. In addition, the assay platform should also be cost-effective and easily adaptable in a clinical setting. In this study we were able to achieve 97–100% sensitivity and 89–91% specificity by combining 2 or 3 genes. We have used MSP, a PCR based semi-quantitative platform which is simple and cost-effective compared to sequencing and mass spectrometric approaches. Pre-analytical parameters such as use of plasma instead of serum [plasma is preferred for cfDNA assays since serum cfDNA has increased genomic DNA levels63,64 and a three-step centrifugation step to avoid lymphocyte contamination was incorporated. Unlike many of the existing studies the specificity of our panel was also tested in benign breast abnormalities and DCIS and validated in an independent sample set. The differential methylation of the said epigenetic markers, between the benign, DCIS and invasive groups suggest that these markers shall be able to differentiate benign lesions from malignant. Although the epigenetic silencing of WIF1, DACT2 and SOSTDC1 in BC has been reported before, we are the first team to evaluate the combined diagnostic efficiency of these genes in circulating cfDNA from BC patients.

The shortcoming of our study is the limited number of DCIS and early-stage BC samples which can contribute to the lack of statistical significance in the differential methylation intensity between controls, benign and early-stage tumours (DCIS and stage I). Hence, validation of these markers in a larger number of DCIS and early-stage breast tumour patients is necessary before it is considered as a complementary tool for BC diagnosis.

Materials and methods

Patients and samples

The study was approved by the Institute Ethical committee, Cancer Institute (WIA). The study population primarily consisted of Indian women. The patient samples were collected after informed consent and all methods were performed in accordance with relevant guidelines and regulations set by the committee. Our Institute’s Tumour bank provided 6 apparent normal (AN) samples; 41 breast tumours (T) with at least 70% tumour cells in the sample provided as confirmed by frozen section and 18 paired normal (PN) samples that were included in the microarray study and for qRT-PCR validation of the microarray data. For the validation of the markers identified, the protein levels were initially estimated in BC tissue lysates comprising of 6 AN, 23 PN and 46 T samples (Supplementary Table S1). The frequency of women opting for reduction mammoplasty or for prophylactic mastectomy is very rare in India, hence it is difficult to procure normal breast tissue samples. Our work hence used samples well away from the benign lesions and confirmed to be morphologically normal, from women with benign breast disease.

Further, an age-matched case–control study (age distribution with 5-year intervals) was performed with 202 BC patients and 203 healthy controls. The inclusion and exclusion criteria are provided in (Supplementary Table S6). Additionally, 16 patients with ductal carcinoma in-situ (DCIS) and 37 patients with benign breast disease were included in the study. The samples were randomly split into training set and test set for validation (Supplementary Fig. S1).

Sample preparation

The tissue lysates were prepared according to Rajkumar et al.65. The concentration of proteins was estimated in 100ul of lysate (1 mg/ml). Post clinical examination, 10 ml blood sample was collected in ethylene diamine tetra acetic acid (EDTA) coated tubes from each individual. The cells and plasma were separated from whole blood by centrifugation at 2500 rpm for 20 min at 37℃ and stored in aliquots in LoBind tubes [Eppendorf, India] at − 80 °C. The plasma samples were thawed and centrifuged at 13,000 rpm and 100 µl (1:1) of diluted cell-free plasma was added to each well of the custom designed Quantibody array slide (Raybiotech, Inc. USA).

Cell lines

The cell lines HBL100, MCF7, MDAMB231, MDAMB468, ZR751, SKBR3 and T47D were used. All the cell lines were purchased from National Centre for Cell Science (Pune, India).

Nucleic acid extraction

Genomic DNA was isolated from cell lines, tissues samples (breast tumour, paired normal and apparent normal) by QIAamp® DNA Mini Kit (Qiagen, Hilden) according to manufacturer’s instructions. Cell free DNA was isolated from 1 ml of plasma using QIAamp Circulating Nucleic Acid Kit (Qiagen, Hilden). DNA and RNA quantification were done using NanoDrop ND1000 (NanoDrop Technologies, USA) spectrophotometer. Cell free DNA was quantitated using Qubit™ dsDNA HS Assay Kit (Invitrogen, USA).

The RNA was extracted from the tissue samples using the RNeasy RNA extraction kit (Qiagen, Hilden; Cat no: 74106) as per the manufacturer’s instructions. The quality of the RNA used for microarray analysis was checked using the Bioanalyzer and samples with RNA Integrity Number (RIN) of 7 or more were included in the study. RNA was quantitated using NanoDrop ND1000 spectrophotometer (NanoDrop Technologies, USA).

Microarray based gene expression profiling

The microarray experiments were done as described previously65. Briefly, 1 µg of total RNA from the tumour/PN/AN sample and universal RNA (Stratagene; Cat no: 740000–41) were reverse transcribed using Array script at 42 °C for 2 h to obtain cDNA using the Amino Allyl MessageAmp II aRNA amplification kit (Ambion, Austin, Tx; Cat no: AM1797). The cDNA was amplified, labelled, hybridized and slides scanned as described earlier. All the raw data files have been submitted to GEO (accession number GSE139038).

Microarray data analysis

The microarray data analysis was done as described previously using BRB-ArrayTools software v 3.7.0 (http://linus.nci.nih.gov/BRB-ArrayTools.html)65. The differentially expressed genes among the 3 classes (tumour/PN/AN) were analysed using the Class comparison module of the BRB-Array Tools software. Univariate F-test was used and the genes were considered statistically significant if their P-value was < 0.001. In addition, a two-fold difference in gene expression was required between the different classes [T vs. PN vs. AN].

Quantitative real time PCR

Validation of the gene expression was done using the TLDA quantitative real time PCR (Applied Biosystems, Foster City, CA). Triplicate cDNA template samples were amplified and analysed on the ABI Prism 7900HT sequence detection system (Applied Biosystems, Foster City, CA). The protocol for validation was adapted from literature65. The raw data from the Prism 7900HT sequence detection system was imported into Microsoft Excel for statistical analysis of the data. Among the endogenous reference genes included on the array (18S ribosomal gene; UBC; GAPDH), UBC and GAPDH were chosen after visualizing the global Ct value distribution, for data normalization. The TLDA assays were run at Lab India Instruments Pvt Ltd laboratories at Gurgaon, New Delhi.

The AN tissue samples were used as calibrators and the relative quantitation values were calculated for all the genes and the samples. Geometric mean was calculated for each of the 67 genes (excluding the 3 endogenous controls). The relative quantitation values for all the samples and genes were imported into BRB Array Tools and Class Comparison analysis was done comparing the different clinicopathological parameters with the gene expression values. Fifteen genes which were secreted or potentially secreted were short-listed, based on our data as well as corroborated with GEO Microarray Dataset (176 primary breast tumors and 10 normal breast samples [from reduction mammoplasties] [GSE22820]). An additional requirement was that the protein could be assayed in the multiplexed sandwich-ELISA platform.

Custom designed protein antibody array (quantibody array)

Quantibody array [multiplexed sandwich ELISA on a glass slide] was used [Ray Biotech, Inc, USA] to study the protein expression in tissue lysates and subsequently in plasma. The assay was done as per the manufacturer’s instructions with modifications65.

Quantibody array data normalization

The data was analysed using the H15S90 Genxbio Q-Analyzer v8 16.4, an array specific, Microsoft Excel based program, supplied with the custom arrays. 2 positive controls for signal normalization, a negative control for background subtraction and an internal control to minimize inter-slide variation were used. A user defined reference array was used, to which the signals of other arrays were normalized. Median + Median absolute deviation was calculated for each protein and the cut-offs which yielded a sensitivity and specificity of more than 75% was chosen.

cDNA conversion and semiquantitative RT-PCR

1 µg of RNA was used to synthesize first strand cDNA (QuantiTect Reverse Transcription Kit, Qiagen, Hilden). The reaction mixture was diluted 50% (v/v) with water and 2 µl of cDNA was used for 25 µl PCR reaction (Eurogentec Takyon™ Low ROX SYBR® 2X Mastermix Blue). RT-PCR was amplified for 35 cycles and GAPDH was used as an internal control for normalization. The products were visualized in 2% agarose gel.

Bisulfite conversion, methylation specific PCR (MSP) and bisulfite sequencing PCR (BSP) analysis

Bisulfite conversion was performed using the EZ DNA Methylation-Gold Kit (Zymo Research, USA), according to manufacturer’s instructions. The amount of input DNA was adjusted to be uniform for each sample. For MSP and BSP analysis of Dapper 2 homolog (DACT2), genomic sequence 1000 bp upstream of transcription start site (TSS) and 200 bp downstream of TSS containing the putative promoter region was retrieved from Database of transcriptional start sites (http://dbtss.hgc.jp/). MethPrimer 2.0 (http://www.urogene.org/methprimer2/) was used for CpG island prediction and generation of primer sets. For Sclerostin domain containing 1 (SOSTDC1)36 and Wnt Inhibitory Factor 1 (WIF1)62, pre-designed primer sets (spanning promoter CpG islands) from previous studies were used (Supplementary Table S8). HotStar Taq® Master Mix (Qiagen, Hilden, Germany) was used for MSP and the conditions were set according to manufacturer’s instructions. The amplicons were resolved in 2% agarose gel.

BSP products of WIF1 and DACT2 were amplified from bisulfite modified MDAMB231 genomic DNA using HotStar Taq® with gene-specific primers (Supplementary Table S8) and cloned using TOPO® TA Cloning® Kit for Sequencing (Invitrogen, Massachusetts USA). The clones were screened using PCR and restriction enzyme digestion. The clones were sequenced using a Big Dye Terminator Cycle Sequencing kit (Applied Biosystems, Foster city, CA, USA) as per the manufacturer’s instructions on ABI 310 genetic analyser. Densitometric analysis was used to quantify (Image Lab 6.0.1.) the intensity of the DNA bands relative to 100 bp DNA ladder (Promega, Madison USA). Ratio of methylated and unmethylated band intensity was calculated for each sample in cases and controls. Median + Median absolute deviation was calculated for each gene and used as cut-offs.

Equipment and settings

The Quantibody array slides were scanned at 5 μm resolution, ratio (635/532) and a PMT of 580 using Molecular Devices GenePix Pro 4100 A scanner. At these scanner settings the signals from the highest standard concentration did not reach saturation.

The MSP amplicon resolved agarose gels were docked using Biorad ChemiDoc MP system. Image Lab version 4.1 software (Bio Rad CA, USA, https://www.bio-rad.com/en-in/product/image-lab-software)66 was used for gel imaging with acquisition settings such as single channel, SyBR Safe mode and gray scale image colour. The image exposure and image area were set as default. For relative quantification of band intensity, Promega 100 kb DNA ladder was used as reference standard.

Statistical analysis

Statistical tests were done using XLSTAT version 2018.1, (Microsoft, WA, USA) https://www.xlstat.com67 and Prism v5.01 (GraphPad software, CA,USA, https://www.graphpad.com). Kruskal Wallis Test was employed to test significance of differential protein levels and methylation intensity. Unpaired t-test and One-way analysis of variance was used to assess association of relative methylation levels of BC patients with clinicopathological parameters. Statistical significance was defined as P < 0.05. The ROC curve analysis and binomial logistic regression was performed using IBM SPSS for Windowsversion 28 (IBM Corp, NY, USA https://www.ibm.com/in-en). The multivariate ROC was generated using Predicted Probabilities module in Binomial Logistic Regression.

Ethics declarations

The study was approved by the Institute’s Ethical committee (IEC), Cancer Institute (WIA) by their letter dated 16–06–2007 (for use of tissue samples from Tumour Bank) and letter dated 18–08–2014 (for case – control study). Patient samples were collected after informed consent and the study was performed in compliance to the National Ethical Guidelines for Biomedical Health Research involving Human Participants and conditions of Indian Council of Medical Research.