Dissecting transcriptional heterogeneity in multiple primary lung cancer by single cell RNA sequencing

Dissecting transcriptional heterogeneity in multiple primary lung cancer by single cell RNA sequencing | Research Square window.SnipcartSettings = { analytics: { enabled: false } }; (function() { var accessVector = localStorage.getItem('access_vector') || ''; window.dataLayer = window.dataLayer || []; if (accessVector) { window.dataLayer.push({ user: { profile: { profileInfo: { snid: accessVector } } } }); } })(); (function(w,d,s,l,i){w[l]=w[l]||[];w[l].push({'gtm.start':new Date().getTime(),event:'gtm.js'});var f=d.getElementsByTagName(s)[0],j=d.createElement(s),dl=l!='dataLayer'?'&l='+l:'';j.async=true;j.src='https://www.googletagmanager.com/gtm.js?id='+i+dl;f.parentNode.insertBefore(j,f);})(window,document,'script','dataLayer','GTM-K279D39R'); Browse Preprints In Review Journals COVID-19 Preprints AJE Video Bytes Research Tools Research Promotion AJE Professional Editing AJE Rubriq About Preprint Platform In Review Editorial Policies Our Team Advisory Board Help Center Sign In Submit a Preprint Cite Share Download PDF Research Article Dissecting transcriptional heterogeneity in multiple primary lung cancer by single cell RNA sequencing Wei Guo, Hao Xia, Qilin Huai, Xuemin Xue, Lei Guo, Fengwei Tan, and 5 more This is a preprint; it has not been peer reviewed by a journal. https://doi.org/ 10.21203/rs.3.rs-1716190/v1 This work is licensed under a CC BY 4.0 License Status: Posted Version 1 posted You are reading this latest preprint version Abstract The progression of multiple primary lung cancer (MPLC) involves complex changes in cell composition and metabolic function. Here, we performed scRNA sequencing of 167,397 cells from six patients with MPLC, combined with bulk whole-exome sequencing. We revealed that both naïve and memory T cells participate in the differentiation of CD8+ T cells. The terminal states of CD8+ T cells are exhausted T cells, which respond to stimuli and positively regulate cell death, and cytotoxic T cells, which are mainly implicated in the regulation of cytokine production. Multiple functional Tregs and naïve T cells contribute to the origin of CD4+ T cells. B cells, with two main functions, also play important roles in the immune response. We also uncovered the cellular metabolic activities that occur during tumor invasion. Positive regulation of blood vessel diameter has been observed in endothelial cells, while angiogenesis has not been found in early glandular neoplasia of the lung. The development of epithelial cells involves two functional states. In one state, cells respond to stress from the immune system, and in the other, some will undergo programmed cell death or enter the cell cycle due to the selective pressure that arises over the course of tumor development in synchronous MPLC. Our study showed the complete landscape of different dynamic cellular changes, which might reveal the key cellular mechanisms of MPLC and therefore provide new clues for the pathogenesis of tumors. multiple primary lung cancer single-cell RNA sequencing heterogeneity tumor microenvironment Figures Figure 1 Figure 2 Figure 3 Figure 4 Figure 5 Figure 6 Figure 7 Introduction Lung cancer is the leading cause of cancer-related mortality worldwide 1 – 3 . Multiple primary lung cancer (MPLC) is defined as the occurrence of two or more primary lung cancers in the same patient. MPLC can be classified as synchronous or metachronous based on the time of occurrence; metachronous MPLC is defined by a diagnostic interval of 6 months or less between tumors 4 , 5 . Port et al. reported that MPLC was observed in approximately 16% of patients with operable stage I, II, and III non-small-cell lung cancer (NSCLC), usually by preoperative imaging analysis 6 . Concomitant with the rising prevalence of lung cancer, the incidence of MPLC is also increasing 4 , 7 . Thus, an in-depth understanding of its mechanism and better treatment options for MPLC are urgently needed. Tumor heterogeneity describes differences between tumors of the same apparent type but with different genetic backgrounds and different carcinogen exposures (intertumoral heterogeneity) or differences between cancer cells within the same tumor (intratumoral heterogeneity) 8 . MPLC, by definition, arises in different areas of the lung parenchyma within a single patient and shares an identical germline genetic background and environmental exposure profile 9 . Previous studies have demonstrated differences in genomic profiles, including somatic point mutations, copy number aberrations, chromosomal structural variations and even mutational spectra, between different MPLC lesions 9 – 11 . The highly complex and heterogeneous tumor ecosystem contains not only malignant cells but also interacting cells from the host, such as endothelial cells, stromal fibroblasts, and a variety of immune cells that control tumor growth and metastasis 12 . However, the heterogeneity of cell types in the tumor microenvironment (TME) in MPLC remains unclear. Thus, it is critically important to decode the complex interplay between tumor cells and the TME in MPLC. In the present study, we performed whole-exome sequencing (WES) and single-cell RNA sequencing (scRNA-seq) on 17 tumor samples and 6 matched normal samples from 6 patients with MPLCs. By comparing tumors at different stages in MPLCs and adjacent normal lung tissue, we comprehensively characterized the transcriptome features of the malignant cells, immune cells, and stromal cells of MPLCs, and we decoded dynamic changes in cell percentage, cell subtype heterogeneity, and intercellular interactions, providing new knowledge regarding the biological basis of MPLC and LUAD development. Methods Patients and sample collection Patients who were diagnosed with pathologically confirmed MPLC with two or three tumors in the ipsilateral or bilateral lungs and underwent surgery (surgery for contralateral lung tumors was completed within 1–2 months after the first surgery, depending on the patients’ recovery status) at the Department of Thoracic Surgery of the National Cancer Center/National Clinical Research Center for Cancer/Cancer Hospital, Chinese Academy of Medical Sciences and Peking Union Medical College, Beijing, China, were enrolled. After MPLC samples were resected, tumor tissues were cut into two pieces along the long axis: one was processed for pathological diagnosis, and the other was processed for scRNA-seq and WES. The clinical information and specimen information of these patients are summarized in supplementary table 1. This study was approved by the National Cancer Center/Cancer Hospital Ethics Committee (Approval number: NCC2208). Written informed consent was obtained from all participants included in this study. Tissue dissociation. For each sample, 100–200 mg of tissue was homogenized by mincing with scissors into smaller pieces (-0.3 mm 2 per piece). Before tissue enzyme digestion, blood and MACS Tissue Storage Solution were washed from the tissue pieces by the addition of 10 ml of ice-cold PBS, followed by further gentle washing and tissue collection using a 70 µm strainer. Then, the tissue pieces were transferred into 5 ml of enzyme mix prepared by combining 4.7 ml of DMEM, 200 µl of Enzyme H, 100 µl of enzyme R and 25 µl of Enzyme A. The samples were incubated at 37°C for 30–60 min in an Eppendorf thermomixer at 700 rpm. After incubation, the cell suspension was filtered through a prewetted 70 µm MACS SmartStrainer, placed in 50 ml tubes and centrifuged at 300 g for 5 min at 4°C. The cell pellet was washed with 5 ml chilled 1X PBS with 0.04% BSA and then centrifuged. One milliliter of chilled 1X red blood cell removal solution was added to the cell pellet, and the cells were resuspended by gentle pipetting with a wide-bore pipette tip and then incubated for 10 min on ice. After incubation, 10 ml chilled 1X PBS with 0.04% BSA buffer was immediately added to the sample tube. The tube was gently inverted several times to mix and centrifuged at 4°C at 300 g for 10 min. The cell pellet was resuspended in 1 ml chilled 1X PBS with 0.04% BSA buffer, and the cell concentration and viability were counted using a Countess II automated cell counter. If the percentage of viable cells was < 70%, dead cells were removed to increase the proportion of viable cells. An appropriate volume of buffer was added to the cell suspension and gently mixed to achieve a target cell concentration of 700–1200 cells/µl. Once the target cell concentration was obtained, GEM generation was immediately performed with 10X Genomics reagents. The remaining cells were frozen at -80°C for WES (no more than 5 × 10 6 cells were used). Whole-exome Sequencing. We used a QIAamp DNA Mini Kit (Catalog# 51304) to process DNA extraction. The cell pellet was resuspended in PBS to a final volume of 200 µl, and 20 µl proteinase K was added to each sample and mixed by vortexing. Then, 4 µl of RNase A stock solution (100 mg/ml) was added to the sample and mixed completely, and 200 µl of Buffer AL was added to the sample and mixed by pulse-vortexing for 15 s. The samples were incubated at 56°C for 10 min, and then 200 µl of ethanol was added. The mixture was carefully transferred to a QIAamp Mini Spin column. The cap was closed, and the tube was centrifuged at 6000 x g for 1 min. The QIAamp Mini Column was placed in a clean 2 ml collection tube, and the tube containing the filtrate was discarded. The column was carefully opened, and 500 µl Buffer AW1 was added. The cap was closed, and the tube was centrifuged at 6000 x g for 1 min. This step was repeated with buffer AW2 and eluted with 100 µl buffer AE. One microliter of DNA sample was taken for quantification on a 1% agarose gel. After passing the QC, the SureSelect XT Library Prep Kit (Catalog#: 5500 − 0132) was used for library construction and QC according to the official protocol. Finally, we used an Illumina Nova 6000 for 150 bp paired end (PE 150) sequencing. Bulk RNA sequencing. We used the RNeasy MinElute Cleanup Kit (Catalog: 74204) to extract RNA. The cell pellet was resuspended in TRIzol reagent to a final volume of 1000 µl. Then, the samples were vortexed, and 200 µl chloroform was added to mix the samples. The lysate was centrifuged for 15 min at 12000 x g at 4°C. The supernatant was pipetted, and one volume of RLT buffer was added and mixed well. Then, one volume of ethanol was added and mixed well. The samples were transferred to a RNeasy MinElute spin column and centrifuged for 15 s at 8000 x g. The flowthrough was discarded, and the same steps were repeated with buffer RW1. Then, 80 µl DNase I solution was added to the samples and incubated for 8 min. After incubation, we added 350 µl of RW1 to the samples and centrifuged them at 8000 x g for 15 s. The flowthrough was discarded, and the wash steps were repeated with 500 µl of RPE buffer and 500 µl of freshly prepared 80% ethanol. Finally, the RNA samples were eluted in 30 µl RNase-free water, and 1 µl of the RNA sample was taken for quantification and evaluation of RNA integrity. After passing the QC, we used the TruSeq RNA Library Prep for Enrichment Kit (Catalog#: 20394426) for library construction according to the official protocol. Finally, we used an Illumina Nova 6000 for PE 150 sequencing. Chromium 10x Genomics library and sequencing. For 10X Chromium Single Cell Immune Profiling, target tissue was prepared into single-cell suspensions by using commercial enzymes from Miltenyi. Then, the cell concentration and viability were determined using a Countess II Automated Cell Counter and recorded. The cells were diluted to 700–1200 cells/µl. Cell loading, GEM generation and library construction of T cells and B cells were performed based on the Chromium Next GEM Single Cell 5' Kit official protocol. After library construction, we used an Agilent 2100 Bioanalyzer to perform a library check. Finally, we used an Illumina Nova 6000 for PE 150 sequencing. Additional sequencing was performed to obtain coverage of at least 50K reads per cell. WES analysis. Raw reads may contain adapter sequences and low-quality bases, which may affect the subsequent analysis. Therefore, it is necessary to review the raw sequencing reads. In this project, read quality was calculated for all samples by using FastQC (version 0.11.2) 13 software with the default parameters. We used the BWA-MEM (Burrows–Wheeler Aligner-maximal exact match) 14 algorithm to compare the sequencing data to the reference genome (GRCh38), obtain the bam file, mark the PCR reads in the bam file, rearrange the regions that may have indel mutations, and recalibrate the quality of each base pair. Gene variants refer to changes in the sequence of the genome, including point mutations caused by single base changes and deletions, duplications, and insertions of multiple bases. In this analysis, the Sentieon (201611.02) 15 workflow was carried out to detect SNVs (single nucleotide variants) and indels (small insertions and deletions). The TMB result was the sum of nonsynonymous mutations and indel numbers in the exon region divided by the length of the total exon. Single-cell RNA sequencing data alignment. The 10X Genomics single-cell raw sequencing paired-end reads were processed using CellRanger software version 5.0.1 and the 10x human genome GRCh38 release as the reference 16 . The output-filtered gene expression matrices were analyzed by R (version 4.0.3) with the Seurat package (version 4.0.0). Low-quality cells with (i) 25% UMIs derived from the mitochondrial genome, or (iii) > 5000 or 0.1%) were removed. Downstream analyses, including normalization, scaling, clustering of cells, and identifying cluster marker genes, were performed by using the function in Seurat. The gene expression matrices were normalized by the NormalizeData function, then the FindVariableFeatures function was used to select 2000 features with high cell-to-cell variation. To reduce the dimensionality of the datasets, the RunPCA function was conducted with default parameters, including the selected 2000 highly variable genes on linear transformation-scaled data generated by the ScaleData function. Finally, we clustered cells using the suitable PCs to calculate the Uniform Manifold Approximation and Projection (UMAP) for dimension reduction with the FindNeighbors and FindClusters functions. Nonlinear dimensional reduction was performed with the RunUMAP function with default settings. Cell type annotation and cluster marker identification. Cells could be clustered together according to common features after nonlinear dimensional reduction. Then, all cells were projected into two-dimensional space by UMAP. We used the FindAllMarkers function in Seurat to find the marker genes with significant differences in expression for each cluster. All clusters were classified and annotated on the basis of these significantly different canonical markers. Clusters that expressed no canonical cell type markers were classified as “Others”. Subclustering of major cell types. The subset function in Seurat was used to extract the subcluster based on the expression of canonical markers of particular cell types. Next, these major cell types were integrated for further subclustering. After integration, the genes were scaled to unit variance. Scaling, PCA, and clustering were performed as described above. Cell number statistics. R was used to calculate the numbers of cells of different cell types. The ggplot2 library was used to generate relevant graphical outputs. The t test was used to test the significance between different groups. Identification of malignant epithelial cells. The R code inferCNV ( https://github.com/broadinstitute/infercnv ) was used to infer the CNVs in epithelial cells from the scRNA-seq data. The epithelial cells (EPCAM+) of the normal sample were considered putative nonmalignant cells, and their CNV estimates were used to define a baseline. The epithelial cells from tumor samples were compared to the normal reference baseline. The CNV R-scores in different regions are readily apparent and are overabundant or less abundant in the regions of the tumor genome compared to that of putative nonmalignant cells. The calculated CNV signal was defined as the mean square of the CNV estimates across all genomic locations. The calculated CNV R-scores were defined as the Pearson correlation coefficient between each cell’s CNV pattern and the average CNV pattern of the top 5% of cells from the same tumor with respect to CNV signal 17 . Epithelial cells (EPCAM+) with CNV R-scores ≥ 0.4 and CNV signal ≥ 0.03 were defined as malignant cells, while those with CNV R-scores ≤ 0.4 and CNV signal ≤ 0.03 were defined as nonmalignant cells. Cells that did not meet either of the above two conditions were defined as unresolved cells. CellPhoneDB. CellPhoneDB was used with default parameters to infer the interactions between the different kinds of cells on a linear transformation scaled data matrix. The method infers potential receptor–ligand interactions based on the expression of a receptor from one cell type and a ligand from another 18 . Only the receptors and ligands that were expressed in more than 30% of the cells in the specific cluster were considered. Network visualization was performed by using Cytoscape (version 3.8.2). The network layout was set to a force-directed layout. Trajectory analysis. Trajectory analysis was performed using Monocle version 3.0.2.3 19 . We ordered the subclusters of cells, such as T cells or epithelial cells, onto a pseudotime trajectory based on the union of highly variable genes obtained from the cells of the subcluster, which were from normal, MIA and IAC samples. All parameters of Monocle3 were set to the defaults (variable genes with the following cutoff criteria: dispersion_empirical > dispersion_fit; and mean expression > 0.001). Dimensional reduction and cell ordering were performed using the DDRTree method and the orderCells function. Data processing of single-cell VDJ libraries. TCR and BCR reads were aligned to the GRCh38 reference genome, and consensus VDJ annotation was performed using Cellranger vdj (version 5.0.1) with default parameters. The R library Immunarch (version 0.6.7) was used to calculate the distribution of the V gene clonotype and calculate statistics for some related indicators, including overlap and correlation. Statistical analyses. The statistical tools, methods, and thresholds used for each analysis are explicitly described within the Results or detailed in the figure legends or Materials and Methods. Results 1. Global cellular landscape in multiple primary lung cancer. To elucidate the cellular development in multiple primary lung cancer (MPLC), we analyzed a total of 167397 cells from 23 sampling sites in 6 patients with multiple foci of early-stage lung cancer from the National Cancer Center/National Clinical Research Center for Cancer/Cancer Hospital (Fig. 1 a-b, Supplementary Fig. 1a, d and Supplementary Table 2). Among these cells, 62714 (23.2%) were from AIS/MIA, 65844 (37.5%) were from IAC, and 38839 (39.3%) were from Normal Lung. To increase the accuracy of cell type designation, we jointly applied principal component analysis (PCA) after correction for read depth and mitochondrial read counts. Through marker-based annotation, 8 major cell types were identified in these unbatched and comparable datasets by graph-based UMAP (Fig. 1 b, Supplementary Fig. 1a, 1d): epithelial cells, alveolar cells, fibroblasts, endothelial cells, myeloid cells, T lymphocytes cells, B lymphocytes cells and mast cells. Global cell type annotations for MPLC patients were categorized according to the molecular subtype of lung tumor tissues. The marker-based annotation correlated well with the cell types obtained via Cell Marker and SingleR v.1.0 20,21 (Fig. 1 b and Supplementary Fig. 1c). T cells and epithelial cells show great differences over the course of development of MPLC. The relative abundance of B cells increased, while that of T cells decreased, from Normal Lung to AIS/MIA and IAC (Fig. 1 d, Supplementary Fig. 1b, and Supplementary 1c). The most frequently mutated genes were EGFR (53%) and RBM10 (35%) (Fig. 1 e and Supplementary Fig. 1g), and the most frequent variant type was missense mutation (Supplementary Fig. 1g). C > T and C > A transversions were enriched in these patients (Supplementary Fig. 1f and Supplementary Fig. 1g). The correlation of the variant allele frequency (VAF) for each individual focus was not significant (Fig. 1 f). The significant mutations and the correlation of the VAF show that the foci in MPLC are independent and do not share an origin. 2. Signatures and metabolic disturbances of malignant cells in MPLC. A total of 13176 normal epithelial cells and 14553 malignant epithelial cells were obtained from Normal Lung, AIS/MIA and IAC samples, and large-scale copy number variations (CNVs) were inferred with Normal Lung cells as references 22 – 24 . Malignant epithelial cells were identified by inferring multipatient copy number variations (CNVs) with normal epithelial cells (mainly immune cells, including T cells and B cells) as references. The CNV patterns showed that malignant cells formed clusters in patients with different origins, indicating a high degree of intertumoral heterogeneity (Fig. 2 a, Supplementary Fig. 2c, 2e). The clusters of malignant cells from different foci were independent and rarely overlapped, while the clusters of normal cells were mixed (Fig. 2 a, Supplementary Fig. 2a, 2c). Alveolar cells are epithelial cells in the lung and play a pivotal role in maintaining the integrity and function of the alveoli. Both alveolar type I cells and alveolar type II cells were positive for EPCAM gene expression 25 , 26 . Other EPCAM + cells, including club, ciliated, basal, goblet, mucous, serous, ionocyte, and neuroendocrine cells, were categorized as epithelial cells. We identified them according to the expression of canonical gene markers (Fig. 2 b). We clustered the normal epithelial cells as alveolar type I cells (AT1; AGER+), alveolar type II cells (AT2; SFTPA1+), ciliated airway epithelial cells (ciliated; TPPP3+), secretory club cells (club; SCGB1A1+) and other cells (no canonical markers to classify) based on the previously described markers (Fig. 2 a, Supplementary Fig. 2a, and Supplementary Fig. 2d) 17 , 27 , independent of patient origin. The relative abundance of malignant epithelial cells, i.e., Ciliated, Club and AT1 cells, was increased in malignant cells compared with Normal Lung cells, and the number of AT2 cells was decreased from Normal Lung to AIS/MIA and IAC. We used lung cancer cells to recapitulate the multilineage differentiation processes of normal epithelia, and subclustering analysis of all epithelial cell datasets revealed divergent differentiation states in MPLC (Fig. 2 d). Intriguingly, the projection of tissue cells and pathological feature classification onto the trajectory assigned the upper right branch and the lower right branch as malignant cells, mainly from IAC tissue. The trajectory indicated that there are two states in the evolution of malignant cells from AIS/MIA to IACs (Fig. 2 c). MSI showed slight changes along the total pseudotime trajectory in MPLC, but TMB did not. In the pseudotime model, this trajectory began from nonmalignant AIS/MIA cells. Then, the trail branched in two directions: some malignant cells became state 1, and some malignant cells became state 2 in IAC. Branch 1 exclusively expressed immune response markers (e.g., CEACAM6, CXCL14, HLA_DPB1, CD74, S100A14, SPINK5, etc.), while branch 2 expressed cell death markers (e.g., PDCD4, AKR1C1, IER3, RGCC, APOH, etc.) (Supplementary Fig. 2g). Functional gene set enrichment analysis showed that the transition from State 1 was involved in the immune response, immune system process and response to transforming growth factor (Fig. 2 e). In addition, the expression of genes implicated in cell death, response to stress and cell cycle increased in State 2 (Fig. 2 f), which shows that some malignant cells will respond to stress from the immune response, some will undergo programmed cell death, and some will enter the cell cycle due to selective pressure during the development of tumors in MPLC. 3. Clonal dynamics and phenotypic transitions of tumor-infiltrating T cells in MPLC. A total of 81,971 T cells were obtained from all samples and clustered into 10 subtypes. Three subtypes of CD4 + T cells were further clustered as Treg cells (TIGIT+, ENTPD1+, FOXP3+), Trm cells (CXCR6+) and CD4 + naïve T cells (CCR7+, IL7R+). Three subtypes of CD8 + T cells were further clustered as activated T cells/exhausted T cells (TNF+, IFNG+, GZMA+, GZMB+, KLRD1+), activated/effector T cells (TNF+, GZMA+, EOMES+), and effector T/Trm cells (ITGAE+, CXCR6+, GZMA+, GZMB+). One subtype of both CD4- and CD8-positive T cells was further clustered as naïve/memory T cells (IL7R+), and 2 subtypes with no CD4 or CD8 expression were further clustered as NK cells/exhausted T cells (HAVCR2+, IFNG+, GZMB+) and cytotoxic T cells/NK cells (XCL1+, XCL2+) (Fig. 3 a, Supplementary Fig. 3a) based on canonical markers (Fig. 3 b, 3 c, Supplementary Fig. 3b), as described previously 28 – 32 . The relative abundances of classified naïve T cells, Treg cells and Trm cells were increased, while those of activated/exhausted T cells were decreased, from Normal Lung to AIS/MIA and IAC (Fig. 3 d). We next analyzed the lineage relationships between phenotypes and clonotypes in T cells. Globally, we found that cells in normal and IAC samples shared more common clonotypes in TCRs than AIS/MIA-group cells (Fig. 3 e, Supplementary Fig. 4a). To investigate whether the development of tumors within clones could inform lineage transitions between T-cell phenotypes, we aggregated all clonotypes and measured the fraction shared in different clusters. There was minimal significant clonotype sharing between T cells in different cell groups. Trm and naïve T cells shared more significant clonotypes with other T cells, and the shared clonotypes increased with the development of infiltration in the tumor tissue, suggesting an increasing transition from Trm and naïve T cells to other activated T-cell states (Supplementary Fig. 4a). The T-cell phenotypes in AIS/MIA samples include more clonotype 1, for which TRB-TRA is “CASSGLAAKPGELF-CAVRRGQNFV”, than other normal and IAC samples (Supplementary Fig. 4b). This sequence may play an important role in antigen presentation during tumor invasion. In normal tissue, there is not enough antigen from tumors to stimulate T cells, and the expression of this antigen presentation factor is low. In IAC tissue, the tumor has experienced stress from T cells in AIS/MIA, and new subclone evolution may not cause the same antigen stimulation. Therefore, clonotype 1, in which TRB-TRA is “CASSGLAAKPGELF-CAVRRGQNFV”, was expressed significantly only in AIS/MIA samples. The TCR clonotype shared more overlap in focus grouped according to the source of the patient compared to those grouped by anatomical pathology info and still had more correlation (Fig. 3 f). 4. Trajectory profiling reveals the branched progression of T cells in MPLC. We next used Monocle3 to analyze the CD8 + and CD4 + T-cell trajectories to further investigate the potential transitions in each of these cell types. The pseudotime trajectory axis derived from Monocle3 indicates that CD8 + naïve and memory T cells transdifferentiate into effector T cells first. Then, some of these cells become exhausted T cells, and others transform into cytotoxic T cells (Fig. 4 a). We compared the expression of genes in the next-stage cells with those in the previous stage (Supplementary Fig. 3c, 3d). Cells on the memory T to effector T branch expressed immune response markers (e.g., TNFRSF14, NFE2L2, TNFRSF1B, IL2RG, etc.) (Supplementary Fig. 3c,3d), while naïve T cells transitioning to effector T cells expressed lymphocyte activation and T-cell differentiation markers (e.g., FOXP1, TGFBR2, TNFRSF4, PRDM1, etc.) (Supplementary Fig. 3c, 3d). Cells on the effector T to exhausted T branch expressed responses to stress and cell death markers (e.g., GNLY, CCL4, and FCGR3A) (Supplementary Fig. 3c, 3d), while effector T cells expressed cytokine production markers (e.g., CXCR4, CRTAM, etc.) (Supplementary Fig. 3c, 3d). Gene Ontology enrichment analysis confirmed that the transition from memory T cells to effector T cells was mainly implicated in the immune response, immune system process and immune response regulation. The transformation from memory to effector in T cells is mainly a process of immune response and regulation (Fig. 4 b (left-upper)). The function of the geneset in the transformation from naïve T to effector T cells was mainly related to lymphocyte activation, T-cell and leukocyte differentiation, T-cell activation and immune system development, which showed that the transformation from naïve to effector T cells was mainly a process of leukocyte and T-cell differentiation (Fig. 4 b (upper right)). Gene Ontology enrichment analysis confirmed that signature genes were enriched for response to stress, positive regulation of cell death and response to stimulus pathway during effector T-cell transformation to exhausted T cells (Fig. 4 b (lower right)). The functional gene set enrichment analysis from effector to cytotoxic T cells showed that the transition was mainly implicated in the regulation of cytokine production, which included the positive regulation of myeloid leukocyte cytokine production involved in the immune response, immune effector processes, macrophage cytokine production and the positive regulation of cytokine production involved in the immune response (Fig. 4 b (lower left)). The pseudotime trajectory axis indicates that CD4 + naïve and Treg cells have four final differentiation states in themselves, which are labeled as Naïve_C1 to Naïve_C4 and Treg_C1 to Treg_C4 (Fig. 4 c, Supplementary Fig. 5a). Pathway analyses showed that Naive_C2 cells were enriched in some effector pathways [Hallmark pi3k_akt_mtor signaling, chemokine signaling pathway, and natural killer cell-mediated cytotoxicity], while Naive_C4 cells did not express some pathways [T-cell receptor signaling pathway, regulation of lymphocyte differentiation, regulation of T-cell differentiation, T-cell differentiation, T-cell proliferation, and T-cell activation]. The function of the T cells in the Naïve_C2 cluster is closer to that of effector or exhausted T cells, while the function of the Naïve_C4 cluster is more closely related to that of other naïve clusters (supplemental Fig. 5b (Left)). The GSVA results showed that Treg_C3 cells were enriched in the regulation of myeloid cell differentiation [myeloid leukocyte differentiation, myeloid cell differentiation, and regulation of myeloid cell differentiation], while Treg_C2 cells mainly participated in the regulation of cell differentiation, activation of the T-cell receptor signaling pathway [regulation of cell differentiation, T-cell receptor signaling pathway, T-cell receptor signaling pathway and T-cell activation]. Treg_C2 cells also had high expression levels of genes related to lymphocyte differentiation, including lymphocyte differentiation, regulation of T-cell differentiation, and regulation of lymphocyte differentiation. Treg_C1 cells mainly participates in the regulation of stem cell differentiation [regulation of stem cell differentiation, stem cell differentiation, and hematopoietic progenitor cell differentiation] (Supplementary Fig. 5b (Right)). 5. B cells have two functional distribution subtypes in MPLC. A total of 4279 B cells were analyzed and classified into five subtypes (Fig. 5 a, Supplementary Fig. 6a, b, supplementary Table 3): Bcells_C1 (MS4A1+, CD19+, CD74+, CD37+, BTG+) (49.7% in IAC, 10.3% in AIS/MIA, 6.4% in Normal); Bcells_C2 (NEIL1+, RGS13+) (0.9% in IAC, 0.4 in AIS/MIA, 0.0% in Normal); Bcells_C3 (CYTOR+, NUCB2+, ANKRD28+, TRIB1+) (7.8% in IAC, 2.9% in AIS/MIA, 1.9% in Normal); Bcells_C4 (6.2% in IAC, 3.8% in AIS/MIA, 0.9% in Normal); and Bcells_C5 (IL7R+, IL32+, CST7+, FYB1+) (5.9% in IAC, 1.5% in AIS/MIA, 1.4% in Normal) (Supplementary Fig. 6e, Supplementary Fig. 7a). The relative abundance of the B-cell subtype showed no uniform pattern in different patient foci and or when grouped by anatomical pathology (Fig. 5 c, Supplementary Fig. 6c, d). The pathway analysis of B cells in MPLC further confirmed that there are two kinds of B cells with different functions. The regulation and activation of the innate immune response was highly activated in the Bcell-C4 and Bcell-C5 clusters (Fig. 5 d, red pathway). However, the Bcell-C1, Bcell-C2, and Bcell-C3 clusters were enriched in pathways related to antigen binding, processing and presentation, including Reactome MHC class II Antigen Presentation (Fig. 5 d, green pathway). Some B cells are involved in immune activation and regulation of the immune response, while others are involved in antigen presentation. We next used Monocle3 to analyze the trajectory performance of B cells to further investigate the potential transition in each of their cell types. The pseudotime trajectory axis derived from Monocle3 indicates that Bcells_C1 (MS4A1+, CD19+, CD74+, CD37+, BTG+) have two differentiation states; one is transdifferentiating into Bcell_C3 (NEIL1+, RGS13+) and Bcell_C5 (IL7R+, IL32+, CST7+, FYB1+), and the other is transitioning to Bcell_C4 (Fig. 5 b, Supplementary Fig. 7b, 7c, 7d). Combined with the pathway enrichment results, these findings show that B cells mainly differentiated into two different types of cells with different functions in MPLC; one is involved in immune activation and regulation of the immune response, while the other is involved in antigen presentation. Different from the TCR clonotypes, which shared more overlap and correlation in foci grouped according to the source within the patient, the BCR clonotypes shared little overlap or correlation in MPLC (Supplementary Fig. 6f). 6. Myeloid cells have two functional distribution subtypes in MPLC. Subclustering of 46836 myeloid cells revealed 10 subsets (Fig. 5 e, Supplementary Fig. 9a.): M1 macrophages (MSR1+, MARCO+, etc.), M2 macrophages (MERTK+, MMP9+, etc.), MDSCs (CEACAM8+, ITGAM+, etc.), monocyte cells (TNFSF10+, FCER2+, CLEC10A+, etc.), dendritic cells (FCER1A+, CD1C+, CLEC4A+, etc.), myeloid dendritic cells (CCR7+, IL15+, CD1C+, etc.), plasma cells (IRF4+, CD36+, CXCR4+, etc.), secretory cells (CXCL10+) (Fig. 5 g, Supplementary Fig. 9d). The percentage of M2 macrophages increased from normal to AIS/MIA to IAC, while the other subtypes were not significantly different in MIA/AIS and IAC (Fig. 5 f, Supplementary Fig. 9b). Macrophages play an important role in maintaining organismal integrity by either directly participating in pathogen elimination or repairing tissue under sterile inflammatory conditions 33 , and Microphage_2 cells support tissue integrity by providing growth factors and healing capacity 34 . From a functional perspective, myeloid cells are divided into two main types: one exhibits high expression of IL1A, FCGR3A and IFN-γ (M1 macrophages), which can positively regulate the innate immune response. The other negatively regulates the adaptive immune response and regulation of stem cell differentiation (M2 macrophages). The increase in M2 macrophages will lead to malignant cells adapting to the immune response and escaping the response from the immune system (Supplementary Fig. 9c). The pseudotime trajectory axis derived from Monocle3 indicates that myeloid cells have four main branch trajectories (Fig. 5 h). M2 macrophages were mainly concentrated in the left branch trajectory, M1 macrophages were mainly in the opposite branch trajectory, monocyte cells were in the middle to left stage, secretory cells were mainly concentrated in the lower left branch trajectory, and other Myeloid2_cells were mainly in the lower right branch trajectory (Fig. 5 h). Next, we evaluated the expression of DEGs in the M1 macrophage (MSR1+; MARCO+) and M2 macrophage (MERTK+; MMP9+) populations in the whole trajectory. The expression of macrophage colony-stimulating factor genes, such as CSF1, MCEMP1, APOE, CD52 and FABP4, increased gradually during the process of transdifferentiating into M1 macrophages, and the expression of genes related to cytomembrane transport, namely, GPR183, CD48, CD93, and CCL5, increased during the process of transdifferentiating into M2 macrophages (Supplementary Fig. 9f). 7. Positive regulation of endothelial cell development is extensive in IAC. We found four distinct subtypes by reclustering 2081 endothelial cells (Fig. 6 a, Supplementary Fig. 8a): Endothelial_C1 (VIPR1+, IL1RL1+, TMEM100+, TNFSF10+, etc.), Endothelial_C2 (EMP1+, CD93+, PLVAP+, etc.), Endothelial_C3 (TNFAIP3+, CXCR4+, CREM+, etc.), and Endothelial_C4 (TFF3+, TIMP2+, STMN1+, etc.) (Fig. 6 c, Supplementary Fig. 8d). The percentage of Endothelial_C1 and Endothelial_C2 in MIA/AIS and IAC were more different than those in normal samples, and the distribution of endothelial subtypes was not significantly different in MIA/AIS and IAC (Fig. 6 b, Supplemental Fig. 8c). Then, we focused on the transcriptomic features and functions of Endothelial_C2 cells. The genes EMP1, TSPAN7, VWF, PLVAP and CD93 were expressed at high levels in these cells, and pathway analyses showed that Endothelial_C2 cells were enriched in positive regulation of endothelial cell development pathways [positive regulation of blood vessel diameter, regulation of endothelial cell differentiation pathway, and abnormal morphology of the great vessels] (supplemental Fig. 8d, Fig. 6 d). Thus, Endothelial_C2 cells are related to the regulation of endothelial development, and their prevalence increases with the development of tumor invasion. 8. Fibroblast subtype distributions are similar in IAC, AIS/MIA and normal samples. Fibroblast cells maintain a very stable state in MPLC (Fig. 6 c, Supplementary Fig. 8c). A total of 2523 fibroblast cells were analyzed, and five subtypes were identified (Fig. 6 c). In MPLC, there was no significant difference in the number of different fibroblast clusters, and there were also no significant changes in the pathway and function of genes expressed in different fibroblast clusters (Fig. 6 f, 6 g, Supplementary Fig. 8f). Fibroblast_C1 cells expressed high levels of genes associated with the Reactome elastic fiber formation annotation, such as COL3A1, COL1A1, COL6A3, COL10A1, and COL8A1, whereas Fibroblast_C5 cells highly expressed genes related to positive regulation of fibroblast proliferation, like RASD1, SFTPB and FOS (Fig. 6 g, Supplementary Fig. 8g). 9. Characterization of cell–cell interactions involved in MPLC. There are many specific signaling network changes in MPLC, which depend on complex interactions between structural and inflammatory cells 35 . One cell expresses a receptor or ligand, and this “ligand–receptor” pair will be defined as an incoming or outgoing interaction for this cell. We used CellPhoneDB to identify potential cell–cell interactions and define their changes in primary lung cancer. Compared with the MIA samples, the cells in tumors exhibited significantly increased interactions (Fig. 7 a). Alveolar and myeloid cells had the most outgoing interactions of all cell types in both MIA and IAC samples. In comparison with MIA, the interaction between T cells and alveolar cells was significantly increased (Fig. 7 b). Myeloid, cytotoxic T, natural killer and B cells are involved in the action of recognizing and presenting the antigens of epithelial cells among the interactions from immune cells to epithelial cells. Analysis of the biological functions revealed that interactions related to antigen recognition and presentation (CD74_COPA, HLA-C_FAM3C, and HLA-DRB1_OGN) were less abundant in IAC than in MIA and normal samples, which means that some escape of immune recognition occurred with deeper invasion in MPLC (Fig. 7 c upper). The interactions related to independent cellular adhesion (NECTIN1_CADM3, FGF1_FGFR4, and CD44_FGFR2) were still less abundant in epithelial cells than in fibroblasts, which means that epithelial cells promoted structural changes in fibroblasts and reduced cell adhesion. The interactions regulating the growth and differentiation of numerous types of cells (ACVR1_BMPR2_BMP6 and ACVR1_BMPR2_BMP5) were more abundant (Fig. 7 c, middle). On the one hand, epithelial cells reduce fibroblast adhesion; on the other hand, they strengthen the differentiation of other mesenchymal cells with deeper invasion in MPLC. The interactions related to endothelial and hematopoietic stem cell differentiation and development (SELL-CD34, SELP-CD34, NPR1-NPPC, and NPR2-NPPC) were more abundant in IAC than in normal and MIA samples (Fig. 7 c, down). Discussion Our knowledge about the dynamics of gene expression and metabolic pathways during cellular changes and evolution in MPLC remains incomplete. In this study, we comprehensively provide important new insights into these changes during the progression from normal AIS to IAC in MPLC. Recent studies have performed similar development trajectory analyses for CD8 + T cells, with naïve T cells as the root and exhausted T cells and cytotoxic T cells as the two end states 17 , 36 , 37 . However, naïve T cells are not the only root for CD8 + T-cell development and evolution in MPLC. We found that memory T cells are also a main source of these cells and play an important role in CD8 + T-cell differentiation. The immune response and regulation of the immune response were significantly enriched during this differentiation process. The transition from CD8 + memory T cells to effector T cells were mainly a process of positive regulation of immune system activities; however, the progression of CD8 + naïve to effector T cells is mainly a process of leukocyte activation and T-cell differentiation. The progress of effector T-cell transformation to exhausted T cells involves responses to stimuli and positive regulation of cell death. The other pathway, of effector T-cell transformation to cytotoxic T cells, mainly involves the regulation of cytokine production, including the positive regulation of myeloid leukocyte cytokine production involved in the immune response and the positive regulation of cytokine production involved in the immune response. CD4 + T cells have two distinct functional statuses: Treg and naïve. We also found that T cells shared more consistent TCR clonotypes than AIS/MIA-group cells in normal and IAC samples. Unfortunately, further in-depth study of this interesting phenomenon is beyond the scope of this study. In contrast to previous studies on the metabolism of malignant cells 17,38−40 , we used trajectory analysis to study the transition of epithelial cells from normal to malignant. We found that although different nodules in the patients were heterogeneous, they all exhibit two states during the evolution of malignant cells from AIS/MIA to IAC. Specifically, some cells respond to stress from immunity, while others undergo programmed cell death or enter the cell cycle due to selective pressure from the development of tumors in MPLC. The proportion of AT2 cells gradually decreases with the development of LUAD from normal to AIS/MIA to IAC, which also provides evidence that AT2 is the most likely origin of malignant cells in MPLC 32 . The functional distribution subtypes of B cells can be distinguished into two different stages. One is the regulation and activation of innate immune response antigen binding, and the other is processing and presentation, including MHC class II antigen presentation. Myeloid cells are divided into two major types, one with high expression of positive regulators of the innate immune response and the other performing negative regulation of the adaptive immune response and regulation of stem cell differentiation presented by M2 macrophages. With the progression of LUAD from normal to AIS/MIA to IAC, the increasing proportion of M2 macrophages causes malignant cells adapt to the immune response and escape the immune response 41 – 43 . We also found that some endothelial cell clusters also participate in the positive regulation of blood vessel development in early glandular neoplasia of the lung. Moreover, the proportion of these endothelial cells gradually increases with the progression of LUAD from normal to AIS/MIA to IAC 44 – 46 . Comprehensive analysis of the cell–cell interactions in MPLC identified that interactions related to antigen recognition and presentation were less abundant in IAC than in AIS/MIA and normal samples, and the interactions related to independent cellular adhesion were still less abundant in epithelial cells than in fibroblasts. The interactions regulating the growth and differentiation of numerous types of cells were more abundant, and the interactions related to endothelial and hematopoietic stem cell differentiation and development were more abundant in the IAC than in the normal and AIS/MIA samples. Despite the striking findings identified in this study, there are still several limitations, including sample collection, technical challenges and software limitations. We were unable to evaluate the continuous development of different nodules in the same patient, and the cells used for single-cell sequencing were from just a part of the tumor tissue, so there might be differences in biological complexity between our samples and the tissues used for pathological diagnosis and real cells in developing tumors. The 10x Genomics single-cell technology cannot cover all the transcripts in cells. There are also still limitations in batch effect removal, clustering and trajectory inference analysis, as these analyses are still maturing, and the complexity of the underlying topology could be underestimated 47 . Declarations Data availability The raw data that support the findings of this study are available from the corresponding author upon reasonable request. Acknowledgments We appreciate the support and participation of the physicians and patients in this study. Competing interests Hao Xia is employee of Mingma Technologies Co., Ltd. All remaining authors have declared no conflicts of interest. Author contributions W.G. and S.G. conceived the study. J.H., J.Y., F.T., Q.X., and S.G. provided project management and supervision. W.G. and L.G. provided or facilitated the accrual of patient samples, pathology, and/or clinical data. W.G. and H.X. performed bioinformatics and genomic analyses. W.G. and H.X. performed statistical analyses. W.G. and Z.B. wrote the original draft, with input from all authors. Funding This study was supported by the National Key R&D Program of China (2021YFC2500900), the National Natural Science Foundation of China (82002451), the CAMS Initiative for Innovative Medicine (2021-1-I2M-015), the Special Research Fund for Central Universities, Peking Union Medical College (3332020024), the Non-profit Central Research Institute Fund of Chinese Academy of Medical Sciences (2018PT32033), and the Beijing Hope Run Special Fund of Cancer Foundation of China (LC2019B15). References Chen, W. et al. Cancer statistics in China, 2015. CA: a cancer journal for clinicians 66, 115–132 (2016). Siegel, R. L., Miller, K. D., Fuchs, H. E. & Jemal, A. Cancer statistics, 2021. CA: a cancer journal for clinicians 71, 7–33 (2021). Sung, H. et al. Global cancer statistics 2020: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA: a cancer journal for clinicians 71, 209–249 (2021). Shintani, Y. et al. Clinical features and outcomes of patients with stage I multiple primary lung cancers. Cancer Sci 112, 1924–1935, doi: 10.1111/cas.14748 (2021). Liu, Y. et al. Genomic heterogeneity of multiple synchronous lung cancer. Nat Commun 7, 13200, doi: 10.1038/ncomms13200 (2016). Port, J. L. et al. Surgical resection for multifocal (T4) non-small cell lung cancer: is the T4 designation valid? The Annals of thoracic surgery 83, 397–400 (2007). Asamura, H. Multiple primary cancers or multiple metastases, that is the question. J Thorac Oncol 5, 930–931, doi: 10.1097/JTO.0b013e3181e15d8c (2010). Meacham, C. E. & Morrison, S. J. Tumour heterogeneity and cancer cell plasticity. Nature 501, 328–337 (2013). Ma, P. et al. Simultaneous evolutionary expansion and constraint of genomic heterogeneity in multifocal lung cancer. Nature communications 8, 1–11 (2017). Liu, Y. et al. Genomic heterogeneity of multiple synchronous lung cancer. Nature communications 7, 1–8 (2016). Murphy, S. J. et al. Identification of independent primary tumors and intrapulmonary metastases using DNA rearrangements in non–small-cell lung cancer. Journal of Clinical Oncology 32, 4050 (2014). Giraldo, N. A. et al. The clinical role of the TME in solid cancer. British journal of cancer 120, 45–53 (2019). S. A., S. B. a. FastQC: a quality control tool for high throughput sequence data. (2010). Li, H. & Durbin, R. Fast and accurate long-read alignment with Burrows-Wheeler transform. Bioinformatics 26, 589–595, doi: 10.1093/bioinformatics/btp698 (2010). Weber JA, A. R., Gallagher BD, Edwards JS. Sentieon DNA pipeline for variant detection - Software-only solution, over 20× faster than GATK 3.3 with identical results. PeerJ PrePrints 4:e1672v2, doi: https://doi.org/10.7287/peerj.preprints.1672v2 (2016). https://github.com/10XGenomics/cellranger. Cellranger. Xing, X. et al. Decoding the multicellular ecosystem of lung adenocarcinoma manifested as pulmonary subsolid nodules by single-cell RNA sequencing. Sci Adv 7, doi: 10.1126/sciadv.abd9738 (2021). Efremova, M., Vento-Tormo, M., Teichmann, S. A. & Vento-Tormo, R. CellPhoneDB: inferring cell-cell communication from combined expression of multi-subunit ligand-receptor complexes. Nat Protoc 15, 1484–1506, doi: 10.1038/s41596-020-0292-x (2020). Cao, J. et al. The single-cell transcriptional landscape of mammalian organogenesis. Nature 566, 496–502, doi: 10.1038/s41586-019-0969-x (2019). Aran, D. et al. Reference-based analysis of lung single-cell sequencing reveals a transitional profibrotic macrophage. Nat Immunol 20, 163–172, doi: 10.1038/s41590-018-0276-y (2019). Zhang, X. et al. CellMarker: a manually curated resource of cell markers in human and mouse. Nucleic Acids Res 47, D721-D728, doi: 10.1093/nar/gky900 (2019). Patel, A. P. et al. Single-cell RNA-seq highlights intratumoral heterogeneity in primary glioblastoma. Science 344, 1396–1401, doi: 10.1126/science.1254257 (2014). Puram, S. V. et al. Single-Cell Transcriptomic Analysis of Primary and Metastatic Tumor Ecosystems in Head and Neck Cancer. Cell 171, 1611–1624 e1624, doi: 10.1016/j.cell.2017.10.044 (2017). Timothy Tickle, I. T., Christophe Georgescu,Maxwell Brown,Brian Haas. Visualizing Large-scale Copy Number Variation in Single-Cell RNA-Seq Expression Data. (2021). Khan, P. et al. Culture of human alveolar epithelial type II cells by sprouting. Respir Res 19, 204, doi: 10.1186/s12931-018-0906-9 (2018). Travaglini, K. J. et al. A molecular cell atlas of the human lung from single-cell RNA sequencing. Nature 587, 619–625, doi: 10.1038/s41586-020-2922-4 (2020). Deprez, M. et al. A Single-Cell Atlas of the Human Healthy Airways. Am J Respir Crit Care Med 202, 1636–1645, doi: 10.1164/rccm.201911-2199OC (2020). Morgan, A. J. et al. Expression of CXCR6 and its ligand CXCL16 in the lung in health and disease. Clin Exp Allergy 35, 1572–1580, doi: 10.1111/j.1365-2222.2005.02383.x (2005). Wein, A. N. et al. CXCR6 regulates localization of tissue-resident memory CD8 T cells to the airways. J Exp Med 216, 2748–2762, doi: 10.1084/jem.20181308 (2019). Ghorani, E. et al. The T cell differentiation landscape is shaped by tumour mutations in lung cancer. Nat Cancer 1, 546–561, doi: 10.1038/s43018-020-0066-y (2020). Delacher, M. et al. Single-cell chromatin accessibility landscape identifies tissue repair program in human regulatory T cells. Immunity 54, 702–720 e717, doi: 10.1016/j.immuni.2021.03.007 (2021). Wang, Z. et al. Deciphering cell lineage specification of human lung adenocarcinoma with single-cell RNA sequencing. Nat Commun 12, 6500, doi: 10.1038/s41467-021-26770-2 (2021). Chavez-Galan, L., Olleros, M. L., Vesin, D. & Garcia, I. Much More than M1 and M2 Macrophages, There are also CD169(+) and TCR(+) Macrophages. Front Immunol 6, 263, doi: 10.3389/fimmu.2015.00263 (2015). Ley, K. M1 Means Kill; M2 Means Heal. J Immunol 199, 2191–2193, doi: 10.4049/jimmunol.1701135 (2017). Holgate, S. T. et al. Asthma. Nat Rev Dis Primers 1, 15025, doi: 10.1038/nrdp.2015.25 (2015). Guo, X. et al. Global characterization of T cells in non-small-cell lung cancer by single-cell sequencing. Nat Med 24, 978–985, doi: 10.1038/s41591-018-0045-3 (2018). Qian, J. et al. A pan-cancer blueprint of the heterogeneous tumor microenvironment revealed by single-cell profiling. Cell Res 30, 745–762, doi: 10.1038/s41422-020-0355-0 (2020). Goveia, J. et al. An Integrated Gene Expression Landscape Profiling Approach to Identify Lung Tumor Endothelial Cell Heterogeneity and Angiogenic Candidates. Cancer Cell 37, 421, doi: 10.1016/j.ccell.2020.03.002 (2020). Ireland, A. S. et al. MYC Drives Temporal Evolution of Small Cell Lung Cancer Subtypes by Reprogramming Neuroendocrine Fate. Cancer Cell 38, 60–78 e12, doi: 10.1016/j.ccell.2020.05.001 (2020). He, Y. et al. Mechanisms of Progression and Heterogeneity in Multiple Nodules of Lung Adenocarcinoma. Small Methods 5, e2100082, doi: 10.1002/smtd.202100082 (2021). Heusinkveld, M. & van der Burg, S. H. Identification and manipulation of tumor associated macrophages in human cancers. J Transl Med 9, 216, doi: 10.1186/1479-5876-9-216 (2011). Conway, E. M. et al. Macrophages, Inflammation, and Lung Cancer. Am J Respir Crit Care Med 193, 116–130, doi: 10.1164/rccm.201508-1545CI (2016). Liu, M. et al. Transcription factor c-Maf is a checkpoint that programs macrophages in lung cancer. J Clin Invest 130, 2081–2096, doi: 10.1172/JCI131335 (2020). Lamalice, L., Le Boeuf, F. & Huot, J. Endothelial cell migration during angiogenesis. Circ Res 100, 782–794, doi: 10.1161/01.RES.0000259593.07661.1e (2007). Betsholtz, C. Cell-cell signaling in blood vessel development and function. EMBO Mol Med 10, doi: 10.15252/emmm.201708610 (2018). Sweeney, M. & Foldes, G. It Takes Two: Endothelial-Perivascular Cell Cross-Talk in Vascular Development and Disease. Front Cardiovasc Med 5, 154, doi: 10.3389/fcvm.2018.00154 (2018). Lahnemann, D. et al. Eleven grand challenges in single-cell data science. Genome Biol 21, 31, doi: 10.1186/s13059-020-1926-6 (2020). Additional Declarations Competing interest reported. Hao Xia is employee of Mingma Technologies Co., Ltd. All remaining authors have declared no conflicts of interest. Supplementary Files Patientsimagingandpathologicalinformation.pdf Supplementalfigures.pdf Supplementarytable1.xlsx Supplementarytable2.xls Supplementarytable3.xls Cite Share Download PDF Status: Posted Version 1 posted You are reading this latest preprint version Research Square lets you share your work early, gain feedback from the community, and start making changes to your manuscript prior to peer review in a journal. As a division of Research Square Company, we’re committed to making research communication faster, fairer, and more useful. We do this by developing innovative software and high quality services for the global research community. Our growing team is made up of researchers and industry professionals working together to solve the most critical problems facing scientific publishing. Also discoverable on Platform About Our Team In Review Editorial Policies Advisory Board Help Center Resources Author Services Accessibility API Access RSS feed Manage Cookie Preferences © Research Square 2026 | ISSN 2693-5015 (online) Privacy Policy Terms of Service Do Not Sell My Personal Information {"props":{"pageProps":{"initialData":{"identity":"rs-1716190","acceptedTermsAndConditions":true,"allowDirectSubmit":true,"archivedVersions":[],"articleType":"Research Article","associatedPublications":[],"authors":[{"id":110513962,"identity":"d4b81cdf-cd90-40de-8d1c-7c1b81494f16","order_by":0,"name":"Wei Guo","email":"","orcid":"","institution":"National Cancer Center, Chinese Academy of Medical Sciences and Peking Union Medical College","correspondingAuthor":false,"prefix":"","firstName":"Wei","middleName":"","lastName":"Guo","suffix":""},{"id":110513964,"identity":"b0d90dd8-2754-4088-b42c-5746fa43c707","order_by":1,"name":"Hao Xia","email":"","orcid":"","institution":"Mingma Technologies Co., Ltd","correspondingAuthor":false,"prefix":"","firstName":"Hao","middleName":"","lastName":"Xia","suffix":""},{"id":110513965,"identity":"4cd182a2-70d8-4653-a68b-227a8eaab663","order_by":2,"name":"Qilin Huai","email":"","orcid":"","institution":"National Cancer Center, Chinese Academy of Medical Sciences and Peking Union Medical College","correspondingAuthor":false,"prefix":"","firstName":"Qilin","middleName":"","lastName":"Huai","suffix":""},{"id":110513966,"identity":"8392cc8c-2ce8-4667-ad59-0f1421a28fc6","order_by":3,"name":"Xuemin Xue","email":"","orcid":"","institution":"National Cancer Center, Chinese Academy of Medical Sciences and Peking Union Medical College","correspondingAuthor":false,"prefix":"","firstName":"Xuemin","middleName":"","lastName":"Xue","suffix":""},{"id":110513967,"identity":"fdb1be76-cfdc-4c4f-85bb-50706260977a","order_by":4,"name":"Lei Guo","email":"","orcid":"","institution":"National Cancer Center, Chinese Academy of Medical Sciences and Peking Union Medical College","correspondingAuthor":false,"prefix":"","firstName":"Lei","middleName":"","lastName":"Guo","suffix":""},{"id":110513968,"identity":"cb392e8f-77ef-4f7f-857a-5f7119037081","order_by":5,"name":"Fengwei Tan","email":"","orcid":"","institution":"National Cancer Center, Chinese Academy of Medical Sciences and Peking Union Medical College","correspondingAuthor":false,"prefix":"","firstName":"Fengwei","middleName":"","lastName":"Tan","suffix":""},{"id":110513969,"identity":"d2149d38-b480-4d05-ace0-9429ff8f773e","order_by":6,"name":"Jun Zhao","email":"","orcid":"","institution":"National Cancer Center, Chinese Academy of Medical Sciences and Peking Union Medical College","correspondingAuthor":false,"prefix":"","firstName":"Jun","middleName":"","lastName":"Zhao","suffix":""},{"id":110513970,"identity":"1b749b05-ede5-4c1d-a190-b303b69d121b","order_by":7,"name":"Qi Xue","email":"","orcid":"","institution":"National Cancer Center, Chinese Academy of Medical Sciences and Peking Union Medical College","correspondingAuthor":false,"prefix":"","firstName":"Qi","middleName":"","lastName":"Xue","suffix":""},{"id":110513971,"identity":"3dce5855-dd4d-4302-b37e-c4a2139d6671","order_by":8,"name":"Shugeng Gao","email":"data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAZAAAAAyAQMAAABI0h/eAAAABlBMVEX///8AAABVwtN+AAAACXBIWXMAAA7EAAAOxAGVKw4bAAAA/klEQVRIiWNgGAWjYJACZih1gCGBRC1sCSRr4TEgTrm8++HDnwvbauXM+dd8/PCgxi5avoH54QeGmjs4tRieSUuTntl23NhyxtvNEgnHknM3HGAzlmA49gy3loYcM2betmOJG26c3caQ2MCcu4GBwYyBseEwbi39b4w/A7XUb7hx5hlQS33u/Ab2b3i1yEvkGEjzttUkGJzvYQNqOZzbcIAHvy0GEs/SpHnOHTDccAPohYRjx3M3HOYpBjLw2NKffPgzT1mdvMH5ww8//qipzp3f3r7xw4caPLYcAFNABRIJUCFQNCVgVw2xpQFM1TEw8B/Ao2wUjIJRMApGNAAA4jtZ9TMhFc4AAAAASUVORK5CYII=","orcid":"","institution":"National Cancer Center, Chinese Academy of Medical Sciences and Peking Union Medical College","correspondingAuthor":true,"prefix":"","firstName":"Shugeng","middleName":"","lastName":"Gao","suffix":""},{"id":110513972,"identity":"0e9d72c8-57d6-4d4d-bc64-1040fb096c69","order_by":9,"name":"Jie He","email":"","orcid":"","institution":"National Cancer Center, Chinese Academy of Medical Sciences and Peking Union Medical College","correspondingAuthor":false,"prefix":"","firstName":"Jie","middleName":"","lastName":"He","suffix":""},{"id":110513973,"identity":"4bb65f26-777b-4d4a-a4d5-e636f24b16cc","order_by":10,"name":"Bolun Zhou","email":"","orcid":"","institution":"National Cancer Center, Chinese Academy of Medical Sciences and Peking Union Medical College","correspondingAuthor":false,"prefix":"","firstName":"Bolun","middleName":"","lastName":"Zhou","suffix":""}],"badges":[],"createdAt":"2022-06-01 14:14:17","currentVersionCode":1,"declarations":"","doi":"10.21203/rs.3.rs-1716190/v1","doiUrl":"https://doi.org/10.21203/rs.3.rs-1716190/v1","draftVersion":[],"editorialEvents":[],"editorialNote":"","failedWorkflow":false,"files":[{"id":22243773,"identity":"dd51278e-2041-4952-9bf0-d925d55b552d","added_by":"auto","created_at":"2022-06-03 22:30:37","extension":"jpg","order_by":1,"title":"Figure 1","display":"","copyAsset":false,"role":"figure","size":694808,"visible":true,"origin":"","legend":"\u003cp\u003eCell type identification in MPLC using a combination of AIS/MIA, IAC and Normal Lung cells. a.Workflow showing the scRNA-seq experimental design and initial data exploration. b.UMAP plot of 167,397 cells in the MPLC dataset. It consists of 12,124 malignant cells, 60032 nonmalignant cells and 9792 unresolved cells. c. Marker genes were used to label clusters by cell identity as represented in the UMAP plot. d. Average proportions of all kinds of cells among all Normal, AIS/MIA, and IAC samples. Proportions of global cell types in MPLC tissues and normal tissues in individual samples (bottom left). Cell number (bottom middle), TMB, MSI, anatomical region and infiltration in each sample (bottom right). e. Mutational landscape of MPLC (n = 17). Genes with mutation frequencies \u0026gt;5% and previously reported as significantly mutated genes in MPLC are shown for each region of the individual patient. Upper, number of driver mutations; middle, heatmap for driver mutations; right, components of mutational signatures. f. The correlation of the variant allele frequency detected by Senteion for the individual focus.\u003c/p\u003e\u003cp\u003e\u003cbr\u003e\u003c/p\u003e","description":"","filename":"Fig1.jpg","url":"https://assets-eu.researchsquare.com/files/rs-1716190/v1/08cc19ecd92caacf87ecac26.jpg"},{"id":22243778,"identity":"9fda63cd-0565-4565-bac8-d50301752d7a","added_by":"auto","created_at":"2022-06-03 22:30:43","extension":"jpg","order_by":2,"title":"Figure 2","display":"","copyAsset":false,"role":"figure","size":676119,"visible":true,"origin":"","legend":"\u003cp\u003eIdentification and characterization of malignant and nonmalignant cells in normal and tumor EPCAM+ cells. a, UMAP plot of all epithelial cells of the patients, color-coded according to malignancy (upper left), cell type (lower left), or patient and sample information (right). b, Canonical marker genes were used to label epithelial subtypes as represented in the UMAP plot. c, Unsupervised pseudotime trajectory analysis of all epithelial cells visualizing information about malignancy, anatomical pathology, cell type and MSI. d, Relative proportions of all epithelial cell subpopulations. e, f, The enriched Gene Ontology terms for all signature genes in epithelial pseudotime analysis state 1 (top) and state 2 (bottom). The significant features are marked with colors.\u003c/p\u003e\u003cp\u003e\u003cbr\u003e\u003c/p\u003e","description":"","filename":"Fig2.jpg","url":"https://assets-eu.researchsquare.com/files/rs-1716190/v1/8c48db7e1fcd9d47dddf9f33.jpg"},{"id":22243767,"identity":"4af08608-69e1-4d9e-ab39-2d3caec798df","added_by":"auto","created_at":"2022-06-03 22:30:29","extension":"jpg","order_by":3,"title":"Figure 3","display":"","copyAsset":false,"role":"figure","size":1270652,"visible":true,"origin":"","legend":"\u003cp\u003eClonal dynamics and phenotypic transitions of tumor-infiltrating T cells in MPLC. a, UMAP plot of a total of 81171 T cells, colored according to cell subtype. b, Canonical marker genes were used to label T cells subtypes as represented in the UMAP plot. c, The balloon plot of marker genes in the subtype of T cells. d, Relative proportions of all T cells subpopulations. e, UMAP of Normal(Upper), AIS/MIA (Middle) and IAC(Down) tumor-infiltrating in T cells colored by selected TCR clones. f. Overlap between T cell clones in different organ and cell types, divided by mouse. Each tile represents the overlap coefficient of clones. Color intensity indicates overlap strength.\u003c/p\u003e\u003cp\u003e\u003cbr\u003e\u003c/p\u003e","description":"","filename":"Fig3.jpg","url":"https://assets-eu.researchsquare.com/files/rs-1716190/v1/d1899952a8c4a83a563a7c91.jpg"},{"id":22243770,"identity":"74d151b5-df5a-4a9b-87d1-ce510535b6cc","added_by":"auto","created_at":"2022-06-03 22:30:35","extension":"jpg","order_by":4,"title":"Figure 4","display":"","copyAsset":false,"role":"figure","size":664806,"visible":true,"origin":"","legend":"\u003cp\u003eTrajectory profiling reveals the branched progression of T cells in MPLC. \u003cstrong\u003ea, \u003c/strong\u003eSingle-cell trajectories of malignant cells in CD8+ T cells. The subplot colored by global umap cluster (Left) and anatomical pathology info (Right). b, The enriched gene ontology terms for all signature genes in CD8+ T cells pseudotime analysis in different states (Immune response (left-up),T cell activation and differentiation(right-up),Response to stress and cell death (left-down) and Cytokine production (right-down) ) and the significantly features are marked with colors. c, Single-cell trajectories of malignant cells in CD4+ T cells colored by global umap cluster and anatomical pathology info (down).\u003c/p\u003e\u003cp\u003e\u003cbr\u003e\u003c/p\u003e","description":"","filename":"Fig4.jpg","url":"https://assets-eu.researchsquare.com/files/rs-1716190/v1/405e3d6c79bbd187d314a570.jpg"},{"id":22243768,"identity":"984c0d74-3277-4693-a72e-25568735865c","added_by":"auto","created_at":"2022-06-03 22:30:33","extension":"jpg","order_by":5,"title":"Figure 5","display":"","copyAsset":false,"role":"figure","size":1061036,"visible":true,"origin":"","legend":"\u003cp\u003eB-cell and myeloid cell distributions in MPLC. a, UMAP projection of 4279 B cells showing the composition of the five main subtypes. b, Unsupervised pseudotime trajectory analysis of all B cells. c, Relative proportions of all B-cell subpopulations. Each bar plot represents different samples, and the last bar plot is grouped by sample type (normal, AIS/MIA or IAC). d, Differentially expressed pathways were scored per cell by GSVA between five B-cell subtypes. e, UMAP projection of 46836 myeloid cells, showing the composition of the ten main subtypes. f, Percentages of each myeloid cell subtype among normal, AIS/MIA and IAC samples. Y axis: Average percentage of samples across the three groups. Groups are separated by blank space. Each bar plot represents different samples, and the last bar plot is grouped by sample type (Normal, AIS/MIA and IAC). g, Balloon plot of marker genes in the myeloid cell subtype. h, Unsupervised pseudotime trajectory analysis of all myeloid cells.\u003c/p\u003e\u003cp\u003e\u003cbr\u003e\u003c/p\u003e","description":"","filename":"Fig5.jpg","url":"https://assets-eu.researchsquare.com/files/rs-1716190/v1/906b5d34a529b70bdde29f7d.jpg"},{"id":22243772,"identity":"f0c910d4-6c86-40ed-8e78-0d9f133dc893","added_by":"auto","created_at":"2022-06-03 22:30:36","extension":"jpg","order_by":6,"title":"Figure 6","display":"","copyAsset":false,"role":"figure","size":1085468,"visible":true,"origin":"","legend":"\u003cp\u003eEndothelial cell and fibroblast cell distributions in MPLC. a, UMAP projection of 2081 endothelial cells showing the composition of the four main subtypes. b, Percentages of each endothelial cell subtype among normal, AIS/MIA and IAC samples. Y axis: Average percentage of samples across the three groups. Groups are separated by blank space. Each bar plot represents different samples, and the last bar plot is grouped by sample type (Normal, AIS/MIA and IAC). c, Violin plot of marker gene expression in Endothelial_C2 cell clusters. d, Differentially expressed pathways were scored per cell by GSVA between the four endothelial cell subtypes. e, UMAP projection of 2525 fibroblast cells showing the composition of the five main subtypes. f, Percentages of each fibroblast cell subtype among normal, AIS/MIA and IAC samples. g, Differentially expressed pathways were scored per cell by GSVA between the five fibroblast cell subtypes.\u003c/p\u003e\u003cp\u003e\u003cbr\u003e\u003c/p\u003e","description":"","filename":"Fig6.jpg","url":"https://assets-eu.researchsquare.com/files/rs-1716190/v1/671601e333949aae73d959c3.jpg"},{"id":22243774,"identity":"e9450507-bab2-4a00-875a-82cf2846af6d","added_by":"auto","created_at":"2022-06-03 22:30:39","extension":"jpg","order_by":7,"title":"Figure 7","display":"","copyAsset":false,"role":"figure","size":942924,"visible":true,"origin":"","legend":"\u003cp\u003eIntercellular interactions in MIA lung tumors and IAC lung tumors. We quantified the predicted cell interactions in MIA and IAC lung cancer. a, Networks depicting cell types as nodes and interactions as edges. The size of the node is proportional to the total number of interactions of the corresponding cell type, and edge thickness is proportional to the number of interactions between the connected types. Left is MIA, and right is IAC. b, Heatmap depicting the number of all possible interactions between the clusters analyzed. Different cell types are identified by different colors. Left is MIA, and right is IAC. c, Dot plot depicting selected immune–epithelial, epithelial–mesenchymal and mesenchymal–epithelial interactions enriched in MPLC (left to right: Normal, MIA, IAC).\u003c/p\u003e\u003cp\u003e\u003cbr\u003e\u003c/p\u003e","description":"","filename":"Fig7.jpg","url":"https://assets-eu.researchsquare.com/files/rs-1716190/v1/03df4270075d253803f42f2b.jpg"},{"id":22244249,"identity":"0567013b-c19c-441d-a5da-1eaca933d5f7","added_by":"auto","created_at":"2022-06-03 22:35:38","extension":"pdf","order_by":0,"title":"","display":"","copyAsset":false,"role":"manuscript-pdf","size":1683422,"visible":true,"origin":"","legend":"","description":"","filename":"manuscript.pdf","url":"https://assets-eu.researchsquare.com/files/rs-1716190/v1/62100fdf-8510-43b2-b5c6-79440383c223.pdf"},{"id":22243777,"identity":"d9a1d412-45fb-45c5-a906-dc06a420a5e4","added_by":"auto","created_at":"2022-06-03 22:30:39","extension":"pdf","order_by":9,"title":"","display":"","copyAsset":false,"role":"supplement","size":62952474,"visible":true,"origin":"","legend":"","description":"","filename":"Patientsimagingandpathologicalinformation.pdf","url":"https://assets-eu.researchsquare.com/files/rs-1716190/v1/c2bd117705923d6c5492c12b.pdf"},{"id":22243775,"identity":"eb53ad1b-caf7-4bfc-a0ea-65d3ea5d5c5b","added_by":"auto","created_at":"2022-06-03 22:30:39","extension":"pdf","order_by":10,"title":"","display":"","copyAsset":false,"role":"supplement","size":2201457,"visible":true,"origin":"","legend":"","description":"","filename":"Supplementalfigures.pdf","url":"https://assets-eu.researchsquare.com/files/rs-1716190/v1/efdc24ea0b2c05a3ee79bed8.pdf"},{"id":22243776,"identity":"313d9882-e557-4711-8469-4e60ee25cbbf","added_by":"auto","created_at":"2022-06-03 22:30:39","extension":"xlsx","order_by":11,"title":"","display":"","copyAsset":false,"role":"supplement","size":10821,"visible":true,"origin":"","legend":"","description":"","filename":"Supplementarytable1.xlsx","url":"https://assets-eu.researchsquare.com/files/rs-1716190/v1/dc92726421117c0915827f8b.xlsx"},{"id":22243769,"identity":"a4d9b523-ca44-427e-a1ca-fe15f4942e48","added_by":"auto","created_at":"2022-06-03 22:30:35","extension":"xls","order_by":12,"title":"","display":"","copyAsset":false,"role":"supplement","size":28160,"visible":true,"origin":"","legend":"","description":"","filename":"Supplementarytable2.xls","url":"https://assets-eu.researchsquare.com/files/rs-1716190/v1/d8e9811c41ae9d053d6b4827.xls"},{"id":22244241,"identity":"e74ffdbe-ce22-449f-8ced-65e432b721f8","added_by":"auto","created_at":"2022-06-03 22:35:35","extension":"xls","order_by":13,"title":"","display":"","copyAsset":false,"role":"supplement","size":19968,"visible":true,"origin":"","legend":"","description":"","filename":"Supplementarytable3.xls","url":"https://assets-eu.researchsquare.com/files/rs-1716190/v1/3b322698dd8f55fee0973b0f.xls"}],"financialInterests":"Competing interest reported. Hao Xia is employee of Mingma Technologies Co., Ltd. All remaining authors have declared no conflicts of interest.","formattedTitle":"Dissecting transcriptional heterogeneity in multiple primary lung cancer by single cell RNA sequencing","fulltext":[{"header":"Introduction","content":"\u003cp\u003eLung cancer is the leading cause of cancer-related mortality worldwide \u003csup\u003e\u003cspan additionalcitationids=\"CR2\" citationid=\"CR1\" class=\"CitationRef\"\u003e1\u003c/span\u003e\u0026ndash;\u003cspan citationid=\"CR3\" class=\"CitationRef\"\u003e3\u003c/span\u003e\u003c/sup\u003e. Multiple primary lung cancer (MPLC) is defined as the occurrence of two or more primary lung cancers in the same patient. MPLC can be classified as synchronous or metachronous based on the time of occurrence; metachronous MPLC is defined by a diagnostic interval of 6 months or less between tumors\u003csup\u003e\u003cspan citationid=\"CR4\" class=\"CitationRef\"\u003e4\u003c/span\u003e,\u003cspan citationid=\"CR5\" class=\"CitationRef\"\u003e5\u003c/span\u003e\u003c/sup\u003e. Port et al. reported that MPLC was observed in approximately 16% of patients with operable stage I, II, and III non-small-cell lung cancer (NSCLC), usually by preoperative imaging analysis\u003csup\u003e\u003cspan citationid=\"CR6\" class=\"CitationRef\"\u003e6\u003c/span\u003e\u003c/sup\u003e. Concomitant with the rising prevalence of lung cancer, the incidence of MPLC is also increasing \u003csup\u003e\u003cspan citationid=\"CR4\" class=\"CitationRef\"\u003e4\u003c/span\u003e,\u003cspan citationid=\"CR7\" class=\"CitationRef\"\u003e7\u003c/span\u003e\u003c/sup\u003e. Thus, an in-depth understanding of its mechanism and better treatment options for MPLC are urgently needed.\u003c/p\u003e \u003cp\u003eTumor heterogeneity describes differences between tumors of the same apparent type but with different genetic backgrounds and different carcinogen exposures (intertumoral heterogeneity) or differences between cancer cells within the same tumor (intratumoral heterogeneity)\u003csup\u003e\u003cspan citationid=\"CR8\" class=\"CitationRef\"\u003e8\u003c/span\u003e\u003c/sup\u003e. MPLC, by definition, arises in different areas of the lung parenchyma within a single patient and shares an identical germline genetic background and environmental exposure profile\u003csup\u003e\u003cspan citationid=\"CR9\" class=\"CitationRef\"\u003e9\u003c/span\u003e\u003c/sup\u003e. Previous studies have demonstrated differences in genomic profiles, including somatic point mutations, copy number aberrations, chromosomal structural variations and even mutational spectra, between different MPLC lesions\u003csup\u003e\u003cspan additionalcitationids=\"CR10\" citationid=\"CR9\" class=\"CitationRef\"\u003e9\u003c/span\u003e\u0026ndash;\u003cspan citationid=\"CR11\" class=\"CitationRef\"\u003e11\u003c/span\u003e\u003c/sup\u003e. The highly complex and heterogeneous tumor ecosystem contains not only malignant cells but also interacting cells from the host, such as endothelial cells, stromal fibroblasts, and a variety of immune cells that control tumor growth and metastasis\u003csup\u003e\u003cspan citationid=\"CR12\" class=\"CitationRef\"\u003e12\u003c/span\u003e\u003c/sup\u003e. However, the heterogeneity of cell types in the tumor microenvironment (TME) in MPLC remains unclear. Thus, it is critically important to decode the complex interplay between tumor cells and the TME in MPLC.\u003c/p\u003e \u003cp\u003eIn the present study, we performed whole-exome sequencing (WES) and single-cell RNA sequencing (scRNA-seq) on 17 tumor samples and 6 matched normal samples from 6 patients with MPLCs. By comparing tumors at different stages in MPLCs and adjacent normal lung tissue, we comprehensively characterized the transcriptome features of the malignant cells, immune cells, and stromal cells of MPLCs, and we decoded dynamic changes in cell percentage, cell subtype heterogeneity, and intercellular interactions, providing new knowledge regarding the biological basis of MPLC and LUAD development.\u003c/p\u003e"},{"header":"Methods","content":"\u003cdiv id=\"Sec3\" class=\"Section2\"\u003e \u003ch2\u003ePatients and sample collection\u003c/h2\u003e \u003cp\u003ePatients who were diagnosed with pathologically confirmed MPLC with two or three tumors in the ipsilateral or bilateral lungs and underwent surgery (surgery for contralateral lung tumors was completed within 1\u0026ndash;2 months after the first surgery, depending on the patients\u0026rsquo; recovery status) at the Department of Thoracic Surgery of the National Cancer Center/National Clinical Research Center for Cancer/Cancer Hospital, Chinese Academy of Medical Sciences and Peking Union Medical College, Beijing, China, were enrolled. After MPLC samples were resected, tumor tissues were cut into two pieces along the long axis: one was processed for pathological diagnosis, and the other was processed for scRNA-seq and WES. The clinical information and specimen information of these patients are summarized in supplementary table 1. This study was approved by the National Cancer Center/Cancer Hospital Ethics Committee (Approval number: NCC2208). Written informed consent was obtained from all participants included in this study.\u003c/p\u003e \u003cp\u003e \u003cb\u003eTissue dissociation.\u003c/b\u003e \u003c/p\u003e \u003cp\u003eFor each sample, 100\u0026ndash;200 mg of tissue was homogenized by mincing with scissors into smaller pieces (-0.3 mm\u003csup\u003e\u003cspan citationid=\"CR2\" class=\"CitationRef\"\u003e2\u003c/span\u003e\u003c/sup\u003e per piece). Before tissue enzyme digestion, blood and MACS Tissue Storage Solution were washed from the tissue pieces by the addition of 10 ml of ice-cold PBS, followed by further gentle washing and tissue collection using a 70 \u0026micro;m strainer. Then, the tissue pieces were transferred into 5 ml of enzyme mix prepared by combining 4.7 ml of DMEM, 200 \u0026micro;l of Enzyme H, 100 \u0026micro;l of enzyme R and 25 \u0026micro;l of Enzyme A. The samples were incubated at 37\u0026deg;C for 30\u0026ndash;60 min in an Eppendorf thermomixer at 700 rpm. After incubation, the cell suspension was filtered through a prewetted 70 \u0026micro;m MACS SmartStrainer, placed in 50 ml tubes and centrifuged at 300 g for 5 min at 4\u0026deg;C. The cell pellet was washed with 5 ml chilled 1X PBS with 0.04% BSA and then centrifuged. One milliliter of chilled 1X red blood cell removal solution was added to the cell pellet, and the cells were resuspended by gentle pipetting with a wide-bore pipette tip and then incubated for 10 min on ice. After incubation, 10 ml chilled 1X PBS with 0.04% BSA buffer was immediately added to the sample tube. The tube was gently inverted several times to mix and centrifuged at 4\u0026deg;C at 300 g for 10 min. The cell pellet was resuspended in 1 ml chilled 1X PBS with 0.04% BSA buffer, and the cell concentration and viability were counted using a Countess II automated cell counter. If the percentage of viable cells was \u0026lt;\u0026thinsp;70%, dead cells were removed to increase the proportion of viable cells. An appropriate volume of buffer was added to the cell suspension and gently mixed to achieve a target cell concentration of 700\u0026ndash;1200 cells/\u0026micro;l. Once the target cell concentration was obtained, GEM generation was immediately performed with 10X Genomics reagents. The remaining cells were frozen at -80\u0026deg;C for WES (no more than 5 \u0026times; 10\u003csup\u003e6\u003c/sup\u003e cells were used).\u003c/p\u003e \u003cp\u003e \u003cb\u003eWhole-exome Sequencing.\u003c/b\u003e \u003c/p\u003e \u003cp\u003eWe used a QIAamp DNA Mini Kit (Catalog# 51304) to process DNA extraction. The cell pellet was resuspended in PBS to a final volume of 200 \u0026micro;l, and 20 \u0026micro;l proteinase K was added to each sample and mixed by vortexing. Then, 4 \u0026micro;l of RNase A stock solution (100 mg/ml) was added to the sample and mixed completely, and 200 \u0026micro;l of Buffer AL was added to the sample and mixed by pulse-vortexing for 15 s. The samples were incubated at 56\u0026deg;C for 10 min, and then 200 \u0026micro;l of ethanol was added. The mixture was carefully transferred to a QIAamp Mini Spin column. The cap was closed, and the tube was centrifuged at 6000 x g for 1 min. The QIAamp Mini Column was placed in a clean 2 ml collection tube, and the tube containing the filtrate was discarded. The column was carefully opened, and 500 \u0026micro;l Buffer AW1 was added. The cap was closed, and the tube was centrifuged at 6000 x g for 1 min. This step was repeated with buffer AW2 and eluted with 100 \u0026micro;l buffer AE. One microliter of DNA sample was taken for quantification on a 1% agarose gel. After passing the QC, the SureSelect XT Library Prep Kit (Catalog#: 5500\u0026thinsp;\u0026minus;\u0026thinsp;0132) was used for library construction and QC according to the official protocol. Finally, we used an Illumina Nova 6000 for 150 bp paired end (PE 150) sequencing.\u003c/p\u003e \u003cp\u003e \u003cb\u003eBulk RNA sequencing.\u003c/b\u003e \u003c/p\u003e \u003cp\u003eWe used the RNeasy MinElute Cleanup Kit (Catalog: 74204) to extract RNA. The cell pellet was resuspended in TRIzol reagent to a final volume of 1000 \u0026micro;l. Then, the samples were vortexed, and 200 \u0026micro;l chloroform was added to mix the samples. The lysate was centrifuged for 15 min at 12000 x g at 4\u0026deg;C. The supernatant was pipetted, and one volume of RLT buffer was added and mixed well. Then, one volume of ethanol was added and mixed well. The samples were transferred to a RNeasy MinElute spin column and centrifuged for 15 s at 8000 x g. The flowthrough was discarded, and the same steps were repeated with buffer RW1. Then, 80 \u0026micro;l DNase I solution was added to the samples and incubated for 8 min. After incubation, we added 350 \u0026micro;l of RW1 to the samples and centrifuged them at 8000 x g for 15 s. The flowthrough was discarded, and the wash steps were repeated with 500 \u0026micro;l of RPE buffer and 500 \u0026micro;l of freshly prepared 80% ethanol. Finally, the RNA samples were eluted in 30 \u0026micro;l RNase-free water, and 1 \u0026micro;l of the RNA sample was taken for quantification and evaluation of RNA integrity. After passing the QC, we used the TruSeq RNA Library Prep for Enrichment Kit (Catalog#: 20394426) for library construction according to the official protocol. Finally, we used an Illumina Nova 6000 for PE 150 sequencing.\u003c/p\u003e \u003cp\u003e \u003cb\u003eChromium 10x Genomics library and sequencing.\u003c/b\u003e \u003c/p\u003e \u003cp\u003eFor 10X Chromium Single Cell Immune Profiling, target tissue was prepared into single-cell suspensions by using commercial enzymes from Miltenyi. Then, the cell concentration and viability were determined using a Countess II Automated Cell Counter and recorded. The cells were diluted to 700\u0026ndash;1200 cells/\u0026micro;l. Cell loading, GEM generation and library construction of T cells and B cells were performed based on the Chromium Next GEM Single Cell 5' Kit official protocol. After library construction, we used an Agilent 2100 Bioanalyzer to perform a library check. Finally, we used an Illumina Nova 6000 for PE 150 sequencing. Additional sequencing was performed to obtain coverage of at least 50K reads per cell.\u003c/p\u003e \u003cp\u003e \u003cb\u003eWES analysis.\u003c/b\u003e \u003c/p\u003e \u003cp\u003eRaw reads may contain adapter sequences and low-quality bases, which may affect the subsequent analysis. Therefore, it is necessary to review the raw sequencing reads. In this project, read quality was calculated for all samples by using FastQC (version 0.11.2) \u003csup\u003e\u003cspan citationid=\"CR13\" class=\"CitationRef\"\u003e13\u003c/span\u003e\u003c/sup\u003e software with the default parameters. We used the BWA-MEM (Burrows\u0026ndash;Wheeler Aligner-maximal exact match)\u003csup\u003e\u003cspan citationid=\"CR14\" class=\"CitationRef\"\u003e14\u003c/span\u003e\u003c/sup\u003e algorithm to compare the sequencing data to the reference genome (GRCh38), obtain the bam file, mark the PCR reads in the bam file, rearrange the regions that may have indel mutations, and recalibrate the quality of each base pair. Gene variants refer to changes in the sequence of the genome, including point mutations caused by single base changes and deletions, duplications, and insertions of multiple bases. In this analysis, the Sentieon (201611.02) \u003csup\u003e\u003cspan citationid=\"CR15\" class=\"CitationRef\"\u003e15\u003c/span\u003e\u003c/sup\u003e workflow was carried out to detect SNVs (single nucleotide variants) and indels (small insertions and deletions). The TMB result was the sum of nonsynonymous mutations and indel numbers in the exon region divided by the length of the total exon.\u003c/p\u003e \u003cp\u003e \u003cb\u003eSingle-cell RNA sequencing data alignment.\u003c/b\u003e \u003c/p\u003e \u003cp\u003eThe 10X Genomics single-cell raw sequencing paired-end reads were processed using CellRanger software version 5.0.1 and the 10x human genome GRCh38 release as the reference\u003csup\u003e\u003cspan citationid=\"CR16\" class=\"CitationRef\"\u003e16\u003c/span\u003e\u003c/sup\u003e. The output-filtered gene expression matrices were analyzed by R (version 4.0.3) with the Seurat package (version 4.0.0). Low-quality cells with (i)\u0026thinsp;\u0026lt;\u0026thinsp;500 UMIs (unique molecular identifiers), (ii)\u0026thinsp;\u0026gt;\u0026thinsp;25% UMIs derived from the mitochondrial genome, or (iii)\u0026thinsp;\u0026gt;\u0026thinsp;5000 or \u0026lt;\u0026thinsp;200 genes detected (genes expressed at a proportion\u0026thinsp;\u0026gt;\u0026thinsp;0.1%) were removed. Downstream analyses, including normalization, scaling, clustering of cells, and identifying cluster marker genes, were performed by using the function in Seurat. The gene expression matrices were normalized by the NormalizeData function, then the FindVariableFeatures function was used to select 2000 features with high cell-to-cell variation. To reduce the dimensionality of the datasets, the RunPCA function was conducted with default parameters, including the selected 2000 highly variable genes on linear transformation-scaled data generated by the ScaleData function. Finally, we clustered cells using the suitable PCs to calculate the Uniform Manifold Approximation and Projection (UMAP) for dimension reduction with the FindNeighbors and FindClusters functions. Nonlinear dimensional reduction was performed with the RunUMAP function with default settings.\u003c/p\u003e \u003cp\u003e \u003cb\u003eCell type annotation and cluster marker identification.\u003c/b\u003e \u003c/p\u003e \u003cp\u003eCells could be clustered together according to common features after nonlinear dimensional reduction. Then, all cells were projected into two-dimensional space by UMAP. We used the FindAllMarkers function in Seurat to find the marker genes with significant differences in expression for each cluster. All clusters were classified and annotated on the basis of these significantly different canonical markers. Clusters that expressed no canonical cell type markers were classified as \u0026ldquo;Others\u0026rdquo;.\u003c/p\u003e \u003cp\u003e \u003cb\u003eSubclustering of major cell types.\u003c/b\u003e \u003c/p\u003e \u003cp\u003eThe subset function in Seurat was used to extract the subcluster based on the expression of canonical markers of particular cell types. Next, these major cell types were integrated for further subclustering. After integration, the genes were scaled to unit variance. Scaling, PCA, and clustering were performed as described above.\u003c/p\u003e \u003cp\u003e \u003cb\u003eCell number statistics.\u003c/b\u003e \u003c/p\u003e \u003cp\u003eR was used to calculate the numbers of cells of different cell types. The ggplot2 library was used to generate relevant graphical outputs. The t test was used to test the significance between different groups.\u003c/p\u003e \u003cp\u003e \u003cb\u003eIdentification of malignant epithelial cells.\u003c/b\u003e \u003c/p\u003e \u003cp\u003eThe R code inferCNV (\u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttps://github.com/broadinstitute/infercnv\u003c/span\u003e\u003cspan address=\"https://github.com/broadinstitute/infercnv\" targettype=\"URL\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e) was used to infer the CNVs in epithelial cells from the scRNA-seq data. The epithelial cells (EPCAM+) of the normal sample were considered putative nonmalignant cells, and their CNV estimates were used to define a baseline. The epithelial cells from tumor samples were compared to the normal reference baseline. The CNV R-scores in different regions are readily apparent and are overabundant or less abundant in the regions of the tumor genome compared to that of putative nonmalignant cells.\u003c/p\u003e \u003cp\u003eThe calculated CNV signal was defined as the mean square of the CNV estimates across all genomic locations. The calculated CNV R-scores were defined as the Pearson correlation coefficient between each cell\u0026rsquo;s CNV pattern and the average CNV pattern of the top 5% of cells from the same tumor with respect to CNV signal \u003csup\u003e\u003cspan citationid=\"CR17\" class=\"CitationRef\"\u003e17\u003c/span\u003e\u003c/sup\u003e. Epithelial cells (EPCAM+) with CNV R-scores\u0026thinsp;\u0026ge;\u0026thinsp;0.4 and CNV signal\u0026thinsp;\u0026ge;\u0026thinsp;0.03 were defined as malignant cells, while those with CNV R-scores\u0026thinsp;\u0026le;\u0026thinsp;0.4 and CNV signal\u0026thinsp;\u0026le;\u0026thinsp;0.03 were defined as nonmalignant cells. Cells that did not meet either of the above two conditions were defined as unresolved cells.\u003c/p\u003e \u003cp\u003e \u003cb\u003eCellPhoneDB.\u003c/b\u003e \u003c/p\u003e \u003cp\u003eCellPhoneDB was used with default parameters to infer the interactions between the different kinds of cells on a linear transformation scaled data matrix. The method infers potential receptor\u0026ndash;ligand interactions based on the expression of a receptor from one cell type and a ligand from another\u003csup\u003e\u003cspan citationid=\"CR18\" class=\"CitationRef\"\u003e18\u003c/span\u003e\u003c/sup\u003e. Only the receptors and ligands that were expressed in more than 30% of the cells in the specific cluster were considered. Network visualization was performed by using Cytoscape (version 3.8.2). The network layout was set to a force-directed layout.\u003c/p\u003e \u003cp\u003e \u003cb\u003eTrajectory analysis.\u003c/b\u003e Trajectory analysis was performed using Monocle version 3.0.2.3 \u003csup\u003e19\u003c/sup\u003e. We ordered the subclusters of cells, such as T cells or epithelial cells, onto a pseudotime trajectory based on the union of highly variable genes obtained from the cells of the subcluster, which were from normal, MIA and IAC samples. All parameters of Monocle3 were set to the defaults (variable genes with the following cutoff criteria: dispersion_empirical\u0026thinsp;\u0026gt;\u0026thinsp;dispersion_fit; and mean expression\u0026thinsp;\u0026gt;\u0026thinsp;0.001). Dimensional reduction and cell ordering were performed using the DDRTree method and the orderCells function.\u003c/p\u003e \u003cp\u003e \u003cb\u003eData processing of single-cell VDJ libraries.\u003c/b\u003e \u003c/p\u003e \u003cp\u003eTCR and BCR reads were aligned to the GRCh38 reference genome, and consensus VDJ annotation was performed using Cellranger vdj (version 5.0.1) with default parameters. The R library Immunarch (version 0.6.7) was used to calculate the distribution of the V gene clonotype and calculate statistics for some related indicators, including overlap and correlation.\u003c/p\u003e \u003cp\u003e \u003cb\u003eStatistical analyses.\u003c/b\u003e \u003c/p\u003e \u003cp\u003eThe statistical tools, methods, and thresholds used for each analysis are explicitly described within the Results or detailed in the figure legends or Materials and Methods.\u003c/p\u003e \u003c/div\u003e"},{"header":"Results","content":"\u003cp\u003e \u003cb\u003e1. Global cellular landscape in multiple primary lung cancer.\u003c/b\u003e \u003c/p\u003e \u003cp\u003eTo elucidate the cellular development in multiple primary lung cancer (MPLC), we analyzed a total of 167397 cells from 23 sampling sites in 6 patients with multiple foci of early-stage lung cancer from the National Cancer Center/National Clinical Research Center for Cancer/Cancer Hospital (Fig.\u0026nbsp;\u003cspan refid=\"Fig1\" class=\"InternalRef\"\u003e1\u003c/span\u003ea-b, Supplementary Fig.\u0026nbsp;1a, d and Supplementary Table\u0026nbsp;2). Among these cells, 62714 (23.2%) were from AIS/MIA, 65844 (37.5%) were from IAC, and 38839 (39.3%) were from Normal Lung. To increase the accuracy of cell type designation, we jointly applied principal component analysis (PCA) after correction for read depth and mitochondrial read counts. Through marker-based annotation, 8 major cell types were identified in these unbatched and comparable datasets by graph-based UMAP (Fig.\u0026nbsp;\u003cspan refid=\"Fig1\" class=\"InternalRef\"\u003e1\u003c/span\u003eb, Supplementary Fig.\u0026nbsp;1a, 1d): epithelial cells, alveolar cells, fibroblasts, endothelial cells, myeloid cells, T lymphocytes cells, B lymphocytes cells and mast cells. Global cell type annotations for MPLC patients were categorized according to the molecular subtype of lung tumor tissues. The marker-based annotation correlated well with the cell types obtained via Cell Marker and SingleR v.1.0 \u003csup\u003e20,21\u003c/sup\u003e (Fig.\u0026nbsp;\u003cspan refid=\"Fig1\" class=\"InternalRef\"\u003e1\u003c/span\u003eb and Supplementary Fig.\u0026nbsp;1c). T cells and epithelial cells show great differences over the course of development of MPLC. The relative abundance of B cells increased, while that of T cells decreased, from Normal Lung to AIS/MIA and IAC (Fig.\u0026nbsp;\u003cspan refid=\"Fig1\" class=\"InternalRef\"\u003e1\u003c/span\u003ed, Supplementary Fig.\u0026nbsp;1b, and Supplementary 1c). The most frequently mutated genes were EGFR (53%) and RBM10 (35%) (Fig.\u0026nbsp;\u003cspan refid=\"Fig1\" class=\"InternalRef\"\u003e1\u003c/span\u003ee and Supplementary Fig.\u0026nbsp;1g), and the most frequent variant type was missense mutation (Supplementary Fig.\u0026nbsp;1g). C\u0026thinsp;\u0026gt;\u0026thinsp;T and C\u0026thinsp;\u0026gt;\u0026thinsp;A transversions were enriched in these patients (Supplementary Fig.\u0026nbsp;1f and Supplementary Fig.\u0026nbsp;1g). The correlation of the variant allele frequency (VAF) for each individual focus was not significant (Fig.\u0026nbsp;\u003cspan refid=\"Fig1\" class=\"InternalRef\"\u003e1\u003c/span\u003ef). The significant mutations and the correlation of the VAF show that the foci in MPLC are independent and do not share an origin.\u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003cp\u003e \u003cb\u003e2. Signatures and metabolic disturbances of malignant cells in MPLC.\u003c/b\u003e \u003c/p\u003e \u003cp\u003eA total of 13176 normal epithelial cells and 14553 malignant epithelial cells were obtained from Normal Lung, AIS/MIA and IAC samples, and large-scale copy number variations (CNVs) were inferred with Normal Lung cells as references\u003csup\u003e\u003cspan additionalcitationids=\"CR23\" citationid=\"CR22\" class=\"CitationRef\"\u003e22\u003c/span\u003e\u0026ndash;\u003cspan citationid=\"CR24\" class=\"CitationRef\"\u003e24\u003c/span\u003e\u003c/sup\u003e. Malignant epithelial cells were identified by inferring multipatient copy number variations (CNVs) with normal epithelial cells (mainly immune cells, including T cells and B cells) as references. The CNV patterns showed that malignant cells formed clusters in patients with different origins, indicating a high degree of intertumoral heterogeneity (Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003ea, Supplementary Fig.\u0026nbsp;2c, 2e). The clusters of malignant cells from different foci were independent and rarely overlapped, while the clusters of normal cells were mixed (Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003ea, Supplementary Fig.\u0026nbsp;2a, 2c). Alveolar cells are epithelial cells in the lung and play a pivotal role in maintaining the integrity and function of the alveoli. Both alveolar type I cells and alveolar type II cells were positive for EPCAM gene expression \u003csup\u003e\u003cspan citationid=\"CR25\" class=\"CitationRef\"\u003e25\u003c/span\u003e,\u003cspan citationid=\"CR26\" class=\"CitationRef\"\u003e26\u003c/span\u003e\u003c/sup\u003e. Other EPCAM\u0026thinsp;+\u0026thinsp;cells, including club, ciliated, basal, goblet, mucous, serous, ionocyte, and neuroendocrine cells, were categorized as epithelial cells. We identified them according to the expression of canonical gene markers (Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003eb). We clustered the normal epithelial cells as alveolar type I cells (AT1; AGER+), alveolar type II cells (AT2; SFTPA1+), ciliated airway epithelial cells (ciliated; TPPP3+), secretory club cells (club; SCGB1A1+) and other cells (no canonical markers to classify) based on the previously described markers (Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003ea, Supplementary Fig.\u0026nbsp;2a, and Supplementary Fig.\u0026nbsp;2d) \u003csup\u003e\u003cspan citationid=\"CR17\" class=\"CitationRef\"\u003e17\u003c/span\u003e,\u003cspan citationid=\"CR27\" class=\"CitationRef\"\u003e27\u003c/span\u003e\u003c/sup\u003e, independent of patient origin. The relative abundance of malignant epithelial cells, i.e., Ciliated, Club and AT1 cells, was increased in malignant cells compared with Normal Lung cells, and the number of AT2 cells was decreased from Normal Lung to AIS/MIA and IAC. We used lung cancer cells to recapitulate the multilineage differentiation processes of normal epithelia, and subclustering analysis of all epithelial cell datasets revealed divergent differentiation states in MPLC (Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003ed). Intriguingly, the projection of tissue cells and pathological feature classification onto the trajectory assigned the upper right branch and the lower right branch as malignant cells, mainly from IAC tissue. The trajectory indicated that there are two states in the evolution of malignant cells from AIS/MIA to IACs (Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003ec). MSI showed slight changes along the total pseudotime trajectory in MPLC, but TMB did not. In the pseudotime model, this trajectory began from nonmalignant AIS/MIA cells. Then, the trail branched in two directions: some malignant cells became state 1, and some malignant cells became state 2 in IAC. Branch 1 exclusively expressed immune response markers (e.g., CEACAM6, CXCL14, HLA_DPB1, CD74, S100A14, SPINK5, etc.), while branch 2 expressed cell death markers (e.g., PDCD4, AKR1C1, IER3, RGCC, APOH, etc.) (Supplementary Fig.\u0026nbsp;2g). Functional gene set enrichment analysis showed that the transition from State 1 was involved in the immune response, immune system process and response to transforming growth factor (Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003ee). In addition, the expression of genes implicated in cell death, response to stress and cell cycle increased in State 2 (Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003ef), which shows that some malignant cells will respond to stress from the immune response, some will undergo programmed cell death, and some will enter the cell cycle due to selective pressure during the development of tumors in MPLC.\u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003cp\u003e \u003cb\u003e3. Clonal dynamics and phenotypic transitions of tumor-infiltrating T cells in MPLC.\u003c/b\u003e \u003c/p\u003e \u003cp\u003eA total of 81,971 T cells were obtained from all samples and clustered into 10 subtypes. Three subtypes of CD4\u0026thinsp;+\u0026thinsp;T cells were further clustered as Treg cells (TIGIT+, ENTPD1+, FOXP3+), Trm cells (CXCR6+) and CD4\u0026thinsp;+\u0026thinsp;na\u0026iuml;ve T cells (CCR7+, IL7R+). Three subtypes of CD8\u0026thinsp;+\u0026thinsp;T cells were further clustered as activated T cells/exhausted T cells (TNF+, IFNG+, GZMA+, GZMB+, KLRD1+), activated/effector T cells (TNF+, GZMA+, EOMES+), and effector T/Trm cells (ITGAE+, CXCR6+, GZMA+, GZMB+). One subtype of both CD4- and CD8-positive T cells was further clustered as na\u0026iuml;ve/memory T cells (IL7R+), and 2 subtypes with no CD4 or CD8 expression were further clustered as NK cells/exhausted T cells (HAVCR2+, IFNG+, GZMB+) and cytotoxic T cells/NK cells (XCL1+, XCL2+) (Fig.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003ea, Supplementary Fig.\u0026nbsp;3a) based on canonical markers (Fig.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003eb,\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003ec, Supplementary Fig.\u0026nbsp;3b), as described previously \u003csup\u003e\u003cspan additionalcitationids=\"CR29 CR30 CR31\" citationid=\"CR28\" class=\"CitationRef\"\u003e28\u003c/span\u003e\u0026ndash;\u003cspan citationid=\"CR32\" class=\"CitationRef\"\u003e32\u003c/span\u003e\u003c/sup\u003e. The relative abundances of classified na\u0026iuml;ve T cells, Treg cells and Trm cells were increased, while those of activated/exhausted T cells were decreased, from Normal Lung to AIS/MIA and IAC (Fig.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003ed). We next analyzed the lineage relationships between phenotypes and clonotypes in T cells. Globally, we found that cells in normal and IAC samples shared more common clonotypes in TCRs than AIS/MIA-group cells (Fig.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003ee, Supplementary Fig.\u0026nbsp;4a). To investigate whether the development of tumors within clones could inform lineage transitions between T-cell phenotypes, we aggregated all clonotypes and measured the fraction shared in different clusters. There was minimal significant clonotype sharing between T cells in different cell groups. Trm and na\u0026iuml;ve T cells shared more significant clonotypes with other T cells, and the shared clonotypes increased with the development of infiltration in the tumor tissue, suggesting an increasing transition from Trm and na\u0026iuml;ve T cells to other activated T-cell states (Supplementary Fig.\u0026nbsp;4a). The T-cell phenotypes in AIS/MIA samples include more clonotype 1, for which TRB-TRA is \u0026ldquo;CASSGLAAKPGELF-CAVRRGQNFV\u0026rdquo;, than other normal and IAC samples (Supplementary Fig.\u0026nbsp;4b). This sequence may play an important role in antigen presentation during tumor invasion. In normal tissue, there is not enough antigen from tumors to stimulate T cells, and the expression of this antigen presentation factor is low. In IAC tissue, the tumor has experienced stress from T cells in AIS/MIA, and new subclone evolution may not cause the same antigen stimulation. Therefore, clonotype 1, in which TRB-TRA is \u0026ldquo;CASSGLAAKPGELF-CAVRRGQNFV\u0026rdquo;, was expressed significantly only in AIS/MIA samples. The TCR clonotype shared more overlap in focus grouped according to the source of the patient compared to those grouped by anatomical pathology info and still had more correlation (Fig.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003ef).\u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003cp\u003e \u003cb\u003e4. Trajectory profiling reveals the branched progression of T cells in MPLC.\u003c/b\u003e \u003c/p\u003e \u003cp\u003eWe next used Monocle3 to analyze the CD8\u0026thinsp;+\u0026thinsp;and CD4\u0026thinsp;+\u0026thinsp;T-cell trajectories to further investigate the potential transitions in each of these cell types. The pseudotime trajectory axis derived from Monocle3 indicates that CD8\u0026thinsp;+\u0026thinsp;na\u0026iuml;ve and memory T cells transdifferentiate into effector T cells first. Then, some of these cells become exhausted T cells, and others transform into cytotoxic T cells (Fig.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e4\u003c/span\u003ea). We compared the expression of genes in the next-stage cells with those in the previous stage (Supplementary Fig.\u0026nbsp;3c, 3d). Cells on the memory T to effector T branch expressed immune response markers (e.g., TNFRSF14, NFE2L2, TNFRSF1B, IL2RG, etc.) (Supplementary Fig.\u0026nbsp;3c,3d), while na\u0026iuml;ve T cells transitioning to effector T cells expressed lymphocyte activation and T-cell differentiation markers (e.g., FOXP1, TGFBR2, TNFRSF4, PRDM1, etc.) (Supplementary Fig.\u0026nbsp;3c, 3d). Cells on the effector T to exhausted T branch expressed responses to stress and cell death markers (e.g., GNLY, CCL4, and FCGR3A) (Supplementary Fig.\u0026nbsp;3c, 3d), while effector T cells expressed cytokine production markers (e.g., CXCR4, CRTAM, etc.) (Supplementary Fig.\u0026nbsp;3c, 3d). Gene Ontology enrichment analysis confirmed that the transition from memory T cells to effector T cells was mainly implicated in the immune response, immune system process and immune response regulation. The transformation from memory to effector in T cells is mainly a process of immune response and regulation (Fig.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e4\u003c/span\u003eb (left-upper)). The function of the geneset in the transformation from na\u0026iuml;ve T to effector T cells was mainly related to lymphocyte activation, T-cell and leukocyte differentiation, T-cell activation and immune system development, which showed that the transformation from na\u0026iuml;ve to effector T cells was mainly a process of leukocyte and T-cell differentiation (Fig.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e4\u003c/span\u003eb (upper right)). Gene Ontology enrichment analysis confirmed that signature genes were enriched for response to stress, positive regulation of cell death and response to stimulus pathway during effector T-cell transformation to exhausted T cells (Fig.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e4\u003c/span\u003eb (lower right)). The functional gene set enrichment analysis from effector to cytotoxic T cells showed that the transition was mainly implicated in the regulation of cytokine production, which included the positive regulation of myeloid leukocyte cytokine production involved in the immune response, immune effector processes, macrophage cytokine production and the positive regulation of cytokine production involved in the immune response (Fig.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e4\u003c/span\u003eb (lower left)). The pseudotime trajectory axis indicates that CD4\u0026thinsp;+\u0026thinsp;na\u0026iuml;ve and Treg cells have four final differentiation states in themselves, which are labeled as Na\u0026iuml;ve_C1 to Na\u0026iuml;ve_C4 and Treg_C1 to Treg_C4 (Fig.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e4\u003c/span\u003ec, Supplementary Fig.\u0026nbsp;5a). Pathway analyses showed that Naive_C2 cells were enriched in some effector pathways [Hallmark pi3k_akt_mtor signaling, chemokine signaling pathway, and natural killer cell-mediated cytotoxicity], while Naive_C4 cells did not express some pathways [T-cell receptor signaling pathway, regulation of lymphocyte differentiation, regulation of T-cell differentiation, T-cell differentiation, T-cell proliferation, and T-cell activation]. The function of the T cells in the Na\u0026iuml;ve_C2 cluster is closer to that of effector or exhausted T cells, while the function of the Na\u0026iuml;ve_C4 cluster is more closely related to that of other na\u0026iuml;ve clusters (supplemental Fig.\u0026nbsp;5b (Left)). The GSVA results showed that Treg_C3 cells were enriched in the regulation of myeloid cell differentiation [myeloid leukocyte differentiation, myeloid cell differentiation, and regulation of myeloid cell differentiation], while Treg_C2 cells mainly participated in the regulation of cell differentiation, activation of the T-cell receptor signaling pathway [regulation of cell differentiation, T-cell receptor signaling pathway, T-cell receptor signaling pathway and T-cell activation]. Treg_C2 cells also had high expression levels of genes related to lymphocyte differentiation, including lymphocyte differentiation, regulation of T-cell differentiation, and regulation of lymphocyte differentiation. Treg_C1 cells mainly participates in the regulation of stem cell differentiation [regulation of stem cell differentiation, stem cell differentiation, and hematopoietic progenitor cell differentiation] (Supplementary Fig.\u0026nbsp;5b (Right)).\u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003cp\u003e \u003cb\u003e5. B cells have two functional distribution subtypes in MPLC.\u003c/b\u003e \u003c/p\u003e \u003cp\u003eA total of 4279 B cells were analyzed and classified into five subtypes (Fig.\u0026nbsp;\u003cspan refid=\"Fig5\" class=\"InternalRef\"\u003e5\u003c/span\u003ea, Supplementary Fig.\u0026nbsp;6a, b, supplementary Table\u0026nbsp;3): Bcells_C1 (MS4A1+, CD19+, CD74+, CD37+, BTG+) (49.7% in IAC, 10.3% in AIS/MIA, 6.4% in Normal); Bcells_C2 (NEIL1+, RGS13+) (0.9% in IAC, 0.4 in AIS/MIA, 0.0% in Normal); Bcells_C3 (CYTOR+, NUCB2+, ANKRD28+, TRIB1+) (7.8% in IAC, 2.9% in AIS/MIA, 1.9% in Normal); Bcells_C4 (6.2% in IAC, 3.8% in AIS/MIA, 0.9% in Normal); and Bcells_C5 (IL7R+, IL32+, CST7+, FYB1+) (5.9% in IAC, 1.5% in AIS/MIA, 1.4% in Normal) (Supplementary Fig.\u0026nbsp;6e, Supplementary Fig.\u0026nbsp;7a). The relative abundance of the B-cell subtype showed no uniform pattern in different patient foci and or when grouped by anatomical pathology (Fig.\u0026nbsp;\u003cspan refid=\"Fig5\" class=\"InternalRef\"\u003e5\u003c/span\u003ec, Supplementary Fig.\u0026nbsp;6c, d). The pathway analysis of B cells in MPLC further confirmed that there are two kinds of B cells with different functions. The regulation and activation of the innate immune response was highly activated in the Bcell-C4 and Bcell-C5 clusters (Fig.\u0026nbsp;\u003cspan refid=\"Fig5\" class=\"InternalRef\"\u003e5\u003c/span\u003ed, red pathway). However, the Bcell-C1, Bcell-C2, and Bcell-C3 clusters were enriched in pathways related to antigen binding, processing and presentation, including Reactome MHC class II Antigen Presentation (Fig.\u0026nbsp;\u003cspan refid=\"Fig5\" class=\"InternalRef\"\u003e5\u003c/span\u003ed, green pathway). Some B cells are involved in immune activation and regulation of the immune response, while others are involved in antigen presentation. We next used Monocle3 to analyze the trajectory performance of B cells to further investigate the potential transition in each of their cell types. The pseudotime trajectory axis derived from Monocle3 indicates that Bcells_C1 (MS4A1+, CD19+, CD74+, CD37+, BTG+) have two differentiation states; one is transdifferentiating into Bcell_C3 (NEIL1+, RGS13+) and Bcell_C5 (IL7R+, IL32+, CST7+, FYB1+), and the other is transitioning to Bcell_C4 (Fig.\u0026nbsp;\u003cspan refid=\"Fig5\" class=\"InternalRef\"\u003e5\u003c/span\u003eb, Supplementary Fig.\u0026nbsp;7b, 7c, 7d). Combined with the pathway enrichment results, these findings show that B cells mainly differentiated into two different types of cells with different functions in MPLC; one is involved in immune activation and regulation of the immune response, while the other is involved in antigen presentation. Different from the TCR clonotypes, which shared more overlap and correlation in foci grouped according to the source within the patient, the BCR clonotypes shared little overlap or correlation in MPLC (Supplementary Fig.\u0026nbsp;6f).\u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003cp\u003e \u003cb\u003e6. Myeloid cells have two functional distribution subtypes in MPLC.\u003c/b\u003e \u003c/p\u003e \u003cp\u003eSubclustering of 46836 myeloid cells revealed 10 subsets (Fig.\u0026nbsp;\u003cspan refid=\"Fig5\" class=\"InternalRef\"\u003e5\u003c/span\u003ee, Supplementary Fig.\u0026nbsp;9a.): M1 macrophages (MSR1+, MARCO+, etc.), M2 macrophages (MERTK+, MMP9+, etc.), MDSCs (CEACAM8+, ITGAM+, etc.), monocyte cells (TNFSF10+, FCER2+, CLEC10A+, etc.), dendritic cells (FCER1A+, CD1C+, CLEC4A+, etc.), myeloid dendritic cells (CCR7+, IL15+, CD1C+, etc.), plasma cells (IRF4+, CD36+, CXCR4+, etc.), secretory cells (CXCL10+) (Fig.\u0026nbsp;\u003cspan refid=\"Fig5\" class=\"InternalRef\"\u003e5\u003c/span\u003eg, Supplementary Fig.\u0026nbsp;9d). The percentage of M2 macrophages increased from normal to AIS/MIA to IAC, while the other subtypes were not significantly different in MIA/AIS and IAC (Fig.\u0026nbsp;\u003cspan refid=\"Fig5\" class=\"InternalRef\"\u003e5\u003c/span\u003ef, Supplementary Fig.\u0026nbsp;9b). Macrophages play an important role in maintaining organismal integrity by either directly participating in pathogen elimination or repairing tissue under sterile inflammatory conditions\u003csup\u003e\u003cspan citationid=\"CR33\" class=\"CitationRef\"\u003e33\u003c/span\u003e,\u003c/sup\u003e and Microphage_2 cells support tissue integrity by providing growth factors and healing capacity\u003csup\u003e\u003cspan citationid=\"CR34\" class=\"CitationRef\"\u003e34\u003c/span\u003e\u003c/sup\u003e. From a functional perspective, myeloid cells are divided into two main types: one exhibits high expression of IL1A, FCGR3A and IFN-γ (M1 macrophages), which can positively regulate the innate immune response. The other negatively regulates the adaptive immune response and regulation of stem cell differentiation (M2 macrophages). The increase in M2 macrophages will lead to malignant cells adapting to the immune response and escaping the response from the immune system (Supplementary Fig.\u0026nbsp;9c). The pseudotime trajectory axis derived from Monocle3 indicates that myeloid cells have four main branch trajectories (Fig.\u0026nbsp;\u003cspan refid=\"Fig5\" class=\"InternalRef\"\u003e5\u003c/span\u003eh). M2 macrophages were mainly concentrated in the left branch trajectory, M1 macrophages were mainly in the opposite branch trajectory, monocyte cells were in the middle to left stage, secretory cells were mainly concentrated in the lower left branch trajectory, and other Myeloid2_cells were mainly in the lower right branch trajectory (Fig.\u0026nbsp;\u003cspan refid=\"Fig5\" class=\"InternalRef\"\u003e5\u003c/span\u003eh). Next, we evaluated the expression of DEGs in the M1 macrophage (MSR1+; MARCO+) and M2 macrophage (MERTK+; MMP9+) populations in the whole trajectory. The expression of macrophage colony-stimulating factor genes, such as CSF1, MCEMP1, APOE, CD52 and FABP4, increased gradually during the process of transdifferentiating into M1 macrophages, and the expression of genes related to cytomembrane transport, namely, GPR183, CD48, CD93, and CCL5, increased during the process of transdifferentiating into M2 macrophages (Supplementary Fig.\u0026nbsp;9f).\u003c/p\u003e \u003cp\u003e \u003cb\u003e7. Positive regulation of endothelial cell development is extensive in IAC.\u003c/b\u003e \u003c/p\u003e \u003cp\u003eWe found four distinct subtypes by reclustering 2081 endothelial cells (Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e6\u003c/span\u003ea, Supplementary Fig.\u0026nbsp;8a): Endothelial_C1 (VIPR1+, IL1RL1+, TMEM100+, TNFSF10+, etc.), Endothelial_C2 (EMP1+, CD93+, PLVAP+, etc.), Endothelial_C3 (TNFAIP3+, CXCR4+, CREM+, etc.), and Endothelial_C4 (TFF3+, TIMP2+, STMN1+, etc.) (Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e6\u003c/span\u003ec, Supplementary Fig.\u0026nbsp;8d). The percentage of Endothelial_C1 and Endothelial_C2 in MIA/AIS and IAC were more different than those in normal samples, and the distribution of endothelial subtypes was not significantly different in MIA/AIS and IAC (Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e6\u003c/span\u003eb, Supplemental Fig.\u0026nbsp;8c). Then, we focused on the transcriptomic features and functions of Endothelial_C2 cells. The genes EMP1, TSPAN7, VWF, PLVAP and CD93 were expressed at high levels in these cells, and pathway analyses showed that Endothelial_C2 cells were enriched in positive regulation of endothelial cell development pathways [positive regulation of blood vessel diameter, regulation of endothelial cell differentiation pathway, and abnormal morphology of the great vessels] (supplemental Fig.\u0026nbsp;8d, Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e6\u003c/span\u003ed). Thus, Endothelial_C2 cells are related to the regulation of endothelial development, and their prevalence increases with the development of tumor invasion.\u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003cp\u003e \u003cb\u003e8. Fibroblast subtype distributions are similar in IAC, AIS/MIA and normal samples.\u003c/b\u003e \u003c/p\u003e \u003cp\u003eFibroblast cells maintain a very stable state in MPLC (Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e6\u003c/span\u003ec, Supplementary Fig.\u0026nbsp;8c). A total of 2523 fibroblast cells were analyzed, and five subtypes were identified (Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e6\u003c/span\u003ec). In MPLC, there was no significant difference in the number of different fibroblast clusters, and there were also no significant changes in the pathway and function of genes expressed in different fibroblast clusters (Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e6\u003c/span\u003ef, \u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e6\u003c/span\u003eg, Supplementary Fig.\u0026nbsp;8f). Fibroblast_C1 cells expressed high levels of genes associated with the Reactome elastic fiber formation annotation, such as COL3A1, COL1A1, COL6A3, COL10A1, and COL8A1, whereas Fibroblast_C5 cells highly expressed genes related to positive regulation of fibroblast proliferation, like RASD1, SFTPB and FOS (Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e6\u003c/span\u003eg, Supplementary Fig.\u0026nbsp;8g).\u003c/p\u003e \u003cp\u003e \u003cb\u003e9. Characterization of cell\u0026ndash;cell interactions involved in MPLC.\u003c/b\u003e \u003c/p\u003e \u003cp\u003eThere are many specific signaling network changes in MPLC, which depend on complex interactions between structural and inflammatory cells\u003csup\u003e\u003cspan citationid=\"CR35\" class=\"CitationRef\"\u003e35\u003c/span\u003e\u003c/sup\u003e. One cell expresses a receptor or ligand, and this \u0026ldquo;ligand\u0026ndash;receptor\u0026rdquo; pair will be defined as an incoming or outgoing interaction for this cell. We used CellPhoneDB to identify potential cell\u0026ndash;cell interactions and define their changes in primary lung cancer. Compared with the MIA samples, the cells in tumors exhibited significantly increased interactions (Fig.\u0026nbsp;\u003cspan refid=\"Fig7\" class=\"InternalRef\"\u003e7\u003c/span\u003ea). Alveolar and myeloid cells had the most outgoing interactions of all cell types in both MIA and IAC samples. In comparison with MIA, the interaction between T cells and alveolar cells was significantly increased (Fig.\u0026nbsp;\u003cspan refid=\"Fig7\" class=\"InternalRef\"\u003e7\u003c/span\u003eb). Myeloid, cytotoxic T, natural killer and B cells are involved in the action of recognizing and presenting the antigens of epithelial cells among the interactions from immune cells to epithelial cells. Analysis of the biological functions revealed that interactions related to antigen recognition and presentation (CD74_COPA, HLA-C_FAM3C, and HLA-DRB1_OGN) were less abundant in IAC than in MIA and normal samples, which means that some escape of immune recognition occurred with deeper invasion in MPLC (Fig.\u0026nbsp;\u003cspan refid=\"Fig7\" class=\"InternalRef\"\u003e7\u003c/span\u003ec upper). The interactions related to independent cellular adhesion (NECTIN1_CADM3, FGF1_FGFR4, and CD44_FGFR2) were still less abundant in epithelial cells than in fibroblasts, which means that epithelial cells promoted structural changes in fibroblasts and reduced cell adhesion. The interactions regulating the growth and differentiation of numerous types of cells (ACVR1_BMPR2_BMP6 and ACVR1_BMPR2_BMP5) were more abundant (Fig.\u0026nbsp;\u003cspan refid=\"Fig7\" class=\"InternalRef\"\u003e7\u003c/span\u003ec, middle). On the one hand, epithelial cells reduce fibroblast adhesion; on the other hand, they strengthen the differentiation of other mesenchymal cells with deeper invasion in MPLC. The interactions related to endothelial and hematopoietic stem cell differentiation and development (SELL-CD34, SELP-CD34, NPR1-NPPC, and NPR2-NPPC) were more abundant in IAC than in normal and MIA samples (Fig.\u0026nbsp;\u003cspan refid=\"Fig7\" class=\"InternalRef\"\u003e7\u003c/span\u003ec, down).\u003c/p\u003e \u003cp\u003e \u003c/p\u003e"},{"header":"Discussion","content":"\u003cp\u003eOur knowledge about the dynamics of gene expression and metabolic pathways during cellular changes and evolution in MPLC remains incomplete. In this study, we comprehensively provide important new insights into these changes during the progression from normal AIS to IAC in MPLC.\u003c/p\u003e \u003cp\u003eRecent studies have performed similar development trajectory analyses for CD8\u0026thinsp;+\u0026thinsp;T cells, with na\u0026iuml;ve T cells as the root and exhausted T cells and cytotoxic T cells as the two end states \u003csup\u003e\u003cspan citationid=\"CR17\" class=\"CitationRef\"\u003e17\u003c/span\u003e,\u003cspan citationid=\"CR36\" class=\"CitationRef\"\u003e36\u003c/span\u003e,\u003cspan citationid=\"CR37\" class=\"CitationRef\"\u003e37\u003c/span\u003e\u003c/sup\u003e. However, na\u0026iuml;ve T cells are not the only root for CD8\u0026thinsp;+\u0026thinsp;T-cell development and evolution in MPLC. We found that memory T cells are also a main source of these cells and play an important role in CD8\u0026thinsp;+\u0026thinsp;T-cell differentiation. The immune response and regulation of the immune response were significantly enriched during this differentiation process. The transition from CD8\u0026thinsp;+\u0026thinsp;memory T cells to effector T cells were mainly a process of positive regulation of immune system activities; however, the progression of CD8\u0026thinsp;+\u0026thinsp;na\u0026iuml;ve to effector T cells is mainly a process of leukocyte activation and T-cell differentiation. The progress of effector T-cell transformation to exhausted T cells involves responses to stimuli and positive regulation of cell death. The other pathway, of effector T-cell transformation to cytotoxic T cells, mainly involves the regulation of cytokine production, including the positive regulation of myeloid leukocyte cytokine production involved in the immune response and the positive regulation of cytokine production involved in the immune response. CD4\u0026thinsp;+\u0026thinsp;T cells have two distinct functional statuses: Treg and na\u0026iuml;ve. We also found that T cells shared more consistent TCR clonotypes than AIS/MIA-group cells in normal and IAC samples. Unfortunately, further in-depth study of this interesting phenomenon is beyond the scope of this study.\u003c/p\u003e \u003cp\u003eIn contrast to previous studies on the metabolism of malignant cells \u003csup\u003e17,38\u0026minus;40\u003c/sup\u003e, we used trajectory analysis to study the transition of epithelial cells from normal to malignant. We found that although different nodules in the patients were heterogeneous, they all exhibit two states during the evolution of malignant cells from AIS/MIA to IAC. Specifically, some cells respond to stress from immunity, while others undergo programmed cell death or enter the cell cycle due to selective pressure from the development of tumors in MPLC. The proportion of AT2 cells gradually decreases with the development of LUAD from normal to AIS/MIA to IAC, which also provides evidence that AT2 is the most likely origin of malignant cells in MPLC \u003csup\u003e\u003cspan citationid=\"CR32\" class=\"CitationRef\"\u003e32\u003c/span\u003e\u003c/sup\u003e.\u003c/p\u003e \u003cp\u003eThe functional distribution subtypes of B cells can be distinguished into two different stages. One is the regulation and activation of innate immune response antigen binding, and the other is processing and presentation, including MHC class II antigen presentation. Myeloid cells are divided into two major types, one with high expression of positive regulators of the innate immune response and the other performing negative regulation of the adaptive immune response and regulation of stem cell differentiation presented by M2 macrophages. With the progression of LUAD from normal to AIS/MIA to IAC, the increasing proportion of M2 macrophages causes malignant cells adapt to the immune response and escape the immune response \u003csup\u003e\u003cspan additionalcitationids=\"CR42\" citationid=\"CR41\" class=\"CitationRef\"\u003e41\u003c/span\u003e\u0026ndash;\u003cspan citationid=\"CR43\" class=\"CitationRef\"\u003e43\u003c/span\u003e\u003c/sup\u003e. We also found that some endothelial cell clusters also participate in the positive regulation of blood vessel development in early glandular neoplasia of the lung. Moreover, the proportion of these endothelial cells gradually increases with the progression of LUAD from normal to AIS/MIA to IAC \u003csup\u003e\u003cspan additionalcitationids=\"CR45\" citationid=\"CR44\" class=\"CitationRef\"\u003e44\u003c/span\u003e\u0026ndash;\u003cspan citationid=\"CR46\" class=\"CitationRef\"\u003e46\u003c/span\u003e\u003c/sup\u003e.\u003c/p\u003e \u003cp\u003eComprehensive analysis of the cell\u0026ndash;cell interactions in MPLC identified that interactions related to antigen recognition and presentation were less abundant in IAC than in AIS/MIA and normal samples, and the interactions related to independent cellular adhesion were still less abundant in epithelial cells than in fibroblasts. The interactions regulating the growth and differentiation of numerous types of cells were more abundant, and the interactions related to endothelial and hematopoietic stem cell differentiation and development were more abundant in the IAC than in the normal and AIS/MIA samples.\u003c/p\u003e \u003cp\u003eDespite the striking findings identified in this study, there are still several limitations, including sample collection, technical challenges and software limitations. We were unable to evaluate the continuous development of different nodules in the same patient, and the cells used for single-cell sequencing were from just a part of the tumor tissue, so there might be differences in biological complexity between our samples and the tissues used for pathological diagnosis and real cells in developing tumors. The 10x Genomics single-cell technology cannot cover all the transcripts in cells. There are also still limitations in batch effect removal, clustering and trajectory inference analysis, as these analyses are still maturing, and the complexity of the underlying topology could be underestimated \u003csup\u003e\u003cspan citationid=\"CR47\" class=\"CitationRef\"\u003e47\u003c/span\u003e\u003c/sup\u003e.\u003c/p\u003e"},{"header":"Declarations","content":"\u003cp\u003e\u003cstrong\u003eData availability\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eThe raw data that support the findings of this study are available from the corresponding author upon reasonable request.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eAcknowledgments\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eWe appreciate the support and participation of the physicians and patients in this study.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eCompeting interests\u0026nbsp;\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eHao Xia is employee of\u0026nbsp;Mingma\u0026nbsp;Technologies Co., Ltd. All remaining authors have declared no conflicts of interest.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eAuthor contributions\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eW.G. and S.G. conceived the study. J.H., J.Y., F.T., Q.X., and S.G. provided project management and supervision. W.G. and L.G. provided or facilitated the accrual of patient samples, pathology, and/or clinical data. W.G. and H.X. performed bioinformatics and genomic analyses. W.G. and H.X. performed statistical analyses. W.G. and Z.B. wrote the original draft, with input from all authors.\u0026nbsp;\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eFunding\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eThis study was supported by the National Key R\u0026amp;D Program of China (2021YFC2500900), the National Natural Science Foundation of China (82002451), the CAMS Initiative for Innovative Medicine (2021-1-I2M-015), the Special Research Fund for Central Universities, Peking Union Medical College (3332020024), the Non-profit Central Research Institute Fund of Chinese Academy of Medical Sciences (2018PT32033), and the Beijing Hope Run Special Fund of Cancer Foundation of China (LC2019B15).\u003c/p\u003e"},{"header":"References","content":"\u003col\u003e\u003cli\u003e\u003cspan\u003eChen, W. et al. Cancer statistics in China, 2015. CA: a cancer journal for clinicians 66, 115\u0026ndash;132 (2016).\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eSiegel, R. L., Miller, K. D., Fuchs, H. E. \u0026amp; Jemal, A. Cancer statistics, 2021. CA: a cancer journal for clinicians 71, 7\u0026ndash;33 (2021).\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eSung, H. et al. Global cancer statistics 2020: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA: a cancer journal for clinicians 71, 209\u0026ndash;249 (2021).\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eShintani, Y. et al. Clinical features and outcomes of patients with stage I multiple primary lung cancers. Cancer Sci 112, 1924\u0026ndash;1935, doi:\u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003e10.1111/cas.14748\u003c/span\u003e\u003cspan address=\"10.1111/cas.14748\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e (2021).\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eLiu, Y. et al. Genomic heterogeneity of multiple synchronous lung cancer. Nat Commun 7, 13200, doi:\u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003e10.1038/ncomms13200\u003c/span\u003e\u003cspan address=\"10.1038/ncomms13200\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e (2016).\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003ePort, J. L. et al. Surgical resection for multifocal (T4) non-small cell lung cancer: is the T4 designation valid? The Annals of thoracic surgery 83, 397\u0026ndash;400 (2007).\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eAsamura, H. Multiple primary cancers or multiple metastases, that is the question. J Thorac Oncol 5, 930\u0026ndash;931, doi:\u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003e10.1097/JTO.0b013e3181e15d8c\u003c/span\u003e\u003cspan address=\"10.1097/JTO.0b013e3181e15d8c\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e (2010).\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eMeacham, C. E. \u0026amp; Morrison, S. J. Tumour heterogeneity and cancer cell plasticity. Nature 501, 328\u0026ndash;337 (2013).\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eMa, P. et al. Simultaneous evolutionary expansion and constraint of genomic heterogeneity in multifocal lung cancer. Nature communications 8, 1\u0026ndash;11 (2017).\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eLiu, Y. et al. Genomic heterogeneity of multiple synchronous lung cancer. Nature communications 7, 1\u0026ndash;8 (2016).\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eMurphy, S. J. et al. Identification of independent primary tumors and intrapulmonary metastases using DNA rearrangements in non\u0026ndash;small-cell lung cancer. Journal of Clinical Oncology 32, 4050 (2014).\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eGiraldo, N. A. et al. The clinical role of the TME in solid cancer. British journal of cancer 120, 45\u0026ndash;53 (2019).\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eS. A., S. B. a. FastQC: a quality control tool for high throughput sequence data. (2010).\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eLi, H. \u0026amp; Durbin, R. Fast and accurate long-read alignment with Burrows-Wheeler transform. Bioinformatics 26, 589\u0026ndash;595, doi:\u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003e10.1093/bioinformatics/btp698\u003c/span\u003e\u003cspan address=\"10.1093/bioinformatics/btp698\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e (2010).\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eWeber JA, A. R., Gallagher BD, Edwards JS. Sentieon DNA pipeline for variant detection - Software-only solution, over 20\u0026times; faster than GATK 3.3 with identical results. PeerJ PrePrints 4:e1672v2, doi:\u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttps://doi.org/10.7287/peerj.preprints.1672v2\u003c/span\u003e\u003cspan address=\"10.7287/peerj.preprints.1672v2\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e (2016).\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003ehttps://github.com/10XGenomics/cellranger. Cellranger.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eXing, X. et al. Decoding the multicellular ecosystem of lung adenocarcinoma manifested as pulmonary subsolid nodules by single-cell RNA sequencing. Sci Adv 7, doi:\u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003e10.1126/sciadv.abd9738\u003c/span\u003e\u003cspan address=\"10.1126/sciadv.abd9738\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e (2021).\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eEfremova, M., Vento-Tormo, M., Teichmann, S. A. \u0026amp; Vento-Tormo, R. CellPhoneDB: inferring cell-cell communication from combined expression of multi-subunit ligand-receptor complexes. Nat Protoc 15, 1484\u0026ndash;1506, doi:\u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003e10.1038/s41596-020-0292-x\u003c/span\u003e\u003cspan address=\"10.1038/s41596-020-0292-x\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e (2020).\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eCao, J. et al. The single-cell transcriptional landscape of mammalian organogenesis. Nature 566, 496\u0026ndash;502, doi:\u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003e10.1038/s41586-019-0969-x\u003c/span\u003e\u003cspan address=\"10.1038/s41586-019-0969-x\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e (2019).\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eAran, D. et al. Reference-based analysis of lung single-cell sequencing reveals a transitional profibrotic macrophage. Nat Immunol 20, 163\u0026ndash;172, doi:\u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003e10.1038/s41590-018-0276-y\u003c/span\u003e\u003cspan address=\"10.1038/s41590-018-0276-y\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e (2019).\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eZhang, X. et al. CellMarker: a manually curated resource of cell markers in human and mouse. Nucleic Acids Res 47, D721-D728, doi:\u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003e10.1093/nar/gky900\u003c/span\u003e\u003cspan address=\"10.1093/nar/gky900\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e (2019).\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003ePatel, A. P. et al. Single-cell RNA-seq highlights intratumoral heterogeneity in primary glioblastoma. Science 344, 1396\u0026ndash;1401, doi:\u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003e10.1126/science.1254257\u003c/span\u003e\u003cspan address=\"10.1126/science.1254257\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e (2014).\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003ePuram, S. V. et al. Single-Cell Transcriptomic Analysis of Primary and Metastatic Tumor Ecosystems in Head and Neck Cancer. Cell 171, 1611\u0026ndash;1624 e1624, doi:\u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003e10.1016/j.cell.2017.10.044\u003c/span\u003e\u003cspan address=\"10.1016/j.cell.2017.10.044\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e (2017).\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eTimothy Tickle, I. T., Christophe Georgescu,Maxwell Brown,Brian Haas. Visualizing Large-scale Copy Number Variation in Single-Cell RNA-Seq Expression Data. (2021).\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eKhan, P. et al. Culture of human alveolar epithelial type II cells by sprouting. Respir Res 19, 204, doi:\u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003e10.1186/s12931-018-0906-9\u003c/span\u003e\u003cspan address=\"10.1186/s12931-018-0906-9\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e (2018).\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eTravaglini, K. J. et al. A molecular cell atlas of the human lung from single-cell RNA sequencing. Nature 587, 619\u0026ndash;625, doi:\u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003e10.1038/s41586-020-2922-4\u003c/span\u003e\u003cspan address=\"10.1038/s41586-020-2922-4\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e (2020).\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eDeprez, M. et al. A Single-Cell Atlas of the Human Healthy Airways. Am J Respir Crit Care Med 202, 1636\u0026ndash;1645, doi:\u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003e10.1164/rccm.201911-2199OC\u003c/span\u003e\u003cspan address=\"10.1164/rccm.201911-2199OC\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e (2020).\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eMorgan, A. J. et al. Expression of CXCR6 and its ligand CXCL16 in the lung in health and disease. Clin Exp Allergy 35, 1572\u0026ndash;1580, doi:\u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003e10.1111/j.1365-2222.2005.02383.x\u003c/span\u003e\u003cspan address=\"10.1111/j.1365-2222.2005.02383.x\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e (2005).\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eWein, A. N. et al. CXCR6 regulates localization of tissue-resident memory CD8 T cells to the airways. J Exp Med 216, 2748\u0026ndash;2762, doi:\u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003e10.1084/jem.20181308\u003c/span\u003e\u003cspan address=\"10.1084/jem.20181308\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e (2019).\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eGhorani, E. et al. The T cell differentiation landscape is shaped by tumour mutations in lung cancer. Nat Cancer 1, 546\u0026ndash;561, doi:\u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003e10.1038/s43018-020-0066-y\u003c/span\u003e\u003cspan address=\"10.1038/s43018-020-0066-y\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e (2020).\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eDelacher, M. et al. Single-cell chromatin accessibility landscape identifies tissue repair program in human regulatory T cells. Immunity 54, 702\u0026ndash;720 e717, doi:\u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003e10.1016/j.immuni.2021.03.007\u003c/span\u003e\u003cspan address=\"10.1016/j.immuni.2021.03.007\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e (2021).\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eWang, Z. et al. Deciphering cell lineage specification of human lung adenocarcinoma with single-cell RNA sequencing. Nat Commun 12, 6500, doi:\u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003e10.1038/s41467-021-26770-2\u003c/span\u003e\u003cspan address=\"10.1038/s41467-021-26770-2\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e (2021).\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eChavez-Galan, L., Olleros, M. L., Vesin, D. \u0026amp; Garcia, I. Much More than M1 and M2 Macrophages, There are also CD169(+) and TCR(+) Macrophages. Front Immunol 6, 263, doi:\u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003e10.3389/fimmu.2015.00263\u003c/span\u003e\u003cspan address=\"10.3389/fimmu.2015.00263\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e (2015).\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eLey, K. M1 Means Kill; M2 Means Heal. J Immunol 199, 2191\u0026ndash;2193, doi:\u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003e10.4049/jimmunol.1701135\u003c/span\u003e\u003cspan address=\"10.4049/jimmunol.1701135\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e (2017).\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eHolgate, S. T. et al. Asthma. Nat Rev Dis Primers 1, 15025, doi:\u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003e10.1038/nrdp.2015.25\u003c/span\u003e\u003cspan address=\"10.1038/nrdp.2015.25\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e (2015).\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eGuo, X. et al. Global characterization of T cells in non-small-cell lung cancer by single-cell sequencing. Nat Med 24, 978\u0026ndash;985, doi:\u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003e10.1038/s41591-018-0045-3\u003c/span\u003e\u003cspan address=\"10.1038/s41591-018-0045-3\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e (2018).\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eQian, J. et al. A pan-cancer blueprint of the heterogeneous tumor microenvironment revealed by single-cell profiling. Cell Res 30, 745\u0026ndash;762, doi:\u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003e10.1038/s41422-020-0355-0\u003c/span\u003e\u003cspan address=\"10.1038/s41422-020-0355-0\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e (2020).\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eGoveia, J. et al. An Integrated Gene Expression Landscape Profiling Approach to Identify Lung Tumor Endothelial Cell Heterogeneity and Angiogenic Candidates. Cancer Cell 37, 421, doi:\u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003e10.1016/j.ccell.2020.03.002\u003c/span\u003e\u003cspan address=\"10.1016/j.ccell.2020.03.002\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e (2020).\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eIreland, A. S. et al. MYC Drives Temporal Evolution of Small Cell Lung Cancer Subtypes by Reprogramming Neuroendocrine Fate. Cancer Cell 38, 60\u0026ndash;78 e12, doi:\u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003e10.1016/j.ccell.2020.05.001\u003c/span\u003e\u003cspan address=\"10.1016/j.ccell.2020.05.001\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e (2020).\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eHe, Y. et al. Mechanisms of Progression and Heterogeneity in Multiple Nodules of Lung Adenocarcinoma. Small Methods 5, e2100082, doi:\u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003e10.1002/smtd.202100082\u003c/span\u003e\u003cspan address=\"10.1002/smtd.202100082\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e (2021).\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eHeusinkveld, M. \u0026amp; van der Burg, S. H. Identification and manipulation of tumor associated macrophages in human cancers. J Transl Med 9, 216, doi:\u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003e10.1186/1479-5876-9-216\u003c/span\u003e\u003cspan address=\"10.1186/1479-5876-9-216\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e (2011).\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eConway, E. M. et al. Macrophages, Inflammation, and Lung Cancer. Am J Respir Crit Care Med 193, 116\u0026ndash;130, doi:\u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003e10.1164/rccm.201508-1545CI\u003c/span\u003e\u003cspan address=\"10.1164/rccm.201508-1545CI\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e (2016).\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eLiu, M. et al. Transcription factor c-Maf is a checkpoint that programs macrophages in lung cancer. J Clin Invest 130, 2081\u0026ndash;2096, doi:\u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003e10.1172/JCI131335\u003c/span\u003e\u003cspan address=\"10.1172/JCI131335\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e (2020).\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eLamalice, L., Le Boeuf, F. \u0026amp; Huot, J. Endothelial cell migration during angiogenesis. Circ Res 100, 782\u0026ndash;794, doi:\u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003e10.1161/01.RES.0000259593.07661.1e\u003c/span\u003e\u003cspan address=\"10.1161/01.RES.0000259593.07661.1e\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e (2007).\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eBetsholtz, C. Cell-cell signaling in blood vessel development and function. EMBO Mol Med 10, doi:\u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003e10.15252/emmm.201708610\u003c/span\u003e\u003cspan address=\"10.15252/emmm.201708610\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e (2018).\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eSweeney, M. \u0026amp; Foldes, G. It Takes Two: Endothelial-Perivascular Cell Cross-Talk in Vascular Development and Disease. Front Cardiovasc Med 5, 154, doi:\u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003e10.3389/fcvm.2018.00154\u003c/span\u003e\u003cspan address=\"10.3389/fcvm.2018.00154\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e (2018).\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eLahnemann, D. et al. Eleven grand challenges in single-cell data science. Genome Biol 21, 31, doi:\u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003e10.1186/s13059-020-1926-6\u003c/span\u003e\u003cspan address=\"10.1186/s13059-020-1926-6\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e (2020).\u003c/span\u003e\u003c/li\u003e\u003c/ol\u003e"}],"fulltextSource":"","fullText":"","funders":[],"hasAdminPriorityOnWorkflow":false,"hasManuscriptDocX":true,"hasOptedInToPreprint":true,"hasPassedJournalQc":"","hasAnyPriority":false,"hideJournal":true,"highlight":"","institution":"","isAcceptedByJournal":false,"isAuthorSuppliedPdf":false,"isDeskRejected":"","isHiddenFromSearch":false,"isInQc":false,"isInWorkflow":false,"isPdf":false,"isPdfUpToDate":true,"isWithdrawnOrRetracted":false,"journal":{"display":true,"email":"[email protected]","identity":"researchsquare","isNatureJournal":false,"hasQc":true,"allowDirectSubmit":true,"externalIdentity":"","sideBox":"","snPcode":"","submissionUrl":"/submission","title":"Research Square","twitterHandle":"researchsquare","acdcEnabled":true,"dfaEnabled":false,"editorialSystem":"","reportingPortfolio":"","inReviewEnabled":false,"inReviewRevisionsEnabled":true},"keywords":"multiple primary lung cancer, single-cell RNA sequencing, heterogeneity, tumor microenvironment","lastPublishedDoi":"10.21203/rs.3.rs-1716190/v1","lastPublishedDoiUrl":"https://doi.org/10.21203/rs.3.rs-1716190/v1","license":{"name":"CC BY 4.0","url":"https://creativecommons.org/licenses/by/4.0/"},"manuscriptAbstract":"\u003cp\u003eThe progression of multiple primary lung cancer (MPLC) involves complex changes in cell composition and metabolic function. Here, we performed scRNA sequencing of 167,397 cells from six patients with MPLC, combined with bulk whole-exome sequencing. We revealed that both naïve and memory T cells participate in the differentiation of CD8+ T cells. The terminal states of CD8+ T cells are exhausted T cells, which respond to stimuli and positively regulate cell death, and cytotoxic T cells, which are mainly implicated in the regulation of cytokine production. Multiple functional Tregs and naïve T cells contribute to the origin of CD4+ T cells. B cells, with two main functions, also play important roles in the immune response. We also uncovered the cellular metabolic activities that occur during tumor invasion. Positive regulation of blood vessel diameter has been observed in endothelial cells, while angiogenesis has not been found in early glandular neoplasia of the lung. The development of epithelial cells involves two functional states. In one state, cells respond to stress from the immune system, and in the other, some will undergo programmed cell death or enter the cell cycle due to the selective pressure that arises over the course of tumor development in synchronous MPLC. Our study showed the complete landscape of different dynamic cellular changes, which might reveal the key cellular mechanisms of MPLC and therefore provide new clues for the pathogenesis of tumors.\u003c/p\u003e","manuscriptTitle":"Dissecting transcriptional heterogeneity in multiple primary lung cancer by single cell RNA sequencing","msid":"","msnumber":"","nonDraftVersions":[{"code":1,"date":"2022-06-03 22:30:16","doi":"10.21203/rs.3.rs-1716190/v1","editorialEvents":[{"type":"communityComments","content":0}],"status":"published","journal":{"display":true,"email":"[email protected]","identity":"researchsquare","isNatureJournal":false,"hasQc":true,"allowDirectSubmit":true,"externalIdentity":"","sideBox":"","snPcode":"","submissionUrl":"/submission","title":"Research Square","twitterHandle":"researchsquare","acdcEnabled":true,"dfaEnabled":false,"editorialSystem":"","reportingPortfolio":"","inReviewEnabled":false,"inReviewRevisionsEnabled":true}}],"origin":"","ownerIdentity":"371b7f57-3d49-42f6-ab52-d50dff4a3830","owner":[],"postedDate":"June 3rd, 2022","published":true,"recentEditorialEvents":[],"rejectedJournal":[],"revision":"","amendment":"","status":"posted","subjectAreas":[],"tags":[],"updatedAt":"2022-06-03T22:30:24+00:00","versionOfRecord":[],"versionCreatedAt":"2022-06-03 22:30:16","video":"","vorDoi":"","vorDoiUrl":"","workflowStages":[]},"version":"v1","identity":"rs-1716190","journalConfig":"researchsquare"},"__N_SSP":true},"page":"/article/[identity]/[[...version]]","query":{"redirect":"/article/rs-1716190","identity":"rs-1716190","version":["v1"]},"buildId":"_2-kVJe1T_tPrBINL-cwx","isFallback":false,"isExperimentalCompile":false,"dynamicIds":[84888],"gssp":true,"scriptLoader":[]}