The study design of the current study is described in Supplementary Material ESM 1. RNA-seq data of 1,098 breast cancer patients were downloaded from The Cancer Genome Atlas (TCGA) Database¹. The mRNA expression levels of FGFR1 gene was estimated using z-score >2.0. The clinicopathological relevance of FGFR1 expression was assessed against multiple features including age, stage, grade, node stage, and metastasis. The details of clinical information are available in Supplementary Material ESM 2. Gene Set Enrichment Analysis (GSEA) based on FGFR1 expression was performed to identify new FGFR1-associated pathways in breast cancer (Subramanian et al., 2005). Clinicopathological association of those pathways in breast cancer was also evaluated. Additionally, another data (GSE20685) of 327 breast cancer patients were also downloaded to validate our findings in an independent cohort (Supplementary Material ESM 3). Similarly, clinical association of FGFR1 expression and FGFR1-associated pathways using GSEA was investigated. Common pathways identified between the two cohorts were selected for further analysis.

The present study was conducted after obtaining approval from institutional ethical review committees at our university and concerned hospital. Tumor biopsies (n = 150) along with matched control tissues were collected for validation of discovery cohort findings after receiving informed consents of participants. Clinicopathological information of these patients was collected through subsequent follow-up from pathology reports. Information regarding clinicopathological features of validation cohort is available in Supplementary Material ESM 4.

Total RNA was extracted from tumors and matched normal samples using TRIzol^® (Invitrogen, Carlsbad, CA, United States). cDNA was generated using RevertAid First-Strand cDNA Synthesis Kit (Thermo Fisher Scientific, Carlsbad, CA, United States) as per manufacturer’s instructions. Conventional PCR was performed with β-actin primers to confirm the cDNA synthesis. Amplified products were electrophoresed on 2% agarose gel and stained with ethidium bromide for further use.

Primers for FGFR1 and GLI1 were designed using Integrated DNA Technology (IDT) software and synthesized for Macrogen, Korea (Supplementary Material ESM 10). Target specificity of these products was confirmed with NCBI Primer Blast to avoid non-specific binding. Quantitative real-time PCR (qRT-PCR) quantitative PCR (qPCR) was performed using VeriQuest SYBR Green qPCR Master Mix (Thermo Fisher Scientific, CA, United States). Expression of target gene was normalized using β-actin as an internal control. The reaction condition included an initial denaturation at 95°C for 15 min, followed by 40 cycles of denaturation at 95°C for 15 s and annealing at 53°C for 1 min in each cycle. Relative mRNA expression and fold change was evaluated using the 2^–ΔΔCt method.

MDA-MB-231 was maintained as per American Type Culture Collection (ATCC). Expression of FGFR1, SHH, and GLI1 were assessed both at RNA and protein level. Similarly, the effect of these aforementioned molecules was evaluated on CRISPR/Cas9-mediated SHH knockout MDA-MB-231 cells (Riaz et al., 2019). Moreover, commercially available small interfering RNAs (siRNAs) targeting SHH (hs.Ri.SHH.13.1) and GLI1 (hs.Ri.GLI1.13.2) were purchased from IDT², and MDA-MB-231 cells were transfected with them for assessing the differential expression of FGFR1, SHH, and GLI1.

A stock of GANT61 (cat: G9048, Sigma) was dissolved in dimethyl sulfoxide (DMSO) maintaining a stock concentration of 1 mM. Briefly, 3 × 10⁵ cells were seeded in six-well plates until confluency and treated with variable concentrations of GANT61. Similarly, AZD4547 (Selleckchem³) was added to the medium at a concentration of 100 nM. Both RNA and protein were quantified from treated and untreated wells for functional assays performed after 48 h of treatment with both inhibitors.

Total protein content from respective cell lines were extracted and quantified using the Pierce BCA protein assay kit (23225, Thermo Fisher Scientific, United States). Extracted proteins were separated on 10% sodium dodecyl sulfate–polyacrylamide gel electrophoresis (SDS-PAGE), transferred to nitrocellulose membrane, and blocked by 5% non-fat milk at room temperature. Membranes were incubated with primary antibodies for FGFR1 (1:1,000), SHH (1:1,000), and GLI1 (1:1,000) (Supplementary Material ESM 11). After overnight incubation at 4°C, membranes were incubated with secondary antibodies at room temperature for 1 h (Supplementary Material ESM 11). Protein signals were visualized using ECL Prime Western Blotting Detection Reagent (GE Healthcare Japan) with β-actin as loading control.

Briefly, 3 × 10⁵ cancer cells from respective cancer cell lines were seeded in six-well plates. The plate was left at 37°C incubation till it reaches confluent monolayer. Once attained, the medium was aspirated, and N-2-hydroxyethylpiperazine-N′-2-ethanesulfonic acid (HEPES) medium was introduced in each well. Wounding measurements were recorded using the previously mentioned protocol (Riaz et al., 2019).

This assay was based on Boyden chamber using inserts (8 μm) placed in 24-well plate. These inserts were precoated with 50 μg/ml Matrigel (BD Biosciences, Berkshire, United Kingdom) prior to cell seeding. A total of 5 × 10⁴ cancer cells were seeded in each insert placed in the respective wells. The plate was left at 37°C incubation for 24 h. After specified time duration, these inserts were fixed with methanol and stained with crystal violet. Cells were counted under light microscope at 40 × magnification as per protocol stated earlier (Malik et al., 2009).

IBM SPSS 21 software (IBM Inc., Armonk, NY, United States) was used for all the statistical analyses. Data were represented as mean ± SD, and Wilcoxon signed-rank test was performed to evaluate difference between tumor and control. Mann–Whitney U test and Kruskal–Wallis ANOVA were applied to explore any probable association of these genes with clinicopathological features. Furthermore, correlation between molecules was observed using Spearmen test. All p-values of 0.05 were considered statistically significant.

Initially, in TCGA dataset of 1,098 patients, clinical relevance of FGFR1 expression with breast cancer patients was assessed. FGFR1 was significantly overexpressed in late tumor stage (p = 0.05) and node-positive patients (p = 0.04) (Table 1).

Next, GSEA was performed to identify common differentially expressed pathways in both breast cancer cohorts. A total of 20 pathways were found to be upregulated along with FGFR1 overexpression in TCGA dataset (Supplementary Material ESM 5). Similarly, a total of 14 pathways were found to be upregulated along with FGFR1 overexpression in the GEO20685 cohort (Supplementary Material ESM 6). Out of all pathways, three pathways including ALK, SHH, and PRION were common in both datasets (Supplementary Material ESM 7). Since SHH pathway dysregulation is reported as an early event in several breast cancer studies, we selected SHH pathway to evaluate the association of FGFR1 and SHH pathway in modulating breast carcinogenesis.

Interestingly, 5 out of 16 SHH pathway genes showed core enrichment in FGFR1-expressed breast cancer patients using leading-edge subset method (Figure 1 and Supplementary Material ESM 8). Of note, all GLI family genes including GLI1, GLI2, and GLI3 were significantly associated with FGFR1 in both datasets. Therefore, to further establish the link of FGFR1 overexpression with GLI genes, clinicopathological association of genes including FGFR1, GLI1, GLI2, and GLI3 were evaluated (Supplementary Material ESM 9). According to the results, GLI1 gene was the most frequently associated gene with the poor prognostic features of breast cancer patients including late stage (p = 0.0001), node positive (p = 0.044), and metastasis (p = 0.022) (Table 1 and Supplementary Material ESM 9). Similarly, in the second dataset of 327 breast cancer patients, GLI1 overexpression was the only gene in the SHH pathway that showed significant associations with late stage (p = 0.047). Therefore, based on all the expression and statistical analyses in discovery cohorts, we further established the potential prognostic association of GLI1 and FGFR1 genes using normal and tumor pairs of 150 breast cancer patients.

The demographics and clinical characteristics of the cohort exhibited that the mean age of breast-cancer-affected patients included in the study was 45 years, ranging from 23 to 75 years. According to Wilcoxon test, both FGFR1 (p = 0.0001) and GLI1 (p = 0.0001) were significantly overexpressed in tumor samples compared to their respective controls. Interestingly, FGFR1 and GLI1 expression showed strong positive correlation (Spearman’s correlation = 0.513, p = 0.0001). Of note, consistent with discovery cohort findings, overexpression of GL1 and FGFR1 genes was significantly associated with high grade (p = 0.05), late stage (p = 0.011, p = 0.004), and metastasis (p = 0.006, p = 0.04), respectively (Table 2). These findings suggest strong biological and prognostic relevance of GLI1 and FGFR1 expression in modulating subtypes of breast cancer progression.

Next, the effect of GANT61 (GLI1 inhibitor) and AZD4547 (FGFR1 inhibitor) on both GLI1 and FGFR1 was studied in MDA-MB-231 cell line. To perform the analysis, multiple genes including SHH, AKT, and MAPK were also selected to compare the effect. According to the results, a significant decrease in the expression of GLI1, SHH, and FGFR1 was observed after treatment with GANT61 at 24 and 48 h. MAPK expression was reduced after 48 h of treatment. However, no change in the expression of AKT gene was observed (Figure 2A). Interestingly, a significant downregulation of FGFR1 expression was observed compared to other genes when cells were treated with AZD4547, GLI1 knockdown (SiGLI1), SHH knockdown (SiSHH), and SHH knockout (SHHKO), suggesting a crosstalk between SHH pathway and FGFR1 activation (Figure 2B).

Next, using Western blot analysis, the impact of GANT61, AZD4547, SiGLI1, and SiSHH was assessed in MDA-MB-231 cells on SHH, GLI1, and FGFR1 expression. First, GANT61 treatment was administered until 48 h with 12-h intervals. A significant decrease in protein expression was observed after 12 h for all three genes including SHH, GLI1, and FGFR1 (Figure 2C). In addition, AZD4547, SiGLI1, SiSHH, and CRISPR/Cas9-mediated SHH knockout in MDA-MB-231 breast cancer cells was also assessed (Figure 2D). Consistently, a similar pattern of decrease in protein expression was observed for SHH, GLI1, and FGFR1, suggesting the synergistic role of FGFR1 and SHH pathway. Furthermore, we also evaluated the effect of SHH and GLI1 expression in mouse embryonic fibroblasts (MEFs) with MEFFGFR WT and MEFFGFR1 KO cells. Significant downregulation of SHH and GLI1 was observed in the FGFR1 knockout mouse MEFs, further indicating crosstalks between FGFR–SHH pathways (Figure 3).

Furthermore, the effect of FGFR1–SHH pathway signaling on the motility of breast cancer cells was evaluated using wound-healing assay. MDA-MB-231 cells were treated with AZD4547 only, GANT61 only, and AZD + GANT61 combined (Figure 4). Breast cancer cell migration was evaluated after 48 h of treatment. Of note, invasion and migration of MDA-MB-231 decreased significantly upon AZD + GANT61 treatment compared to AZD4547 or GANT61 treatment alone. Therefore, the data suggest that targeting FGFR1–GLI1 simultaneously significantly reduce cell invasion and metastasis.