Single‐Cell RNA Sequencing Reveals the Critical Role of SEC16B in Lung Metastasis of Osteosarcoma

Abstract

ABSTRACT

Osteosarcoma (OS) is highly malignant and easily prone to lung metastasis. The mechanisms of lung metastasis in OS remain unclear. The single‐cell RNA sequencing (scRNA‐seq) samples in this study included six primary osteosarcoma samples (published in‐house data), two lung metastasis samples (GSE152048), and four normal bone tissue samples (GSE169396). To identify potential targets for metastasis, bulk RNA sequencing data from four primary tumors and four lung metastases (in‐house data) were also analyzed. scRNA‐seq identified five tumor cell subpopulations. CytoTRACE and lung metastasis scores indicated that the C1 subpopulation was most closely associated with lung metastasis. By intersecting lung metastasis‐related genes identified via hdWGCNA analysis with differentially expressed genes from bulk RNA sequencing, SEC16B was identified as the key gene influencing lung metastasis. qRT‐PCR results revealed that SEC16B expression was significantly downregulated in OS cell lines. Transwell assay demonstrated that overexpression of SEC16B significantly inhibited the invasion and migration capabilities of OS cells. Additionally, analyses using Scissor, CellphoneDB, and CSOmap suggested that fibroblasts, endothelial cells, and OS cells in the tumor microenvironment formed a pre‐metastatic niche through mechanisms involving angiogenesis and extracellular matrix remodeling. Overall, this study identifies a new population that may promote lung metastasis by downregulating SEC16B in OS. Moreover, fibroblasts and endothelial cells in the tumor microenvironment play a critical role in OS lung metastasis.

Abbreviations

• bulkRNA-seq
• bulk RNA sequencing
• CNV
• copy number variation
• DEG
• differentially expressed gene
• ECM
• extracellular matrix
• ECs
• endothelial cells
• EMT
• epithelial-mesenchymal transition
• ER
• endoplasmic reticulum
• GMP
• granulocyte-macrophage progenitors
• GSEA
• Gene Set Enrichment Analysis
• GSVA
• Gene Set Variation Analysis
• hdWGCNA
• high-dimensional weighted gene co-expression network analysis
• ML
• lung metastasis
• MSigDB
• Molecular Signatures Database
• NM
• normal bone tissue
• OB
• osteoblasts
• OCs
• osteoclasts
• OS
• osteosarcoma
• PCA
• principal component analysis
• PT
• primary osteosarcoma
• scRNA-seq
• single-cell RNA sequencing
• TME
• tumor microenvironment
• UMAP
• uniform manifold approximation and projection
• UMI
• unique molecular identifier
• VEGF
• vascular endothelial growth factor

Background

Osteosarcoma (OS) is the most prevalent primary bone tumor in adolescents, characterized by its high aggressiveness [1]. In previous reports, osteosarcoma has been identified as the second leading cause of death in adolescents [2]. Currently, the standard treatment for osteosarcoma is the neoadjuvant chemotherapy-surgery-conservative chemotherapy treatment regimen, which has shown significant benefits in patients with localized osteosarcoma, with an overall survival rate of 60%–70% [3, 4]. However, the lungs are the most common site for metastatic colonization in OS [5]. Patients with lung metastases from osteosarcoma have a 5-year survival rate of only 20%, indicating that the occurrence of lung metastasis is associated with poor prognosis [6, 7]. Meanwhile, traditional treatments such as surgery, chemotherapy, and radiotherapy have seen some development over the past 40 years, but the prognosis for metastatic osteosarcoma remains unsatisfactory. Given the limitations of traditional treatment regimens, emerging therapies such as targeted therapy and immunotherapy are being explored with hope. While it has been found that emerging therapies do have some effect on metastatic osteosarcoma, they also reveal their own limitations. On one hand, targeted therapy is constrained by the short progression-free survival time and insufficient objective response rate of targeted drugs, while on the other hand, immunotherapy is not yet approved for osteosarcoma patients due to the low immunogenicity of osteosarcoma and the immunosuppressive microenvironment [8, 9]. Therefore, there is an urgent need to utilize new therapeutic strategies to improve the treatment outcomes for patients with lung metastases from osteosarcoma. In-depth research into the mechanisms of OS lung metastasis is crucial for providing theoretical support for new treatments.

The mechanisms of osteosarcoma lung metastasis are complex, involving multiple molecules and pathways [10–13]. Research by Adam A. Sabile et al. has demonstrated that Cyr61 induces lung metastasis of OS by activating the integrin-PI3K/Akt/GSK3β signaling cascade [7]. Additionally, Yu Han et al. have reported that tumor-associated macrophages promote OS metastasis and invasion by activating the COX-2/STAT3 axis and inducing epithelial-mesenchymal transition (EMT) [14]. Other studies have suggested that the tumor microenvironment(TME), composed of immune cells, stromal cells, blood vessels, and extracellular matrix(ECM), serves as a “fertile ground” for the occurrence and development of OS, providing favorable conditions for the colonization of OS cells in the lungs [15–17]. Angiogenesis and ECM remodeling within the TME are key factors contributing to lung metastasis and are closely associated with the formation of a pre-metastatic niche [18–21]. Additionally, single-cell RNA sequencing (scRNA-seq) and bulk RNA sequencing (bulk RNA-seq) technologies have become essential tools for investigating disease mechanisms [22–24].

In this study, scRNA-seq data were utilized to uncover the cellular heterogeneity within the TME, while bulk RNA-seq data were used to validate the gene expression trends. Based on scRNA-seq data, a cell subpopulation with lung metastatic characteristics was identified. Subsequently, hub genes associated with metastasis were identified through high-dimensional weighted gene co-expression network analysis (hdWGCNA) and differential expression analysis of bulk RNA-seq data [25]. To further investigate the mechanisms by which the TME influences lung metastasis, CellphoneDB and CSOmap were utilized [26].

Methods

Data Source

This study incorporated data from both single-cell RNA sequencing (scRNA-seq) and bulk RNA sequencing (bulk RNA-seq). The scRNA-seq data consisted of six primary osteosarcoma (OS) samples from our previous research (published as internal data) [17], along with two OS lung metastasis samples from GSE152048 and four normal bone tissue samples from the GEO database (GSE169396). To observe the overall gene expression patterns across cells, bulk RNA sequencing was also performed on four primary OS samples and four lung metastasis samples (new internal data). Additionally, phenotypic information from 18 samples in the GSE14359 dataset was integrated into the scRNA-seq data, and the GSE220538 dataset was used for validating the expression patterns of target genes.

Initial Processing of Raw Data

The clinical information of the patients is provided in Table S1. For the processing of scRNA-seq data, tumor tissue samples were obtained from six OS patients who had not undergone chemotherapy. The samples were cut into 1 mm³ pieces and digested into tissue suspensions. ScRNA-seq was performed in paired-end sequencing mode, with 150 base pairs read on both the Read 1 and Read 2 ends. The CellRanger software (version 4.0.0) provided by 10× Genomics () was downloaded and used to complete the upstream analysis of the scRNA-seq data according to the user manual. The scRNA-seq data were aligned to the human reference genome GRCh38, filtered, and then processed to generate count barcodes and unique molecular identifier (UMI) matrices, including the standard triplet files: barcodes.tsv, genes.tsv, and matrix.mtx.

The raw bulk RNA-seq data were obtained from four primary OS lesions and four OS lung metastasis lesions. Samples were collected through clinical surgery, and total RNA was extracted using the TRIzol reagent according to the manufacturer's instructions. RNA purity and quantification were assessed using a NanoDrop 2000 spectrophotometer (Thermo Scientific, USA), and RNA integrity was evaluated with an Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA). Transcriptome libraries were constructed using the VAHTS Universal V6 RNA-seq Library Prep Kit following the manufacturer's protocol. Sequencing was performed on the Illumina NovaSeq 6000 platform(), generating 150 bp paired-end reads. The raw reads in fastq format were processed using fastp software (version 0.20.0) to remove low-quality reads, resulting in clean reads for subsequent data analysis. Alignment to the reference genome and calculation of gene expression levels were performed using HISAT2 software (version 2.2.0). Gene expression profiles, including read counts for each gene, were obtained using the HTSeq-count tool.

Quality Control and Normalization of scRNA-Seq Data

The standard triplet files generated from upstream analysis were processed using the Seurat package (version 4.3.0) in R software (version 4.1.0) to extract data from each sample, which were then converted into individual Seurat objects. Subsequently, all Seurat objects were merged into a single comprehensive Seurat object using the merge() function in Seurat. To filter out low-quality cells, strict filtering criteria were applied: cells were retained if their gene expression levels were between 500 and 4500 genes, with mitochondrial gene content not exceeding 10% and erythrocyte gene content not exceeding 1%. Cells with RNA counts above the 97th percentile or below 1000 were excluded. The Subset() function in Seurat was used to remove cells that did not meet these standards.

Following quality control, normalization of the scRNA-seq data was performed to enhance the suitability of the data for further analysis. The NormalizeData() function in Seurat, using the “LogNormalize” parameter, was employed to normalize the sparse gene expression matrix within the Seurat object. Finally, the FindVariableFeatures() function in Seurat was used to identify the 2000 genes with the highest variability from the expression matrix, and the ScaleData() function was applied to scale and center the gene expression matrix.

Dimensionality Reduction, Clustering, and Cell Type Identification of scRNA-Seq Data

Dimensionality reduction transforms the complex, multidimensional information within single-cell samples into a more interpretable form, thereby clarifying the biological significance of the analysis results. Principal component analysis (PCA) was performed using the RunPCA() function in the Seurat package to identify the most variable directions in the data. To eliminate potential batch effects, the RunHarmony() function in the Harmony package (version 0.1.1) was applied. To further explore cellular heterogeneity, unsupervised clustering of the cells was conducted using the FindNeighbors() and FindClusters() functions in Seurat. The high-dimensional data were then projected into lower-dimensional space using the RunTSNE() and RunUMAP() functions. Marker genes for the different cell clusters were identified using the FindAllMarkers() function and were matched with marker genes from the Human Primary Cell Atlas to define the cell types.

Construction of Lung Metastasis Score

A total of 24 genes positively regulating OS lung metastasis were collected from the PubMed database () and are listed in Table S2. Based on this lung metastasis gene set, the AUCell package (version 1.16.0) was used to evaluate the association strength between individual cells and the lung metastasis phenotype by constructing a gene regulatory network (AUCell). Violin plots of the lung metastasis scores were generated using the ggplot2 package (version 3.3.4).

Identification of Malignant Cells and hdWGCNA

Single-cell copy number variation (CNV) analysis is commonly used to distinguish between benign and malignant cells, with increases or deletions in copy number suggesting a tendency toward malignancy. A CNV object was constructed using the CreateInfercnvObject() function in the infercnv package (version 1.10.1), and the degree of malignancy was visualized using heatmaps and violin plots. The hdWGCNA is suitable for scRNA-seq data, allowing for the construction of co-expression networks across multiple scales of cellular and spatial hierarchies to decipher specific gene expression patterns involved in various biological processes. Following the official tutorial (), the Seurat object was converted into a Metacells object by the hdWGCNA package (version 0.2.19), and a co-expression network with a soft threshold of 10 was constructed for standard downstream analysis.

Trajectory Analysis

Trajectory analysis in single-cell transcriptomics can infer the temporal sequence and state transitions during cellular development or differentiation. For this purpose, the Monocle2 package (version 2.22.0) and the slingshot package (version 2.2.1) were utilized. These algorithms perform trajectory inference on dimensionally reduced data and arrange cells along pseudotime to reconstruct the developmental trajectories of target cells. The CytoTRACE package (version 0.3.3) was used to assess the stemness characteristics of each cell, aiding in the identification of the starting points of cellular developmental trajectories.

Additionally, scRNA-seq velocity analysis was conducted to reveal the dynamic changes in gene expression over time, providing insights into the molecular mechanisms underlying cell differentiation. The Velocyto package (version 0.17.17) and the scvelo package (version 0.3.1) were employed within the Python software environment (version 3.8.19) to construct temporal gene expression maps of cell clusters.

Differential Analysis and Scissor Algorithm

Based on bulk RNA-seq data, differentially expressed gene (DEG) between primary OS and lung metastasis was identified using the limma package (version 3.50.3), and the results were visualized as a volcano plot. Additionally, the Scissor algorithm can automatically select cell subpopulations from scRNA-seq data that are associated with the phenotypic information from bulk RNA-seq data. The Scissor package (version 2.0.0) was utilized to map the lung metastasis phenotype from bulk RNA-seq data onto the corresponding cell subpopulations in the scRNA-seq data.

Cell–Cell Communication

Complex cell–cell communication plays a crucial role in the interactions between cells. The CellphoneDB package (version 5.0.0) in Python (version 3.9.0) was used to infer and analyze cell–cell communication by identifying potential mechanisms of interaction through ligand-receptor pairs. Additionally, the CSOmap tool was employed to simulate direct cell-to-cell communication in spatial contexts based on single-cell sequencing data. In this study, the CSOmapR package (version 1.0) was utilized to identify cell communication networks within the TME, providing deeper insights into the nature of intercellular interactions.

Enrichment Analysis

Gene Set Variation Analysis (GSVA) is an enrichment analysis method used to assess the activity levels of gene sets across different groups. Gene Set Enrichment Analysis (GSEA) identifies significant differences in the expression of entire gene sets between two biological states. GSVA and GSEA were conducted based on the single-cell gene expression matrix to determine the pathways and functions significantly enriched in each cell population. The gene sets from the Molecular Signatures Database (MSigDB)—including GO v2023.1, KEGG v2023.1, and HALLMARK v2023.1—were used as reference gene sets for GSVA and GSEA enrichment analyses (MSigDB: ). The GSVA package (version 1.42.0) and the GSEA package (version 1.2) are both suitable for this study.

Cell Culture

In this study, the hFOB 1.19 cell line was used as human normal osteoblasts, while the 143B, HOS, MG63, U2OS, and SAOS2 cell lines were used as human OS cells. All cell lines were purchased from Procell Life Science & Technology Co. Ltd. (Wuhan, China). The 143B, HOS, and MG63 cells were cultured in DMEM medium, whereas U2OS and SAOS2 cells were cultured in McCoy's 5A medium. Both media were supplemented with 1% penicillin/streptomycin (Solarbio, Beijing, China), of which the SAOS2 cell line was supplemented with 15% fetal bovine serum (FBS, Gibco), and the rest of the cell lines were supplemented with 10% FBS. OB cells were cultured in DMEM medium. The OS cell lines were maintained in a humidified incubator at 37°C with 5% CO₂, while hFOB 1.19 was cultured at 33.5°C with 5% CO₂.

Cell Transfection

According to the instructions of Lipofectamine 3000 reagent (Invitrogen, USA), the overexpression plasmid containing SEC16B was mixed with the transfection reagent and introduced into 143B and HOS cells to achieve transient transfection within 24–48 h. After transfection, SEC16B expression levels were detected by RT-qPCR to verify the overexpression efficiency.

RT-qPCR Detection

Total RNA was extracted using the RNA fast 200 Kit (Fastagen Biotech, Shanghai, China) according to the manufacturer's instructions. RNA was reverse transcribed into complementary DNA (cDNA) using a cDNA synthesis kit (Takara, Japan). RT-qPCR was performed on a StepOnePlus Real-Time PCR System (ABI7500) using SYBR Green (FastStart Universal SYBR Green Master ROX, Germany) as the detection dye. The PCR protocol was as follows: an initial denaturation at 95°C for 10 min, followed by 40 cycles of 95°C for 10 s and 60°C for 1 min. With GAPDH as a control and gene expression levels in the cell lines were presented as relative expression and calculated using the 2^−ΔΔCt method. Primer sequences can be found in Table S3.

Cell Migration and Invation Assay

HOS and 143B cells were cultured in 24-well plates. An 8 μm pore Transwell chamber was placed in each well, and 50,000 cells were added to each chamber. Three replicates were set up for each group according to the experimental conditions and cell types. For the invasion assay, a layer of Matrigel was also applied to the chamber. The lower chamber contained 600 μL of DMEM medium supplemented with 10% FBS as a chemoattractant. After 48 h of incubation, the cells were fixed with 4% paraformaldehyde and stained with crystal violet. After wiping off non-migrated cells inside the chamber, the migrated and invasive cells were observed and photographed using an inverted phase-contrast microscope (Olympus BX53). The number of cells in each group was quantified using Image J software (version Fiji, Image-J Software, USA).

Statistical Analysis

Data analysis was performed using R (version 4.1.0) and GraphPad Prism (version 9.0; San Diego, CA, USA). Unpaired Student's t-test, Wilcoxon rank-sum test, Kruskal–Wallis rank-sum test, and Spearman correlation analysis were employed in this study. A p value of < 0.05 was considered statistically significant.

Results

This study employed multiple bioinformatics analysis methods to explore the mechanisms underlying lung metastasis in OS and validated the findings using RT-qPCR and Transwell assays. The overall workflow is detailed in Figure 1A.

[IMAGE OMITTED. SEE PDF]

ScRNA-Seq Landscape of the Tumor-Infiltrating Microenvironment

In this study, cellular information was extracted from scRNA-seq data of 12 samples, resulting in a total of 57,810 high-quality cells after quality control and data normalization. To further explore cellular heterogeneity, the high-dimensional cell data were reduced to two-dimensional space and clustered into 11 cell clusters using unsupervised clustering. Based on classic marker genes from the Human Primary Cell Atlas, nine distinct cell types were identified: OB (IBSP, RUNX2, ALPL), macrophages (C1QB, C1QC, C1QA), NKT cells (NKG7, CD3D, CD3E), granulocyte-macrophage progenitors(GMP) (MPO, AZU1, MKI67), neutrophils (CAMP, S100A12, S100A8), B cells (MZB1, CD79A, CD79B), cancer-associated fibroblasts (CAFs) (TAGLN, ACTA2, LUM), endothelial cells (ECs) (PLVAP, VWF, ENG), and osteoclasts(OCs) (ACP5, MMP9, CTSK) (Figure 1B–D).

By comparing the differences in cell proportions among the normal bone tissue(NM), primary OS (PT), and lung metastasis (ML) groups, it was observed that the OB cells increased progressively with disease progression (Figure 1E). To differentiate between OS cells and normal bone tissue cells, OB were categorized as either OS cells or OCs based on the sample origin. The results of the CNV analysis revealed that OS cells exhibited higher CNV levels, confirming the malignant characteristics of OS cells (Figure S1A). Thus, the scRNA-seq landscape reveals the heterogeneity and diversity of the TME, highlighting the close association between OS lung metastasis and poor prognosis. These findings prompt further investigation into potential biological insights.

Identification of the Lung Metastasis Subpopulation in OS Cells

OS is a heterogeneous tumor, and therefore, the subpopulations of OS cells within the TME exhibit variability. Using multiple algorithms and scoring methods, we identified the OS cell subpopulations associated with lung metastasis at the genetic level. Initially, the OS cells were categorized into five groups: C0 (7031 cells), C1 (5421 cells), C2 (793 cells), C3 (704 cells), and C4 (439 cells) (Figure 2A). To further pinpoint the OS cell subpopulation responsible for lung metastasis, a lung metastasis score was constructed based on a lung metastasis gene set, revealing that C1 cells were most strongly associated with lung metastasis (Figure 2B). Trajectory analysis performed using the slingshot and scvelo algorithms indicated that C1 cells are the origin of OS cell development and differentiation, suggesting that C1 cells possess higher stemness characteristics (Figure 2C,D, Figure S1B). Furthermore, by calculating the CNV scores of OS cell subpopulations with normal OB as a reference, C1 cells were found to exhibit the highest degree of malignancy compared to the other subpopulations (Figure 2E). In subsequent GSVA analysis, the functions of C1 cells were primarily enriched in pathways related to EMT and angiogenesis (Figure 2F). Previous studies have indicated that these pathways are closely linked to OS cell metastasis [10]. Therefore, this study identifies C1 cells as the subpopulation involved in OS lung metastasis.

[IMAGE OMITTED. SEE PDF]

Experimental Identification of Target Genes for OS Lung Metastasis

The hdWGCNA was employed to associate C1 cells with lung metastasis characteristics, with an optimal soft threshold set at 10 (Figure 3A). Genes within the C1 subpopulation were categorized into 13 modules (M1 ~ M13) (Figure 3B). The M1 module, consisting of 1137 genes, was significantly downregulated in the ML samples, while the M6 module, containing 267 genes, exhibited a significant positive correlation with lung metastasis (Figure 3C).

[IMAGE OMITTED. SEE PDF]

To compare the overall gene expression patterns between the PT and the ML groups, differential analysis was conducted using bulk RNA-seq data from four primary tumors and four lung metastasis samples. The results of this differential analysis were visualized as a volcano plot (Figure 3D). To further identify potential target genes influencing the lung metastasis of C1 cells, the downregulated genes from the M1 module and the upregulated genes from the M6 module were intersected with the DEG identified in the bulk RNA-seq analysis. The Venn diagram revealed three common downregulated genes: SEC16B, APELA, and FAM50B (Figure 3E). Finally, the density plot confirmed that these genes were significantly downregulated in C1 cells (Figure 3F).

Experimental Validation of Target Gene Functions

To validate the functions of the identified genes, RT-qPCR and Transwell assays were conducted for in vitro functional verification. The RT-qPCR results showed that SEC16B expression was significantly reduced in OS cell lines (143B, HOS, MG63, U2OS, SAOS2) compared to OB cells, particularly in 143B and HOS cells (Figure 4A). To further investigate the role of SEC16B, overexpression experiments were performed in the 143B and HOS cell lines. The transfection results demonstrated that in the overexpression group, SEC16B expression was significantly higher in 143B and HOS cells compared to the control group, indicating a significant increase in overexpression efficiency (Figure S2B,C). Transwell assays demonstrated that SEC16B overexpression resulted in a marked inhibition of invasive and migratory abilities (Figure 4B,C). External BULK-seq datasets also confirm that the expression level of SEC16B is relatively lower in the lung metastasis group (Figure S2A). These in vitro experiments suggest that low SEC16B expression exhibits enhanced invasive and migratory capabilities.

[IMAGE OMITTED. SEE PDF]

Cell Communication Patterns Mediating OS Lung Metastasis in the TME

The TME exhibits complex cell–cell communication patterns, including molecular mechanisms involved in OS lung metastasis. The Scissor algorithm was utilized to correlate lung metastasis characteristics from bulk RNA-seq data with the cell clusters identified in scRNA-seq data. As shown in Figure 5A, macrophages, neutrophils, ECs, and fibroblasts were identified as Scissor+ cells, which are likely to promote the lung metastasis of OS. Conversely, OCs and B cells were classified as Scissor− cells, potentially inhibiting the metastasis of OS to the lungs. The study further delineated subpopulations within the TME and presented these findings using UMAP visualization (Figure 5B). Subsequently, cellphoneDB was utilized to analyze the cell communication network between C1 cells and stromal cells. The heatmap of the cell communication network revealed strong interactions between C1 cells and both LUM_CAFs and ECs. Additionally, a notable interaction was observed between C1 cells and MCAM_CAFs (Figure 5C). Ligand-receptor interactions were further explored to infer the potential mechanisms of OS metastasis to the lungs. Remarkably, we found that the PPIA/BSG signaling axis was activated between C1 cells, LUM_CAFs, and MCAM_CAFs, which may enhance the migratory capacity of OS cells [27]. Moreover, C1 cells activated angiogenesis-related signaling pathways, including VEGFB/NRP1, VEGFA/NRP2, and SDC1/ADGRA2, leading to the proliferation and migration of ECs and MCAM_CAFs to mediate angiogenesis (Figure 5D). CSOmap analysis was employed to simulate the spatial distribution of cells, thereby indirectly inferring cell communication patterns and addressing the absence of spatial information in cellphoneDB analysis. The CSOmap results demonstrated that C1 cells were in direct contact with the majority of LUM_CAFs, with the spatial proximity between these cells being statistically significant for direct cell–cell interactions (q < 0.05) (Figure 5E,F, Figure S1C). Similarly, the results suggested that communication between ECs and C1 cells may be mediated through the secretion of signaling factors.

[IMAGE OMITTED. SEE PDF]

Mechanisms by Which CAFs and ECs Establish the Pre-Metastatic Niche

The establishment of a pre-metastatic niche is a critical process in the distant metastasis of tumor cells. OS lung metastasis similarly depends on the formation of such a niche. CAFs and ECs exhibit close communication with C1 cells within the cell communication network and are key players in the TME's establishment of the pre-metastatic niche. This study investigated the molecular mechanisms by which these two cell types contribute to the formation of the pre-metastatic niche, based on their functional characteristics.

CAFs were classified into two functional subtypes, LUM_CAFs and MCAM_CAFs, based on marker genes. During OS lung metastasis, the proportion of LUM_CAFs increased (Figure 6A,B, Figure S1D). Furthermore, monocle trajectory analysis and CytoTRACE scoring revealed that the MCAM_CAFs subtype possesses higher stemness properties and transitions into the LUM_CAFs subtype as OS lung metastasis progresses (Figure S1E,F). It was observed that LUM_CAFs and MCAM_CAFs influence EC proliferation and migration through the VEGFB/NRP1 and VEGFA/NRP2 signaling axes, thereby remodeling the vasculature. They also enhance the migratory capacity of C1 cells through the THY1/ADGRE5 and EFNB1/EPHA2 signaling axes (Figure 6C,D). Functional enrichment analysis of CAFs indicated that MCAM_CAFs are associated with angiogenesis-related functions, such as vascular development, blood vessel morphogenesis, circulatory system development, and tubular structure morphogenesis. In contrast, LUM_CAFs are enriched in functions related to ECM remodeling and cell migration, including integrin binding, the non-canonical Wnt signaling pathway, regulation of locomotion, cell motility, and negative regulation of cell adhesion (Figure 6E,F). Furthermore, ECs were found to activate multiple molecular signaling pathways that promote the establishment of the pre-metastatic niche. Within the cell communication network, ECs promote CAF proliferation and angiogenesis through the PDGFB/PDGFRB and JAM2/JAM3 signaling axes and maintain the stemness of C1 cells by activating the JAG2/NOTCH2 signaling axis (Figure 6G). The PPIA-BSG signaling axis was significantly activated between C1 cells, ECs, MCAM_CAFs, and LUM_CAFs, suggesting its critical role in osteosarcoma lung metastasis (Figure 6H).

[IMAGE OMITTED. SEE PDF]

Discussion

Osteosarcoma (OS) is highly prone to lung metastasis, posing a significant clinical challenge. In other cancers, Jiahao Guo and colleagues have utilized single-cell RNA sequencing (scRNA-seq) and bulk RNA sequencing (bulk RNA-seq) data to identify RAB13 as a potential therapeutic target for ovarian cancer metastasis [25]. These studies highlight the importance of integrating scRNA-seq and bulk RNA-seq in recent cancer research, as they provide both single-cell resolution and a comprehensive view of gene expression, respectively [28–30]. scRNA-seq enables the precise identification of cellular heterogeneity and the localization of cells associated with metastasis, while bulk RNA-seq offers a broad overview of gene expression patterns. Building on this foundation, our study aims to delve deeper into the complex mechanisms underlying lung metastasis in osteosarcoma and to propose novel therapeutic strategies for this aggressive form of cancer. By combining the strengths of both sequencing techniques, we hope to provide valuable insights into the molecular landscape of osteosarcoma metastasis and to guide the development of more effective treatment approaches.

This study successfully constructed single-cell atlases of normal bone tissue, primary osteosarcoma, and osteosarcoma lung metastasis using scRNA-seq technology, revealing significant heterogeneity in the tumor microenvironment (TME) across different groups. The lung metastasis score indicates that C1 cells exhibit the highest potential for lung metastasis, and their stem cell-like characteristics have been confirmed through trajectory analysis, consistent with previous findings [21]. In CNV analysis, C1 cells display higher malignancy scores compared to other cell subpopulations, aligning with existing understanding of cancer biology and indicating that their identity is not that of migratory mesenchymal cells within the tumor microenvironment [31]. Functional enrichment analysis of this study shows that C1 cells are primarily activated in pathways including EMT (epithelial-to-mesenchymal transition), angiogenesis, Hedgehog signaling, and Wnt signaling. The EMT process promotes tumor invasion and metastasis, and the close relationship between EMT and osteosarcoma lung metastasis has been reported in previous studies [32, 33]. Notably, the Wnt signaling pathway, as a key activator of EMT, upregulates Twist, Slug, and N-cadherin, while downregulating E-cadherin, facilitating this transition, making Wnt signaling pathway activation a critical factor in driving osteosarcoma lung metastasis [32, 34]. Additionally, angiogenesis and Hedgehog signaling pathways are also implicated in osteosarcoma lung metastasis [35, 36]. Based on these findings, we have strong evidence that C1 cells are a subpopulation of osteosarcoma with significant metastatic potential.

The molecular alterations observed in C1 cells, which possess characteristics of lung metastasis in osteosarcoma (OS), are crucial factors influencing metastatic progression. By analyzing gene expression differences across different groups in scRNA-seq and bulk RNA-seq data, SEC16B (SEC16 homolog B) was identified as a key hub gene. Previous studies have shown that SEC16B, an endoplasmic reticulum (ER) export factor associated with human obesity, induces ER stress, thereby inhibiting cell proliferation [37, 38]. Furthermore, overexpression of SEC16B has been found to inhibit cell proliferation and growth in breast cancer with bone metastasis, as well as suppress the invasion and migration abilities of cancer cells [39]. In this study, we also found downregulation of SEC16B expression in both primary OS tumors and lung metastases through RT-qPCR experiments, with similar results confirmed by Transwell assays. These findings suggest that the reduced expression of SEC16B may contribute to the ability of OS cells to evade apoptosis and facilitate distant metastasis. This study is the first to identify SEC16B as a target gene regulating lung metastasis in OS, providing a solid theoretical basis for clinical treatment. However, the specific regulatory mechanisms of SEC16B in OS remain to be further explored in future studies.

TME represents a complex ecosystem composed of interactions between tumor cells and various host-derived cells [39]. The interaction networks within TME cells significantly influence tumor behavior, particularly their invasiveness [39]. Therefore, the lung metastatic phenotype was integrated with cell types from the single-cell atlas, revealing that myeloid cells, CAFs, and ECs promote lung metastasis in OS. Further analysis demonstrated close communication between CAFs, ECs, and C1 cells. Previous studies have shown that activated tumor-infiltrating fibroblasts promote lung metastasis in OS [40, 41]. Notably, this study is the first to discover that C1 cells, together with CAFs, activate the PPIA/BSG signaling axis. Overactivation of the PPIA/BSG signaling axis promotes metastasis and stem-cell-like properties in cancer cells, facilitating ECM remodeling by downstream activation of the MMP gene family through LUM_CAFs and MCAM_CAFs [27, 42]. Additionally, LUM_CAFs and MCAM_CAFs were found to influence ECs through the VEGFB/NRP1 and VEGFA/NRP2 signaling axes, reshaping the vasculature. Previous research has demonstrated that the expression of BSG (CD147) can induce angiogenesis by increasing vascular endothelial growth factor (VEGF) production, which promotes malignant tumor invasion and metastasis, leading to poor prognosis [27]. Furthermore, CAFs and ECs positively regulate the JAG2/NOTCH2 signaling axis, maintaining the stem-cell-like properties of OS cells, which facilitates lung metastasis [43]. Based on these findings, we propose that OS cells with low SEC16B expression utilize the cell interaction network to excessively activate pathways related to ECM remodeling and angiogenesis, thereby establishing a pre-metastatic niche for osteosarcoma metastasis.

Although this study identified a group of osteosarcoma cells that promote lung metastasis and explored their molecular characteristics, certain limitations remain. First, the samples used in this study were obtained from the GEO database and prior research by our team. While the validity of the results was supported by various bioinformatics analyses, the relatively small sample size may limit the generalizability of the findings. Second, due to time and research cost constraints, comprehensive functional mechanism experiments have not yet been conducted. The in vitro experiments only validated the function of the target gene SEC16B, but its specific mechanism in influencing lung metastasis in osteosarcoma has not been confirmed experimentally.

Conclusions

In this study, a subset of OS cells with lung metastatic potential was identified, characterized by low expression of the SEC16B gene. Additionally, it was demonstrated that CAFs, ECs, and OS cells contribute to the remodeling of the TME and the formation of a pre-metastatic niche, thereby promoting lung metastasis.

Author Contributions

All authors contributed to the article and approved the submitted version.

Acknowledgments

We thank the patients, their families, and the investigators and research staff involved.

Ethics Statement

The studies involving human participants were reviewed and approved by the Ethics Committee and Institutional Review Board of the First Affiliated Hospital of Guangxi Medical University (2019KY-E-097).

Consent

Written informed consent for participation in this study was provided by the participants' legal guardians/next of kin. Written informed consent was also obtained from the legal guardians/next of kin of the individual(s) for the publication of any potentially identifiable images or data included in this article.

Conflicts of Interest

The authors declare no conflicts of interest.

Data Availability Statement

The original data presented in this study are included in the article/Supporting Information. Further inquiries can be directed to the corresponding authors.

References

Z. Si, Z. Shen, F. Luan, and J. Yan, “PINK1 Regulates Apoptosis of Osteosarcoma as the Target Gene of Cisplatin,” Journal of Orthopaedic Surgery and Research 18 (2023): 132.

H. C. Beird, S. S. Bielack, A. M. Flanagan, et al., “Osteosarcoma,” Nature Reviews Disease Primers 8 (2022): 77.

P. Pilavaki, A. Gahanbani Ardakani, P. Gikas, and A. Constantinidou, “Osteosarcoma: Current Concepts and Evolutions in Management Principles,” Journal of Clinical Medicine 12 (2023): 2785.

J. Wang, M. Li, P. Guo, and D. He, “Correction: Survival Benefits and Challenges of Adjuvant Chemotherapy for High‐Grade Osteosarcoma: A Population‐Based Study,” Journal of Orthopaedic Surgery and Research 18 (2023): 834.

I. Singh, N. Rainusso, L. Kurenbekova, et al., “Intrinsic Epigenetic State of Primary Osteosarcoma Drives Metastasis,” bioRxiv [Preprint], November 13, 2023.

Y. C. Huang, W. C. Chen, C. L. Yu, et al., “FGF2 Drives Osteosarcoma Metastasis Through Activating FGFR1‐4 Receptor Pathway‐Mediated ICAM‐1 Expression,” Biochemical Pharmacology 218 (2023): 115853.

A. A. Sabile, M. J. Arlt, R. Muff, et al., “Cyr61 Expression in Osteosarcoma Indicates Poor Prognosis and Promotes Intratibial Growth and Lung Metastasis in Mice,” Journal of Bone and Mineral Research 27 (2012): 58–67.

F. Duffaud, O. Mir, P. Boudou‐Rouquette, et al., “Efficacy and Safety of Regorafenib in Adult Patients With Metastatic Osteosarcoma: A Non‐Comparative, Randomised, Double‐Blind, Placebo‐Controlled, Phase 2 Study,” Lancet Oncology 20 (2019): 120–133.

S. Yu and X. Yao, “Advances on Immunotherapy for Osteosarcoma,” Molecular Cancer 23 (2024): 192.

H. Kong, W. Yu, Z. Chen, et al., “CCR9 Initiates Epithelial‐Mesenchymal Transition by Activating Wnt/Beta‐Catenin Pathways to Promote Osteosarcoma Metastasis,” Cancer Cell International 21 (2021): 648.

C. W. Lee, Y. C. Chiang, P. A. Yu, et al., “A Role of CXCL1 Drives Osteosarcoma Lung Metastasis via VCAM‐1 Production,” Frontiers in Oncology 11 (2021): 735277.

B. Moukengue, H. K. Brown, C. Charrier, et al., “TH1579, MTH1 Inhibitor, Delays Tumour Growth and Inhibits Metastases Development in Osteosarcoma Model,” eBioMedicine 53 (2020): 102704.

S. Shao, L. Piao, J. Wang, et al., “Tspan9 Induces EMT and Promotes Osteosarcoma Metastasis via Activating FAK‐Ras‐ERK1/2 Pathway,” Frontiers in Oncology 12 (2022): 774988.

Y. Han, W. Guo, T. Ren, et al., “Tumor‐Associated Macrophages Promote Lung Metastasis and Induce Epithelial‐Mesenchymal Transition in Osteosarcoma by Activating the COX‐2/STAT3 Axis,” Cancer Letters 440–441 (2019): 116–125.

C. Brachelente, F. Torrigiani, I. Porcellato, et al., “Tumor Immune Microenvironment and Its Clinicopathological and Prognostic Associations in Canine Splenic Hemangiosarcoma,” Animals (Basel) 14 (2024): 1224.

C. Hou, M. Lu, Z. Lei, et al., “HMGB1 Positive Feedback Loop Between Cancer Cells and Tumor‐Associated Macrophages Promotes Osteosarcoma Migration and Invasion,” Laboratory Investigation 103 (2023): 100054.

H. Tang, S. Liu, X. Luo, et al., “A Novel Molecular Signature for Predicting Prognosis and Immunotherapy Response in Osteosarcoma Based on Tumor‐Infiltrating Cell Marker Genes,” Frontiers in Immunology 14 (2023): 1150588.

C. Deng, Y. Xu, H. Chen, et al., “Extracellular‐Vesicle‐Packaged S100A11 From Osteosarcoma Cells Mediates Lung Premetastatic Niche Formation by Recruiting gMDSCs,” Cell Reports 43 (2024): 113751.

Y. Liu and X. Cao, “Characteristics and Significance of the Pre‐Metastatic Niche,” Cancer Cell 30 (2016): 668–681.

H. Peinado, H. Zhang, I. R. Matei, et al., “Pre‐Metastatic Niches: Organ‐Specific Homes for Metastases,” Nature Reviews. Cancer 17 (2017): 302–317.

Y. Xu, C. Deng, H. Chen, et al., “Osteosarcoma Cells Secrete CXCL14 That Activates Integrin alpha11beta1 on Fibroblasts to Form a Lung Metastatic Niche,” Cancer Research 84 (2024): 994–1012.

H. He, Y. Yang, L. Wang, et al., “Combined Analysis of Single‐Cell and Bulk RNA Sequencing Reveals the Expression Patterns of Circadian Rhythm Disruption in the Immune Microenvironment of Alzheimer's Disease,” Frontiers in Immunology 14 (2023): 1182307.

Y. Su, X. Zhang, Y. Liang, J. Sun, C. Lu, and Z. Huang, “Integrated Analysis of Single‐Cell RNA‐Seq and Bulk RNA‐Seq to Unravel the Molecular Mechanisms Underlying the Immune Microenvironment in the Development of Intestinal‐Type Gastric Cancer,” Biochimica et Biophysica Acta‐Molecular Basis of Disease 1870 (2024): 166849.

T. Song, Y. Yang, Y. Wang, Y. Ni, Y. Yang, and L. Zhang, “Bulk and Single‐Cell RNA Sequencing Reveal the Contribution of Laminin gamma2 −CD44 to the Immune Resistance in Lymphocyte‐Infiltrated Squamous Lung Cancer Subtype,” Heliyon 10 (2024): e31299.

J. Guo, X. Han, J. Li, et al., “Single‐Cell Transcriptomics in Ovarian Cancer Identify a Metastasis‐Associated Cell Cluster Overexpressed RAB13,” Journal of Translational Medicine 21 (2023): 254.

X. Ren, G. Zhong, Q. Zhang, L. Zhang, Y. Sun, and Z. Zhang, “Reconstruction of Cell Spatial Organization From Single‐Cell RNA Sequencing Data Based on Ligand‐Receptor Mediated Self‐Assembly,” Cell Research 30 (2020): 763–778.

J. M. Han and H. J. Jung, “Cyclophilin A/CD147 Interaction: A Promising Target for Anticancer Therapy,” International Journal of Molecular Sciences 23 (2022): 9341.

G. Huang, X. Zhang, Y. Xu, et al., “Prognostic and Predictive Value of Super‐Enhancer‐Derived Signatures for Survival and Lung Metastasis in Osteosarcoma,” Journal of Translational Medicine 22 (2024): 88.

K. Liu, H. Han, K. Xiong, et al., “Single‐Cell Landscape of Intratumoral Heterogeneity and Tumor Microenvironment Remolding in Pre‐Nodal Metastases of Breast Cancer,” Journal of Translational Medicine 22 (2024): 804.

Z. Wang, Y. Wang, M. Chang, et al., “Single‐Cell Transcriptomic Analyses Provide Insights Into the Cellular Origins and Drivers of Brain Metastasis From Lung Adenocarcinoma,” Neuro‐Oncology 25 (2023): 1262–1274.

J. Rosai, E. A. Saxen, and L. Woolner, “Undifferentiated and Poorly Differentiated Carcinoma,” Seminars in Diagnostic Pathology 2 (1985): 123–136.

S. Brabletz, H. Schuhwerk, T. Brabletz, and M. P. Stemmler, “Dynamic EMT: A Multi‐Tool for Tumor Progression,” EMBO Journal 40 (2021): e108647.

C. H. Hou, F. L. Lin, S. M. Hou, and J. F. Liu, “Cyr61 Promotes Epithelial‐Mesenchymal Transition and Tumor Metastasis of Osteosarcoma by Raf‐1/MEK/ERK/Elk‐1/TWIST‐1 Signaling Pathway,” Molecular Cancer 13 (2014): 236.

Y. Wu, C. Ginther, J. Kim, et al., “Expression of Wnt3 Activates Wnt/Beta‐Catenin Pathway and Promotes EMT‐Like Phenotype in Trastuzumab‐Resistant HER2‐Overexpressing Breast Cancer Cells,” Molecular Cancer Research 10 (2012): 1597–1606.

B. K. Nirala, T. Yamamichi, and J. T. Yustein, “Deciphering the Signaling Mechanisms of Osteosarcoma Tumorigenesis,” International Journal of Molecular Sciences 24 (2023): 11367.

L. Yu, H. Zhu, Z. Wang, et al., “Circular RNA circFIRRE Drives Osteosarcoma Progression and Metastasis Through Tumorigenic‐Angiogenic Coupling,” Molecular Cancer 21 (2022): 167.

R. Shi, W. Lu, Y. Tian, and B. Wang, “Intestinal SEC16B Modulates Obesity by Regulating Chylomicron Metabolism,” Molecular Metabolism 70 (2023): 101693.

M. Yamaguchi, N. Z. Ghanem, K. Hashimoto, J. W. Ramos, and T. Murata, “The Overexpressed Transcription Factor RGPR‐p117 Suppresses the Proliferation of Normal Rat Kidney Proximal Tubular Epithelial NRK‐52E Cells: Involvement of Diverse Signaling Pathways,” Life Sciences 306 (2022): 120795.

M. Yamaguchi, T. Murata, and N. Shimokawa, “Overexpression of RGPR‐p117 Reveals Anticancer Effects by Regulating Multiple Signaling Pathways in Bone Metastatic Human Breast Cancer MDA‐MB‐231 Cells,” IUBMB Life 77 (2025): e2939.

Y. Liu, J. Zhang, Z. Wang, et al., “Multi‐Omics Characterization Reveals the Pathogenesis of Liver Focal Nodular Hyperplasia,” iScience 25 (2022): 104921.

Y. Zhang, Z. Liu, X. Yang, et al., “H3K27 Acetylation Activated‐COL6A1 Promotes Osteosarcoma Lung Metastasis by Repressing STAT1 and Activating Pulmonary Cancer‐Associated Fibroblasts,” Theranostics 11 (2021): 1473–1492.

H. Fan, S. Lu, S. Wang, and S. Zhang, “Identification of Critical Genes Associated With Human Osteosarcoma Metastasis Based on Integrated Gene Expression Profiling,” Molecular Medicine Reports 20 (2019): 915–930.

H. Jin, S. Luo, Y. Wang, et al., “miR‐135b Stimulates Osteosarcoma Recurrence and Lung Metastasis via Notch and Wnt/Beta‐Catenin Signaling,” Molecular Therapy—Nucleic Acids 8 (2017): 111–122.