1. Introduction

Lung adenocarcinoma (LUAD), the most common histological subtype of non-small cell lung cancer (NSCLC), remains a leading cause of cancer-related mortality worldwide [1]. The genetic landscape of this disease is complex, involving multiple driver mutations such as EGFR, KRAS, and genomic instability alterations, which collectively contribute to tumor progression and metastasis [24]. The tumor microenvironment (TME) plays a critical role in shaping tumor progression, immune evasion, and therapeutic response, particularly in the context of genetic mutations affecting DNA repair mechanisms [5, 6]. Homologous recombination repair (HRR) deficiency, characterized by mutations in key DNA repair genes, has been increasingly recognized as a determinant of tumor genomic instability in breast cancer, ovarian cancers, and lung cancer [79].

BReast CAncer gene 1 (BRCA1) and gene 2 (BRCA2) are crucial mediators of genomic stability, and their mutations often lead to increased DNA damage, heightened immune surveillance, and altered tumor-immune interactions [10, 11]. Previous studies revealed that mutations in Brca1 and Brca2 mediate distinct immune microenvironment effects in murine models of breast cancer, and showed that Brca2 deficiency improved response to immune checkpoint inhibitors [10, 12]. Among the key genetic mutations that drive LUAD, BRCA1 and BRCA2 mutations have been increasingly implicated in the pathogenesis and progression of lung cancer [9, 13]. Recent studies suggest that HRR-deficient tumors exhibit distinct immunogenic profiles, influencing immune cell infiltration, antigen presentation, and responses to immune checkpoint blockade (ICB) therapy in NSCLC [14, 15]. BRAC1 mRNA can predict the efficacy of second-line cisplatin chemotherapy in patients with metastatic NSCLC and has a good predictive value for patients with advanced NSCLC receiving platinum-based chemotherapy [16, 17]. One report showed a case of BRCA2-positive lung adenocarcinoma, in which the disease was stable for about 2 years after treatment with olaparib [18]. Despite this, the molecular architecture of the TME in LUAD with BRCA1 and BRCA2 mutations remains incompletely understood.

Understanding the molecular mechanisms through which BRCA1 and BRCA2 mutations alter the TME in LUAD is crucial for identifying potential therapeutic targets and improving treatment strategies. To address these challenges, we leveraged single-cell RNA sequencing (scRNA-seq) and single-cell T cell receptor sequencing (scTCR-seq) to comprehensively dissect the BRCA1/2-driven TME in LUAD. scRNA-seq enables transcriptomic profiling at the single-cell level, providing unparalleled insights into cellular heterogeneity, lineage differentiation, and gene expression dynamics within tumors [19, 20]. Unlike bulk RNA sequencing, which averages gene expression across heterogeneous cell populations, scRNA-seq allows for the identification of distinct malignant and immune cell subsets, as well as their functional states [21]. Furthermore, scTCR-seq offers a powerful approach for understanding T cell diversity, clonal expansion, and antigen-specific immune responses [22].

In this study, we leveraged large-scale genomic and transcriptomic datasets to investigate the prognostic implications of BRCA1/2 mutations in LUAD patients (Nearly 2,000 samples). By integrating scRNA-seq and scTCR-seq, this study provides a high-resolution characterization of the immune microenvironment and antigen-specific T cell responses in BRCA1/2-mutant LUAD. Our findings reveal how BRCA1/2 mutations shape the immune landscape, influence CD8+ and CD4+ T cell populations, and modulate immunotherapy outcomes. These insights offer comprehensive view of the molecular mechanisms underlying BRCA1 and BRCA2 mutations in LUAD and highlights new opportunities for personalized therapeutic strategies.

2. Materials and methods

2.1. Human specimens

Patient-derived samples for this study were collected at Affiliated Cancer Hospital & Institute of Guangzhou Medical University, Guangdong, China. All patients were treatment-naïve and had invasive lung adenocarcinoma and PD-L1 negative confirmed based on immunohistochemistry. Fresh tumor and adjacent normal tissues of patients were collected through surgery. This study was approved by the Ethics Committee of Affiliated Cancer Hospital and Institute of Guangzhou Medical University (Approval No. KY-2025003) in compliance with the international guidelines. All patients have signed informed consent. The samples were divided into two parts, one for single-cell sequencing and the other for whole exome sequencing. In addition, we collected and obtained public datasets of multiple cohorts. The genomic variation, gene expression, and clinical data of LUAD patients (n=564) of the cancer genome atlas (TCGA) program were downloaded from cBioPortal data resource (https://www.cbioportal.org/). The genomic variation and survival data of LUAD in individuals of East Asian ancestry (OncoSG cohort, n=305) [23] also from cBioPortal. The mutation and survival data of LUAD patients receiving ICB treatment from SU2C MARK cohort (n=309) [24]. The gene expression, survival, and response data of NSCLC patients receiving ICB treatment from OAK cohort (n=699) [25]. Tumor mutation burden (TMB), homologous recombination deficiency (HRD) score, and neoantigen load data for TCGA-LUAD patients were obtained from a previous study [26]. The HRD score was determined bysumming specific genomic alterations, including loss of heterozygosity (LOH), large-scale state transitions (LST), and telomeric allelic imbalances (ntAI).

2.2. DNA extraction and library construction of WES

Genomic DNA were extracted from lung adenocarcinoma tissues and blood using CretMagTM Multi Sample DNA Kit. Genomic DNA were captured using Agilent SureSelect Human All Exon v6 library following the manufacturer’s protocol (Agilent Technologies, USA). DNA libraries were constructed following the protocols provided by Illumina. Then these libraries were sequenced on the Illumina sequencing platform (Illumina NovaSeq X plus) and 150 bp paired-end reads were generated. The whole exome sequencing was conducted by OE Biotech Co., Ltd. (Shanghai, China).

2.3. Somatic variant identification

The raw reads (fastq) were pre-processed and trimmed with fastp (Version: 0.23.4) based on default parameters. Clean reads were aligned to the reference human genome (GRCh38) utilizing the BWA (version 0.7.17). The mapped reads were sorted and indexed by using SAMtools (Version 1.5). The GATK (Version 4.0.5.1) was used for recalibration of the base quality score and for single nucleotide polymorphism (SNP) and insertion/deletion (INDEL) realignment, and used for marking duplicate reads, to obtain analysis-ready BAM files. The final BAM files were used as input files for variant calling. Identification of somatic mutations using the GATK Mutect2 method based on matched blood sample and population germline resource from gnomAD (https://gnomad.broadinstitute.org/). ANNOVAR was used to annotated variations during SNP&INDEL calling.

2.4. Single-cell suspension preparation

The solid tissues (including tumor and adjacent tissues) were processed to prepare single-cell suspensions. The tissue samples were cut into small pieces, and the cells were washed twice with precooled RPMI 1640 + 0.04% BSA medium under aseptic conditions. The tissues were cut into small pieces of approximately 0.5 mm3 using surgical scissors, placed in freshly prepared hydrolase and digested at a constant temperature of 37°C for 60 minutes, mixing every 10 minutes. The tissue was filtered 2 times using a cell sieve and centrifuged (300g) for 5 minutes at 4°C. After resuspending the sediment with an appropriate amount of medium, an equal volume of erythrocyte lysis buffer (MACS, catalog number: 130-094-183) was added. After mixing, the cells were incubated at 4°C for 10 minutes and centrifuged for 5 minutes (300g), and the number of viable cells, as determined by Trypan blue staining, was determined with a hemocytometer. Human fresh peripheral blood was collected in an anticoagulant tube and diluted with PBS at 1:1. Then, the sample was added to the Ficoll separation medium (Merck) and centrifugated 20 minutes (75g). After centrifugation, a white-grayish layer consisting of mononuclear cells will be found on top of the Ficoll-Paque. We carefully collected the mononuclear cells with a Pasteur pipette, transferred them to a new centrifuge tube, and washed them by centrifuging for 10 minutes (300g). The supernatant was discarded and resuspended in RPMI-1640 supplemented with 0.04% BSA. The number of viable cells, as determined by Trypan blue staining, was determined with a hemocytometer.

2.5. Single-cell RNA and V(D)J sequencing

The cell concentration of the freshly prepared single cell suspension was adjusted to 700-1200 cell/μl. The cDNA library amplification was performed according to the operating instructions of 10× Genomics Chromium Next GEM Single Cell 5ʹ Reagent Kits v2.0 (PN-1000263). DNA library construction was performed using the 10x Genomics Chromium Single Cell 5’ Kit (PN-1000190) following the manufacturer’s guidelines. Purified cDNA libraries were sequenced using the paired-end 150 mode on an Illumina platform (Illumina NovaSeq X plus). For T-cell receptor (TCR), T cell V(D)J enrichment was performed using the 10× Genomics single cell V(D)J Enrichment Kit (TCR, PN-1000005) according to the human T cell operating instructions. The TCR library was amplified using the Chromium Single Cell TCR amplification Kit (TCR, PN-1000252), and the experimental operation was performed according to the product instructions. The constructed libraries were performed high-throughput sequencing (Illumina NovaSeq X plus) using the paired-end 150 mode at OE Biotech Co. Ltd (Shanghai, China).

2.6. scRNA-seq and scTCR-seq data processing

Raw read files were processed with Cell Ranger 7.1.0 and mapped to the GRCh38 genome assembly and counted with GRCh38 annotations. Unique molecular identifier (UMIs) was calculated by using the “cellranger count” function. As a result, a cellular gene expression matrix containing the number of gene UMIs detected in each cell was generated. In addition, we took several steps to obtain high-quality data. First, cells with high mitochondrial gene expression were removed, as dead cells often exhibit mitochondrial contamination [27]. Specifically, we followed the previous method and used the median-centered median absolute deviation (MAD)-variance normal distribution to fit the expression levels of mitochondrial genes, and removed cells with expression levels significantly higher than expected (Benjamini-Hochberg-corrected FDR < 0.01) [28]. Second, we removed cells with less than 500 or more than 6000 detected genes. Third, we used DoubletFinder [29] and scDblFinder [30] methods to identify and remove potential doublets. The 90th percentile of the doublet score of DoubletFinder was used as the cutoff value. And the expected doublet rate is 10% in scDblFinder method. Raw TCR reads were aligned to the GRCh38 reference genome and consensus TCR annotation was performed using “cellranger vdj” function (Cell Ranger 7.1.0). TCR annotation was performed using the default parameters of the 10X cellranger vdj pipeline. The R package scRepertoire was employed for clonotype assignment and dynamics analysis.

2.7. Normalization, clustering and visualization of scRNA-seq data

The gene expression matrix of the selected high-quality cells was subjected to downstream analysis using the R package Seurat (v5.0.3) [31]. The main steps include the following. First, only genes expressed in more than 5 cells were retained. Second, the raw UMI count matrix was log-normalized and scaled to 10,000 using the “LogNormalize” function. Third, 2,000 highly variable genes were identified using the “FindVariableFeatures” function based on the normalized matrix. Fourth, the gene expression matrix was scaled and centered using the “ScaleData” function. Fifth, unsupervised clustering was performed by constructing a shared nearest neighbor (SNN) graph using the “FindNeighbors” function and the Louvain algorithm. The first 30 principal components were considered and the clustering resolution was 1.6. For uniform manifold approximation and projection (UMAP) visualization, the first two dimensions of UMAP was calculated using the “RunUMAP” function. Sixth, differential expression analysis between clusters was performed using the “FindAllMarkers” function from Seurat (using default parameters but setting min.pct to 0.5). Use the “FindMarkers” function to find differential expression genes for specific clusters. Finally, two complementary approaches are used to annotate cell clusters: (1) canonical markers; (2) markers were in the top rank of differential expressed genes.

2.8. Signature score calculation at single-cell and bulk level

Signature (module) score at single-cell level was calculated for selected gene-sets (including pathways and signatures) using “AddModuleScore_UCell” function from UCell R package. The pathways or signatures involved in this study include: immune pathways and receptors from ImmPort database (https://immport.org/shared/home), tissue-specific T cell signatures [32] and known cancer meta-programs [33] from previous studies, as were as several specific signatures from the molecular signatures database (MSigDB) human collections. In addition, tumor tissue upregulated and downregulated gene signatures are the top 60 up-regulated genes and the bottom 60 down-regulated genes, respectively, calculated using the TCGA-LUAD cohort. CD8+ Trm and CD4+ Trm signatures were identified based on “FindMarkers” function at single-cell level. Signature score at bulk level was calculated using single sample gene set enrichment analysis (ssGSEA) algorithm in R package GSVA.

2.9. Functional enrichment analysis

Differential expression analysis between malignant cells and normal cells was performed, and the up-regulated genes of malignant cells were screened according to the thresholds avg_log2FC≥0.5 and Padj≤0.01. According to the descending sorting of avg_log2FC, the top 200 genes of malignant cells were selected. The GO BP functional enrichment analysis was performed using the R package clusterProfiler (the threshold used were P≤0.001 and Padj≤0.01). Padj calculations were based on Benjamini-Hochberg correction. The same approach was used for the comparison between the two types of malignant cells and the comparison between malignant cells in BRCA1 and BRCA2 mutation samples. In addition, we obtained immune signatures and gene modules of immune cell type from MSigDB hallmarks, ImmPort database, and previous studies [34, 35]. Gene set enrichment analysis (GESA) was employed to compare the immune response profile between BRCA1 and BRCA2 mutations using the R package enrichplot.

For all malignant cells, we performed pathway enrichment analysis using ActivePathways (v.2.0.3) R package for the terms from MSigDB (human collection C5) based on the gene sets of two meta-programs. The terms were removed once the number of genes was less than 10 or more than 500. Terms with BH-corrected significance adjusted P≤0.05 are considered to be significantly enriched by the gene sets of two meta-programs and will be retained. The significantly enrichment results were used as input to the EnrichmentMap plugin in cytoscape (v.3.9.1) software to draw a network diagram.

2.10. Single-cell copy number variation inference

The inferCNV (v.1.19.1) R package was used to distinguish malignant cells by inferring chromosomal CNVs based on the single-cell expression data. The cells from CD8+ T cells as normal reference were used to estimate CNVs for the population of malignant cells. We prepared gene order files containing the chromosomal start and end positions of each gene from the GRCh38 assembly as input to the “gene_order_file” parameter in the “CreateInfercnvObject” function. To perform the infercnv operation, the count matrix and annotation file are input to create the infercnv object, and then inferCNV is run with cutoff=0.1 in the “run” function.

2.11. Identification of potential drugs targeting BRCA1-mutant genes

This study identified potential drug small molecules based on BRCA1 mutation-related cancer-promoting factors. Gene expression profiles of lung adenocarcinoma cell lines before and after small molecule drug treatment were obtained from the NIH library of integrated network-based cellular signatures (LINCS) data resource. The connectivity Map (cMap) analysis was employed to identify potential drugs associated with molecular perturbations of four BRCA1 mutation prognostic risk factors (4-R genes) [36]. To demonstrate the inhibitory effect of HDAC inhibitors on the growth of lung adenocarcinoma cell lines, the replication levels at different HDAC inhibitor concentrations were obtained from the cancer cell line encyclopedia (CCLE) data resource. However, due to limited data, we only obtained data on two HDAC inhibitors, entinostat and vorinostat, in 43 lung adenocarcinoma cell lines. Since drug concentrations (C) were divergent, we categorized them into 4 groups: 0 µM<C≤2 µM ([0–2]), 2 µM <C≤6 µM ([26]), 6 µM <C≤10 µM ([610]) and vehicle control (C=0 µM). The mean of each group is the average replication level of the cell line.

2.12. Cell culture and reagents

The human lung cancer cell line A549 (#SCSP-503) and HEK293T (#SCSP-502) were purchased from the Type Culture Collection of the Chinese Academy of Sciences, China. Short Tandem Repeat (STR) analyses were performed to authenticate the identity of each Cell line used in this article. The A549 cells were cultured in RPMI-1640 medium (#C11875500BT; Gibco). The HEK293T cells were cultured in Dulbecco’s Modified Eagle’s medium (#C11995500BT; Gibco). All two cell lines were supplemented with 10% fetal bovine serum (#S712-012S; Lonsera), and grown in a humidified 5% CO2 atmosphere at 37 °C.

2.13. Vectors and lentiviral transfection

All the short hairpin RNAs (shRNA) were cloned into pLKO.1 vector. Target sequences are as follows:

LDHA#1: GCCACAGATTTACCCGTGGAT,

LDHA#2: GCCAACAACTTGTGTCTCAAT,

S100A10#1: TAAGGAGCCAAATACCTTGCG,

S100A10#2: CATGAAACACAAACGGCAAAT,

GAPDH#1: GTGCGGAGTGTAATCAGTATT,

GAPDH#2: TCAGGTTGTACGGGATCAAAT.

Lentiviral particles were produced by transfecting the shRNA plasmid, packaging vectors psPAX2 and pMD2.G in HEK293T cells at a ratio of 4:3:1 by JetPRIME In Vitro DNA Transfection Reagent (#101000046 (114-15); Polyplus). The media was changed after 6 h. Then 48 h after media change, the media containing the virus were collected and concentrated. Viral supernatants were centrifuged at 1500 x g for 45 min and viral pellets were resuspended with DMEM medium. The lentivirus was frozen at –80 °C for further use. The efficiency of lentiviral shRNA clones was determined by western blotting.

2.14. Lentiviral infection

1 × 105 A549 cells were seeded into each well of 6well plates, at the same time, added with 400 μL of lentivirus and 8 μg/mL polybrene (#28728-55-4; Selleck). Plates were incubated at 37 °C for 12 h. The cells were then passaged and seeded into new culture dishes. 24 hours after cell passage, infected lentiviral cells were selected with Puromycin (1 μg/mL, #P8833-10mg; Sigma-Aldirch) for 3 days.

2.15. Cell proliferation assay

Cell proliferation was measured over a 4-day time course using the CellTiter-Glo® Luminescent Cell Viability Assay kit (#G9242; Promega). For this assay, 1000 cells were seeded into each well of a 96-well plate containing 100 μL of medium per well and cultured in normal conditions. At the time point of detection, plates were cooled down to room temperature, then 100 μL of the CellTiter-Glo reagent was added to each well. Mix the contents for 2 minutes on an orbital shaker to induce cell lysis. Allow the plate to incubate at room temperature for 10 minutes to stabilize the luminescent signal. Finally, plates were read using a microplate reader (#SLXFA-SN; Gene Company). All data were normalized to day 1 and presented as mean ± SD.

2.16. Drug treatment

For Cell viability assay, 5000 cells of A549 were seeded into 96-well plates containing 100 μL of medium per well. After 24 h of culture at 5% CO2 at 37 °C, cells were treated with a gradient concentration of corresponding agents (Vorinostat, S1047; Belinostat, S1085; Entinostat, S1053; Pracinostat, S1515; Selleckchem) for 48 h. Cell viability was determined using the CellTiter-Glo® Luminescent Cell Viability Assay kit (#G9242; Promega) as described above. For detecting the protein level of the targeted gene induced by drug treatment, 4 × 105 A549 cells were seeded into each 60 mm cell culture dish. After 24 h of culture at 5% CO2 at 37 °C, cells were treated with 1 and 3 μM Vorinostat or Belinostat for 2 days. Then cells were harvested for detecting the protein level of the targeted genes.

2.17. Western blotting

In the case of the efficiency of lentiviral shRNA clones and drug treatment studies. Cells were washed in PBS twice and collected by trypsinization and centrifugation (1500 x g for 5 min). Cell pellets were washed in PBS twice, resuspended in RIPA (#P0013B; Beyotime), supplemented with protease inhibitor cocktail (#78440, Invitrogen), incubated on ice and boiled for 10 min at 100°C. Protein lysates were resolved on a 4-12% YoungPAGE™ Precast Gels (#PK0931, GenScript) and transferred to 0.45 μm Immobilon-P PVDF membrane (#IPVH00010, Merck millipore). Membranes were incubated for 1 hr in TBST (#BL602A, Biosharp) containing 5% non-fat dry milk and then incubated for 1 hr at RT or overnight at 4°C with primary antibodies. Detection was achieved using appropriate horseradish peroxidase (HRP)-conjugated secondary antibodies. Finally, the membrane was incubated with an enhanced chemiluminescence ECL (#CW0049M; CWBIO), and the images were obtained by using MiniChemi®580 + Imaging System (MiniChemi®580; SINSAGE). The following antibodies were used: GAPDH Rabbit mAb (#A19056, 1:100,000; Abclonal), S100A10 Rabbit mAb (#A13634, 1:1000; Abclonal), LDHA Rabbit mAb (#A0861, 1:1000; Abclonal), β-Tubulin Rabbit mAb (#A12289, 1:10000; Abclonal) and HRP Goat anti-Rabbit IgG(H+L) (#111-035-003, 1:10000; Jackson).

2.18. RNA isolation and RT-qPCR

Total cellular RNA was isolated using TRIzol Reagent (#15596018; Thermo Fisher Scientific). The RNA concentration was measured using Nanodrop 2000 (Thermo Fisher Scientific, USA). Total RNA was reverse transcribed into cDNA using HiScript III 1st Strand cDNA Synthesis Kit (+gDNA wiper) (#R312-02; Vazyme) according to the manufacturer’s instructions. Quantitative reverse transcription PCR was performed per the manufacturer’s protocol on the SLAN-96P Real-time PCR system using ChamQ SYBR qPCR Master Mix (#Q311-02; Vazyme) and gene-specific primers. The qPCR primers sequences are as follows:

LDHA-fwd: 5′-ACCGTGTTATTGGAAGCGGT-3′,

LDHA-rev: 5′-CTCCATGTTCCCCAAGGACC-3′,

GAPDH-fwd: 5′-GTCAAGGCTGAGAACGGGAA-3′,

GAPDH-rev: 5′-AAATGAGCCCCAGCCTTCTC-3′,

β-actin-fwd: 5′-CACCATTGGCAATGAGCGGTTC-3′,

β-actin-rev: 5′-AGGTCTTTGCGGATGTCCACGT-3′.

2.19. Statistics and survival analysis

All statistical analysis were conducted using R software (version 4.2.3, http://www.r-project.org). The two-sided Wilcoxon rank-sum test was employed to count and compare the differences between the two groups. Kruskal-Wallis test was used to compare the differences among the three groups. The significance level of discrete variables was calculated by Chi-squared test. The identification of prognostic factors associated with BRCA1 and BRCA2 mutations was based on the up-regulated genes of the corresponding mutations using univariate Cox regression analysis. Survival analysis was performed by the log-rank test and compared using Kaplan-Meier curves. Statistical significance for pairwise comparisons across the multiple groups was assessed using the Benjamini-Hochberg false discovery rate (FDR) correction. For all statistical tests, *P≤0.05; **P≤0.01; ***P≤0.001; ****P≤0.0001.

3. Results

3.1. BRCA mutations are associated with genomic instability and poor prognosis in lung adenocarcinoma

Using the TCGA-LUAD cohort data, we examined the effect of somatic mutations of BRCA1 and BRCA2 on genomic instability. The results showed that the mutations in both genes could significantly increase the homologous recombination repair deficiency (HRD) score and present a higher tumor mutation burden (TMB) and fraction of genome altered (FGA; Fig. 1A and Fig. S1A), indicating elevated genomic instability in LUAD tumors with BRCA1/2 mutations. This phenomenon was also observed in BRCA1/2 mutated patients with NSCLC (Fig. S1B-C). Furthermore, BRCA1 and BRCA2 mutations increased the neoantigen load in LUADs (Fig. 1A), suggesting that BRCA1/2-mutated tumors might be more immunogenic. By analyzing the proportion of immune subtypes [26], we found that patients with BRCA1 mutations had a higher proportion of C1 subtype (wound healing; P=5.80e-04, Chi-squared test), while patients with BRCA2 mutations and wild-type had a higher proportion of C3 subtype (inflammatory; P=6.61e-03, Fig. 1B). BRCA1-mutated patients also showed activated DNA damage repair pathways and exhibited higher DNA replication and cell cycle activity compared to wild-type patients (Fig. 1C-D), suggesting that BRCA1 mutations may promote tumor growth through compensatory DNA repair activation—a pattern that may reflect partial HRD rather than complete loss of function.

BRCA1 and BRCA2 mutations are associated with genomic instability and poor prognosis of LUAD.

A. Distribution of HRD score, TMB and neoantigen load in patients with wild-type, BRCA1- and BRCA2-mutated lung adenocarcinoma from the TCGA-LUAD cohort (n=566). B. Proportion of immune subtypes. The significance level was calculated by Chi-squared test. C1: wound healing, C2: IFN-g dominant, C3: inflammatory, C4: lymphocyte depleted, C6: TGF-b dominant. C-D. Activity of DNA damage repair and cell cycle-related pathways in wild-type, BRCA1- and BRCA2-mutated patients (TCGA-LUAD cohort). E-F. Survival analysis of BRCA mutations in TCGA-LUAD (E) and OncoSG-LUAD (F) cohorts. G. Survival analysis of BRCA mutations in LUAD patients after ICB treatment. H. Survival analysis based on BRCA1 (left) and BRCA2 (right) expression in patients with ICB treatment. I. Lymphocyte infiltrates in patients with wild-type, BRCA1 and BRCA2 mutations (TCGA-LUAD cohort). J. Mutation frequencies of the oncogenic signaling pathway genes (in-house data).

Next, we explored the impact of BRCA1/2 mutations on clinical outcomes of LUAD patients and found that these mutations were associated with worse prognosis for disease-free survival (DFS, P=8.46e-03, Log-rank test; Fig. 1E) and overall survival (OS, P=0.03; Fig. 1F) in multiple LUAD cohorts. In particular, BRCA1 and BRCA2 mutations as risk factors shortened the survival of LUAD patients (P=0.028; Fig. 1E). Interestingly, LUAD patients with BRCA mutations had significantly longer progression-free survival after ICB treatment (Fig. 1G). Through transcriptional analysis, we found that upregulation of both BRCA1 and BRCA2 was associated with worse prognosis after ICB treatment (P<0.05, Fig. 1H), consistent with BRCA1/2 expression as a risk factor for lung cancers [9]. The apparent paradox between favorable ICB outcomes linked to BRCA1/2 mutations and poorer outcomes associated with high BRCA1/2 expression may be explained by distinct biological roles—where mutations induce genomic instability and immunogenicity, while high expression reflects proficient DNA repair and tumor cell fitness. These results suggest that BRCA abnormalities may have potential to predict ICB treatment outcomes. In addition, we found the inconsistencies in lymphocyte infiltration between BRCA1 and BRCA2 mutations, suggesting heterogeneity in the immune microenvironment (Fig. 1I).

3.2. Single-cell landscape of BRCA mutations in LUAD patients

Exome sequencing of this study showed that the two types of tumor tissues carried non-synonymous mutations of driver genes TP53, EGFR, and FAT1 (Table S1), and had similar mutation frequencies in oncogenic signaling pathways, including p53 and Notch pathways, as well as significant mutation gene set (MutSig) scrutinized by MutSigCV [37] (Fig. 1J). To preliminarily explore the impact of BRCA1/2 mutations on the immune microenvironment, scRNA-seq was performed on LUAD patients carrying BRCA1 or BRCA12 somatic mutations (a total of 6 samples, including 2 tumors, 2 adjacent normal tissues, and 2 peripheral blood samples from one BRCA1-mutant and one BRCA2-mutant patient). A total of 69,180 high-quality cells were obtained after QC, of which 27,851 (40.3%) were from tumor tissues, 25,035 (36.2%) were from adjacent normal tissues, and 16,294 (23.5%) were from blood (Fig. 2A and Fig. S2A-B). Of these cells, 35,450 cells (51.2%) and 33,730 cells (48.8%) were from the BRCA1- and BRCA2-mutant samples, respectively. After batch correction (Harmony), we performed cell clustering and cataloged all cells into 19 main cell types annotated with canonical marker genes (Fig. 2A-C and Fig. S2A-C), identifying 4 epithelial compartments (alveolar type 1 [AT1], alveolar type 2 [AT2], ciliated, and club cells), 2 stromal compartments (fibroblasts and endothelial cells), 5 lymphocyte compartments (CD4+ T, CD8+ T, natural killer-NK/NKT, B cells, and plasma) and 6 myeloid compartments (macrophages, dendritic cells-DCs, classical/nonclassical monocytes, neutrophil, and mast cells). Among immune cells, the most abundant cell types were T lymphocytes and macrophages (Fig. 2B).

Single-cell transcriptome analysis and LUAD malignant cells identification.

A. UMAP visualization of cell types in patients with BRCA1/2 mutations (in-house data). B. The number of each cell type and its proportion in different samples, including tumor tissue (T), adjacent normal tissue (N), and peripheral blood (B). C. The bubble plot shows the expression of the markers used for cell type identification. D. The heat map inferred CNVs in epithelial cells, malignant cells, and reference cells. E. Distribution of CNV scores across different cell types on specific chromosomes. F. UMAP visualization of CNV score in epithelial cells and malignant cells. G. UMAP visualization (left) and box plot (right) of tumor upregulated signature score. H. Expression of CEACAM6 in epithelial and malignant cells.

To distinguish malignant cells, we first inferred copy number variation (CNV) of all epithelial cell types in solid tissues using the inferCNV pipeline [38] (Fig. S2D-E). Using immune cells as a control and comparing the CNV scores of each epithelial cell subset, we identified two malignant cell subsets that were significantly clustered together using the K-means clustering algorithm and exhibited higher scores of oncogenic CNV in LUAD (Fig. 2D-E). For example, both chromosome 1q21 region and chromosome 7—frequently amplified in LUAD tumors [39, 40]—showed significantly higher CNV scores in the malignant cells of this study (Fig. 2D-E). The previously reported loss of chromosome 3q [41] was also observed in malignant cells (Fig. 2D-E). To demonstrate the properties of malignant cells, we used bulk RNA-seq data from the TCGA-LUAD cohort to identify tumor upregulated and downregulated gene signatures, and calculated module scores in different epithelial cell types. The malignant cell subsets had significantly higher tumor upregulated signature scores and lower downregulated signature scores (Fig. 2F-G and Fig. S2F). In addition, tumor marker genes, such as CEACAM6, CLDN4, and WFDC2, showed higher expression in malignant cells (Fig. 2H and Fig. S2C).

3.3. Distinct malignant transcriptional programs associated with BRCA1 and BRCA2 mutations in LUAD

We have identified two subsets of malignant cells (13,694 cells in total) from the BRCA1/2-mutant samples. We then asked which key molecules were specifically activated in malignant cell subsets and whether there were differences in the oncogenic programs between the BRCA1- and BRCA2-mutant cases. To address this, we first performed differential expression analysis comparing malignant cell subsets versus other epithelial cells. The results showed that tumor markers WFDC2 and CLDN4 were significantly upregulated in the malignant1 subset, while CEACAM6 was significantly upregulated in both malignant cell types (Fig. S3A and Table S2). Previous studies show that DUSP1 could promote angiogenesis, invasion, and metastasis in NSCLC, and FTH1 expression protects cells from cell death induced by the GPX4 inhibitor RSL3 [42, 43]. These genes upregulated in malignant cells are often also highly expressed in tumor tissues relative to paracancerous tissue (Fig. S3B). The upregulated genes in malignant1 subset were mainly enriched in ribonucleotide biosynthetic and nucleoside-related metabolic process (Fig. 3A and Table S3), suggesting enhanced DNA replication and growth abilities. Differently, the upregulated genes of malignant2 were mainly enriched in cell adhesion, antigen processing and presentation of peptide antigen, and MHC protein complex assembly (Fig. 3A and Table S3), suggesting involvement in lymphocyte recruitment. These patterns were also observed in direct comparison between the two malignant subsets (Fig. S4A and Table S4-5).

Identification of key programs in BRCA1/2 mutated tumors.

A. Functional enrichment analysis of genes that upregulated in malignant cells relative to normal epithelial cells. B. The heat map of the similarity between pairs of programs is identified based on the NMF algorithm. C. MP score between BRCA1 and BRCA2 mutated malignant cells. D. K-M curves for MP1 score survival analysis. E. Enrichment map of the functions enriched by MP1 and MP2 genes. F. The overlap of in-house MPs and known cancer MPs. G. The heat map shows the activity of known cancer MPs in BRCA1 and BRCA2 mutated malignant cells. H. Volcano diagram of differentially expressed genes between BRCA1 and BRCA2 mutated malignant cells.

We then applied nonnegative matrix factorization (NMF) algorithm to identify transcriptional programs consisting of co-expressed genes across all malignant cells. Using this method, we extracted a total of 4 programs and identified two meta-programs (MP1, MP2) based on program similarities (Fig. 3B). The MP1 module contains 156 genes (Table S6), whose activity was significantly higher in BRCA2-mutant tumor cells (Fig. 3C) and was associated with better prognosis in LUAD patients (P=0.031, Fig. 3D). The favorable prognosis linked to MP1 was related to its involvement in lymphocyte migration, T-cell chemotaxis, and myeloid leukocyte differentiation (Fig. 3E and Table S7). The MP2 module contains 200 genes involved in functions such as ATP synthesis/metabolic process and mitochondrial translation, and its activity was higher in BRCA1-mutant tumors (Fig. 3C, 3E, Table S6-7). To further contextualize MP1 and MP2, we obtained 41 known cancer MPs defined by Gavish et al.[33]. Among these, 15 MPs irrelevant to LUAD were excluded and 26 MPs were included. Comparing our MPs with known MPs, we found that MP1 was mainly enriched in alveolar, stress and secretion-related MPs, while MP2 overlapped with DNA damage cell cycle (HMG-rich) and MYC-related MPs (Fig. 3F). Furthermore, BRCA1 mutations were associated with hypoxia and oncogenic MYC-related MPs, while BRCA2 mutations were linked to alveolar, stress and secretion-related MPs (Fig. 3G and Fig. S4B). Comparing key molecules, BRCA1 mutation was associated with upregulation of type I IFN genes (e.g., IFI27, IFITM3), angiogenesis genes (S100A4, S100A10) and LDHA (Fig. 3H, Fig. S4C, and Table S8). BRCA2 mutation was associated with upregulation of stress-related genes (e.g., GADD45B, DUSP1, and SOD2), immune cell-related genes (NFKBIA, ICAM1), inflammatory-related genes (TNFAIP3, CEBPB, and CXCL2), and apoptosis gene CFLAR (Fig. 3H, Fig. S4C, and Table S8).

3.4. Inferred transcriptional trajectories from epithelial to malignant cells in BRCA1- and BRCA2-mutant tissues

Tumor evolution is a process of gradual accumulation of somatic mutations from normal cells [44]. To infer potential transcriptional transitions during malignant transformation, we analyzed the transcriptional trajectories using Monocle 2 [45], revealing a dynamic transitional spectrum from AT1/2 cells to malignant cells (Fig. 4A and Fig. S5A). The early pseudotime was dominated by AT1 and AT2 cells, primarily from tumor-adjacent and BRCA2-mutant tumor tissues (Fig. 4A). The trajectory then extended into two branches: one dominated by malignant1 subset (87.6% from BRCA1-mutant tumor), and the other almost all malignant2 subset, allowing us to explore expression gradients associated with BRCA1 and BRCA2 mutations (Fig. 4A). Branch fate determination gene analysis showed that the malignant1 branch upregulated genes related to oxidative phosphorylation, ATP synthesis and metabolism, and ROS response (Fig. S5B). The malignant2 branch upregulated immune-related genes such as IL1B, TNFAIP3, and CD74, enriched in cytokine-mediated signaling, antigen processing and presentation, and T cell activation and differentiation (Fig. S5B). This suggests heterogeneous transcriptional fates during the transition from normal lung epithelial cells to malignant cells in these two cases. We further investigated the expression differences along the inferred trajectories. The BRCA1 mutation branch showed activation of innate immune response-related pathways such as type I IFN production, response to IFN-α/β, and upregulated genes related to T cell homeostasis/migration and T cell receptor signaling (Fig. 4B), suggesting an association with innate immune response and T cell activation. The BRCA2 mutation branch showed activation of cell growth and development, response to wounding, and wound healing-related functions, along with promotion of MHC protein complex assembly and lymphocyte-mediated immunity (Fig. 4B), suggesting that BRCA2 mutation may be associated with immune engagement alongside tumor growth.

Tumor evolutionary analysis of BRCA1 and BRCA2 mutations.

A. Pseudotime analysis of epithelial and malignant cells. The branched trajectory was colored by pseudotime (left) and cell types (right). The pie chart shows the proportion of sample groups in each branch. B. The top 50 cell fate genes and their enriched GO BP terms of BRCA1/2 mutations in malignant2 cell type. C-E. Pseudotime was broken down into 10 bins to smooth gene expression patterns. Average gene module score between BRCA1 and BRCA2 mutation groups for tumor down-regulation (C), MP1 cell cycle G2/M (D), and MP5 stress (E) gene modules. F, J. Cell density distribution based on cell cycle HMG rich versus stress gene modules (F) and MHC class I versus MHC class II gene modules (J) in the same cells between BRCA1 and BRCA2 mutations. Color represents the density of cells. The dotted line represents the median value of the corresponding module score. G. Average gene module score for IFN-α response (top) and MHC class II (bottom) gene signatures. H-I. The activities of representative signatures between BRCA1 and BRCA2 mutations in bin 10 of malignant cells according to pseudotime.

Based on pseudotime order, divided trajectories into 10 bins and analyze the activity changes of related features. The module scores of tumor downregulated genes gradually decreased along pseudotime (Fig. 4C), while tumor upregulated genes increased (Fig. S5C). Cancer MPs, such as cell cycle G2/M (MP1) and stress (MP5) increased along pseudotime in both BRCA1- and BRCA2-mutant tissues (Fig. 4D-E). Density analysis suggested that stress activation was more prominent in BRCA2-mutant tumor cells, while cell cycle HMG-rich program was mainly activated in the BRCA1-mutant group (Fig. 4F). Differences in expression programs between BRCA1 and BRCA2 mutations were mainly apparent in mid-to-late pseudotime. For example, BRCA1 mutation-associated fate-determination functions such as response to IFN-α and type I IFN showed higher activity in mid-to-late pseudotime bins (Fig. 4G-H and Fig. S5D). This pattern was also observed in cGAS-STING innate immune sensing pathways, such as cGAS and IFN-γ response (Fig. 4H and Fig. S5E). In contrast, TGF-β signaling and stress modules were more active in late pseudotime in BRCA2-mutant cells (Fig. 4I). MHC class I and II molecules were more active in late pseudotime in BRCA1- and BRCA2-mutant cells, respectively (Fig. 4G-I), a pattern also reflected in cell density analysis (Fig. 4J). CD8+ Tcm and Th1 signatures showed higher activity in late pseudotime in BRCA1- and BRCA2-mutant cells, respectively (Fig. S5F-G), suggesting differential association with CD8+ vs. CD4+ T cell engagement.

3.5. BRCA1 and BRCA2 mutations reshape distinct immune microenvironment in LUAD tumors

We next explored how BRCA1/2 mutations reshape tumor immune microenvironment of. BRCA2 mutation-upregulated genes were enriched in MHC class II-related antigen presentation and lymphocyte-mediated immunity functions (Fig. 5A). GSEA showed that MHC class II molecules were upregulated in BRCA2 mutation group, while MHC class I molecules were upregulated in BRCA1 mutation group (Fig. 5B-C), consistent with cell density distribution (Fig. 4J). Similarly, response to IFN-α/γ and type I IFN-related genes, as well as TNF-α family members and receptors, were enriched in BRCA1-mutant tumor cells (Fig. 5B-C and Fig. S6A-B). Similar patterns were observed for DNA repair, fatty acid metabolism, and CD8+ T cell signatures (Fig. 5B, 5D, Fig. S6A). DNA double-strand breaks (DSBs) can lead to release of cytoplasmic dsDNA, which can be sensed by cGAS-STING signaling to trigger innate immune responses and type I IFN production, promoting anti-tumor immunity [4648]. Based on these observations, we hypothesized that BRCA1 mutation triggers the activation of cGAS-STING signaling and CD8+ T cell infiltration. Indeed, BRCA1 mutations were upregulated DNA damage-related checkpoints, significantly activated DSBs-related DNA repair pathways (e.g., DNA double-strand break repair, homology-directed repair, and HRR), and had higher activities of cGAS-STING signaling and STING mediated induction of host immune responses compared to BRCA2 mutations (Fig. 5E-G and Fig. S6C-F). cGAS, STING1, and downstream factors STAT1 and CCL5 were upregulated in BRCA1-mutant tumor cells (Fig. 5H). Immune activity analysis suggested that BRCA1 mutation associated with enhanced activity of CD8+ T cells and myeloid cells (Fig. 5I), consistent with GSEA results. Notably, BRCA1 mutation increased infiltration of CD8+ T cells was also confirmed in lymphocyte subsets (Fig. 5J).

BRCA1/2 mutations are associated with tumor lymphocyte activation.

A. Functional enrichment of differentially expressed genes in malignant cells between BRCA1 and BRCA2 mutation groups. B. GESA of immune response profile between BRCA1 and BRCA2 mutations. Red (blue) indicates the gene set enriched in BRCA2 (BRCA1)-mutated malignant cells. C-D. GESA plot of immune signature (C) and gene module of immune cell type (D) in malignant2 subset. NES: normalized enrichment score. E. Density ridge plot of representative pathways in two types of malignant subset with BRCA1 and BRCA2 mutations. HRR: homologous recombination repair. F. UMAP visualization of cGAS pathway score. G. Mean of pathway activities related to DNA damage repair across two types of malignant cells with BRCA1 and BRCA2 mutations. H. The expression of the representative markers. I. Immune cell activities between BRCA1 and BRCA2 mutations in bin 10 of malignant cells according to pseudotime. J. Lymphocyte activity between BRCA1 and BRCA2 mutations in lymphocyte cells.

Unlike BRCA1 mutation, upregulated genes in BRCA2-mutant tumor cells were enriched in inflammatory response, interleukins receptor, and chemokines (Fig. 5B and Fig. S6A), which was confirmed by cell activity analysis. For example, the IL2/IL17 pathway, 4-1BB pathway, and cytokine pathway were more active in BRCA2-mutant tumor cells (Fig. 5E and Fig. S6C), consistent with the higher proportion of inflammatory subtypes in BRCA2-mutant patients (Fig. 1B). In addition, BRCA2 mutation was also associated with activation and differentiation of CD4+ T cells (Fig. 5B and 5J). For example, BRCA2-upregulated genes were enriched in CD4+ Tem signature and showed higher activity of Th1/Th2 pathways, Th1 toxicity and MHC class II modules (Fig. 5B-E and Fig. S6C). Indeed, immune activity analysis also suggested that BRCA2-mutant tumor cells were more associated with activation of Th1, Th2, and Treg cells compared to BRCA1 mutations (Fig. 5I). Together, these observations suggested that BRCA1 mutations are associated with cGAS-STING-mediated host immune responses and CD8+ T cell infiltration, while BRCA2 mutations are linked to MHC II molecule-mediated antigen presentation, inflammatory responses and CD4+ T cell differentiation.

3.6. BRCA1 and BRCA2 mutations facilitate the clonal expansion of tissue-resident memory T cells

To examine the effect of BRCA1/2 mutations on T lymphocytes, we obtained T cell data from BRCA wild-type in LUAD as controls from a previous study (LME) [49] and combined with our T cell data for integrated analysis. Nine T cell subsets (67,823 cells) were identified (Fig. 6A and Fig. S7A-C). Naïve T cells were mainly distributed in blood, while Teff cells were more abundant in blood and paracancerous tissues than in tumors (Fig. 6B). Treg cells were decreased in BRCA1/2 mutant tumors compared to wild-type (Fig. S7D). We identified two tissue-resident memory T cell (Trm) subsets, CD8+ Trm and CD4+ Trm, both predominantly derived from tumor tissues (Fig. 6B). CD8+ Trm cells were more abundant in BRCA1-mutant sample, whereas CD4+ Trm cells were higher in BRCA2-mutant sample (Fig. 6B and Fig. S7D). CD8+ Trm cells highly expressed TNF family members and showed high activity for interferons and chemokines, similar to Teff and NK/NKT cells, whereas CD4+ Trm cells upregulated chemokine receptors and cytokines receptors (Fig. S7E). Trm score of both ubsets could predict the overall survival of NSCLC patients after ICB treatment (P<0.01, Fig. S7G-H) and were associated with ICB response (P<0.05, Fig. S7I), suggesting Trm cells as potential biomarkers for ICB treatment in lung cancer.

T lymphocyte infiltration and TCR clonal expansion analysis.

A-B. UMAP visualization (A) and proportion (B) of T lymphocyte subsets in combined dataset. C-D. Box plot (C) and UMAP plot (D) of representative immune pathways from ImmPort database. E-F. The activity of tissue-specific T cell signatures from Judith Wienke et al. G. Terminal exhaustion marker expression levels. H. UMAP of lymphocyte subsets in in-house dataset. I. Quantification of clonal size across sample types. Clonotypes are ranked by expansion level, including: single (1 cell), small (>1 and <5 cells), medium (>5 and <20 cells), large (>20 and <100 cells), and hyperexpanded (>100 and <500 cells). J. Expanded clonotypes distribution in different samples. K. Proportional and dynamic changes in shared clonotypes (top5 per group) between different samples within the same patient. The colors represent different clonotypes. L. The proportion of clonally expanded cells (≥5 clone size) in BRCA1/2 mutated tumor tissues for specific cell types. M. The sample and cell proportions of the 50 most abundant TCR motifs.

Consistent with malignant cell results, BRCA1 mutation was associated with upregulation of TNF family members in T lymphocytes compared to BRCA2 mutation and wild type (Fig. S7E-F and Fig. S8). Chemokine receptor and TGF-β were more active in BRCA2-mutant tumor tissues (Fig. 6C-D, Fig. S7J, and Fig. S8A). BRCA1/2 mutations were associated with enhanced antigen presentation and cytokine signaling compared to wild type (Fig. S7K and Fig. S8A). Using tissue-specific T-cell signatures [32], we calculated activity at the single-cell level (Fig. 6E and Fig. S9A), and found that BRCA1/2 mutations were associated with reduced CD28 family co-stimulation and decreased cytotoxic T lymphocytes (CTLs) and NKT cell infiltration compared to wild type (Fig. 6E-F), which explains why patients with BRCA1/2 mutations showed worse prognosis (Fig. 1E-F). BRCA1 mutations were associated with reduced TIL dysfunction and exhaustion dysfunction in tumors (Fig. 6E-F), consistent with lower expression of T cell exhaustion markers (Fig. 6G). In blood samples, BRCA1 mutations were associated with enhanced TCR and its downstream signaling, while BRCA2 mutations were linked to enhanced T and NK-mediated cytotoxicity and NK cell activation receptors (Fig. 6E).

Furthermore, we performed scTCR-seq to analyze the TCR clonal expansion. Teff, CD8+ Trm, and CD4+Trm showed the highest hyperexpanded clonal types (Fig. 6H-I and Fig. S9B-C), suggesting antigen-driven activation. BRCA1-mutant samples had a higher proportion of expanded TCR clonotypes compared to BRCA2, across sample types (Fig. 6J and Fig. S9D). In addition, BRCA1 mutations were associated with higher TCR diversity in tumor and blood, while BRCA2 mutations increased TCR richness in blood (Fig. S9E). In BRCA1-mutant patient, shared clonotypes tended to expand in paracancerous tissue and were more shared between blood and tumor tissues (Fig. 6K and Fig. S9F-G). However, in BRCA2-mutant patient, shared clonotypes tended to expand in blood and paracancerous tissue (Fig. 6K and Fig. S9F). BRCA1 mutations were associated expansion of Teff, CD8+ Trm, and Tem in solid tissues, while CD4+ Trm expansion was more prominent in BRCA2 mutations (Fig. 6L and Fig. S9H). This was especially evident among the top 50 expanded TCR clonotypes (Fig. 6M). Teff was the most expansion subset in blood, regardless of mutation type (Fig. 6M and Fig. S9I). Together, these results suggest that BRCA1 mutations are associated with increase antigen exposure and CD8+ Trm expansion, while BRCA2 mutations are linked to CD4+ Trm expansion in tumor tissue and enhanced T/NK toxicity in blood.

3.7. Targeting a BRCA1 mutation-associated transcriptional program inhibits LUAD tumor growth in vitro

Our analyses suggested that BRCA1/2 mutations are associated with distinct transcriptional programs. Then, we identified BRCA mutation-related prognostic factors using Cox regression analysis. The factors associated with BRCA2 mutation were linked to longer patient survival (HR<1, Table S9), and a BRCA2 signature score predicted better survival (P=0.0015, Fig. S10A), consistent with upregulation of MHC complex assembly genes. This score can also predicted better survival in TCGA-LUAD (P<0.0001; Fig. S10B and Fig. S11A) and OAK (P<0.0001; Fig. S10C-D) cohorts. In contrast, a BRCA1 mutation prognostic factor was associated with shorter survival (HR>1, Table S9), and its signature score predicted worse survival in TCGA-LUAD (P=0.0012; Fig. 7A and Fig. S11B) and OAK (P=0.0013; Fig. 7B and Fig. S10D, suggesting that BRCA1 mutation activates a tumor-promoting program. Based on this, we identified four BRCA1 mutation-associated risk genes, S100A10, LDHA, MYL12A, and GAPDH (4-R genes) for subsequent analysis. Although the expression of 4-R genes was heterogeneous in LUAD patients, it was associated with worse overall survival (Fig. 7C and Fig. S10E-F). Previous studies have shown that S100A10 can promote cancer metastasis by recruiting MDSC cells, and increased LDHA activity contributes to tumor immune escape [50, 51]. Knockdown of S100A10, LDHA, and GAPDH reduced LUAD cell proliferation in vitro (Fig. 7D-E).

Targeting BRCA1 mutation-related risk genes inhibits tumor growth.

A-B. Survival analysis of BRCA1-mutant signature score in TCGA-LUAD (A) and OAK (B) cohorts. C. Survival analysis of S100A10 (left) and LDHA (right) in TCGA-LUAD cohort according to the expression median. D. The effects of sh-LDHA, S100A10, GAPDH, and vector on cell proliferation were determined by a cell proliferation assay in A549 LUAD cell line. E. The effects of sh-LDHA, S100A10, and GAPDH on protein levels in A549 cell line (The original image is shown in Fig. S12A). F. The molecular perturbations of 4-R genes in A549 cell line at 10 µM concentration (24 h) from the LINCS data resource. NCS: Normalized connectivity score. G. The effects of HDAC inhibitors concentration on cell proliferation in A549 cell line. H. Cell replication levels after treatment with different concentrations of entinostat (left) and vorinostat (right) in LUAD cell lines. I. The protein levels of LDHA, S100A10, and GAPDH after treatment with vorinostat (left) and belinostat (right) in A549 cell line were measured using western blotting (The original image is shown in Fig. S12B). J. The mRNA levels of LDHA and GAPDH after treatment (3 µM) with vorinostat (left) and belinostat (right) in A549 cell line were measured using qPCR. Significance levels were determined by Student’s t-test. For all statistical tests, *P≤0.05; **P≤0.01; ***P≤0.001.

We next asked whether small molecules could reduce the expression of these risk genes and inhibit LUAD growth. Using gene expression profiles of A549 cells before and after drug treatment from LINCS data resource, we performed cMap analysis to identify drugs that perturb 4-R gene expression. Interestingly, histone deacetylase (HDAC) inhibitors consistently downregulated 4-R genes in A549 cells (Fig. 7F and Table S10). Cell growth decreased significantly with increasing HDAC inhibitor concentrations (Fig. 7 G-H, Fig. S10G-H, and Table S11). Vorinostat (FDA-approved for cutaneous T cell lymphoma) significantly inhibited 4-R gene expression (NCS<-1, P<0.01; Fig. 7F and Table S10). Furthermore, vorinostat and belinostat reduced expression of LDHA, S100A10, and GAPDH in A549 cells in a dose-dependent manner (Fig. 7I-J and Fig. S10I). Other HDAC inhibitors (mocetinostat, entinostat, and SB-939) also showed similar effects on the perturbation of 4-R genes (Fig. 7F-H and Table S10). These findings suggest that HDAC inhibitors may target BRCA1 mutation-associated risk genes and inhibit LUAD tumor growth in vitro, offering a preliminary therapeutic hypothesis for further testing.

4. Discussion

Our study provides a comprehensive molecular and cellular characterization of the TME in LUAD with BRCA1 and BRCA2 mutations, using large-scale multi-omics data and single-cell data. We describe transcriptional and immune patterns associated with these mutations, which may contribute to genomic instability, immune cell infiltration, and immune modulation within the TME. These observations offer hypotheses for understanding HRR deficiency in LUAD and for developing potential therapeutic strategies, including ICB and targeted therapies. HRR deficiency is a key contributor to tumorigenesis and treatment response in multiple cancers [7]. Although HRR deficiency has been correlated with enhanced sensitivity to DNA-damaging agents and ICB in cancers such as breast, ovarian, and NSCLC, our study suggests that in LUAD, BRCA1/2 mutations may be associated with distinct TME features that could influence immunotherapy effectiveness. This underscores the complexity of BRCA-driven biology in LUAD.

A key aspect of this study is the single-cell transcriptomic profiling of BRCA1/2-mutant LUAD tumors, paracancerous tissues, and peripheral blood patients. This enabled high-resolution examination of immune signatures, though the sample size limits generalizability. BRCA1 mutations were associated with the upregulation of cGAS-STING and type I IFN response gene, which can enhance innate immune responses, antigen presentation, and cytotoxic T-cell activity [52]. BRCA1-mutant tumors showed higher CD8+ T cell activation and MHC-I expression, possibly due to increased neoantigen load. This suggests that BRCA1-mutant tumors may be more immunogenic and potentially responsive to immunotherapy. In contrast, BRCA2 mutation was associated with upregulated MHC-II expression, which plays a role in antigen presentation and CD4+ T cell recruitment [53], suggesting a shift toward adaptive immunity and antigen-presenting cell involvement. BRCA2-mutant tumors showed higher CD4+ T cell activation, particularly inflammatory subsets, and were associated with increased NK cell activity and T cell cytotoxicity in blood, suggesting systemic immune engagement.

Despite these immune signatures, BRCA1/2 mutations were associated with downregulation of CD28 family co-stimulatory molecules, essential for effective T cell activation [54]. This may indicate adaptive resistance mechanisms. ICB may overcome such resistance (e.g., blocking PD-L1), unleashing the pre-existing, infiltration-potentiated immune response (Fig. 1G). We identified two Trm subsets (CD8+ Trm and CD4+ Trm) as potential biomarkers for ICB response in NSCLC. BRCA1 mutation was associated with CD8+ Trm TCR expansion, while BRCA2 mutation was linked to CD4+ Trm expansion, which may contribute to maintain a chronic inflammatory TME [55, 56]. Given the observed T cell infiltration and IFN-driven immune activation, BRCA1-mutant tumors might initially respond well to anti-PD-1/PD-L1 therapies. However, BRCA1 mutations also activated risk genes (S100A10, LDHA, GAPDH, and MYL12A). We identified HDAC inhibitors as potential agents targeting this program. HDAC inhibitors can modulate chromatin accessibility, reduce immune suppression, and enhance tumor immunogenicity, suggesting that combining them with ICB might enhance immune response and overcome resistance mechanisms in BRCA1-driven tumors.

Despite the comprehensive multi-omics and single-cell analyses conducted in this study, several limitations should be acknowledged. Our study primarily relied on retrospective data from public cohorts and a limited number of in-house samples. Although we included both tumor tissues and matched paracancerous and blood samples, the sample size remains modest, particularly for patients with BRCA1 or BRCA2 mutations. While our single-cell RNA-seq and TCR-seq analyses provide high-resolution insights into the cellular and clonal dynamics of the TME, the functional validation of key mechanisms—such as the causal role of BRCA1-driven cGAS-STING activation in CD8+ T cell recruitment—remains largely correlative. Experimental models, such as BRCA1/2-knockout cell lines or patient-derived organoids, coupled with in vivo studies, are required to establish causal relationships and elucidate the underlying molecular pathways. Finally, the therapeutic implications of HDAC inhibitors in BRCA1-mutant LUAD, although supported by in vitro data, warrant further preclinical and clinical investigation. The efficacy and safety of combining HDAC inhibitors with ICB in BRCA1-mutant LUAD patients remain unexplored and represent an important avenue for future research.

5. Conclusions

Overall, this study provides a high-resolution molecular map of the TME in BRCA1/2-mutant LUAD, highlighting distinct immune and transcriptional signatures associated with each mutation in the profiled cases. By uncovering these molecular distinctions, our study offers a hypothesis-generating framework for understanding HRR deficiency in LUAD and lays groundwork for future validation and precision immunotherapy strategies tailored to BRCA1/2 mutation status.

Data availability

The raw data of single-cell sequencing (including scRNA-seq and scTCR-seq) and the processed data generated in this study have been deposited at the Gene Expression Omnibus (GEO) repository, with the accession code GSE292700. Open-source R packages and software, as well as standard workflows were used in this study. No previously unreported custom code was developed for the analyses presented. All codes and other processed data are available from the corresponding author upon reasonable request.

Acknowledgements

The authors wish to thank the staff members from Affiliated Cancer Hospital and Institute of Guangzhou Medical University. High-throughput sequencing was performed by the OE Biotech Co., Ltd. (Shanghai, China). The contribution of the bioinformatics core facility at Guangzhou medical university is gratefully acknowledged.

Additional information

CRediT authorship contribution statement

Gaoming Liao: Writing - original draft, Visualization, Formal analysis, Data curation, Conceptualization, Software, Methodology, Funding acquisition. Xinbin Yang: Writing - review & editing, Investigation, Methodology. Qi Liu: Writing - review & editing, Investigation, Validation. Shufeng Nan: Validation, Software. Yan Liu: Software. Jinwei Li: Validation. Si Huang: Validation. Wang Ning: Validation. Xionghai Qin: Writing - review & editing, Supervision, Conceptualization. Gang Xu: Writing - review & editing, Project administration, Supervision, Conceptualization.

Funding

This study was supported by the National Natural Science Foundation of China (Grant No. 82403816), the Postdoctoral Fellowship Program of CPSF (Grant No. GZB20230180), the Guangdong Basic and Applied Basic Research Foundation (Grant No. 2023A1515110297), the China Postdoctoral Foundation (Grant No. 2023M740846), and the Fundamental Research Funds for the Provincial Universities (2023-KYYWF-0212).

Abbreviations

  • NSCLC: Non-small cell lung cancer

  • LUAD: Lung adenocarcinoma

  • HRR: Homologous recombination repair

  • HRD: Homologous recombination repair deficiency

  • TMB: Tumor mutation burden

  • TME: Tumor microenvironment

  • ICB: Immune checkpoint blockade

  • Trm: Tissue-resident memory T cell

  • HDAC: Histone deacetylase

  • BRCA1/2: BReast CAncer gene 1/2

  • TCGA: The cancer genome atlas

  • UMI: Unique molecular identifier

  • UMAP: Uniform manifold approximation and projection

  • TCR: T-cell receptor

  • MSigDB: Molecular signatures database

  • ssGSEA: Single sample gene set enrichment analysis

  • NMF: Nonnegative matrix factorization

  • cMap: Connectivity Map

  • LINCS: Library of integrated network-based cellular signatures

  • CCLE: Cancer cell line encyclopedia

Funding

MOST | National Natural Science Foundation of China (NSFC) (82403816)

  • Gaoming Liao

China Postdoctoral Science Foundation (中国博士后科学基金会职) (2023M740846)

  • Gaoming Liao

Fundamental Research Funds for the Provincial Universities (2023-KYYWF-0212)

  • Xionghai Qin

Postdoctoral Fellowship Program of CPSF (GZB20230180)

  • Gaoming Liao

Guangdong Basic and Applied Basic Research Foundation (2023A1515110297)

  • Gaoming Liao

Additional files

Supplementary figures

Uncropped immunoblot images

Supplemental tables