Single-cell and spatial transcriptomics reveal a potential role of ATF3 in brain metastasis of lung adenocarcinoma

Introduction

Lung cancer is the leading cause of cancer-related death worldwide, wherein lung adenocarcinoma (LUAD) is the most common type, comprising approximately 40% of lung cancers (1). About 30–50% of patients with advanced LUAD develop brain metastasis (BrM) (2). However, treatment of LUAD BrM has been a clinical challenge. Surgical resection is only applicable to localized and small lesions (3). While systemic treatments, such as chemotherapy and targeted therapy, have achieved certain efficacy, the prognosis of LUAD patients with BrM remains poor (4,5). Therefore, exploring the mechanisms driving LUAD BrM is of great significance.

The epithelial-mesenchymal transition (EMT) is known to play crucial roles in tumor metastasis (6). Notably, the mesenchymal-epithelial transition (MET), a biological process in which cells switch from a motile, mesenchymal phenotype to a stationary, epithelial phenotype, is also important for tumor metastasis (7). Current research has highlighted that MET enables disseminated tumor cells at metastatic sites to regain epithelial traits and initiate secondary tumor growth (8). Thus, targeting the MET program holds promise in impeding the colonization of metastatic cells (7). However, studies on the MET in LUAD BrM are limited, and key regulators driving this process remain largely unknown.

The tumor microenvironment (TME) comprises immune and stromal cells surrounding tumor cells, along with extracellular matrix and a plethora of molecular signals, thus serving as a pivotal factor in tumor development and progression (9). The TME also plays vital roles in BrM, providing a sanctuary for metastatic cells through various mechanisms (10). Recently, single-cell RNA sequencing (scRNA-seq) has been widely used to explore the tumor ecosystem. Some studies have dissected LUAD BrM and its TME at single-cell resolution. However, these current studies focus on comparing BrM to primary lung tumors, with limited exploration of tumor cell evolution and their interactions with the TME (11,12). Additionally, scRNA-seq causes the inevitable loss of spatial architecture, which limits further research. Spatial transcriptomics (ST), as a complementary technology, can be integrated with scRNA-seq to gain a deeper understanding of the spatial distribution and interactions of cell components within the tumor ecosystem (13).

Thus, we conducted scRNA-seq and ST analyses to dissect the mechanisms underlying BrM of LUAD. Our findings reveal the MET process of tumor cells within BrM. Notably, the tumor cell-intrinsic ATF3 may promote BrM through inducing the MET of tumor cells. Additionally, macrophages may also promote the MET process through LGALS3/ANXA2-mediated cell-cell interaction. Collectively, this study sheds light on the LUAD BrM and identifies ATF3 and LGALS3/ANXA2 as potential therapeutic targets. We present this article in accordance with the MDAR reporting checklist (available at https://tlcr.amegroups.com/article/view/10.21037/tlcr-24-784/rc).

Methods

Tumor specimens

This study was approved by the Ethics Committee of Shanghai General Hospital (approval #2021SQ054) and was conducted in accordance with the Declaration of Helsinki (as revised in 2013). Written informed consent was obtained from all participants in the study. Fresh BrM tissues were obtained by surgical resection from patients diagnosed with LUAD BrM (Table S1).

Tissue processing for scRNA-seq

The tissue was prepared into single-cell suspensions within 2 hours after surgical resection. Firstly, resected tissue was rinsed with PBS, cut into small pieces of 1–2 mm³, and transferred to a centrifuge tube containing Dulbecco’s modified Eagle’s medium (DMEM) supplemented with 1 mg/mL collagenase (C5138, Sigma-Aldrich). The tissue was dissociated at 37 ℃ with shaking at 100 rpm for 45 min. Subsequently, the digested tissue suspension was sequentially filtered through 70 and 40 µm cell strainer. After centrifugation at 400 g for 5 min, the cell pellet was resuspended in red blood cell lysis buffer (#420301, BioLegend) at 4 ℃ for 2 min. After washing with PBS, the cell pellet was resuspended in PBS containing 1% serum. Cell viability was measured using an automated cell counter, and cell suspensions with more than 85% viability were further used for scRNA-seq.

Single-cell RNA library preparation and sequencing

The single-cell RNA libraries were constructed using Chromium Single Cell 3' v3 Kit (10× Genomics) following the experimental protocol provided by the manufacturer. Library quality control was performed using the Agilent 2100 Bioanalyzer system. The libraries were sequenced on the Illumina NovaSeq 6000 PE150 platform and mapped to the GRCh38 human reference genome by using Cell Ranger (version 3.0.1).

The scRNA-seq dataset from open database

To verify the findings from our data, scRNA-seq data of 10 cases of LUAD BrM originated from GSE131907 were utilized as a validation set. The raw UMI matrix was downloaded from the GEO database (https://www.ncbi.nlm.nih.gov/geo/).

Preprocessing of scRNA-seq data and cell type annotation

The raw gene expression matrices of each sample were transformed into Seurat objects using the Seurat R package (v5.0.1) (14). The DoubletFinder (v2.0.4) was utilized to identify doublets, and the decontX (v1.0.0) was applied to remove ambient RNA (15). After quality control and normalization, the FindVariableFeatures function was used to select highly variable genes. The data from different samples were integrated via the anchor-based approach by Seurat’s FindIntegrationAnchors and IntegrateData function. Subsequently, principal component analysis (PCA) was conducted by RunPCA, and dimensionality reduction and clustering were performed by the FindNeighbors and FindClusters functions with the first 30 principal components (PCs). Finally, the single-cell data were visualized using two-dimensional Uniform Manifold Approximation and Projection (UMAP). To annotate the cell clusters, the FindAllMarkers function was used to compute highly expressed genes for each cluster with parameters logfc.threshold =2 and min.pct =0.5. In the subclustering analysis of myeloid cells, the FindAllMarkers was also employed to obtain marker genes for each cluster with parameters adjusted to logfc.threshold =1 and min.pct =0.2. Manual annotation was then performed based on canonical markers obtained from the CellMarker database.

ST tissue processing, library preparation and sequencing

For one tumor (LUAD_BrM-3) with sufficient tissue, some of it was embedded in optimal cutting temperature compound (OCT, #4583, Sakura) and frozen rapidly for ST sequencing. The OCT-embedded tissue was sliced at −12 ℃ into a 10 µm-thick tissue section to match the capture region. Following the manufacturer’s protocol, the specimen was mounted on a Visium Spatial slide (10× Genomics) containing RNA-binding capture probes. After fixation and H&E staining, the section was imaged under a microscope to confirm coverage of the target region. After permeabilization, the mRNAs released from cells in the tissue were captured by the probes, and subsequently used as templates for reverse transcription. Then, the Visium Spatial Library Construction Kit (10× Genomics) was applied to construct the spatial libraries. The libraries were sequenced on the Illumina NovaSeq 6000 PE150 platform. Space Ranger (version 1.1.0) was used to process the raw sequencing reads and map the spatial libraries to the GRCh38 human reference genome.

ST data processing

The Seurat R package (v5.0.1) (14) was used to analyze the output files from Space Ranger. The SCTransform function was employed for normalization, scaling, and detecting highly variable features in ST data. Subsequently, PCA and dimensionality reduction was also performed using RunPCA and FindNeighbors with the first 20 PCs. The FindClusters was executed with the resolution =1.5 for optimal clustering.

Geneset scoring

The gene signature of each major cell type in the scRNA-seq data was obtained using the FindMarkers function with parameters logfc.threshold =2 and min.pct =0.5. The gene signature of each macrophage subpopulations (Mac_CCL18, Mac_FCGBP and Mac_SLC2A1) and tumor cluster (t2 and t6) in the scRNA-seq data was also defined by FindMarkers with parameters adjusted to logfc.threshold=1 and min.pct=0.2. Moreover, previous reported macrophage-related genesets identified from scRNA-seq by Xing et al. (16) and Cheng et al. (17) were used to reflect the functional characteristics of macrophage subpopulations. The lists of epithelial and mesenchymal markers predefined by Taube et al. (18), as well as the hallmark EMT geneset downloaded from MSigDB database were utilized to evaluate the EMT status of tumor clusters. The BMAECs marker genes (ZFAS1, IFI6, C15orf48, SAA1, S100A9, IFI27) defined by Wang et al. (11) were used to assess the presence of BrM-associated BMAEC cells in LUAD BrM tumors. The scoring of all the genesets mentioned above in the scRNA-seq and ST data was conducted by the AddModuleScore function in Seurat with default parameters.

Gene set variation analysis (GSVA)

GSVA, an unsupervised, non-parametric analysis method, was used to evaluate the enrichment of HALLMARK genesets in each macrophage subpopulations. The genesets were downloaded from the MSigDB database on the GSEA website (https://www.gsea-msigdb.org/gsea/msigdb/index.jsp). The AverageExpression function in Seurat was used to calculate the average expression value of genes in each macrophage subpopulations. Then the gsva function in the GSVA R package (v1.50.2) (19) was used to estimate enrichment scores. The result was visualized by the pheatmap function.

Survival analysis

The Cancer Genome Atlas (TCGA) database was used to analyze the prognostic value of specific genes or cell clusters. The gene expression and survival data of TCGA low-grade glioma (LGG) cohort were retrieved from UCSC Xena (https://xenabrowser.net/datapages/). After calculating counts per million mapped reads (CPM) from read counts and integrating with survival data, a total of 467 tumor RNA-seq profiles with survival information were obtained. For macrophage subpopulations or tumor clusters, the gsva function was utilized to calculate the gene signature scores for each sample, which were used to assess the levels of these cell populations in bulk samples. Finally, the survival and survminer R packages were used to calculate and generate Kaplan-Meier survival curves.

Copy number variation analysis

To demonstrate the malignancy of epithelial cells in our scRNA-seq data, copy number variation analysis was performed by the infercnv R package (v1.18.1) (20) on the epithelial clusters. Myeloid, B, and T/NK cell clusters were utilized as references as non-malignant cells with normal copy numbers. Key parameter settings employed during the infercnv analysis included: cutoff=0.1, cluster_by_groups=TRUE, denoise=TRUE, HMM=FALSE.

RNA velocity analysis

RNA velocity analysis was performed using the scVelo python package (v0.3.2) (21). A total of 2,000 highly variable genes were extracted for analysis by the scv.pp.filter_and_normalize function, and moments of each cells were calculated using the scv.pp.moments function with parameters n_pcs=30 and n_neighbors=30. RNA velocity was subsequently calculated by scv.tl.velocity in stochastic mode, and velocity graph was constructed by scv.tl.velocity_graph. The RNA velocity driving genes for each tumor clusters were identified by the scv.tl.rank_velocity_genes function with parameters min_corr=0.4.

Inference of highly active TFs

To identify TFs involved in phenotypic transition of tumor cells, the regulon activity was analyzed using the pySCENIC python package (v0.12.1) (22). Files for the analysis were all downloaded from the cisTarget website (https://resources.aertslab.org/cistarget/), and the program was run with default parameters. The output loom file was imported into the R environment, and the calcRSS function in the SCENIC R package (v1.3.1) was used to calculate regulon specificity scores for TFs in each tumor cluster, which reflected the regulatory activity of TFs.

Cell-cell interaction analysis

NicheNet was applied using nichenetr R package (v2.1.0) (23). The nichenet_seuratobj_cluster_de function was employed to calculate potential ligands in macrophages that could drive the MET process of tumor cells. The top 30 highly active ligands ranked by the area under the precision-recall curve (AUPR) value were displayed by heatmap and used for further analysis. Potential receptors matching these ligands were also visualized in heatmap by the ligand_receptor_heatmap command.

Cell culture and plasmids

The human LUAD cell lines H2009 (RRID: CVCL_1514), H441 (RRID: CVCL_1561), A549 (RRID: CVCL_0023) and H1755 (RRID: CVCL_1492) were obtained from ATCC (Manassas, VA). All of them were cultured in DMEM added with 10% FBS and 1% penicillin/streptomycin, in a humidified 5% CO₂ environment at 37 ℃. The lentiviral vector pLKO.1-puro (#8453, Addgene) was used to construct the ATF3 shRNA plasmids. The ATF3 shRNA sequences designed are as follows: shATF3-1 (CTTCATCGGCCCACGTGTATT), shATF3-2 (AGATGAGAGAAACCTCTTTAT). Lentivirus production was performed by co-transfecting the shRNA plasmids with packaging plasmids psPAX2 and pMD2.G into 293T cells using HighGene transfection reagent (ABclonal, China). After 48 hours, supernatants with lentiviral particles were collected, filtered through 0.45 µm PVDF filters, and used for cell infection with 8 µg/mL polybrene. Stably transduced cells were then selected with puromycin.

Western blotting

In this study, proteins in cell pellets were extracted using T-PER reagent (ThermoFisher) containing protease inhibitor cocktail, and the Pierce BCA protein assay kit (Thermo Scientific) was used to determine the protein concentration. The western blotting assay was conducted using a standard protocol. A 10% sodium dodecyl-sulfate polyacrylamide gel was used for protein electrophoresis. The total amount of protein per well was around 20 µg. The 0.45 µm nitrocellulose membrane (Cytiva) was used for protein transfer at a constant current of 260 mA for 1.5 hours. Blocking was conducted with 5% non-fat milk at room temperature for 1 hour. The primary antibodies and dilutions ratios used in this study were as follows: ATF3 (1:1,000, #18665, CST, RRID:AB_2827506), ZO-1 (1:1,000, #8193, CST, RRID: AB_10898025), E-cadherin (1:1,000, #3195, CST, RRID:AB_2291471), N-cadherin (1:1,000, #13116, CST, RRID:AB_2687616), Vimentin (1:1,000, #5741, CST, RRID:AB_10695459), and β-actin (1:20,000, A3854, Sigma, RRID:AB_262011). Protein bands were visualized by the FluorChem E Imaging System (ProteinSimple, USA).

Quantitative real-time polymerase chain reaction (qRT-PCR)

Cellular total RNA was extracted using Trizol reagent (Invitrogen) and subsequently reverse-transcribed to cDNA with HiScript III All-in-one RT SuperMix (Vazyme, Nanjing, China), following the manufacturer’s instructions. qRT-PCR was then performed using ChamQ Universal SYBR qPCR Master Mix (Vazyme, Nanjing, China) in the QuantStudio 7 System (Applied Biosystems) with ACTB as an internal control. The primer sequences used are as follows: ATF3 (forward AGGAAAAAGAGGCGACGAGAA; reverse TTTCTGCAGGCACTCCGTC), TJP1 (forward TCGCGTCCCCATGTATGAAG; reverse GAGGACTCTTCGGGATGCTG), CDH1 (forward AGTCACTGACACCAACGATAAT; reverse ATCGTTGTTCACTGGATTTGTG), CDH2 (forward CGATAAGGATCAACCCCATACA; reverse TTCAAAGTCGATTGGTTTGACC), VIM (forward GGACCAGCTAACCAACGACA; reverse GGGTGTTTTCGGCTTCCTCT), ACTB (forward GCTGTGCTATCCCTGTACG; reverse GGGCATACCCCTCGTAGAT).

Transwell assay

In the transwell migration assay, LUAD cells transfected with control or shATF3 plasmids were suspended in serum-free medium at a concentration of 1.5×10⁵/mL. Then, 200 µL of cell suspension was added to the upper chamber, achieving a total cell count of 3.0×10⁴ per well, and 500 µL of medium containing 10% FBS was added to the lower chamber. After incubating at 37 ℃ for 20 hours, the migrated cells were fixed with 4% paraformaldehyde for 30 minutes and stained with 2.5% crystal violet for 20 minutes. Cell number was counted in three random fields. The entire experiment was repeated three times.

Cell proliferation assay

In Cell Counting Kit-8 (CCK-8) assay, cells were seeded at 1,000/well in a 96-well plate and measured every 24 hours after cell attachment. For detection, supernatant was removed and replaced with 100 µL of medium and 10 µL of CCK-8 reagent (TargetMol, USA). After incubating at 37 ℃ for 2 hours, absorbance at 450 nm was measured by the microplate reader. The data are presented as the mean ± standard error of the mean (SEM) of three replicates.

Statistical analysis

For statistical comparisons between two groups, the Student’s t-test was used for Gaussian distributed data, while the Mann-Whitney test was used for non-Gaussian distributed data. For Kaplan-Meier survival analysis, the log-rank test was employed for comparison of two groups. Pearson’s correlation tests were utilized for analyzing the expression correlation between ATF3 and both CDH1 and VIM in 123 non-small cell lung cancer (NSCLC) cell lines. All statistical analyses were conducted in the R environment (v4.3.3) and GraphPad Prism 8.0 software. A P value <0.05 was regarded to be statistically significant.

Results

ScRNA-seq analyses of human LUAD BrMs

To identify mechanisms driving BrM, we conducted scRNA-seq on three BrM tissues from LUAD patients (Figure 1A). After quality control and sample integration, a total of 19,614 cells were identified (Figure 1B). Following dimensionality reduction and clustering, these cells were classified into 16 clusters (Figure S1A) and annotated into 6 major cell types based on the expression of canonical cell markers (Figure 1C). The non-immune cells included epithelial cells, fibroblasts and endothelial cells, while the immune cells comprised myeloid cells, T/NK cells and B cells (Figure 1C,1D and Figure S1B). The cellular composition exhibited inter-tumoral heterogeneity among these samples (Figure 1E and Figure S1C). Epithelial cells constituted the highest percentage (80.25%), followed by myeloid cells (7.06%), with the remaining cell types being less represented (<5%) (Figure 1F). Therefore, we primarily focused on epithelial and myeloid cells within the tumor ecosystems in this study.

Figure 1 ScRNA-seq analyses of LUAD BrMs. (A) Schematic workflow depicting the experimental design. The schematic was created by Figdraw (https://www.figdraw.com/#). (B) UMAP plot showing the sample origin of all single cells from three LUAD BrM tissues. (C) UMAP visualization of all cells (n=19,614) after quality control and dataset integration. Cells are color-coded by major cell types. (D) Bubble chart displaying the expression of representative markers for each major cell type. (E) Histogram indicating the proportion of major cell type within each LUAD BrM sample. (F) Pie chart summarizing the percentage of each cell type within three LUAD BrMs. scRNA-seq, single-cell RNA sequencing; LUAD, lung adenocarcinoma; BrM, brain metastasis; UMAP, uniform manifold approximation and projection; NK, natural killer.

ST of human LUAD BrM

To analyze the spatial architecture of LUAD BrM, we further performed ST analysis on one LUAD BrM specimen (Figure 1A). Accordingly, a total of 3,609 spots in the LUAD_BrM-3 sample were divided into 14 ST_Clusters (Figure 2A and Figure S2A). Given that marker-based analysis was insufficient to identify the predominant cell types distributed within each ST_Cluster (Figure S2B), the gene signature of each cell type defined by scRNA-seq above was applied to score (Figure 2B and Figure S2C). We found high epithelial scores in most ST_Clusters, particularly S7/5/4/1/0/6, which showed tumor-like histopathological features in H&E staining (Figure 2A-2C). The myeloid signature was highly enriched in S11 (Figure 2B,2C). The spatial proximity between S11 and S4/5 suggested close communication between myeloid and tumor cells in LUAD BrM (Figure 2D). Moreover, B cells and stromal cells (fibroblasts and endothelial) were mainly distributed in S8 and S13, respectively (Figure 2B and Figure S2D). The classical marker GFAP for astrocytes and MBP for oligodendrocytes were highly expressed in S10/3/12/2, indicating these regions tended to be the brain tissue (Figure S2E,S2F). Thereby, we depicted the spatial architecture of LUAD BrM by combining scRNA-seq and ST data (Figure 2E).

Figure 2 ST analysis of LUAD BrM. (A) Representative image of H&E staining of the LUAD_BrM-3 ST sample (left), and spatial distribution of all spots (n=3,609) within this sample color-coded by ST_Clusters (right). Scale bar: 500 µm. (B) Violin plots showing the gene signature score across ST_Clusters. The gene signature was defined by highly expressed genes of each major cell type derived from scRNA-seq analysis. ST_Clusters are ranked by the epithelial cell scores. (C) ST plots showing the epithelial and myeloid scores within the ST sample, respectively. (D) ST plots showing the spatial localization of ST_Cluster 4/5 (high epithelial score) and ST_Cluster 11 (high myeloid score), respectively. (E) ST plot showing the spatial distribution of annotated regions based on cell type scoring and histopathological features. ST, spatial transcriptomics; LUAD, lung adenocarcinoma; BrM, brain metastasis; H&E, hematoxylin and eosin; scRNA-seq, single-cell RNA sequencing.

Tumor cell evolution in LUAD BrMs

To analyze the tumor cells, we performed inferCNV on the epithelial clusters (0/1/2/3/4/6/9) (Figure 1C and Figure S1A) and detected obvious copy number aberrations, confirming their malignancy in BrM (Figure S3A). The scRNA-seq data of 15,741 tumor cells were then segregated into 11 tumor clusters through subclustering (Figure 3A). These clusters exhibited notable differences in proportions across samples, underscoring intra- and inter-tumoral heterogeneity of tumor cells (Figure S3B). We focused on clusters 0/2/3/4/6 present in all samples to uncover common regulatory mechanisms (Figure 3B). They co-expressed epithelial and mesenchymal markers (Figure 3C), suggesting a partial EMT state associated with increased plasticity for tumor progression and metastasis (24). Within these clusters, we noticed significantly higher CDH1 and TJP1 and relatively lower VIM expression in t6 compared to t2 (Figure 3C). Geneset scoring further revealed that t6 showed significantly higher epithelial scores and lower mesenchymal and HALLMARK_EMT scores compared to t2 (Figure 3D and Figure S3C). Remarkably, unsupervised RNA velocity analysis illustrated an evolution trajectory starting from t2 and advancing to t6 (Figure 3E), suggesting a switch from mesenchymal to epithelial phenotypes. Wang et al. identified a group of BrM-associated tumor cells named BMAECs, present in both primary lung tumors and BrM tumors (11). The high BMAECs scores further supported that t2 represented early metastatic cells originating from the primary site (Figure S3D). Additionally, spatial localization analysis shown that t2 and t6 were co-localized in myeloid S11 and tumor S5 (Figure S3E,S3F). Thus, these results revealed the MET process of tumor cells within LUAD BrM. To further validate this discovery, we integrated previously published scRNA-seq data of 10 LUAD BrMs and extracted tumor clusters (Figure S4A-S4C). Cluster 0 exhibited the highest t6 scores. In contrast, cluster 7 had the highest t2 scores, along with significantly elevated VIM expression and HALLMARK_EMT scores (Figure S4D-S4G). Considering the proximity of cluster 7 and 0 in the UMAP and high BMAECs scores in cluster 7 (Figure S4C,S4H), these results further supported the MET process of tumor cells within LUAD BrM.

Figure 3 Tumor cell evolution within LUAD BrMs. (A) UMAP visualization of all tumor cells (n=15,741) color-coded by clusters. (B) Box plots displaying the proportion of each tumor cluster among all tumor cells in different LUAD BrM samples. Tumor clusters are ranked by median value. (C) Violin plots showing the CDH1, TJP1 and VIM expression levels among tumor clusters ranked by mean value. Statistical analysis was performed using an unpaired Mann-Whitney test. (D) Violin plots showing the hallmark EMT gene signature score among tumor clusters ranked by mean value. The genesets was obtained from MSigDB database. Statistical significance was computed with the unpaired Mann-Whitney test. (E) RNA velocity stream plot depicting the tumor cell evolution trajectory from cluster t2 to t6. UMAP, uniform manifold approximation and projection; LUAD, lung adenocarcinoma; BrM, brain metastasis; EMT, epithelial-mesenchymal transition; MSigDB, Molecular Signatures Database.

ATF3 may drive MET program in tumor cells during LUAD BrM

To identify key regulators associated with the MET process in LUAD BrM tumor cells, we computed RNA velocity driving genes by scVelo and analyzed highly active TFs by pySCENIC in each tumor cluster (available online: https://cdn.amegroups.cn/static/public/10.21037tlcr-24-784-1.xlsx; https://cdn.amegroups.cn/static/public/10.21037tlcr-24-784-2.xlsx and Figure S5A). We overlapped the velocity driving genes with the top 20 TFs in t2 and t6, revealing ATF3 as the only gene in the intersection (Figure 4A), supported by its relatively high RNA velocity, regulon activity and expression level in t2 and t6 (Figure S5B-S5D). ATF3 also displayed high expression in MET-related clusters 0 and 7 in the validation scRNA-seq dataset (Figure S5E,S5F). These results suggested that ATF3 may play a crucial role in the MET process. To validate this finding, we assessed expression of ATF3 and EMT markers in human LUAD cell lines with high metastatic potential. Western blotting revealed that H2009 and H441 cells with high ATF3 levels exhibited elevated E-cadherin and ZO-1, and minimal vimentin and N-cadherin expression compared to A549 and H1755 cells with low ATF3 levels (Figure 4B). Notably, ATF3 expression in 123 NSCLC cell lines from the CCLE database positively correlated with CDH1 and negatively correlated with VIM expression (Figure 4C,4D). Subsequent knockdown of ATF3 using shRNA in H2009 and A549 cells demonstrated that ATF3 suppression significantly reduced epithelial markers and increased mesenchymal markers expression (Figure 4E,4F and Figure S5G,S5H). These results supported ATF3’s role in promoting the MET program. Functional experiments further revealed that ATF3 knockdown significantly enhanced the migratory capacity of H2009 and A549 cells (Figure 4G and Figure S5I), while suppressing tumor cell proliferation (Figure 4H and Figure S5J), consistent with the finding that tumor cells with an epithelial phenotype typically possess enhanced colonization abilities (24). Collectively, our analyses and experiments demonstrated that ATF3 may drive the MET process in tumor cells, thereby fostering their survival and colonization in LUAD BrM.

Figure 4 ATF3 may promote MET program in LUAD. (A) Venn diagram showing the overlapping genes between RNA velocity driving genes and top 20 highly active TFs in clusters t2 and t6. (B) Western blot analyses of indicated protein expression in human LUAD cell lines, with β-actin used as a loading control. (C,D) Scatter plots indicating the association between the expression levels of ATF3 with CDH1 (C) and VIM (D) in NSCLC cell lines (n=123) from the CCLE database. Cells are color-coded by predefined EMT status. Statistical analysis was conducted using Pearson’s correlation tests. The expression and annotation data were downloaded from the CellMinerCDB database (https://discover.nci.nih.gov/rsconnect/cellminercdb/). (E) qRT-PCR analyses of indicated gene expression in H2009 cells upon ATF3 knockdown. Data are shown as the mean ± SD. (F) Western blot analyses of indicated protein expression in H2009 cells upon ATF3 knockdown, with β-actin used as a loading control. (G) Transwell migration assays of H2009 (upper) and A549 (lower) cells upon ATF3 knockdown. Migrated cells were fixed with 4% paraformaldehyde and stained with 2.5% crystal violet for 20 min. Scale bar: 250 µm. (H) CCK-8 assay of the in vitro proliferation ability of H2009 cells upon ATF3 knockdown. Data are shown as the mean ± SEM of three replicates. Statistical analysis was performed using an unpaired two-sided Student’s t-test. *, P<0.05; **, P<0.01; ***, P<0.001. TF, transcription factor; NSCLC, non-small cell lung cancer; MET, mesenchymal-epithelial transition; shNC, negative control shRNA; LUAD, lung adenocarcinoma; CCLE, Cancer Cell Line Encyclopedia; EMT, epithelial-mesenchymal transition; qRT-PCR, quantitative real-time polymerase chain reaction; SD, standard deviation; CCK-8, Cell Counting Kit-8; SEM, standard error of the mean.

Exploration of macrophages in LUAD BrMs

Next, we re-clustered all myeloid cells and identified three major cell types, including monocytes, macrophages and dendritic cells (DCs) in LUAD BrM ecosystems (Figure S6A,S6B). The macrophages were further divided into Mac_CCL18, Mac_FCGBP, Mac_SLC2A1, and proliferative macrophages based on the specific highly expressed genes (Figure S6A,S6B). These macrophage subpopulations showed inter-tumoral heterogeneity across three samples (Figure S6C). We found that Mac_CCL18, Mac_FCGBP and Mac_SLC2A1 all expressed classical M2 markers, while M1 markers were expressed only in a small proportion of macrophages (Figure S6D). Moreover, using previous macrophage-associated genesets derived from scRNA-seq (16,17), we observed that all macrophage subpopulations exhibited high M2-like, anti-inflammatory, and tumor-associated macrophage scores (Figure S6E). Consistently, GSVA revealed the enrichment of Notch and Hedgehog signaling pathways in Mac_FCGBP, TGFβ signaling pathway in Mac_CCL18, and PI3K-AKT and Wnt signaling pathways in Mac_SLC2A1 (Figure S6F), which are all related to M2 macrophage polarization (25,26). Given the absence of a public LUAD BrM cohort with survival data and the similarity in TME between BrM and gliomas, the prognostic value of these macrophages was evaluated in the TCGA-LGG cohort. Notably, we observed a significant association of these macrophages with poor prognosis of LGG patients (Figure S6G). Additionally, spatial localization analysis revealed that Mac_FCGBP and Mac_SLC2A1 were primarily distributed in regions of myeloid S11 and tumor S4/5, suggesting potential interactions between these macrophages and tumor cells in LUAD BrM (Figure S6H,S6I). To further validate the existence of these macrophage subpopulations, we found that Mac_CCL18, Mac_FCGBP and Mac_SLC2A1 showed relatively high scores in the myeloid cluster of the validation scRNA-seq dataset (Figures S4B,S7A,S7B). Together, our findings indicate that most macrophages exhibit a tumor-promoting M2 phenotype in LUAD BrM.

Macrophages-mediated LGALS3/ANXA2 interaction may promote the MET program

The spatial co-localization of macrophages and tumor cells undergoing MET implies a potential role of macrophage-tumor interactions during the MET process (Figure S3E,S3F and Figure S6H,S6I). To investigate how macrophages regulate the transition of tumor cells from cluster t2 to t6, we further employed Nichenet to identify potential ligands with high regulatory activity in macrophages (Figure 5A and available online: https://cdn.amegroups.cn/static/public/10.21037tlcr-24-784-3.xlsx). These ligands were expressed at different levels and proportions across distinct macrophage subpopulations (Figure 5B). Since macrophages were primarily localized in S11, we focused on highly expressed ligands LGALS3, FGG, and TYROBP in S11 (Figure 5C,5D and Figure S8A,S8B). Nichenet predicted the receptors associated with these three ligands (Figure S8C). ANXA2, the receptor for LGALS3, exhibited elevated expression in tumor clusters t2 and t6 (Figure 5E) and spatial clusters S5 and S11 (Figure 5C,5D and Figure S8D), indicating a possible LGALS3/ANXA2 interaction between macrophages and tumor cells. In the validation scRNA-seq data, we also analyzed ligands in macrophages driving tumor cells from cluster 7 to 0. Notably, LGALS3 also exhibited high regulatory activity in this transition (Figure S8E). Notably, high expression of LGALS3 and ANXA2 correlated with unfavorable patient prognosis of patients (Figure 5F). Collectively, these findings suggest that macrophages may facilitate the MET program through LGALS3/ANXA2-mediated cell-cell interaction.

Figure 5 Macrophage-tumor cell interactions in LUAD BrM. (A) Heatmap showing the top 30 active ligands in macrophages ranked by AUPR. (B) Bubble chart displaying the expression of top 30 active ligands among three macrophage subpopulations. (C) ST plots showing the expression of LGALS3 (left) and ANXA2 (right) within the LUAD BrM tissue. (D) Violin plots showing the expression of LGALS3 and ANXA2 across ST_Clusters ranked by mean value, respectively. (E) Bubble chart displaying the expression levels of predicted receptors of LGALS3, FGG and TYROBP across tumor clusters. (F) Kaplan-Meier survival analyses of the expression of LGALS3 and ANXA2 in the TCGA-LGG cohort, respectively. P value was calculated by log-rank test. AUPR, the area under the precision-recall curve; ST, spatial transcriptomics; LUAD, lung adenocarcinoma; BrM, brain metastasis; TCGA, The Cancer Genome Atlas; LGG, low-grade glioma.

Discussion

In this study, we unveiled the MET process in LUAD BrM at single-cell resolution. It has been found that EMT is actually a spectrum where tumor cells reside between epithelial and mesenchymal phenotypes, known as partial EMT (27). Unlike complete EMT, cells undergoing partial EMT display increased plasticity and migrate more aggressively in the form of cell clusters (28). We identified that tumor cells in LUAD BrM indeed manifested a partial EMT state and revealed the presence of MET program in these cells. We found that tumor cell transition from t2 to t6 was accompanied by a significant increase in CDH1 expression and epithelial scores. Interestingly, apart from t6, tumor clusters like t0 and t3 also displayed epithelial traits (Figure 3C). Although clear evolution trajectories from t2 to these clusters were not identified, their transcriptional similarities suggest potential evolution from t2 or a rare tumor cluster that were not detected in this study.

Currently, the mechanisms underlying the MET program remain elusive (7). Our study identified ATF3 as a potential regulator of the MET program. As a stress-response gene, ATF3 participates in regulating tumor cell proliferation, apoptosis, ferroptosis, and metastasis through diverse mechanisms (29,30). ATF3 has been found to exert dual roles in various cancer types, including lung cancer (31). For example, ATF3 induces DR5-dependent apoptosis in NSCLC cells, yet its overexpression has been found to promote cell malignancy (32,33). Our study revealed elevated ATF3 levels in NSCLC cells with epithelial traits. ATF3 knockdown in LUAD cells inhibited proliferation but enhanced EMT and migration. These findings suggest that ATF3 may suppress metastasis in primary lung lesions and facilitate metastatic cells to regain the epithelial phenotype, thereby promoting their colonization and survival in BrM. Moreover, the pro-metastatic factor adenylate kinase-4 (AK4) inhibits ATF3 expression in lung cancer, and low ATF3 expression associates with poor prognosis in clinical cohorts, further supporting our findings (34). However, our study only conducted two-dimensional transwell assay in vitro. Further research utilizing three-dimensional models, such as tumor spheroid invasion assay, as well as in vivo mouse metastasis models with bioluminescence or magnetic resonance imaging techniques, would provide a more comprehensive understanding of the actual role of ATF3 in LUAD BrM. Once its function is clearly defined, targeting ATF3 could reverse the epithelial phenotype of tumor cells and potentially impede disease progression in LUAD BrM patients.

Consistent with previous scRNA-seq studies on LUAD BrM (12,35), we discovered myeloid cells as the predominant non-tumor components, suggesting their pivotal roles in the LUAD BrM TME. We identified three macrophage subpopulations with pro-tumoral M2-like traits. Interestingly, our findings suggest that macrophages might induce tumor cell MET through LGALS3/ANXA2-mediated cell-cell interaction. LGALS3 (galectin-3), a glycoprotein, is a critical target linked to tumor growth, metastasis, and immune suppression, notably high expressed in M2 macrophages (36,37). It also serves as a cell-surface adhesion molecule mediating cell-cell interactions (38). ANXA2, a cell membrane phospholipid-binding protein, participates in various cancer physiological processes and is overexpressed in multiple tumor types, including NSCLC (39). Notably, the LGALS3/ANXA2 interaction observed in breast cancer cells has been shown to promote cell survival (38). However, the precise roles of LGALS3/ANXA2 in LUAD BrM warrant further experimental validation. Techniques like Co-immunoprecipitation or surface plasmon resonance could be employed to assess the binding affinity between LGALS3 and ANXA2. Knockout of LGALS3 in macrophages or ANXA2 in tumor cells would clarify whether this interaction plays an inductive role in the MET program. Given the promising efficacy of galectin-3 inhibitor and ANXA2 antibody in preclinical tumor models (36,39), this interaction also presents a potential therapeutic target for LUAD BrM patients.

There are several limitations in this study. First, the study is limited by a small sample size, with only three tumors were examined, which may affect the generalizability and reproducibility of our findings. Future research with more samples will strength our findings. Second, while bioinformatics analyses have provided several insights into LUAD BrM, further experimental investigations are needed to validate these hypotheses.

Conclusions

Our study reveals potential roles of ATF3 and macrophage-mediated LGALS3/ANXA2 interaction during BrM of LUAD, shedding light on the mechanisms of LUAD BrM and providing potential therapeutic targets.

Acknowledgments

Funding: This work was supported by the National Natural Science Foundation of China (grant/award Nos. 82072571, 82273008, 82103012), Nurture projects for basic research of Shanghai Chest Hospital (grant/award No. 2021YNJCM09), Natural Science Foundation of Shanghai (grant/award No. 24ZR1464100) and Science and Technology Commission of Shanghai Municipality Experimental Animal Research Fund (grant/award No. 24141900500).

Footnote

Reporting Checklist: The authors have completed the MDAR reporting checklist. Available at https://tlcr.amegroups.com/article/view/10.21037/tlcr-24-784/rc

Data Sharing Statement: Available at https://tlcr.amegroups.com/article/view/10.21037/tlcr-24-784/dss

Peer Review File: Available at https://tlcr.amegroups.com/article/view/10.21037/tlcr-24-784/prf

Conflicts of Interest: All authors have completed the ICMJE uniform disclosure form (available at https://tlcr.amegroups.com/article/view/10.21037/tlcr-24-784/coif). The authors have no conflicts of interest to declare.

Ethical Statement: The authors are accountable for all aspects of the work in ensuring that questions related to the accuracy or integrity of any part of the work are appropriately investigated and resolved. The study was conducted in accordance with the Declaration of Helsinki (as revised in 2013). The study was approved by the Ethics Committee of Shanghai General Hospital (approval #2021SQ054) and informed consent was obtained from all individual participants.

Open Access Statement: This is an Open Access article distributed in accordance with the Creative Commons Attribution-NonCommercial-NoDerivs 4.0 International License (CC BY-NC-ND 4.0), which permits the non-commercial replication and distribution of the article with the strict proviso that no changes or edits are made and the original work is properly cited (including links to both the formal publication through the relevant DOI and the license). See: https://creativecommons.org/licenses/by-nc-nd/4.0/.

References

1. Thai AA, Solomon BJ, Sequist LV, et al. Lung cancer. Lancet 2021;398:535-54. [Crossref] [PubMed]
2. Zhao Y, Gu S, Li L, et al. A novel risk signature for predicting brain metastasis in patients with lung adenocarcinoma. Neuro Oncol 2023;25:2207-20. [Crossref] [PubMed]
3. Yen CT, Wu WJ, Chen YT, et al. Surgical resection of brain metastases prolongs overall survival in non-small-cell lung cancer. Am J Cancer Res 2021;11:6160-72. [PubMed]
4. Rybarczyk-Kasiuchnicz A, Ramlau R, Stencel K. Treatment of Brain Metastases of Non-Small Cell Lung Carcinoma. Int J Mol Sci 2021;22:593. [Crossref] [PubMed]
5. Yuan J, Cheng Z, Feng J, et al. Prognosis of lung cancer with simple brain metastasis patients and establishment of survival prediction models: a study based on real events. BMC Pulm Med 2022;22:162. [Crossref] [PubMed]
6. Bakir B, Chiarella AM, Pitarresi JR, et al. EMT, MET, Plasticity, and Tumor Metastasis. Trends Cell Biol 2020;30:764-76. [Crossref] [PubMed]
7. Yao D, Dai C, Peng S. Mechanism of the mesenchymal-epithelial transition and its relationship with metastatic tumor formation. Mol Cancer Res 2011;9:1608-20. [Crossref] [PubMed]
8. Łazarczyk M, Mickael ME, Skiba D, et al. The Journey of Cancer Cells to the Brain: Challenges and Opportunities. Int J Mol Sci 2023;24:3854. [Crossref] [PubMed]
9. de Visser KE, Joyce JA. The evolving tumor microenvironment: From cancer initiation to metastatic outgrowth. Cancer Cell 2023;41:374-403. [Crossref] [PubMed]
10. Aili Y, Maimaitiming N, Qin H, et al. Tumor microenvironment and exosomes in brain metastasis: Molecular mechanisms and clinical application. Front Oncol 2022;12:983878. [Crossref] [PubMed]
11. Wang Z, Wang Y, Chang M, et al. Single-cell transcriptomic analyses provide insights into the cellular origins and drivers of brain metastasis from lung adenocarcinoma. Neuro Oncol 2023;25:1262-74. [Crossref] [PubMed]
12. Wu Y, Yang F, Luo S, et al. Single-cell RNA sequencing reveals epithelial cells driving brain metastasis in lung adenocarcinoma. iScience 2024;27:109258. [Crossref] [PubMed]
13. Marx V. Method of the Year: spatially resolved transcriptomics. Nat Methods 2021;18:9-14. [Crossref] [PubMed]
14. Hao Y, Stuart T, Kowalski MH, et al. Dictionary learning for integrative, multimodal and scalable single-cell analysis. Nat Biotechnol 2024;42:293-304. [Crossref] [PubMed]
15. Yang S, Corbett SE, Koga Y, et al. Decontamination of ambient RNA in single-cell RNA-seq with DecontX. Genome Biol 2020;21:57. [Crossref] [PubMed]
16. Xing X, Yang F, Huang Q, et al. Decoding the multicellular ecosystem of lung adenocarcinoma manifested as pulmonary subsolid nodules by single-cell RNA sequencing. Sci Adv 2021;7:eabd9738. [Crossref] [PubMed]
17. Cheng S, Li Z, Gao R, et al. A pan-cancer single-cell transcriptional atlas of tumor infiltrating myeloid cells. Cell 2021;184:792-809.e23. [Crossref] [PubMed]
18. Taube JH, Herschkowitz JI, Komurov K, et al. Core epithelial-to-mesenchymal transition interactome gene-expression signature is associated with claudin-low and metaplastic breast cancer subtypes. Proc Natl Acad Sci U S A 2010;107:15449-54. [Crossref] [PubMed]
19. Hänzelmann S, Castelo R, Guinney J. GSVA: gene set variation analysis for microarray and RNA-seq data. BMC Bioinformatics 2013;14:7. [Crossref] [PubMed]
20. Patel AP, Tirosh I, Trombetta JJ, et al. Single-cell RNA-seq highlights intratumoral heterogeneity in primary glioblastoma. Science 2014;344:1396-401. [Crossref] [PubMed]
21. Gayoso A, Weiler P, Lotfollahi M, et al. Deep generative modeling of transcriptional dynamics for RNA velocity analysis in single cells. Nat Methods 2024;21:50-9. [Crossref] [PubMed]
22. Aibar S, González-Blas CB, Moerman T, et al. SCENIC: single-cell regulatory network inference and clustering. Nat Methods 2017;14:1083-6. [Crossref] [PubMed]
23. Browaeys R, Saelens W, Saeys Y. NicheNet: modeling intercellular communication by linking ligands to target genes. Nat Methods 2020;17:159-62. [Crossref] [PubMed]
24. Huang Y, Hong W, Wei X. The molecular mechanisms and therapeutic strategies of EMT in tumor progression and metastasis. J Hematol Oncol 2022;15:129. [Crossref] [PubMed]
25. Xia T, Zhang M, Lei W, et al. Advances in the role of STAT3 in macrophage polarization. Front Immunol 2023;14:1160719. [Crossref] [PubMed]
26. Zhou J, Lyu N, Wang Q, et al. A novel role of TGFBI in macrophage polarization and macrophage-induced pancreatic cancer growth and therapeutic resistance. Cancer Lett 2023;578:216457. [Crossref] [PubMed]
27. Nieto MA, Huang RY, Jackson RA, et al. EMT: 2016. Cell 2016;166:21-45. [Crossref] [PubMed]
28. Aiello NM, Maddipati R, Norgard RJ, et al. EMT Subtype Influences Epithelial Plasticity and Mode of Cell Migration. Dev Cell 2018;45:681-695.e4. [Crossref] [PubMed]
29. Ku HC, Cheng CF. Master Regulator Activating Transcription Factor 3 (ATF3) in Metabolic Homeostasis and Cancer. Front Endocrinol (Lausanne) 2020;11:556. [Crossref] [PubMed]
30. Qian X, Zhu L, Xu M, et al. Shikonin suppresses small cell lung cancer growth via inducing ATF3-mediated ferroptosis to promote ROS accumulation. Chem Biol Interact 2023;382:110588. [Crossref] [PubMed]
31. Thompson MR, Xu D, Williams BR. ATF3 transcription factor and its emerging roles in immunity and cancer. J Mol Med (Berl) 2009;87:1053-60. [Crossref] [PubMed]
32. Liu G, Su L, Hao X, et al. Salermide up-regulates death receptor 5 expression through the ATF4-ATF3-CHOP axis and leads to apoptosis in human cancer cells. J Cell Mol Med 2012;16:1618-28. [Crossref] [PubMed]
33. Li X, Zhou X, Li Y, et al. Activating transcription factor 3 promotes malignance of lung cancer cells in vitro. Thorac Cancer 2017;8:181-91. [Crossref] [PubMed]
34. Jan YH, Tsai HY, Yang CJ, et al. Adenylate kinase-4 is a marker of poor clinical outcomes that promotes metastasis of lung cancer by downregulating the transcription factor ATF3. Cancer Res 2012;72:5119-29. [Crossref] [PubMed]
35. Kim N, Kim HK, Lee K, et al. Single-cell RNA sequencing demonstrates the molecular and cellular reprogramming of metastatic lung adenocarcinoma. Nat Commun 2020;11:2285. [Crossref] [PubMed]
36. Farhad M, Rolig AS, Redmond WL. The role of Galectin-3 in modulating tumor growth and immunosuppression within the tumor microenvironment. Oncoimmunology 2018;7:e1434467. [Crossref] [PubMed]
37. Shanshiashvili L, Tsitsilashvili E, Dabrundashvili N, et al. Metabotropic glutamate receptor 5 may be involved in macrophage plasticity. Biol Res 2017;50:4. [Crossref] [PubMed]
38. Shetty P, Bargale A, Patil BR, et al. Cell surface interaction of annexin A2 and galectin-3 modulates epidermal growth factor receptor signaling in Her-2 negative breast cancer cells. Mol Cell Biochem 2016;411:221-33. [Crossref] [PubMed]
39. Wang T, Wang Z, Niu R, et al. Crucial role of Anxa2 in cancer progression: highlights on its novel regulatory mechanism. Cancer Biol Med 2019;16:671-87. [Crossref] [PubMed]