HIF2A mediates lineage transition to aggressive phenotype of cancer-associated fibroblasts in lung cancer brain metastasis

ABSTRACT

Brain metastasis is the most devasting form of lung cancer. Recent studies highlight significant differences in the tumor microenvironment (TME) between lung cancer brain metastasis (LCBM) and primary lung cancer, which contribute significantly to tumor progression and drug resistance. Cancer-associated fibroblasts (CAFs) are the major component of pro-tumor TME with high plasticity. However, the lineage composition and function of CAFs in LCBM remain elusive. By reanalyzing single-cell RNA sequencing (scRNA-seq) data (GSE131907) from lung cancer patients with different stages of metastasis comprising primary lesions and brain metastasis, we found that CAFs undergo distinctive lineage transition during LCBM under a hypoxic situation, which is directly driven by hypoxia-induced HIF-2α activation. Transited CAFs enhance angiogenesis through VEGF pathways, trigger metabolic reprogramming, and promote the growth of tumor cells. Bulk RNA sequencing data was utilized as validation cohorts. Multiplex immunohistochemistry (mIHC) assay was performed on four paired samples of brain metastasis and their primary lung cancer counterparts to validate the findings. Our study revealed a novel mechanism of lung cancer brain metastasis featuring HIF-2α-induced lineage transition and functional alteration of CAFs, which offers potential therapeutic targets.

KEYWORDS:

• Brain metastasis
• Cancer-Associated Fibroblasts
• HIF
• lung cancer
• tumor microenvironment

Backgrounds

Lung cancer is one of the most common malignant cancers in adults.¹ Although improved systemic therapies have led to better control of primary lung cancer at early and advanced stages,¹ 40–50% of the patients will eventually develop brain metastasis.^2,3 The prognosis of patients with lung cancer brain metastasis (LCBM) is extremely poor. Even with the latest integrative treatment, the 5-year survival rate remains less than 5%.^4–6 Drug resistance occurs within a year after treatment⁷，suggesting the brain may provide a sanctuary for metastatic tumor cells. Taken together, their divergent treatment response implies that primary lung cancer and LCBM have different molecular and cellular mechanisms during progression.

Recent studies showed that tumor microenvironment (TME) in LCBM is different from that in primary lung cancer.^8,9 LCBM has an immunosuppressive TME with an accumulation of TGF-β secreting cancer-associated fibroblasts (CAFs),¹⁰ and anti-proinflammatory macrophages.⁹ Alterations in stromal components of TME are also observed, including elevated expression of extracellular matrix and the adhesion molecule ALCAM.^11,12 The unique TME in LCBM may explain tumor cell settlement and drug resistance¹³ and highlight the importance of TME in LCBM.

CAFs, one of the most abundant cell populations in TME, can remodel the TME landscape by adjusting their secretion patterns and matrix deposition levels, thus promoting tumor progression and metastasis.^14,15 Sun et al. identified active interaction and spatial colocalization of CAFs and endothelial cells in both LCBM and gliomas, indicating a potential role of CAFs in promoting endothelial cell function, likely through hypoxia-responsive pathways.¹⁶ Apart from that, a pan-cancer analysis of CAFs revealed their diverse origins and robust plasticity in cellular function and lineage composition.¹⁷ Mucciolo et al. revealed that TGF-β-induced EGFR activation of a subset of myofibroblasts in pancreatic ductal adenocarcinoma, which promotes metastasis. Blockade of EGFR signaling completely depletes the subset of myofibroblasts. Ma et al. found diverse metabolic reprogramming mechanisms within distinct CAF subpopulations across six solid tumors.¹⁸ However, the transitions in lineage composition and functions of CAFs in LCBM remain elusive.

The advancement of single-cell RNA sequencing (scRNA-seq) has revolutionized our ability to comprehensively profile the tumor microenvironment at the single-cell resolution, facilitating a nuanced understanding of cellular lineage transitions in LCBM.^19,20 In this study, we performed an unbiased analysis of bulk RNA and single-cell RNA sequencing data of both primary lung cancer and LCBM. We categorized CAFs and characterized their transitions and validated them using samples from LCBM and their primary counterparts. Our results indicate that HIF-2α drives lineage transition of CAFs to malignant phenotypes during LCBM, which enhances angiogenesis through VEGF pathway, trigger metabolic reprogramming, and promote the growth of tumor cells. Collectively, we elucidated the HIF-2α mediated lineage transition mechanisms of CAFs in LCBM, thereby providing valuable insights into potential therapeutic avenues.

Materials and methods

Data collection and processing

Bulk-RNA sequence data (GSE121968, GSE112996, and GSE141685) was assessed from Gene Expression Omnibus (GEO) and cBioPortal dataset Metastatic Melanoma (DFCI, Nature Medicine 2019).²¹ Data were cleaned to remove duplicate rows and NA values and were normalized for the following analysis. Data from the GSE112996 and GSE141685 datasets was merged, and the R package Combat was utilized to reduce batch effect.²² The scRNA-seq data of different metastatic stages of lung adenocarcinoma were accessed from the GSE131907 dataset.⁸ Sequencing data from normal lung tissue in the dataset were excluded from this analysis. Clinical information of the patients and metadata were extracted simultaneously. R package Seurat (Version 4.2.0)²³ was used to process scRNA-seq data. QC control of raw data was the same as the source literature. Principal component analysis (PCA) was utilized and incorporated with the top 2000 highly variable features by the FindVariableFeatures function in the Seurat package to reduce the dimensionality of this dataset, and the first 20 PCs were identified for analysis. Uniform Manifold Approximation and Projection (UMAP) reduction was then adopted for further dimensionality reduction and visualization. Cell cluster annotation was conducted by the expression of known marker genes.

Differential analysis and GO enrichment analysis

Differential analysis of GSE121968 Bulk-RNA sequence data was performed by R package DESeq2 (Version 1.38.1)²⁴ by original counts on 14,280 genes across six paired samples. Differential genes (DEG) threshold value setting was |Log2FoldChange| >1 and p < 0.05. 703 genes passed filtration and were identified as DEGs. Functional annotation and GO enrichment analysis were conducted on up-regulated DEGs using the Gene Ontology database by R package clusterProfiler.²⁵

GSEA enrichment analysis

Gene Set Enrichment Analysis (GSEA) was utilized to gain a holistic understanding of the functional landscape by using specific gene sets.²⁶ All 14,280 genes were pre-ranked based on their Log2FoldChange value and are used to assess whether predefined gene sets representing biological pathways or functional categories are significantly enriched. In our analysis, four gene sets were used and downloaded from the GSEA MSigDB database (https://www.gsea-msigdb.org), and one gene set was defined according to literatures (Table S3). The enrichment score and p value were calculated using the R package clusterProfiler (Version 4.6.2)²⁵ and R package irGSEA (Version 3.2.2). Data visualization was conducted using the R package GseaVis (Version 0.0.9).

TCGA RNA‑seq data processing and scRNA-Seq data deconvolution

RNA-seq data and phenotype profiles of the lung adenocarcinoma (LUAD) cohort from the TCGA database were downloaded from the UCSC Xena database. Raw count data was directly utilized after removing NAs. CIBERSORTx website (https://cibersortx.stanford.edu) was used to deconvolve the cell fraction details from the TCGA LUAD bulk-seq dataset by the reference matrix generated by scRNA-seq data.²⁷ Cell composition data were linked to survival data and compared between high/low groups using Kaplan–Meier curves.

Fibroblast cell extraction and re-clustering

All fibroblasts were extracted for subsequent analysis. R package Harmony (Version 0.1.1)²⁸ was utilized to ease the batch effect in fibroblast subgroups across all samples. Data then again underwent the same dimensionality reduction and cell annotation process as described above.

Cell–cell interaction analysis

We used the R package Cellchat (Version 1.6.1) to project receptor-ligand pairs onto our cell subsets within tissues of each origin to identify cell–cell interactions.²⁹ This method takes gene expression data from cells as input, and the probability of cell–cell communication is simulated by combining gene expression with prior knowledge of the interaction between signaling ligands, receptors, and their co-factors. Heatmap functions were used to display the whole landscape of outgoing and incoming signal patterns in different statuses. Only significant ligand-receptor pairs (p < 0.05) were extracted for illustration. Visualization of the potential interaction strength between ligand and receptor in specific pathways is performed using Cellchat build-in function netVisual_aggregate. netVisual_bubble function was used to demonstrate significant legend-receptor pairs change trend in brain metastases.

Trajectory analysis

Naive-CAF, Matrix-CAF, and Myo-CAF were subsets from Fibroblasts subgroups and were investigated by utilizing the Monocle2 algorithm (R package Version 2.24.0).³⁰ Insight into dynamic biological processes, such as interconversion and evolutionary trajectories of these three cell types are obtained through this analysis. The NewCellDataSet function was utilized to create a new object for Monocle, incorporating transcript count data from the included cell populations. Differences in mRNA across different cells were normalized using estimateSizeFactors and estimateDispersions functions. Signature genes expressed in at least 10% of the dataset cells, with a p-value <0.01 were calculated using the differentialGeneTest function and were incorporated to define the trajectory progress. After dimension reduction, orderCells function ordered the cells according to gene expression. Branched Expression Analysis Modeling (BEAM) analysis was utilized by function BEAM and plot_genes_branched_heatmap to analyze specified nodes and time-ordered cell data to find and cluster differential genes and deduce gene modules related to the branches.

Single-cell regulatory network inference and clustering (SCENIC) analysis

We employed the pySCENIC Python package (Version 0.12.0) for SCENIC analysis.³¹ Initially, the GRNBoost2 function was utilized to infer co-expression modules between transcription factors (TFs) and candidate target genes. Subsequently, the RcisTarget function identified enriched motifs in each co-expression module, defining a module as a TF and its potential direct target genes. Finally, the AUCell function assessed the activity of each module in individual cells. Cell-type specificity of a module was evaluated by the Module Specificity Score (RSS) method developed by Suo et al.³²

Acquisition of tumor samples

Paraffin-embedded blocks of two patient cohorts were utilized in this study, all of them underwent brain tumor resection surgery at Huashan Hospital (Shanghai, China) and are pathologically confirmed brain metastasis of lung cancer. Cohort 1 is comprised of paired resected tumor samples at the lung primary site and brain metastasis site from four patients. Cohort 2 is made up of samples from 10 patients. Five of them have lung cancer history and are on Tyrosine Kinase Inhibitor (TKI) medication but eventually developed LCBM (TKI-resistant group). The other five patients were diagnosed with LCBM as the initial diagnosis (TKI naive group). Tumors were examined by pathologists and confirmed to be lung cancer or LCBM. Written consents were obtained from all patients.

H&E staining and immunohistochemistry staining

All involved surgical tissue samples were fixed in formalin for 48 h, dehydrated, and embedded in paraffin. Paraffin blocks were then cut into 4-μm-thick slices. Hematoxylin and Eosin Staining (H&E) staining was performed using a H&E Kit (Beyotime, Cat#C0105S). Briefly, following dehydration and clearing, the specimens underwent hematoxylin staining to highlight the nuclei. Samples were subsequently stained with eosin for differentiation. Immunohistochemistry was conducted according to the manufacturer’s protocol. After deparaffinization and antigen retrieval, sections were blocked and incubated with the primary antibody, followed by incubation with secondary antibodies conjugated with horseradish peroxidase (HRP) (Vector Laboratories, Cat#MP-7500). An immunohistochemical reaction was developed using DAB substrate (Servicebio, Cat#C0105S), and counterstaining with hematoxylin was performed. The stained slides were scanned with an Olympus VS200 Slide Scanner.

Multiplex immunofluorescence staining

Multiplex immunofluorescence staining was performed using a 4-plex IHC Kit according to the manufacturer’s protocol (Ruchuang, Cat#RC0086Plus-34 R). Briefly, after deparaffinization and antigen retrieval, specimens were blocked and incubated with the primary antibody, followed by incubation with a secondary HRP-conjugated antibody (Vector Laboratories, Cat#MP-7500). The relative primary antibody markers were α-SMA (14395–1-AP, proteintech, 1:200), COL1A1 (67288–1-Ig, proteintech, 1:200), ADH1B (66939–1-Ig, proteintech, 1:200), PanCK (26411–1-AP, proteintech, 1:200), CD31 (66065–2-Ig, proteintech, 1:200), HIF-1α (66730–1-Ig, proteintech, 1:200) and HIF-2α (66731–1-Ig, proteintech, 1:200). Visualization was conducted by Opal TSA working solution for another 10 min, nuclei were stained with DAPI. The stained slides were scanned with an Olympus VS200 Slide Scanner.

Results

The TME of LCBM comprises of distinct CAF populations

To delineate the heterogeneity of TME during LCBM, we performed a comprehensive analysis of the scRNA-seq dataset (GSE131907) in primary and metastatic lung cancers (Figure 1(a)).⁸ We identified 106,771 cells across 11 distinct cell clusters by unsupervised clustering of cells from 32 patients (Figures 1(b) and S1a). Using the canonical markers of different cell types (Figure 1(c)), we further categorized these clusters into eight cell populations, i.e., epithelial cells, CAFs, endothelial cells, T lymphocytes, B lymphocytes, myeloid cells, MAST cells, and oligodendrocytes (Figure 1(b,c)). We calculated and compared the proportion of each cell type in different stages of lung cancers (Figure 1(d)). Compared with primary lung cancers, the proportion of nonmalignant cells in LCBM notably decreases, while the proportion of tumor cells increases (Figures 1(e) and S1c). Interestingly, the proportion of CAFs remains relatively unchanged (Figure 1(e)), suggesting that the function of CAFs in LCBM may not be determined by their overall quantity.

Figure 1. Landscape of the CAF composition in LCBM.

(a) Schematic diagram of this study. (b) Uniform Manifold Approximation and Projection (UMAP) plots of cells from 32 patients across all metastatic stages, which are grouped into eight major cell types. (c) Dot plot of canonical marker genes for nine major cell types by mean expression. (d) Cell composition proportions of patients with different metastasis status (i.e., tLung, tL/B, mLN, mBrian). (e) Proportion of nonmalignant cells (left) and CAFs (right) relative to the total cell count in different patients, stratified by metastatic status, p-values are calculated by t-test. *p value < 0.05, ns: not significant. (f) UMAP plots of cells from the “CAF” cluster, which are reclustered into eight major cell types. (g) Dot plot of canonical marker genes for four major CAF types by mean expression. (h) Illustration of canonical marker genes and enrichment analysis of all marker genes in each CAF subtype with Metascape database. (i) Representative mIF images of different kinds of CAFs, arrows highlight cells of interest. Scale bars, 50 μm.

As CAFs are essential for the TME of metastasis, we then performed a detailed analysis of CAFs in LCBM to elucidate their association with lung cancer metastasis, we used modularity-optimized clustering to cluster CAFs into five distinct subpopulations (Figure 1(f)) including previously defined CAF types.^17,33 Matrix-CAFs (POSTN⁺, COL5A1⁺) account for 38.3% of all CAFs. Naive-CAFs (ADH1B⁺) account for 30.8%, followed by Myo-CAFs (Myofibroblast-like CAFs) (21.6%, marked by ACTA2 and RGS5) and Ap-CAF (Antigen-presenting CAF) (9.2%, HLA-DRB1⁺, HLA-DPB1⁺). Cluster 5 did not manifest distinctive features of any known fibroblast subgroups and was therefore labeled “ND (Not Determined)” and excluded from subsequent analysis (Figure 1(g) and Table S1). To assess the functions of the aforementioned CAF subpopulations, we conducted a functional comparison analysis using the Metascape database. Pathway analysis revealed differential enrichment of upstream regulators and pathways (Figure 1(h) and Table S2). Specifically, collagen degradation is enriched in Matrix-CAFs, smooth muscle contraction in Myo-CAFs, morphogenesis of a branching epithelium pathway in Naive-CAFs, and leukocyte chemotaxis pathway in Ap-CAFs. Finally, we validated the existence of Naive-CAF, Myo-CAF, and Matrix-CAF in LCBM by mIF analysis of the samples of primary lung cancer and paired brain metastasis. Indeed, we observed distinct compositions of CAFs in LCBM using α-SMA, COL1A1, and ADH1B as mIF markers for Myo-CAF, Matrix-CAF, and Naive-CAF, respectively (Figure 1(i)). Together, we found that the TME of LCBM has distinct CAF populations, although the overall proportion of CAFs remains the same as primary lung cancers.

Lineage transition of CAFs occurs in LCBM

CAFs exhibit significant plasticity, whose composition is highly variable during tumor progression and strongly associated with the malignant advancement of tumors.¹⁹ We compared CAF compositions in the TME of primary and metastatic lung cancer and found distinct CAF composition in LCBM (Figure 2(a)). Notably, LCBM has significantly fewer Naive-CAFs (1.1%) than primary lung cancers (39.2%), while it has significantly more Myo-CAFs and Matrix-CAFs (Figure 2(b)). The number of Ap-CAFs remained unchanged (Fig. S1D). This was validated by mIF images of paired samples (Figure 2(c–e)). Consistently, we observed a significant decrease in Naive-CAFs in LCBM (p < 0.001). The increase of Matrix-CAFs (p = 0.097) and Myo-CAF (p = 0.34) was not significant probably due to the small sample size. After analyzing the lung adenocarcinoma (LUAD) cohort in TCGA, we found that the Matrix-CAF to Naive-CAF ratio directly correlates with poorer prognosis (HR = 1.43, p = 0.013, n = 585) (Figures 2(f) and S1e,f).

Figure 2. mIF images of paired lung cancer primary/brain metastasis samples.

(a) Left: UMAP plots of cells at primary state (tLung) and brain metastatic state (mBrain). Right: Cell composition proportions are compared between tLung and mBrain. (b) Proportion of Naive-CAF (left), Matrix-CAF (middle), and Myo-CAF (right) relative to the total CAF count in LCBM, p-values are calculated by t-test. *p value < 0.05, **p value < 0.01, ns: not significant. (c) Representative mIF images of nuclei (blue), α-SMA (red), COL1A1 (green), and ADH1B (purple) in lung cancer primary site (tLung), Scale bars: both 100 μm. (d) Representative mIF images of nuclei (blue), α-SMA (red), COL1A1 (green), and ADH1B (purple) in paired brain metastasis site (mBrain), Scale bars: both 100 μm. (e) Violin plot showing different cell proportions of Matrix-CAF, Myo-CAF, and Naive-CAF between lung cancer primary site and brain metastasis, the p values are calculated by the unpaired t-test. mBrain: lung cancer brain metastasis. (f) Kaplan – Meier curves for overall survival in TCGA LUAD patients according to Matrix-CAF/Naive-CAF ratio. HR, the p-value is calculated by the log-rank test, HR, hazard ratio.

We observed distinct composition of CAFs in LCBM, which may play different roles during cancer progression and metastasis. We then asked whether there are lineage transitions of CAFs in LCBM. We performed pseudotime trajectory by Monocle2, which revealed a gradual transition of the Naive-CAFs, leading to the emergence of two distinct branches A and B (Figure 3(a,b)). Myo-CAFs are enriched in the terminus of branch A, while Matrix-CAFs are in branch B (Figure 3(a,b)). The expression of key markers (i.e., ADH1B, ACTA2, and COL1A1) was coupled with the transition of CAFs in LCBM.

Figure 3. The lineage transition of CAFs in LCBM is associated with hypoxia.

(a) Top: Unsupervised trajectory of Naive-CAF, Matrix-CAF, and Myo-CAF, colored by pseudotime. Bottom: Expression level of ADH1B (a marker of Naive-CAF), ACTA2 (a marker of Myo-CAF), and COL1A1 (a marker of Matrix-CAF) along pseudotime. (b) Trajectory colored by cell type (up), EMT score (middle), and HIF1A-Glycolysis score (bottom). (c) Branched heatmap of differentially expressed genes between Branch A and Branch B with GO BP annotation of enriched terms of all function modules. (d) UMAP plots of cells from the “CAF” cluster, which are colored in accordance with the hypoxia hallmarks enrichment score. (e) Box plot of hypoxia hallmarks enrichment score between lung cancer primary site and LCBM, p-value is calculated by t-test. *p value < 0.05, **p value < 0.01, ***p value < 0.001, ****p value < 0.0001, ns: not significant. (f) Box plot of hypoxia hallmarks enrichment score between different kinds of CAFs, p values are calculated by one-way anova test. (g) Expression level of HIF1A and HIF2A between primary lung cancer and brain metastasis site among all CAF types. (h) Box plot of HIF-1α (left) and HIF-2α (right) pathway enrichment scores between primary lung cancer and brain metastasis site, p values are calculated by t-test. (i) Expression level of HIF1A and HIF2A along pseudotime on Branch A and B. (j) Box plot of HIF1A and HIF2A expression level between different LCBM status of merged datasets comprising GSE112996 and GSE141685 (left) and DFCI metastatic melanoma cohort (right). nonBM: lesions outside brain parenchymal, BM: brain metastasis. p values are calculated by the Wilcox test. (k) (Left) Representative H&E and IHC images showing paired samples of lung cancer primary lesion and brain metastasis lesion from the same patient, tLung: primary lung cancer, mBrain: lung cancer brain metastasis, scale bar: 100 μm. (Right) Violin plot showing IHC scores of HIF-1α and HIF-2α in lung cancer primary site and brain metastasis, the p values are calculated by the unpaired t-test.

Four distinct functional modules of the differentially expressed genes in CAFs were identified during the transitions from Naive-CAF into Matrix-CAF or Myo-CAF (Figure 3(c)). In addition to the enrichment of epithelial–mesenchymal transition (EMT) and angiogenesis, a shared pathway in responding to low oxygen levels emerged as a common feature in both branches of differentiation (Figure 3(c)). In line with the functional annotation of Myo-CAFs and Matrix-CAFs, genes associated with angiogenesis and vasculature development were enriched in the transition toward Myo-CAF, while those related to extracellular matrix (ECM) remodeling in the transition toward Matrix-CAF. Together, CAFs undergo extensive lineage transitions during LCBM, which is associated with their divergent functions.

Hypoxia promotes the lineage transition of CAFs

Hypoxia is a significant hallmark of tumors and plays a pivotal role in tumor metastasis.^34–36 Therefore, we hypothesized that the lineage transition of CAFs is driven by hypoxia. We first analyzed the bulk RNA-sequencing data and observed that hypoxia-related pathways are highly enriched in the progression and metastasis of lung cancers (Fig. S2A-S2C). In scRNA-seq data, the overall enrichment of hypoxia pathways in CAFs of LCBM is significantly higher than that in the primary sites (Figure 3(d,e)). After examining each CAF subgroup (Figure 3(f)), we found that hypoxia-related pathways are statistically enriched in Matrix-CAFs and Myo-CAFs.

Hypoxia-induced factors (HIF) mediate the cellular responses to hypoxia.³⁴ In LCBM, both HIF1A and HIF2A are significantly upregulated in CAFs (Figure 3(g)). To assess the overall activity of these two HIF pathways, we defined key genes downstream of HIF1A and HIF2A (Table S3) and compared pathway scores between primary and brain metastatic CAFs. HIF-1α and HIF-2α pathway activity are both significantly upregulated in LCBM (Figure 3(h)). In the pseudotime analysis, HIF2A, but not HIF1A, is notably upregulated in branch A (Figure 3(i)), suggesting that HIF-2α may play a more predominant role in the lineage changes of CAFs.

To confirm the aforementioned findings, we integrated datasets of GSE112996 and GSE141685 which comprise primary lung cancers and LCBM (Fig. S2D). Consistently, we observed the significant upregulation of HIF2A (p = 0.018) but not HIF1A (p = 0.75) (Figure 3(j)). This was also validated by IHC staining (Figures 3(k) and S2e). Furthermore, this was also observed in paired samples of peripheral and brain metastases lesions of melanoma from the DFCI cohort (Figure 3(j)). Taken together, our results imply that hypoxia pathways could regulate the lineage transition of CAFs in solid tumor brain metastasis.

HIF-2α activation drives the lineage transition of CAFs

To investigate the mechanisms of HIF2A mediated pathways in LCBM, we first examined the expression pattern in various cell types of LCBM. We found that HIF1A and HIF2A are predominantly distributed in CAFs, whose expression levels are significantly increased in LCBM (Figure 4(a)). This was validated by the protein levels of HIF-1α and HIF-2α by mIF staining of paired lung cancer primary/brain metastasis samples (Figure 4(b)).

Figure 4. Hypoxia-induced HIF-2α activation drives lineage transition of CAFs.

(a) Violin plot shows the expression level of HIF1A and HIF2A between primary lung cancer (tLung) and brain metastasis (mBrain), comparison was also made between groups by t-test. *p value < 0.05, **p value < 0.01, ***p value < 0.001, ****p value < 0.0001, ns: not significant. (b) Representative multiplex immunofluorescence (mIF) images of nuclei (blue), α-SMA (green), HIF-1α (purple), and HIF-2α (red) in lung cancer primary (up) and brain metastasis (bottom) sites, Scale bars: 100 μm. (c) Dot plot of most relevant modules to primary lung cancer and brain metastasis. (d) Box plot of HIF-2α pathway enrichment scores between primary lung cancer and brain metastasis site among all CAF types, p values are calculated by t-test. (e) RSS plot of most relevant modules to different CAF types. (f) Feature plot of HIF-2α and SIX1 modules in primary lung cancer and brain metastasis. (g) Core regulation network for HIF-2α and SIX1 modules related to Vessel development and ECM organization.

We then employed SCENIC analysis to identify key transcription factors (TF) modules which include the TFs and their downstream targets associated with CAF plasticity and lineage transition. Indeed, the HIF-2α module is highly relevant to LCBM (Figure 4(c)). The pathway scores of HIF-2α downstream genes are significantly increased in Myo-CAF and Matrix-CAF in LCBM (Figure 4(d)). In addition, the HIF-2α module is the most specific marker for Myo-CAF, suggesting its pivotal role in reshaping the CAF landscape (Figure 4(e)). Of note, the HIF-1α module is not identified in SCENIC analysis. It should be noted that the SIX1 module, a hypoxia-responding aerobic glycolysis regulating module,³⁷ is also enriched in LCBM (Figure 4(f)). Therefore, we analyzed the metabolic characteristics of the three CAFs and found that Myo-CAFs are metabolically inactive, primarily relying on oxidative phosphorylation (OXPHOs). In contrast, Matrix-CAFs exhibit enhanced glycolytic and OXPHO activity compared to Naive-CAFs, with increased activity in the pentose phosphate pathway and amino sugar and nucleotide sugar metabolism, reflecting their rich intracellular biosynthetic processes (Fig. S3a). In LCBM, the overall aerobic glycolysis is significantly upregulated, and Matrix-CAF has the highest aerobic glycolysis level (Fig. S3b-c). Analysis of the regulatory network of HIF-2α and SIX1 revealed the significant roles of these two transcription factors in angiogenesis and ECM remodeling (Figure 4(f)), shedding light on potential regulatory mechanisms of hypoxia on CAFs. Taken together, we specified HIF-2α and SIX1 as core TFs which are associated with the lineage transition and metabolic reprogramming of CAFs in response to hypoxia in LCBM.

CAFs in LCBM promote tumorigenesis by interacting with TME

As stated above, HIF-2α and SIX1 are associated with angiogenesis and ECM remodeling of CAFs in LCBM. So, we speculated that the signaling transduction patterns between CAFs and other TME components may be reshaped in LCBM. We examined the cell–cell communication and found that the interaction numbers between CAFs and endothelial cells are significantly increased (Figure 5(a)). Signal analysis shows that the communication probability of PIGF-VEGFR1 is significantly upregulated in the interplay between Myo-CAF and endothelial cells (Figure 5(b,c)). PIGF (Placental Growth Factor) is a direct downstream regulatory target of HIF-2α (Figure 4(g)). In LCBM, Myo-CAFs actively interact with endothelial cells and stimulate vessel growth, which was validated by increased colocalization of Myo-CAFs (marked by α-SMA) and endothelial cells (marked by CD31) (Figure 5(d)). Our results provided evidence that Myo-CAFs promote angiogenesis by engaging with endothelial cells through the PIGF-VEGFR1 pathway in LCBM.

Figure 5. Transited CAFs in LCBM interact with TME and malignant cells.

(a) Changes of interaction numbers during brain metastasis among CAFs, endothelial cells, and epithelial cells. Specific numbers of changes are also presented. (b) Bubble plot of L-R pairs with increased communication probability during NSCLC brain metastasis. (c) Violin plot shows the expression level of PIGF and VEGFR1 between primary lung cancer (tLung) and brain metastasis (mBrain) among different cell types. (d) left: Representative mIF images of nuclei (blue), α-SMA (green), PanCK (white), and CD31 (red) in lung cancer primary site (up) and paired brain metastasis site (bottom), Scale bars: 100 μm. right: Violin plot showing the percentage of Myo-CAF (green)/Endothelial cells (red) overlap between lung cancer primary site and brain metastasis, the p values are calculated by the unpaired t-test. (e) Relative information flow from CAFs to macrophages in LCBM. (f) String diagram of SPP1, PTN, CCL signal pathway in LCBM. The most significant L-R pairs are labeled below. (g) Branched heatmap of differentially expressed genes between Branch A and Branch B with Metascape annotation of enriched terms of all function modules. (h) Violin plot shows the expression level of CCL3, BMP1, and FAP between primary lung cancer (tLung) and brain metastasis (mBrain) among different cell types. (i) left: UMAP plots of CytoTRACE stemness score for all malignant cells. middle: UMAP plots of identified CSCs at primary state (tLung) and brain metastatic state (mBrain). right: Comparison of CSC proportion at primary state (tLung) and brain metastatic state (mBrain), *p value < 0.05. (j) Metascape enrichment of upregulated pathways of CSCs in LCBM. (k) SCAP analysis of metabolic pathways in Matrix-CAF between lung primary site and brain metastasis. (l) Representative mIF images of nuclei (blue), COL1A1 (green), and α-SMA (red) on samples of TKI treatment-naive and TKI treatment-resistant patients, Scale bars: 100 μm. (m) Violin plot showing the percentage of Myo-CAF in TKI naïve and resistant cases, **p value < 0.01.

Several studies have uncovered a close interaction between CAFs and the immune microenvironment,^38,39 highlighting a particularly tight interplay between CAFs and macrophages. We studied the changes in the interaction between CAFs and macrophages during LCBM and observed a significant increase in the SPP1, PTN, and CCL pathways (Figure 5(e,f)). We also traced the lineage changes of macrophages in LCBM and identified two distinctly different fate branches, one of which is to SPP1⁺ macrophages (Figure 4(d)), which is defined as an adverse prognostic factor in multiple studies.^38,40 FAP⁺ CAFs were discovered to regulate the chemotaxis of SPP1⁺ macrophages through CCL3 in colorectal cancer, leading to poor prognosis.⁴¹ Our research aligned with this finding and observed an upregulation of CCL3 and FAP expression in Matrix-CAFs within LCBM (Figure 5(f)). Also, we clustered relevant genes related to macrophage lineage differentiation to SPP1⁺ phenotype, and found that receptor tyrosine kinase (RTK) pathway and NF-κB mediated pathways may be driving factors (Figure 5(g)). BMP is a crucial activator of the RTK protein and NF-κB pathways. We observed significant upregulation of BMP1 and other RTK ligands in Matrix-CAF in LCBM (Figures 5(h) and S4e). Our study reveals that the hypoxia and HIF-2α regulated lineage transition of CAFs promotes macrophage differentiation toward the pro-tumor SPP1⁺ phenotype through upregulation of CCL3 and BMP1.

Collectively, our results reveal that under hypoxia and HIF-2α regulation, Myo-CAFs promote angiogenesis by actively interacting with endothelial cells through the PIGF-VEGFR1 pathway in LCBM. On the other hand, Matrix-CAFs may drive macrophage lineage transition toward the SPP1⁺ phenotype.

CAFs in LCBM promote tumorigenesis by nurturing cancer stem cells

Cancer stem cells (CSC) are a core component in sustaining tumor growth, heterogeneity, and therapy resistance, serving as the initiating factor for malignant progression and metastasis of tumors. Tumor cells with multidirectional differentiation potential and high stemness settle in the brain, leading to LCBM. Several studies have reported a close interconnection between CAFs and tumor cells in terms of cellular communication and metabolism.⁴² The lineage alterations of CAFs may play a crucial role in supporting the invasiveness of CSC. We defined the genuine malignant tumor cells from epithelial cell subgroups through copy number analysis (CNA) and evaluated the stemness of all tumor cells using the CytoTRACE algorithm (Figures 5(g) and 4(f)). In LCBM, the proportion of CSC (defined as the top 10% of tumor cells by CytoTRACE stemness score) significantly increased (Figure 5(i)), indicating their crucial role in LCBM development and drug resistance. Pathway enrichment analysis of upregulated genes related to CSCs in LCBM revealed a prominent elevation in amino acid metabolism pathways. These features resemble amino acid metabolism pathways observed in Triple-Negative Breast Cancer (TNBC), a subtype of breast cancer known for its poor prognosis and high propensity for brain metastasis (Figure 5(j)). We identified a significant increase in amino acid synthesis pathways in Matrix-CAFs (Figure 5(k)). Several studies have reported metabolic reprogramming in CAFs induced by tumor cells, promoting tumor growth by synthesizing and secreting the amino acids required by tumor cells.^42,43 We speculate that similar mechanisms may also be present in LCBM.

The issue of targeted therapy resistance in LCBM has always been a hot topic. Tyrosine kinase inhibitors (TKI) are now the most common drugs for treating lung cancer, but ultimately all treated LCBM patients will develop resistance to TKI.⁷ We collected samples from five TKI-naive and five TKI-resistant LCBM patients for mIHC staining. We found a significant increase in Myo-CAF in samples of LCBM that recurred after treatment (Figure 5(l,m)). Combined with the promoting role of Myo-CAF in angiogenesis, we speculate that CAF lineage transformation enhances angiogenesis and tumor vascular mimicry to support tumor cell resistance to targeted drugs. Collectively, our results revealed that the adaption of CAFs to hypoxia under HIF-2α regulation is pro-tumorigenic, by directly facilitating the growth of CSCs and enhancing drug resistance.

Discussion

Our investigation reveals the presence of diverse CAF populations during the metastatic progression of lung cancer. Although the stromal composition and immunosuppressive characteristics of CAFs in primary cancers have been intensively studied and unveiled,^42,44 the CAF signatures at different stages of metastasis remain poorly understood. By the integrated analysis of scRNA-seq and bulk data, we provided a comprehensive exploration of CAF dynamics during lung cancer metastasis, illustrating a unique CAF landscape of LCBMs compared with their primary counterparts. CAFs undergo substantial lineage changes and functional reprogramming during brain metastasis, primarily mediated through the hypoxia-inducible factor (HIF) pathway (Figure 6).

Figure 6. Lineage transition of CAFs in LCBM.

Hypoxia-induced HIF-2α activation drives distinctive lineage transition from Naive-CAF to Myo-CAF and Matrix-CAF during LCBM. Transited CAFs could boost tumor cell proliferation, induce metabolic reprogramming, and promote angiogenesis.

Concordantly, our present work reinforces the previous findings of the plasticity of fibroblasts and CAFs in response to different environmental stimuli, particularly hypoxic TME.^17,45 Fibroblasts, normal in their context, can be “educated” by the pro-tumor TME, differentiating into myofibroblast-like CAF and Matrix-CAF, thus promoting tumor and facilitating metastasis. The connection between hypoxic environments and fibroblast plasticity has been observed in various non-tumor pathological settings, emphasizing the role of signaling pathways such as VEGF and PDGF.⁴⁶ In the context of tumors, the intricate interplay between hypoxia-induced pathways and CAF plasticity remains largely unexplored. Our results emphasize the role of the HIF-2α pathway in mediating the plasticity of CAFs, thereby promoting directed differentiation and functional specialization. This process ultimately facilitates the brain metastasis of lung cancer. HIF signaling is a major pathway responding to hypoxia.³⁴ Given that HIF-1α is associated with short-term acute hypoxia pathways and HIF-2α is associated with long-term chronic hypoxia pathways,^47,48 we propose that this observation could be attributed to the tumor gradually adapting to prolonged hypoxia during metastasis and “switching on” to long-term hypoxia regulatory processes.

Targeting TME has emerged as a promising avenue in cancer therapy. Successful applications of immune checkpoint inhibitors underscore the importance of reversing pro-tumorigenic immunosuppressive microenvironments. Our study suggests that targeting HIF-2α, a major player in the transition of CAF lineages to a malignant phenotype, holds therapeutic potential in LCBM. The HIF-2α pathway plays a significant role in multiple biological processes and diseases,^49–52 and the blockade of this pathway has recently been proven to be an effective treatment. The FDA-approved HIF-2α antagonist, Belzutifan (MK-6482), initially designed for VHL syndrome-associated renal clear cell carcinoma and other tumors, may offer a viable option for the prevention and treatment of LCBM.⁵³ Recently, Belzutifan has been approved by the FDA for previously treated advanced renal clear cell carcinoma, indicating that the HIF-2α pathway also plays a significant role in somatically mutated tumors.⁵⁴ The observed changes in the CAF population and function during the progression of lung cancer, driven by the HIF-2α pathway, present unique insights into the underlying mechanisms of LCBM and propose a potential treatment target. A future clinical trial to verify the efficacy of HIF-2α blockade on LCBM is hence needed.

The manuscript still bears some inadequacies. Owing to the inherent heterogeneity of CAFs, a unified classification for CAFs is currently lacking. Hence, the definition of CAF grouping and the analysis of its evolutionary patterns as outlined in this manuscript necessitate validation on a larger scale. The mechanisms we proposed shall be strengthened by further in-depth experimental validation.

In summary, our study sheds light on the distinctive changes in CAF cell population and cellular function driven by the HIF-2α pathway during the progression of lung cancer, particularly in the context of brain metastasis. Our work may serve as a foundation for further exploration of the role of CAF plasticity in tumor progression and the development of effective strategies for the prevention and treatment of lung cancer brain metastasis.

Abbreviations

LCBM	=	lung cancer brain metastasis
CAFs	=	cancer-associated fibroblasts
TME	=	tumor microenvironment
CNS	=	central nerve system
bulk-seq	=	bulk sequencing
scRNA-seq	=	single-cell RNA sequencing
H&E	=	hematoxylin-eosin staining
IHC	=	immunohistochemistry
mIF	=	multiplex immunofluorescence
CSC	=	cancer stem cell
TKI	=	tyrosine kinase inhibitor

Authors’ contributions

Conceptualization, W.H., Z.Z. and Y.M.; methodology, M.F., M.Y., L.Z. and Y.F.; formal analysis, M. F. and L.Z.; investigation, M.F., Z.S., W.H. and M.Y.; resources, M.Y., Z.S., H.W. and M.Y.; data curation, M. F., L.Z. and F.Y.; writing – original draft, M.Y., M.F.; writing – review & editing, X.Z., W.H. and Y.M.; visualization, M. F., L.Z. and M.Y.; supervision, W.H. and Y.M.; funding acquisition, W.H. and Y.M.

Code availability

The code used for all processing and analysis is available upon request.

Ethics approval

The ethics committee of Huashan Hospital approved the present study. Written consents were obtained from all patients.

Acknowledgments

We sincerely thank Medical Science Data Center in Shanghai Medical College of Fudan University for providing data analysis platforms. We thank Dr. Fengyun Zheng and Core Facility of Shanghai Medical College of Fudan University for support with slide image scanning and data analysis.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Data availability statement

Raw data of all datasets analyzed in this paper are publicly accessible and can be assessed according to the method given in the paper.

Supplementary material

Supplemental data for this article can be accessed online at https://doi.org/10.1080/2162402X.2024.2356942

Correction Statement

This article has been republished with minor changes. These changes do not impact the academic content of the article.

Additional information

Funding

The study was funded by the National Key R&D Program of China [2022YFC2406703, 2022YFC3401600], CAMS Innovation Fund for Medical Sciences (CIFMS)2022-I2M-C&T-B-112, and Fuqing Scholar Program of Shanghai Medical College.