Comprehensive genomic and transcriptomic profiling of pulmonary nodules in synchronous multiple primary lung cancer

Abstract

Background

Synchronous multiple primary lung cancer (sMPLC) exhibits distinct histopathological characteristics among pulmonary nodules. However, a comprehensive understanding of the somatic mutation landscape and transcriptome heterogeneity is lacking. Therefore, our study aims to meticulously investigate genomic distinctions among multiple pulmonary nodules within individual patients.

Methods

We performed targeted DNA sequencing on tumor specimens and conducted bulk RNA transcriptome analysis on 53 multiple nodules originating from 26 lung cancer patients. The multiple nodules from the same patient was determined as major nodule and minor nodule. Additionally, the tumor tissues underwent histopathological evaluation through H&E staining, complemented by a comprehensive series of immunohistochemistry (IHC) analyses to detect protein expression. The detected protein markers encompassed PD-L1, Ki67, and others.

Results

For the 53 nodule samples from 26 MPLCs patients, EGFR was the mostly mutated genes, and the TP53 mutation frequency was notably different between major and minor nodules. Furthermore, pathway enrichment analysis based on the differentially expressed genes (DEGs) between major and minor nodules revealed the significantly active cell cycle and p53 pathways in the major nodules. Additionally, both major and minor nodules demonstrated mostly similar immune microenvironment and PD-L1 protein expression, and a significantly higher expression of Ki67. A noteworthy suppression was observed in the immune microenvironment in nodules, revealed by the expression of macrophage, neutrophils, and NK cells. Furthermore, minor nodules exhibited a modestly elevated expression of macrophages compared to major nodules. Additionally, among the significantly up-regulated cell cycle-related genes in the major nodules when compared with minor nodules, CCNE1 mRNA expression demonstrated significant correlation with poor prognosis in the lung cancer. Furthermore, the MYC inhibitor demonstrated more sensitivity for the major nodules than minor nodules.

Conclusions

This study validated molecular distinctions between samples from major and minor nodules in patients with sMPLC at both genomic and transcriptomic levels. The major nodules exhibited heightened activity in tumor cell proliferation pathways and demonstrated malignancy-related biological characteristics, which correlated with pathological assessment results.

Background

Lung cancer ranks second among global malignant tumors, comprising 11.4% of new cancer cases and standing as the primary cause of cancer-related deaths. It accounts for a global mortality rate of 1.8 million annually [1]. The utilization of high-resolution chest imaging systems and lung cancer screening programs has led to a rising incidence of multiple pulmonary nodules, particularly multiple ground-glass opacities (GGOs), in clinical practice worldwide [2]. MPLCs refer to the presence of two or more primary pulmonary cancers within a single patient, and can be classified as synchronous multiple primary lung cancer (sMPLC) or metachronous multiple primary lung cancer (mMPLC) according to the time when the carcinoma were confirmed. When primary lung cancers coexist, they are termed synchronous multiple primary lung cancer (sMPLC) [3].

The Fleischer Society's position and the consensus of the International Association for the Study of Lung Cancer classify multiple GGOs as MPLCs rather than intrapulmonary metastasis. However, there is insufficient evidence to definitively ascertain whether multiple GGO lesions are clonal or not [4]. Traditionally, the classification of synchronous multiple primary lung cancers (sMPLCs) requires histological distinctions, anatomical separation, or a clear origin from known carcinoma in situ for coexisting pulmonary nodules [5]. Indeed, molecular profiling has shown the capability to discern distinctions in multiple primary lung cancers (MPLCs) that may elude histological assessment [6]. In general, employing broader mutation panels for clonality assessment in non-small cell lung cancer (NSCLC) enhances accuracy by capturing mutations on a larger scale. A more refined interpretation of the clonal origin and carcinogenesis of sMPLC within the existing literature can aid in treatment strategies and contribute to improved patient survival [7]. In concordance with my research findings, driver mutations with equivalent functions in tumorigenesis exhibit mutual exclusivity. Meanwhile, sMPLCs may also demonstrate molecular similarities due to parallel evolution resulting from shared carcinogen exposures or chance occurrences. Both of these factors introduce bias in clonality analysis and necessitate corresponding correction strategies [8]. Furthermore, exploring alternative programs that integrate genomic variations across different levels proves to be a stimulating avenue for assessing clonal relatedness.

Since the sMPLCs reveal different molecular characteristics, it implies that individual nodules may respond differently to a uniform therapeutic regimen. For instance, Zhang et al. observed that only one of the three nodules in an MPLC patient responded to three cycles of neoadjuvant pembrolizumab [9]. Due to the existence of the heterogeneous molecular clonal driving, the different nodules in the same patients reveal different clinical risk and therapeutic response to the same drugs. Consequently, there is an urgent need to delineate the heterogeneous molecular characteristics of sMPLCs to facilitate personalized and precise treatment strategies, thereby enhancing survival outcomes.

While extensive literature exists on gene mutations in solitary lung cancers, there is a notable dearth of research investigating clinicopathological and demographic characteristics associated with gene mutations in multiple primary lung cancers, particularly through comprehensive large sample studies. To bridge this knowledge gap, we conducted a systematic exploration of genomic profiles in sMPLC and delved into immune profiles within the tumor microenvironment. Through comprehensive profiling encompassing DNA, RNA-seq, TME, and I-TED among sMPLC cases, our study offers novel insights into the biology and potential carcinogenesis of sMPLC. These findings have implications for disease monitoring and inform future therapeutic interventions.

Material and methods

Patients and study design

From March 2020 to April 2023, the study enrolled 26 patients with sMPLC who had been undergone curative resection without involvement of mediastinal nodes and metastatic disease as identified by ACCP modified guidelines in our center [10]. Patients with a history of preoperative induction therapy, including chemotherapy or radiotherapy, or those with a prior malignancy were excluded from the study. The exclusion criteria for patients with prior malignancies or metastatic disease included: (1) any other malignancy diagnosed or relapsed at any time, which is currently being treated; (2) any other current malignancy or malignancy diagnosed or relapsed within the past 3 years. Written informed consent for genetic research was obtained from all participants, in accordance with protocols approved by the institutional review board at XuanWu Hospital (Ethic approval number: 2020002; approval data: between March 25, 2020 and March 25, 2021). Inclusion criteria for multiple tumors within the same lobe required that they be situated in distinct segments, originate from carcinoma in situ, or exhibit identical histologic features but different subtyping. Gene analysis was conducted on all excised nodules. Basic patient demographic information (i.e., age, gender, height, weight, smoking status, tumor markers, pulmonary function, tumor location, nodule size, and CT features), operative information (i.e., operative time, hemorrhage and operative type), and clinical outcomes (i.e., hospital stay, drainage and pathological type) from a multidisciplinary tumor board.

Additionally, multiple publicly available GEO datasets were utilized for exploring and validating the prognostic significance of genes in lung cancer, including GSE72094 [11], GSE13213 [12], GSE42127 [13], and GSE41271 [14].

CT evaluation

All patients underwent preoperative high-resolution CT scans. The lesions were classified based on established criteria from previous studies, including pure ground-glass opacification (pGGO), mixed ground-glass opacification (mGGO) and solid nodule [15]. Subsequent evaluation of all nodules involved estimating the extent of the GGO lesion using a thin-section CT scan with 1.5-mm collimation. Nodule sizes were preoperatively determined through careful analysis of thin-section CT scans. The two nodules are categorized as major and minor nodules by the multidisciplinary team based on the following criteria: (1) the major nodule was identified as the lesion posing the highest risk to patient survival, while the secondary nodule was deemed suitable for resection at the same time. (2) In patients (not suspected of having a second focus of cancer) who are found intraoperatively to have a second lesion in a lobe, the second lesion thought to be resected was the minor nodule.

Histopathology

Archived formalin-fixed, paraffin-embedded (FFPE) slides of sMPLCs were subjected to simultaneous reevaluation by two independent, board-certified pathologists and classified as AIS, MIA and IVA. Histological typing was performed by paired lesions according to 8th TNM classification guidelines [16]. In subsequent subgroup analysis, MIA and IVA were categorized as invasive pathological types, while AIS was classified as a non-invasive pathological type. Paired of lesions in the same patient were belonged to the same pathology subgroup or the different pathology subgroup according to their pathological type.

Nucleic acid extraction and NGS (next-generation sequencing) experiment

Professional pathologists assessed the tumor cell content in all FFPE samples. Samples with a tumor cell content exceeding 20% were chosen for DNA and RNA extraction in this study. Genomic DNA was extracted from tumor tissues employing the QIAamp DNA FFPE Tissue Kit (Qiagen, Hilden, Germany), utilizing FFPE samples. Peripheral blood samples served as the normal control of the DNA mutation calling, and DNA extraction was carried out using the MagPure Tissue&Blood DNA LQ Kit (Magen, Guangzhou, China). Nucleic acid concentrations were determined using either the Qubit HS dsDNA kit or the Qubit RNA HS Assay Kit (Invitrogen, Carlsbad, USA). The DNA and RNA samples were undergoing the next library construction of the next-generation sequencing (NGS), including targeted DNA sequencing and RNA sequencing.

The targeted DNA sequencing were performed following the standard operation procedure of 654-gene panel from Kanghui biotech Co., Ltd., Liaoning, China. The somatic mutations of SNVs, indels, CNVs, and fusions were called using a targeted DNA sequencing panel including 654 cancer-related genes 654 genes, and the germline mutations were excluded using normal blood samples. The 654 genes include tumor oncogene and suppressor genes, drug targets genes, therapeutic sensitive and resistance genes, etc., in lung cancer. Also, this large gene panel could provide the sufficient estimation of total mutational burden. The DNA libraries were prepared using the CS2.0 Tissue DNA Library Prep Kit (Kanghui biotech Co., Ltd., Liaoning, China). The genomic DNA (gDNA) was fragmented, and underwent end-repair, A-tailing, and ligation with adapters containing unique molecular identifiers (UMIs). Next, the ligation products were purified and PCR-amplified. Lastly, the enrichment were conducted for regions of interest through hybridization reaction using the CS2.0 DNA Hybridization and Wash Kit (Kanghui biotech Co., Ltd., Liaoning, China). And then, the DNA library was sequenced on the NovaSeq 6000 (Illumina, San Diego, USA) in paired-end 150 bp mode, with an average depth exceeding 1,000X for tissue samples. Data analysis was conducted using dedicated pipeline.

Total RNA of tumor tissues from FFPE samples and adjacent paracancerous tissues were extracted using TRIzol reagent (Invitrogen, San Diego, CA, USA). The total RNA concentrations were also quantified before the next experiment. RNA-sequencing (RNA-seq) analysis was conducted on an Illumina Hiseq 2500 platform. The library was prepared using rRNA-depleted RNA by NEBNext® Ultra™ Directional RNA Library Prep Kit for Illumina® (NEB, Ipswich, MA, USA) in a 20 ng-RNA reaction system. All sequencing procedures and analyses were carried out by Novogene Company (Beijing, China).

Targeted DNA sequencing data collecting and processing

All the raw data of NGS sequencing had passed the quality control, and in an average depth of 1000 × for targeted DNA sequencing. Next, these sequencing results were filtered following the strict criteria to ensure the accuracy of the data. Furthermore, duplicate reads and soft-clipped reads were removed before data analysis. Then, the sequences were aligned to the human reference genome (GRCh38) by Burrows-Wheeler Aligner (version 0.1.22), and the somatic mutations of SNVs and indels were called using VarScan software and IndelRealigner tool from Genome Analysis Toolkit (GATK, version 3.6). Functional annotations were established through ANNOVAR.

Somatic mutation profiles in nodule samples

The mutational landscape of 26 sMPLC patients was analyzed and visualized using 'ComplexHeatmap' R package, and the mutation landscape was plotted by 'Heatmap' function. The mutation waterfall plot displayed only nonsynonymous mutations, which encompassed missense mutations, nonsense mutations, frame-shift deletions, frame-shift insertions, in-frame deletions, and in-frame insertions.

Differentially expressed genes (DEGs) and pathway enrichment analysis

DEGs between major and minor nodules were identified using the 'DESeq2' R package. DEGs were screened with criteria of an absolute log2-fold change (|logFC|) > 1 and p < 0.05. The significantly up-regulated and down-regulated DEGs were subsequently subjected to Gene Ontology (GO) annotation and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis was performed by Metascape (http://metascape.org), aiming to obtain a comprehensive understanding of their biological distinctions. Additionally, the Gene Set Enrichment Analysis (GSEA) was analyzed on WebGestalt toolkit (www.webgestalt.org) to explore specific pathways associated with major and minor nodules, and the normalized enrichment scores (NES) represented the gene sets with significant enrichment (adjusted p < 0.05) after GSEA analysis.

Additionally, the single-sample gene set enrichment analysis (ssGSEA) score of gene signatures from previously experimentally validated and MsigDB was calculated using the GSVA R package.

The different immune cells infiltrating analysis

The fraction of immune infiltrating cells were evaluated using six deconvolution algorithm, including CIBERSORT [17], TIMER [18], quanTIseq [19], xCell [20], EPIC [21], and MCP-counter [22]. Additionally, the p-values were transformed using logarithm base 10 (log10). Due to the differences of assessed immune cell types in these methods, we mainly focused on the eight common immune cells, including B cell, T cell, CD4 T cell, CD8 T cell, neutrophil cell, macrophage monocyte cell, and myeloid dendritic cell.

The intratumour expression distance (I-TED) analysis

The correlation distance was calculated using the function 'dcor()' in the R package 'energy' (v.1.7-6). For each nodule, 1 minus the correlation distance between gene expression to all other nodules in the same patient was calculated for the top 500 most variable genes.

Drug sensitivity prediction

The Genomics of Drug Sensitivity in Cancer (GDSC) (https://www.cancerrxgene.org/) and the Cancer Therapeutics Response Portal (CTRP) database were used to download the drug sensitivity information and gene expression data. The R package 'oncoPredict' was used to predict the half-maximal inhibitory concentration (IC50) of drug sensitivity.

Statistical analysis

The primary statistical analyses were performed using R software (version 4.2.1), and the ‘ggplot2’ R packages was used for data visualization. In addition to the aforementioned software and R packages, we employed additional fundamental statistical analysis methods. Fisher's exact test was employed to assess differential mutation patterns among distinct histological subtypes. The Wilcoxon test was utilized to elucidate differences in gene expression across various histological subtypes. Survival analysis was performed utilizing the 'survival' and 'survminer' R packages to assess differences in survival rates between groups. Statistical significance was determined based on a two-sided FDR-adjusted p-value (or q-value) or a p < 0.05.

Results

Patients characteristics and operative outcomes

The workflow of this study is shown in Fig. 1. A total of 26 MPLCs patients with 53 independent lesions, as determined based on both histological and clinical assessments, were scheduled for curative thoracico-surgery. Patients had a mean age of 60.3 years. Eighteen of 26 patients were female (18/26, 69.2%), and eight were male (8/26, 30.8%). There were 21 never smokers (21/26, 80.8%), and 5 active smokers at the time of surgery (5/26, 19.2%). Only 5 patients (5/26, 19.2%) had a family history of lung cancer. All patients underwent video-assisted thoracoscopic resection. Tumors were found in all pulmonary lobes but were predominantly in the right upper lobe of the major nodules (10/26 nodules, 38.4%) and right upper lobe of the minor nodules (9/27 nodules, 33.4%). The number of nodules was 2 in 25 patients (25/26, 96.2%) and the three nodules we found was one patient (1/26, 4.8%). Other demographic characteristics are summarized in Table 1. As seen in Table 2, ten patients received lobectomy (10/26, 38.5%), and 16 patents received partial lobectomy (16/26, 61.5%). The mean operative time and hospital stay were 164.2 min and 9.3 days, respectively.

Fig. 1

The workflow of the study

Table 1 Clinicopathologic Characteristics of MPLC Cases (N = 26) and Tumors (N = 53)

Table 2 Operative characteristics of MPLC cases (N = 26) and tumors (N = 53)

Histologic review

All lesions, as confirmed by pathological reports, were limited to a maximum dimension of 3.0 cm. The results of the histologic review are documented for each nodule for every patient in Table 2. Among the 26 major nodules, the pathologic diagnosis of resected nodules was AIS in 4 of 26 specimens (4/26, 15.4%), MIA in 2 of 26 specimens (2/26, 7.7%), and IVA in 20 of 26 specimens (20/26, 76.9%). Among the 27 minor nodule, the pathologic diagnosis of resected specimens was AIS in 10 of 27 specimens (10/27, 37.0%), MIA in 8 of 27 specimens (8/27, 29.6%), IVA in 9 of 27 specimens (9/27, 33.4%). Consequently, among the 26 patients, 18 had the same pathology type between major and minor nodules, whereas eight exhibited different pathological type between the two nodule types (Table 3).

Table 3 The subgroups of major and minor nodules based on the historical pathology types

Somatic mutation profiles between major and minor nodules

For the all 53 nodules samples, the most frequently mutated genes were EGFR (69.81%), followed by RBM10 (22.64%), TP53 (15.09%), BRAF (7.55%), ERBB2 (5.66%), KRAS (5.66%), MED12 (5.66%), and MUC16 (5.66%) (Fig. 2A). Of which, the mutation frequency of TP53 significantly differed between major and minor nodules (p = 0.02), predominantly manifesting within the major nodules (Fig. 2B). Additionally, the single-base substitutions (SBS) also demonstrated no significant difference between major and minor nodule samples (Fig. 2C). The alterations in oncogene pathway revealed that checkpoint factors (CPF) (p = 0.024) and p53 (p = 0.024) pathways were obviously higher in the major nodules samples when compared with minor nodules (Fig. 2D). In addition, there were no obvious difference in mutation counts (p = 0.66, Fig. 2E), mutant-allele tumor heterogeneity (MATH) score (p = 0.31, Fig. 2F), and PD-L1 expression (p = 0.38, Fig. 2H) between the two nodule types, while the significantly higher Ki67 expression was shown in major nodules than minor nodules (p = 0.048, Fig. 2G).

Fig. 2

Genomic disparities between major and minor nodules. A The heatmap displays the somatic profiles and clinical pathology information in patients with major nodules, indicated by yellow boxes, and minor nodules, represented by purple boxes. B The SNV variations in mutation frequency between major and minor nodules. C The frequency differences of single base substitutions between major and minor nodules. D The pathway alterations in the major and minor nodules. The differences between major and minor nodules in E mutation counts, F Math score, G ki67 protein expression, and H PD-L1 protein expression

The differences in RNA expression profile between nodules and adjacent paracancerous tissues in the MPLC patients

To identify the specific characteristics of major and minor nodules, differentially expressed genes (DEGs) were analyzed among the groups. A total of 2006 DEGs were found between major nodules and adjacent paracancerous tissue (MAP-DEGs), consisting of 1207 up-regulated and 799 down-regulated genes. Additionally, 1120 DEGs were detected between minor nodules and adjacent paracancerous tissue (MIP-DEGs), comprising 756 up-regulated and 264 down-regulated genes.

The KEGG enrichment analysis revealed that the up-regulated MAP-DEGs in major nodule samples mainly involved in the biological functions of cell cycle, p53 signaling pathway (Fig. 3A). Conversely, the significantly down-regulated MAP-DEG in major nodule samples related to calcium signaling pathway, cAMP signaling pathway, and ECM–receptor interaction in adjacent paracancerous tissues (Fig. 3A). The GSEA analysis further confirmed these findings, and also demonstrated the active drug metabolism pathway in the major nodules and a suppressed natural killer cell mediated cytotoxicity in the adjacent paracancerous tissues (Fig. 3B).

Fig. 3

Identification of differentially expressed genes (DEGs) and pathway enrichment analysis between nodules and adjacent paracancerous tissue. A The pathway enrichment based on the significantly DEGs in the major nodules and adjacent paracancerous tissue. B GSEA analysis unveiled the different function between major and adjacent paracancerous tissue. C The pathway enrichment based on the significantly DEGs in the minor nodules and adjacent paracancerous tissue. D GSEA analysis unveiled the different function between minor and adjacent paracancerous tissue. E The pathway enrichment based on the overlapping DEGs in the nodules and adjacent paracancerous tissue. In the GSEA plot, a positive NES indicates that the gene set is associated with major nodules or minor nodules, while a negative NES indicates that the gene set is associated adjacent paracancerous tissue

In addition, the up-regulated MIP-DEGs revealed significant correlations with immune response pathways, such as cytokine–cytokine receptor interaction, IL-17 signaling pathway, and intestinal immune network for IgA production pathways (Fig. 3C). The GSEA analysis revealed an active taurine and hypotaurine metabolism, oxocarboxylic acid metabolism, glycosphingolipid biosynthesis, and central carbon metabolism in cancer in minor nodules (Fig. 3D).

Subsequently, we identified 903 overlapping DEGs (MM-OGs) of MAP-DEGs and MIP-DEGs, of which 674 were up-regulated and 229 were down-regulated. We found the up-regulated MM-OGs was associated with cytokine–cytokine receptor interaction, salivary secretion, cell adhesion molecules, drug metabolism, B cell receptor signaling pathway, and complement and coagulation cascades pathways. This observation indicates that within the major and minor nodular samples, the body counteracted the proliferation of tumor cells through the regulation of metabolism transduction and immune stress responses (Fig. 3E).

The functional difference between major nodule and minor nodule in the MPLC

Furthermore, to elucidate the molecular and functional characteristics of major and minor nodules, a total of 213 DEGs were identified between major and minor nodules (MM-DEGs) for the function enrichment analysis, with 184 significantly up-regulated genes and 29 significantly down-regulated gene in major nodules.

Furthermore, the up-regulated MM-DEGs in major nodules mostly involved in cell cycle pathways (Fig. 4A).

Fig. 4

Identification of differentially expressed genes (DEGs) and pathway enrichment analysis between major nodules and minor nodules samples. A The pathway enrichment based on the significantly up-regulated and down-regulated between major and minor nodules. B GSEA analysis unveiled the different function between major and minor nodules. In the GSEA plot, a positive NES indicates that the gene set is associated with major nodules, while a negative NES indicates that the gene set is associated minor nodules. C The different ssGSEA scores between major and minor nodules in cell proliferation and development related pathways

These findings unveiled a more active cell cycle pathways in the major nodule samples. Furthermore, the GSEA enrichment analysis also highlighted the different biology function in cell cycle, p53 signaling pathway, DNA replication, homologous recombination, proteasome, and base excision repair, which were active in the major nodule samples and suppressed in the minor nodules (Fig. 4B).

To compare active proliferation characteristics and progression risk, we analyzed the gene set scores between major and minor nodules. The major nodules exhibited higher signature scores in the MYC, epithelial–mesenchymal transition (EMT), and cell cycle-related signatures, including the E2F target, G2M, cell cycle, and mitotic spindle signatures (Fig. 4C). These findings suggest significant differences in tumor cell proliferation and development, indicating that major nodules present a higher risk of metastasis than minor nodules.

The subgroup analysis unveiled the influences of pathology types on the functions of major and minor nodules

Given the notable differences in gene expression and pathway enrichment observed between major and minor nodules, our investigation aimed to determine whether consistent disparities existed among subgroups.

We compared the DEGs between major nodule and minor nodules in subgroups that the nodules were with the same pathology type (SaMM-DEGs) or the different pathology type (DiMM-DEGs), and the subgroups that the major nodule was invasive pathology type (ISMM-DEGs) or non-invasive pathology type (nISMM-DEGs) in one patient.

Based on the 68 significantly down-regulated SaMM-DEGs, the GSEA analysis showed an active of cell cycle, DNA replication, p53 signaling pathway, and homologous recombination pathway, and a suppressed cGMP-PKG signaling pathway and cholesterol metabolism in the major nodules (Figure S1A). In addition, a total of 214 up-regulated and 98 down-regulated DiMM-DEGs were identified. The pathways of cell cycle, p53 signaling, and chemical carcinogenesis pathways were significantly enriched in the major nodules, revealing the regulation of tumor cell proliferation and the relaxation of cell cycle arrest. Additionally, the PPAR signaling and neuroactive ligand–receptor interaction pathways was suppressed in the major nodules (Figure S1B).

Regardless of whether the pathological types of the major and minor nodules were the same or not, the function-related signaling pathways identified by GSEA enrichment analysis were consistent (Figure S1C). The active cell cycle-related pathways were the typical character of the major nodules when compared with minor nodules.

Additionally, in instances where patients exhibited invasive pathology type of major nodules, we identified 145 up-regulated and 36 down-regulated ISMM-DEGs. The up-regulated ISMM-DEGs demonstrated significant enrichment in functions related to cell cycle and epithelial cell differentiation pathway. The down-regulated ISMM-DEGs indicated significant enrichment in pathways related to positive regulation of peptide hormone secretion biological function (Figure S2A). The GSEA analysis results further supported these findings, and the major nodules exhibited obvious active of cell cycle, mismatch repair, and p53 signaling pathway in major nodules (Figure S2B).

Conversely, in cases where major nodules presented non-invasive pathology type, we identified 157 up-regulated and 100 down-regulated nISMM-DEGs. The up-regulated DEGs were primarily associated with pathways related to cell cycle. Conversely, the down-regulated DEGs were predominantly linked to pathways involving C-type lectin receptor signaling, cAMP signaling, and calcium signaling pathways (Figure S2C).

GSEA analysis revealed a significant enrichment of the cell cycle-related pathways in major nodules when major nodules exhibited non-invasive pathology (Figure S2D).

The results of these subgroup analyses further substantiated the biological functional differences between major and minor nodules, indicating that such distinctions were not influenced by the concordance of pathological types between major and minor nodules or the specific pathological type of major nodules.

Identification of survival-related cell cycle genes in lung cancer

Due to the significant enrichment of cell cycle-related pathways in major nodules, we further identified DEGs associated with the cell cycle function. Notably, CCNE1, a key oncogene in the cell cycle pathway, exhibited higher expression in major nodules compared to minor nodules and normal tissues within our cohort (p = 0.026, Fig. 5A).

Fig. 5

Prognostic and survival implications of CCNE1. A The different expression of CCNE1 between major and minor nodule samples. The validation of prognosis value of CCNE1 in multiple GEO datasets, including B GSE72094 (stage IA LUAD), C GSE42127 (stage IA NSCLC), D GSE42127 (stage I NSCLC), E GSE41271 (stage I NSCLC from Japan), F GSE13213 (stage I–IV LUADs from Japan)

Given the absence of follow-up survival information, we explored the prognostic implications of CCNE1 using multiple public GEO datasets. High CCNE1 expression was significantly associated with poor prognosis in stage IA lung adenocarcinoma (LUAD) patients in the GSE72094 cohort (p = 0.048, Fig. 5B). Furthermore, consistent prognostic value was observed for stage IA NSCLC patients (p = 0.028, Fig. 5C) and stage I NSCLC patients (p = 0.0011, Fig. 5D) in the GSE42127 cohort. Additionally, we validated the prognostic value of CCNE1 expression across the stage I of the GSE41271 cohort (stage I LUADs), revealing a significant association between high CCNE1 expression and poor prognosis in NSCLC patients (p = 0.0013, Fig. 5E). Furthermore, the prognosis value of CCNE1 mRNA expression was also confirmed in stage I LUAD patients (p = 0.0018, Fig. 5F) in the entire population of GSE13213 (stage I–IV LUADs) cohort from Asian population cohort, comprising patients from Japan.

These findings underscore the potential significance of CCNE1 as a prognostic marker in lung cancer.

Immune cells infiltration difference between major and minor nodules

For a comprehensive analysis of immune cell infiltration, we employed six RNA-seq data-based immune cell quantification methods, namely CIBERSORT, TIMER, quanTIseq, xCell, EPIC, and MCP-counter. The immune cell expression profiles of nodules and corresponding normal tissue samples were thoroughly examined (Fig. 6A). Subsequently, we compared the distinctions among major nodules, minor nodules, and normal samples. The enrichment score of immune cell sub-population in the major and minor nodules showed no significant differences, except for macrophage which expressed a little higher in the minor nodules. Both major and minor nodules exhibited a marked up-regulation of B cells. Conversely, the immune microenvironment in both major and minor nodules demonstrated a notable suppression in the expression of macrophage, neutrophils, and NK cells.

Fig. 6

The comparison of immune microenvironment and I-TED between major and minor nodules in different groups. A The immune cell infiltration analysis of major nodule, minor nodule, and normal tissues using different methods. B The enrichment proportions of immunologically 'hot' and 'cold' tumors among major nodules, minor nodules, and their respective paired adjacent paracancerous tissue samples. The different fraction of immune 'hot' and 'cold' tumors in the major and minor nodules when the pathology type of major and minor nodules were C identical and D different. E The different proportions of 'cold' tumors, 'hot' tumors, and mixed 'cold-hot' tumors in the patients when the pathology type of major and minor nodules were identical and different. The different fraction of immune 'hot' and 'cold' tumors in the major and minor nodules when the pathology type of major nodule was F invasive and G non-invasive. H The different proportions of 'cold' tumors, 'hot' tumors, and mixed 'cold-hot' tumors between the patients with pathology type of major nodule was invasive and non-invasive. I The I-TED difference when the pathological types of major and minor nodules were either identical or different. J The I-TED difference when the major nodule exhibited either invasive or non-invasive pathology types. I-TED: intratumour expression distance

Additional, based on the identification method for immunologically 'hot' and 'cold' tumors in Jia et al.[23]. We assessed the enrichment proportions of immunologically 'hot' and 'cold' tumors among major nodules, minor nodules, and their respective paired adjacent paracancerous tissue samples (Fig. 6B). Among the 26 major nodules, eight (30.8%) were classified as 'hot' tumors, while 18 (69.28%) were categorized as 'cold' tumors. In contrast, among the 27 minor nodules, 13 (48.1%) were characterized as 'hot' tumors, and 14 (51.9%) as 'cold' tumors. We found no significant difference in the fraction of immunologically 'hot' and 'cold' tumors between major and minor nodules.

When the pathology types of major and minor nodules were either identical or different, the ratio of 'cold' and 'hot' tumors exhibited no notable variance (Fig. 6C, D). Regardless of whether the pathological types of major and minor nodules were the same in the patients, the proportions of 'cold' tumors, 'hot' tumors, and mixed 'cold-hot' tumors exhibit remarkable similarity (Fig. 6E).

In instances where the major nodule exhibited either invasive or non-invasive pathology types, the distribution ratio of 'cold' and 'hot' tumors between major and minor nodules showed no significant differences (Fig. 6F, G). Remarkably, a statistically significant difference in the distribution ratio of cold and hot tumors was observed when contrasting patients with invasive and non-invasive types of major nodules (p = 0.05, Fig. 6H). This underscores that the immune cell infiltration of major and minor nodules remain unaffected by the pathological type.

When the major nodule was 'hot', regardless of whether the minor nodule was 'hot' or 'cold', the genes significantly up-regulated in the major nodules were associated with allograft rejection and the NFKB/TNF regulatory pathway. In contrast, when the major nodule was 'cold' it exhibited distinct high-risk features. The up-regulated genes in these 'cold' major nodules were significantly related to metastasis pathways, while the down-regulated genes were associated with hypoxia, IL-6 response, inflammatory response, and apoptosis biological processes. These findings reveal that the heterogeneous immune microenvironment of major and minor nodules influences biological function, with a more pronounced effect observed in the major nodules (Figure S3).

These findings implied the pathology type of major nodules was the primary factors affected the immunologically “hot” and “cold” tumors distribution in the patients, and the immune infiltration of major nodules affected more on the biological function of major and minor nodules.

Expression diversity in MPLC patient

In accordance with the TRACERx series studies, the intratumour expression distance (I-TED) metric is determined as the mean normalized gene expression correlation distance within a specific region paired with every other region from the same tumor. A higher I-TED value is indicative of increased intratumoral heterogeneity (ITH) in gene expression [24].

The I-TED value was not distinctly different when the pathological types of major and minor nodules were either identical or different (Fig. 6I). Furthermore, when the major nodule exhibited either invasive or non-invasive pathology types, the I-TED value demonstrated no significantly difference (Fig. 6J).

Due to the diversity of major and minor nodules, we aimed to predict the sensitivity of commonly used drugs for lung cancer. As is well-known, a lower IC50 value indicates higher drug sensitivity. We identified several drugs with significantly different sensitivities between major and minor nodules. Major nodules exhibited a notably better response to the PI3K inhibitor (AZD8186), MYC inhibitor (sepantronium bromide), and BCL-2 inhibitor (ABT-199), suggesting more active proliferation characteristics. Sepantronium bromide may be a better choice for combination therapy to improve the response of major nodules. Conversely, minor nodules showed a better response to staurosporine, the IQGAP3 inhibitor (BRD-K88742110), and the ARFGAP1 inhibitor (QS-11). Additionally, both types of nodules demonstrated similar sensitivity to the EGFR inhibitor, paclitaxel, carboplatin, and docetaxel. Furthermore, these drugs revealed lower IC50 values than other drugs in lung cancer. Afatinib and osimertinib demonstrated lower IC50 values than other EGFR inhibitors. In summary, these drug sensitivity predictions provide valuable insights for the clinical combination of drugs to effectively treat major and minor nodules (Fig. 7).

Fig. 7

The predicted IC50 of different antitumor drugs in major and minor nodules based on the GDSC and CTRP database

Discussion

Several conventional technologies, including polymerase chain reaction (PCR), array-based comparative genomic hybridization (aCGH), and fluorescence in situ hybridization (FISH), have been utilized for detecting genetic alterations. However, when assessing numerous genes across multiple lesions, these traditional techniques often prove to be time-consuming and have limitation. Simultaneously, data collected from various methods may exhibit contradictions. Fortunately, with the development of bioinformatic analyses methods, it tends to be easy to explore the genome characteristics of MPLCs through targeted DNA sequencing. Leveraging bioinformatic tools, NGS has the capacity to explore sequencing data beyond somatic mutations, capturing information at the chromosomal level or beyond. Moreover, advanced innovations, such as RNA-seq, coupled with subclone reconstruction, metabolomics, and multiomics, offer significant potential beyond conventional techniques to assist scientists in unraveling the intricate complexities of MPLC.

MPLC can arise from both intrinsic and non-intrinsic cancer risk factors. Intrinsic risk factors encompass genetic mutations induced by DNA replication errors, including mutations in EGFR, KRAS, TP53, or PARP1 [25]. In many previous study, paired nodules in the single patient with sMPLC had the same gene mutation profiles [26]. Our study analyzed mutation profiling of MPLCs by NGS technology, which revealed that mutations in EGFR, RBM10, TP53, BRAF, ERBB2, KRAS, MED12, and MUC16 mutations were the most prevalent in MPLCs and randomly allocated to both groups. However, TP53 were significantly more frequent in the major nodule than the minor group. The tumor suppressor gene TP53 is the most frequently mutated gene in human tumors [1]. In TP53, many of the observed mutations in cancer are found to be single nucleotide missense variants [27]. These variants are broadly distributed throughout the gene, but with the majority localizing in the DNA binding domain. To fulfill its proper biological function, four TP53 polypeptides must form a tetramer which functions as a transcription factor [28]. While a large proportion of cancer genomics research is focused on somatic variants, TP53 is also of note in the germline [28]. Germline TP53 mutations are the hallmark of Li–Fraumeni syndrome [29], and many (both germline and somatic) variants have been found to have a prognostic impact on patient outcomes and might be commonly associated with intrinsic factors in major nodules of sMPLCs.

Due to rapid advancements in modern molecular biology technology, the treatment paradigm for lung cancer is centered around targeting aberrant molecules within specific signaling pathways [30]. Ki-67 antigen is a marker of cellular cycle and proliferation, usually used to estimate the cell’s population proliferation, also indicating the cell growth ratio [31]. Martin et al. suggested the existence of an association between elevated Ki-67 expression and the reduced survival in patients with lung cancer [32]. Motohiro Yamauchi and his coworkers indicate that checkpoint factors plays a pivotal role in G1 arrest, which amplifies G1 checkpoint signals sufficiently for phosphorylating p53 in cells leading to tumor development [33]. The cell cycle pathway is one of the most important cellular signal pathways that determine whether cells will survive or die when encountering any DNA-damaging factors. Cell cycle checkpoints occur at the G1/S and G2/M transitions as well as at the intra-S-phase [34]. One of the hallmarks of human cancers is the alteration of many signaling pathways, leading to the loss of basic cell cycle [35]; this results in unrestrained cell proliferation, cell cycle deregulation and ultimately, cancer development [34]. The most frequent mutational events in lung cancer occur to TP53, which are critical players in cell cycle control [36].

Cyclin E1 (CCNE1) is a protein-coding gene that belongs to the cyclin family of genes which controls the G1/S phase transition of the cell cycle. Previously, its abnormal expression pattern has been examined and found to be correlated with ovarian and breast cancer progression [37, 38]. In the study of Md Asad Ullah, the prognostic and therapeutic values of the CCNE1 gene in LUAD and LUSC have been explored using comprehensive bioinformatics analysis and a database mining approach [39]. In our study, the CCNE1, as significantly up-regulated cell cycle associated genes in the major nodules, are negatively associated with the OS of the early-stage LUAD patients in multiple public GEO datasets, which postulated that the survival of patients with sMPLCs is mainly determined by the stage of the major nodules. Because CCNE1 elevated expression was enriched in TP53-mutant tumors, it is possible that TP53 mutation is important to prevent cell death or cell-cycle arrest upon increased Cyclin E/CDK2 activity [40]. Therefore, TP53, Ki-67, and CCNE1 are thought to be the intrinsic factor leading to the formation of the major nodule of sMPLCs via the cell signal pathway of cell cycle.

To enhance comprehension of the sMPLCs, immune microenvironmental characteristics were also investigated. Prior investigations have examined the feasibility and safety of neoadjuvant therapy involving immune checkpoint inhibitors (ICIs) for early-stage NSCLC [41]. Tumor mutational burden (TMB) has been established as a valuable predictive biomarker in these studies [42]. In both major and minor nodules, a noteworthy suppression was observed in the immune microenvironment, affecting the expression of macrophage, neutrophils, and NK cells, and PD-L1 protein expression. Furthermore, minor nodules exhibited a modestly elevated expression of macrophages compared to major nodules. Notably, traditional TME markers may not apply to the description of all sMPLCs, particularly those with indolent microenvironment such as multiple ground-glass nodules (GGNs). Therefore, we conducted a more comprehensive analyses of TME immunophenotypes and found that hot and cold tumors allocated randomly in major and minor groups. Like previous study, multi-parameter and multi-locus analysis will aid efforts to assess the territorial and heterogeneous features of immune niches among nodules, resulting in the development of more precise and personalized immune therapies to treat sMPLCs [23]. Liu et al. postulated that ICI may be recommended for patients with MPLC [43], and this hypothesis was further confirmed by our data in which the paired nodules in a same patient own the similar TME markers. A subsequent clinical study is planned to investigate adjuvant therapy with ICIs in sMPLC patients.

In a research of NSCLC based on mouse models, it has been demonstrated that transcriptomic plasticity underlies ITH [44]. Cancers feature the accumulation of genetic abnormalities and dynamic development of subclones with distinct molecular variations. Here, we assessed the I-TED referring the TRACERx study, which is calculated by the mean normalized gene expression correlation distance for a given major nodule paired with minor nodule from the same patient. I-TED was not associated with the heterogeneity of subclonal mutations nor the number of regions sampled per tumor [24]. Our results revealed no significant I-TED value difference between the paired two nodules when the pathology type differs between major and minor nodules, which demonstrated that both genomic and transcriptomic data within these samples are able to be reality and have clinical significance between these two nodules type.

Additionally, we conducted preliminary analysis on the sensitivity of major and minor nodules to various antitumor drugs. Our predictive findings revealed that both major and minor nodules exhibit high sensitivity to EGFR inhibitors and anti-angiogenic drugs. However, major nodules showed significantly greater sensitivity to targeted therapies involving cell proliferation pathways such as MYC, PI3K, and BCL2, compared to minor nodules. The PI3K and MYC inhibitor could be promising as targeted or combined therapeutic drugs in patients with sMPLC. These predictions suggest a new insight for further investigation of clinical trials into the differential drug sensitivities of major versus minor nodules, which could aid clinicians in choosing effective combination therapy strategies to improve the survival benefit for sMPLC patients. Given that these results are predicted solely on public drug sensitivity experimental data, further exploration is required through clinical trials.

This study has several limitations. First, only a single arm of patients with sMPLCs was included without a control group. Second, the sample size was relatively small, which may introduce a bias in the results. The results in this study will be more robust if more samples were included and validated in external validation cohorts. We will try to improve the study in further exploration. Third, despite adhering strictly to the diagnostic criteria for sMPLCs in patient selection, the possibility of mixed metastatic cancer cases in the lung cannot be entirely rule out. This may lead to a certain degree of bias in the results. Thus, a larger sMPLC cohort is needed for further studies to validate more solid genetic features of MPLC.

Conclusions

Our study revealed the significantly different genome and transcriptome characteristics of major and minor nodules, underlying the molecular foundation for the nodule classification in the sMPLC patients. These findings also revealed potential distinctive biomarkers and pathway function enrichment, suggesting the importance of additional investigation for a comprehensive understanding of major and minor nodule development and the formulation of therapeutic strategies.