Collecting and preprocessing data

RNA sequencing data of 111 LSCC samples and 12 normal laryngeal tissues were obtained from the TCGA (The Cancer Genome Atlas) database (https://portal.gdc.cancer.gov). Clinical information and RNA sequencing data of LSCC samples were obtained from the GSE27020 and GSE51985 Gene Expression Datasets (https://www.ncbi.nlm.nih.gov/geo/). Gene expression data were normalized using the "SangerBox" tool (http://sangerbox.com/) before further analysis.

Identification of Differentially Expressed Genes (DEGs)

In total, 334 VM-related genes were obtained from the GeneCards Database (https://www.genecards.org). Using the "DESeq2" package, differentially expressed genes (DEGs) related to apoptosis were identified in the 111 LSCC and 12 normal samples. The threshold was set at P < 0.05, and log2 fold change (FC) > 0.585. Using the R "survival" and "survminer" packages, univariate analysis was conducted on DEGs related to overall survival (OS). The correlation of prognostic genes was further analyzed using the R packages "igraph," "psych," "reshape2," and "RColorBrewer."

Construction and validation of a VM-related genes(VMRGs) prognostic model

Based on 334 VMRGs, we screened 141 VM-related differential genes (VMDEGs) from normal and LSCC samples and used these 141 differential genes to construct a prognostic model. LASSO Cox regression analysis and least absolute shrinkage were applied to 141 VMDEGs to construct a prognostic model. Cross-validation (CV) algorithm was used for model selection criteria and validation. A minimum mean cross-validated error rule and a 10-fold cross-validation approach were used to select the penalization parameter λ in LASSO. We calculated the risk score using the following formula:

$$\mathrm\;\mathrm=\sum_^ncoef\left(^K\right)\;\times\;\exp\;r\left(^K\right)$$

where n is the number of risk genes, coef is the regression coefficient in the multivariate Cox regression analysis (VMRGs), and expr is the expression of the risk-associated genes (BMP2, EPO, and AGPS). The median risk score was used to divide patients with LSCC into high- and low-risk groups, followed by OS assessment. The R package "timeROC" was used to plot time-dependent receiver operating characteristic (ROC) curves and areas under the curves (AUCs). We validated the model using an LSCC cohort from the GEO database to enhance its credibility. We then determined whether the risk score contributed to overall survival (OS) and progression-free survival (PFS) in the cohort using multivariate and univariate Cox regression analyses. Covariates included age, sex, and TNM stage.

Functional enrichment and Tumor Mutational Burden (TMB) analyses

Functional enrichment analyses (GO, KEGG, and GSEA) were performed using the "clusterProfiler" package as previously reported [26]. In addition, Gene Set Variation Analysis (GSVA) was conducted using the "GSVA" R package to further identify differentially enriched KEGG pathways. Using the R package "MAfTools," we calculated the TMB of each tumor based on the number of mutations per million bases calculated from the somatic mutation data.

Cell culture

HaCaT, TU212, and TU686 cells were purchased from iCell Biotechnology Co., LTD. (Shanghai, China). All three cells had STR testing certificates (Fig. S1). According to Mycoplasma Detection Kit (#G238; abm, Canada), all cell lines were negative for mycoplasma (Fig. S2). TU212 is a laryngeal squamous cell carcinoma cell line and TU686 is a hypopharyngeal carcinoma cell line. TU212 and TU686 cells were cultured in DMEM (Gibco, Carlsbad, CA, USA) supplemented with 10% fetal bovine serum (Gibco) and 1% penicillin/streptomycin (Sigma-Aldrich) at 37 °C in a humidified atmosphere containing 5% CO2. Prior to subsequent experiments, cells were treated with MK2206 (5um), a highly selective Akt inhibitor, for 12 h.

Western Blotting (WB) and immunohistochemistry (IHC)

Total protein was extracted using a protein extraction buffer (#GK10023; GLPBIO, Montclair, CA, USA) containing protease (#GK10014; GLPBIO) and phosphatase (#GK10013; GLPBIO) inhibitors. Protein samples were then separated by sodium dodecyl sulfate–polyacrylamide gel electrophoresis and transferred to PVDF membranes (Millipore, Billerica, MA, USA), which were then blocked by incubation with 5% BSA for 1 h. Subsequently, membranes were incubated with a blocking buffer containing the primary antibody overnight at 4 °C, washed in TBST, and then incubated in TBST containing the HRP-conjugated secondary antibodies (M21008, Abmart, Shanghai, China). A Bio-Rad ECL Western Blotting Detection System (Bio-Rad, Hercules, CA, USA) was used for visualization. For measuring molecular weight, PierceTM unstained protein MW Marker (Thermo Fisher, 26610; UElandy, P6110, Suzhou, China) was used. The antibodies used were as follows: BMP2 (dilution 1:2000, #66383–1-Ig; ProteinTech, Wuhan, China), ACTB (dilution 1:3000, #81115–1-RR; ProteinTech), RUNX1 (dilution 1:1000, #25315–1-AP; ProteinTech), VE-cadherin (dilution 1:500, #A22659; Abclonal, Wuhan, China), N-cadherin (dilution 1:500, #A19083; Abclonal), E-cadherin (dilution 1:500, #A3044; Abclonal), AKT (dilution 1:1000, #4691S; CST, Danvers, MA, USA), p-AKT (dilution 1:1000, #4060S; CST), p-PTEN (dilution 1:1000, #9551S; CST), and p-PDK1(dilution 1:1000, #3438S; CST).

An LSCC tissue microarray (HN049La01) was purchased from Zhongke Guanghua Intelligent Biotechnology Co., Ltd. (Xian, China). IHC was performed as previously described [27], using an immunohistochemistry kit (Zhongshan Goldenbridge Biotechnology, Beijing, China). Tumor tissue sections were incubated with the primary antibody overnight at 4 °C. Next, sections were rinsed using phosphate-buffered saline (PBS) and then incubated with horseradish peroxidase (HRP)-labeled IgG for 60 min, followed by the addition of DAB (Zhongshan Goldenbridge). Finally, sections were counterstained with hematoxylin and eosin (H&E). The staining intensities of sections were independently assessed by two pathologists. The following antibodies were used: BMP2 (dilution 1:50), RUNX1 (dilution 1:150), VE-ca (dilution 1:150), N-ca (dilution 1:150), and E-ca (dilution 1:150).

Plasmids and shRNA constructs

The BMP2-overexpressing plasmid (pCDNA3.1-BMP2), shRNAs targeting BMP2 or RUNX1, and their respective negative control RNAs were purchased from Tsingke Biotech Co., Ltd. (Beijing, China). Transfection was performed using Lipofectamine 2000 (Invitrogen, Carisbad, CA) with a plasmid-to-transfection reagent ratio of 1:2. Using the pLKO.1-CMV-LUC-PURO vector, Lv-shBMP2 and Lv-NC cells were constructed. Lentiviral transfection was performed according to the manufacturer’s instructions. Fresh complete medium was added to cultures 2 d after infection. Cells were selected using puromycin (3 μg/mL, Invitrogen) for 7 d to obtain stably transfected cells.

The shRNA sequences used were as follows: shBMP2-1: 5′-CCGGCCGGAGATTCTTCTTTAATTTCTCGAGAAATTAAAGAAGAATCTCCGGTTTTTT-3′; shBMP2-2: 5′-CCGGGATCATCTGAACTCCACTAATCTCGAGATTAGTGGAGTTCAGATGATCTTTTTT-3′; and shRUNX1: 5′-CCGGCCTCGAAGACATCGGCAGAAACTCGAGTTTCTGCCGATGTCTTCGAGGTTTTTT-3′.

Quantitative real-time PCR (qRT-PCR)

TRIzol reagent (Invitrogen) was used to extract total RNA, which was then reverse transcribed into cDNA using a reverse transcription system (Takara, Tokyo, Japan). qRT-PCR was performed using the SYBR Green Mastermix (Takara) with the CFX96 Touch qRT-PCR System (Bio-Rad). For relative quantification, the mRNA levels of target genes were normalized to those of β-actin. The primer sequences used were as follows: β-actin: forward sequence, 5′-CTCCATCCTGGCCTCGCTGT-3′; reverse sequence, 3′-GCTGTCACCTTCACCGTTCC-5′. BMP2: forward sequence, 5′-CGCTGCCCAGAGGACTTC-3′; reverse sequence, 3′-GACTTTAGCGGTCTCGGAGC-5′. RUNX1: forward sequence, 5′-CTGCCCATCGCTTTCAAGGT-3′; reverse sequence, 3′-GCCGAGTAGTTTTCATCATTGCC-5′.

Wound healing and Transwell assays

The migratory ability of cells was evaluated using a wound healing assay. A 6-well plate was seeded with shBMP2-transfected TU686 and TU212 cells. After reaching 100% confluence, the cell monolayer was wounded with a 10 μL pipette tip and then cultured in a basic medium. At specified times, the wound area was quantified using the ImageJ software(Invitrogen (Carisbad, CA). Each experiment was repeated thrice. Two types of Transwell assays were used: migration and invasion. Using Transwell inserts (#3422; Corning, NY, USA) with an 8-µm-pore size, we added 5 × 104 cells resuspended in basic medium to the upper chamber, whereas 10% fetal bovine serum was added to the lower chamber. Transwell invasion assays involved the use of 60 µL Matrigel (#356234; Corning) precoated onto the upper membrane. Following an 18 h migration or 36 h invasion assay, cells were fixed with 4% paraformaldehyde(P8430, Solarbio, Beijing, China) and stained with 0.1% crystal violet (G1061, Solarbio). Digital images of membranes were obtained by capturing three random fields in each chamber.

Spheroid invasion assay

For the spheroid invasion assay, spheroids were formed using a previously described method [26]. After 3d, they were embedded in a gel prepared using a 1:1 mixture of Matrigel and DMEM/F-12 medium. After 2 d, the total invasion area was divided by the central spheroid area, and measurements were taken using ImageJ.

Matrigel tube formation and plug assays

The Matrigel tube formation assay was conducted to evaluate the ability of TU212 and TU686 cells to form vessel-like structures, which is a key event in VM. After dissolving it at 4 °C, 60 μL Matrigel was added to each well in a 96-well plate and allowed to polymerize at 37 °C for 1 h. Subsequently, 4 × 104 LSCC cells were resuspended in 100 μL DMEM/F12 medium and incubated at 37 °C in 5% CO2. The tube formation ability of cells was observed every 3 h. The Cytation5 Cell Imaging Multimode Reader (BioTek, Winooski, VT, USA) was used for image acquisition.

Nude mice (4-week-old) were used for in vivo Matrigel plug assays. Nude mice were subcutaneously injected in the axilla of the forelimb with TU212 cells (5 × 106) mixed with Matrigel (0.2 mL). On the seventh day, mice were euthanized, and the Matrigel plugs were excised. The density of newly formed VM was measured by H&E staining of histological sections.

BALB/c nude mice and zebrafish tumor models

TU212 cells were transfected with Lv-luc-shBMP2 or empty vector, and stable transfectants were selected. Next, 2 × 106 luciferase-labeled TU212 cells in 100 µL DMEM medium were injected into the tail veins of male nude mice aged 5 weeks (n = 4 per group). Over a 5-week course, D-luciferin (150 mg/kg, #GC43496; GLPBIO, Montclair, CA, USA) was intraperitoneally injected into mice, and whole-body photon flux was measured weekly using a Xenogen IVIS imaging system (Caliper Life Sciences, Mountain View, CA, USA) to monitor lung metastasis. To verify the inhibitory effect of the BMP2 knockdown on LSCC growth in vivo, TU212 cells (5 × 107) were resuspended in 200 μL DMEM medium and injected subcutaneously into the right axilla of nude mice. Tumor volume was measured every 3 d. Tumor volume was calculated using the following equation: tumor volume (mm3) = (tumor width)2 × tumor length / 2. At 21 d after tumor implantation, mice were euthanized, and xenografts were removed, fixed, weighed, photographed, and preserved.

The zebrafish tumor model was developed as previously described [22]. Animals were housed at the Zebrafish Centre of the First Affiliated Laboratory of the Air Force Medical University. Using a microinjection system, TU212 cells (200 cells/5 nL) were injected into the yolk sac cavity of zebrafish embryos after being labelled with 2 g/mL DiI (V228885, ThermoFisher, Waltham, MA, USA).

Immunofluorescent (IF) staining

Cells cultured in 10-mm glass bottom dishes were fixed with 4% paraformaldehyde and blocked with 5% BSA at 25 °C for 30 min. Fixed cells were incubated overnight at 4 °C with primary antibodies against E-cadherin (1:200); N-cadherin (1:200); VE-cadherin (1:200, #ab313632; Abcam); BMP2 (1:50); and RUNX1 (1:100). After washing, cells were incubated with fluorescence-labeled secondary antibodies (1:200; Invitrogen) at 37 °C for 30 min. Nuclei were counterstained with Hoechst solution (1:1000, #33258; Beyotime). Images were captured using a FV3000 confocal fluorescence microscope (Olympus, Center Valley, PA, USA).

Chromatin immunoprecipitation (ChIP)

A ChIP assay kit (17–295; Merck, Darmstadt, Germany) was used according to the manufacturer’s instructions. Briefly, cells were crosslinked with 1% formaldehyde and sonicated using the Covaris M220 system (Covaris, Woburn, USA) to generate DNA fragments of 100–500 bp in length. After preclearing, the supernatant was incubated with antibodies against RUNX1 (10 g rabbit IgG, ab272456; Abcam) or isotype control antibody (2 g rabbit IgG, ab172730; Abcam). PCR was conducted using primers targeting the BMP2 promoter binding sites (RUNX1 set1, forward sequence: 5′-ATTTCCAGCCTGCTGTTTTCTT-3′; reverse sequence: 3′-CCACTCCCTGCTCTCAAAGGA-5′ and RUNX1 set2, forward sequence: 5′-ACATATTAACCGAAATGTGGCCC-3′; reverse sequence: 3′-GGAAAATTAAAAGAAAACAGCAGGC-5′). DNA isolated from the total nuclear extract was used as a control. PCR products were run on 2% agarose gel.

Dual-luciferase reporter assay.

The dual-luciferase reporter system (Dual-Luciferase Reporter Assay; Promega, Madison, WI, USA) was used. The sequence of the BMP2 promoter region (from the transcription start site: -2000 bp to + 100 bp) was obtained from the UCSC database (Table S1). Transcription factors associated with BMP2 were predicted using the AnimalTFDB (AnimalTFDB4 (hust.edu.cn)), GTRD (http://gtrd20-06.biouml.org/), and JASPAR (https://jaspar.genereg.net/) databases. Potential binding sites for RUNX1 within the BMP2 promoter region were predicted using JASPAR (https://jaspar.genereg.net/). The pcDNA3.1-RUNX1 and pcDNA3.1-NC vectors were constructed, and the wild type (WT) and mutant (MT) sequences (Table S2) of the promoter region were cloned into the pGL3 vector (TsingKe Biotech Co., Ltd.). The promoter-specific and pcDNA3.1 constructs were cotransfected into TU212 cells with the control Renilla luciferase construct PRL-TK, using the Hieff TransTM reagent (YEASEN, Shanghai, China). Luciferase and Renilla luciferase assays were performed for normalization.

Statistical analysis

All data are expressed as the mean ± standard deviation. Student’s t-test was conducted using the GraphPad Prism (GraphPad Prism v9.0; GraphPad Software, USA) and R software(R version 4.1.1, R Core Team Vienna, Austria). One-way analysis of variance was used to compare data among multiple groups. A P < 0.05 was considered statistically significant. Statistical significance was defined as follows: *P < 0.05, **P < 0.01, and ***P < 0.001.