The tumor microenvironment fosters tumor cell survival and is a key player in tumor progression (1, 2). The recent development of single-cell RNA sequencing (scRNA-seq) has allowed a systemic view of the tumor microenvironment. Comparisons between tumor and normal tissues have revealed alterations in immune cells such as a decrease in NK cells and an increase in suppressive or regulatory T cells (3, 4). Tumor associated macrophages demonstrate high levels of heterogeneity (5), and have emerged as a target for cancer treatment (6). In most solid tumors, cancer associated fibroblasts, including myofibroblasts, populate the tumor microenvironment (7, 8). However, the timing and sequence of the events remain unclear due to the long period of tumor progression before diagnosis and sampling.

Mouse models serve as useful tools for understanding the cellular and molecular dynamics of human cancer (9). Even in mouse models, capturing the early stages of endogenous cancer can be challenging due to low penetrance and long latency (10). Transplant models using cell lines may offer a practical alternative. In our study, we used a lung orthotopic tumor transplant model using the Lewis lung carcinoma (LLC) mouse lung cancer cell line (11). The LLC cell line harbors a heterozygous Kras G12C mutation (12), a prevalent driver mutation in both human and mouse lung cancer (13). KRAS mutations cause tumor-promoting inflammation and immune escape (14), which allows this model to reflect immunosuppressive lung adenocarcinoma that requires novel therapeutic approaches. Moreover, we used a method of intracardiac injection to simulate a multi-organ metastatic scenario (15). Our focus is on lung transplants, aiming to characterize the array of cellular and molecular changes that occur during the early stages of tumor engraftment.

Overview of scRNA-seq on LLC orthotopic mouse transplant model

To characterize cellular alterations following LLC engraftment in the lung, we established an orthotopic transplant model using intracardiac injection (Fig. 1A) (16). We assessed the presence of LLC cells using an in vivo imaging system (IVIS) to measure luciferase activity and detected multi-organ signals (Supplementary Fig. 1A). The hematoxylin and eosin (H&E) staining of tissue sections demonstrated the formation of tumor masses in the lung and liver (Fig. 1B). Notably, LLC tumor cells fail to grow in brain parenchyma (17). Subsequently, we performed scRNA-seq on the lungs (Fig. 1A). In the IVIS images, we distinguished between tumor-negative (without luminescence) and tumor-positive (with luminescence) lungs. In flow cytometry data, we observed immune (CD45⁺), epithelial (CD326⁺) and endothelial (CD31⁺) cells, but failed to identify LLC tumor cells (Supplementary Fig. 1B). We anticipated that the scRNA-seq data would enable us to annotate diverse cell types and provide a more detailed characterization.

Cellular landscape analysis identifies LLC tumor cells with Myc and Kras G12C expression

As illustrated in Fig. 1A, scRNA-seq was conducted on a cell suspension obtained by the enzymatic dissociation of lung tissues. After quality control filtering based on the number of genes and the percentage of mitochondrial genes (Supplementary Fig. 2A), 5,001 cells from tumor-negative and 3,007 cells from tumor-positive samples were subjected to further analysis (Fig. 2A and Supplementary Fig. 2B). Data integration and clustering using the Seurat pipeline (18), resulted in the identification of 12 cell types within the uniform manifold approximation and projection (UMAP) space (Fig. 2A). These cell types included tentative LLC tumor cells, mesothelial cells, fibroblasts, endothelial cells, normal epithelial cells (including AT1, AT2, goblet/club, and ciliated cells), macrophages, monocytes, T/NK cells, and B cells, identified by canonical marker gene expression (Fig. 2B). The LLC tumor cell clusters were verified by the KRAS mutation in sequencing BAM files, as marked on the feature plot (Fig. 2C, upper panel). Sequence comparisons in the integrated genomics viewer (IGV) (Supplementary Fig. 2C) reflected the heterozygous Kras G12C mutation in the LLC cell line. Notably, the Myc gene was highly expressed in LLC tumor clusters (Fig. 2C, lower). The tumor-positive lung contained most of the LLC tumor and mesothelial cells (Fig. 2A, right). Other cell types, including normal epithelial cells, were more abundant in the tumor-negative lung. To estimate overall changes in the cellular landscape, we compared differentially expressed genes (DEGs) between the total cells from tumor-negative and tumor-positive lung samples (Fig. 2D). In the tumor-positive lung, top DEGs including Myc and Cxcl1 were those most enriched in the LLC tumor cell clusters. By comparison, Sftpc in the top DEGs of tumor-negative lung is a AT2-specific gene. These results are consistent with the clustering analysis, showing tumor or AT-2 cell enrichment in the tumor-positive or tumor-negative lungs, respectively. We validated expression levels of the Myc, Cxcl1, and Sftpc genes by qRT-PCR (Fig. 2E).

Sub-clustering analyses identify cell type-specific genes modulated during tumor engraftment

Next, we performed sub-clustering analyses for fibroblast, myeloid (monocyte and macrophage clusters), and T/NK cell populations, which were present abundantly in both tumor-negative and positive lungs. First, fibroblasts were subdivided into ten clusters (18 PCs, res 0.6) representing six distinct subsets: adventitial fibroblasts, lipo-fibroblasts, peribronchiolar fibroblasts, Adh7+ fibroblasts, pericytes, and smooth muscle cells (SMCs) based on canonical marker expression (Fig. 3A and Supplementary Fig. 3A) (19). Notably, adventitial fibroblasts (Pi16, Serfinf1) exhibited elevated Fbln2 and Timp1 gene expression in the tumor-positive lung (Fig. 3A, lower). Lipo-fibroblasts displayed more pronounced gene expression (Inmt, Limch1) in the tumor-negative lung, whereas SMCs expressed the Acta2 and Myh11 genes in the tumor-positive lung. In the proportion, lipo-fibroblasts comprised the largest portion, with a slightly higher ratio of adventitial fibroblasts and a lower ratio of pericytes in the tumor-positive lung compared to the tumor-negative lung (Supplementary Fig. 3A, right).

Second, the analysis of myeloid cells revealed the presence of seven distinct clusters (14 PCs, res 0.3) representing four types of macrophages (monocyte-derived, inflammatory, interstitial, and alveolar), two types of dendritic cells (CD209a+ and plasmacytoid DCs), and monocytes/neutrophils (Fig. 3B and Supplementary Fig. 3B) (20). Interstitial macrophages expressed higher levels of C1qa, C1qb, and C1qc genes in the tumor-negative lung and showed downregulation of MHC class II genes (H2-Aa, H2-Eb1) in the tumor positive lung (Fig. 3B, lower). In the cell counts per cluster, macrophages were the most abundant in both tumor-positive and tumor-negative lungs (Supplementary Fig. 3B, right).

Finally, the analysis of the T/NK subset revealed eight clusters (15 PCs, res 0.6) representing 6 subpopulations: CD4 T cells, including naïve-like and regulatory T cells (Tregs), CD8 T cells with naïve-like or cytotoxic signatures, proliferating cells mostly with CD8+ T cell identity, NK cells, innate lymphoid cells (ILC2), and mucosal-associated invariant T cells (MAIT) with canonical marker expression (Fig. 3C and Supplementary Fig. 3C) (21, 22). In the proliferating cell cluster, the expression of Ctla4 gene was more pronounced in the tumor-positive lung, whereas the expression of proliferating genes was eminent only in the tumor-negative lung (Fig. 3C, lower). Cytotoxic CD8 T cells and proliferating cells were found at lower proportions in the tumor-positive lung (Supplementary Fig. 3C, right).

Analysis of cellular interactions reveals molecular connections between the LLC tumor and lung microenvironment

To gain insight into the cellular and molecular interactions driving tumor engraftment, we performed a correlation analysis between tumor cell clusters and the constituents of the tumor microenvironment (Fig. 4A). The analysis revealed a high gene expression correlation between LLC tumor cell clusters 4 and 6. Furthermore, copy number variation (CNV) inference confirmed similar genomic patterns in these tumor clusters (Fig. 4B). Then, we conducted a CellPhoneDB analysis (23) to infer molecular interactions between various cell types. Visualized as a heatmap, the interaction count network suggested vigorous cellular interactions (Fig. 4C). With a specific focus on LLC tumor cells, we identified receptor-ligand pairs that could mediate interactions with the microenvironment (Fig. 4D). Among these potential molecular interactions, SPP1-CD44, SPP1-integrin-complex, FN1-integrin-complex, and COL-integrin-complex interactions were prominent between LLC tumor cells and other cell types. ICAM1-AREG and LAMP1-FAM3C interactions were significant between mesothelial, fibroblast, endothelial, or epithelial cells and LLC tumor cells. FAM3C-CLEC2D interaction was notable among fibroblast, endothelial, and immune cells. TYROBP-CD44, C5AR1-RPS19, and CD74-MIF interactions were most prominent between myeloid cells and LLC tumor cells (Fig. 4D). We confirmed MIF protein expression in the LLC tumor mass and CD74 in the surrounding cells (Fig. 4E), as well as reduced MHC class II expression in the macrophages from the tumor-positive sample (Fig. 4F). Collectively, these molecular alterations suggest strong tissue remodeling and immunosuppression activities, which play a crucial role in promoting tumor engraftment within the lung microenvironment. 

In this study, we used scRNA-seq to explore the genomic features of the tumor transplant and the associated lung microenvironment. While previous studies have described the lung orthotopic model (24, 25), ours is the first to report systemic assessment of the total cell populations covering the tumor-stromal-immune axis.

In the analysis, we identified the LLC cells as a distinct cluster in the tumor-positive sample. The cluster exhibited high expression levels of Myc, Cxcl1, and Rhox5 genes, which are associated with tumor progression (26-28), yet lacked the typical marker of lung adenocarcinoma NK2-related homeobox transcription factor Nkx2-1 (29). In both mouse lung cancer models and human patients, Nkx2-1 is often lost in metastatic stages (2, 12). Notably, the LLC cluster showed high expression levels of mesenchymal genes such as Vim, Col3a1, and Mif corroborating the high metastatic potential of LLC cells (30).

To investigate the stromal-immune axis facilitating LLC engraftment, we employed two approaches, DEG and cellular network analysis. In the DEG analysis, a subset of fibroblasts displayed upregulation of Fbln2 and Timp1, which are genes involved in tissue remodeling and repair (31, 32). On the other hand, myeloid subsets showed a downregulation of MHC class II genes, suggesting dysfunctional antigen presentation (33, 34). Moreover, T/NK cell subsets exhibited a more regulatory phenotype but reduced proliferation-associated gene expression in CD8+ T cells. These findings indicate that tissue remodeling and immune suppression play crucial roles in creating a tumor-permissive environment. In the cellular network analysis, we found that mesenchymal genes in the LLC tumor probably engage in robust molecular interaction with various stromal cells, via pathways such as SPP1-CD44, SPP1-integrin-complex, FN1-integrin-complex, and COL-integrin-complex. Many of these alterations and interactions align with reports in human cancer data (4, 7), underscoring their contributions to tissue remodeling and cancer progression (35-37). Additional validation of these cellular and molecular interactions will be achieved through further analyses using spatial transcriptomics (38).

In summary, our study has revealed that an LLC transplant triggers significant changes in the tumor microenvironment. These changes are marked by the upregulation of tissue remodeling and repair genes in the stroma and the downregulation of MHC class II genes in macrophages. Notably, there is a more regulatory, less proliferative phenotype in the T/NK cell compartment, suggesting impaired antigen presentation and reduced T cell infiltration. These insights are crucial for understanding and potentially treating low-immunogenic tumors, particularly in the early stages.

Cell lines

The LL/2-Luc2 (CRL-1642-LUC2) Lewis lung carcinoma cell line was purchased from the American Type Culture Collection (ATCC, Manassas, VA). Cells were cultured in DMEM containing 10% FBS and 2 μg/ml blasticidin S hydrochloride in a 5% CO₂ incubator at 37°C.

Orthotopic mouse tumor model

Female C57BL/6 mice were purchased from Orient Bio Inc (Seongnam, Republic of Korea). LL/2-Luc2 cells (1 × 10⁴ to 1 × 10⁵ in 100 μl phosphate-buffered saline) were injected into the left ventricle of 8- to 9-week-old mice (16). Mice were monitored 2 to 3 times a week starting on the 5th day. Each mouse received an intraperitoneal injection of 100 μl of IVISbrite D-luciferin (Perkin Elmer, Waltham, MA), and bioluminescence images were acquired 10 minutes later using IVIS Lumina XRMS (Perkin Elmer).

All procedures involving animal research were in accordance with the Laboratory Animals Welfare Act, the Guide for the Care and Use of Laboratory Animals, and the Guidelines and Policies for Rodent experiment provided by the Institutional Animal Care and Use Committee (IACUC) in the School of Medicine, The Catholic University of Korea (Approval number: CUMS-2021-0363-06).

Sample collection for scRNA-seq

The lung tissue was minced, dissociated for 40 minutes at 37°C using the enzymes D, R, and A in RPMI1640 medium of the Tumor Dissociation Kit, mouse (Miltenyi Biotec, Bergisch Gladbach, Germany), and filtered through a 70-μm filter. Red blood cells were removed using RBC lysis solution (QIAGEN, Hilden, Germany).

Flow cytometry

Cells were incubated with TruStain FcX anti-mouse CD16/32 antibody (BioLegend, San Diego, CA) for 10 min and then stained with primary antibodies (Biolegend) for 40 minutes at 4°C before data analysis by FACS-Canto (BD Biosciences, San Jose, CA) and FlowJo software v.10.10.0.

Antibodies for total cell profiling: APC anti-CD45 (clone 30-F11); PE anti-CD90.2 (53-2.1); PerCP-Cy5.5 anti-CD31 (390); BV421 anti-CD326 (G8.8).

Antibodies for MHC class II expression by macrophages: PE Anti-CD45 (30-F11); FITC Anti-F4/80 (BM8); APC Anti-MHC II (M5/114.15.2).

Single-cell RNA sequencing

Single cell suspensions were processed using Chromium Next GEM Single Cell V(D)J Reagent Kits v2 (10× Genomics, Pleasanton, CA) with a cell recovery target of 7,000 per library according to the manufacturer’s instructions. Libraries were sequenced on an Illumina Hiseq X for 5’ GEX library, with 100PE.

scRNA-seq data analysis

The data was processed according to the Cell Ranger pipeline (version 5.0.0, 10× Genomics) and aligned to the mouse reference genome (mm10-2020-A). Gene expression profiles and the annotation information of single cells were compiled using the Seurat package (version 4) (18). Cells with less than 1,000 and more than 5,000 unique features or more than 10% of mitochondrial genes were removed from the merged data. Filtered cells were normalized with the LogNormalize function and highly variable features were identified using the FindVariableFeatures function. For batch correction, we used the Harmony package and performed principal component analysis (PCA) with the RunPCA function. Dimensionality reduction and clustering using UMAP were performed using the RunUMAP, Findneighbors, and FindClusters functions with 9 PCs, resolution = 0.3. Subset analyses were performed using a similar method with PC and resolution specified in the main text.

Kras G12C mutation detection

For the detection of the Kras G12C mutation, we used a method recommended by 10X Genomics, which involves creating a BAM file filtered by a specific barcode list (39). The nucleotide substitution from C to A marking the amino acid codon change from glycine (GGT) to cysteine (TGT) in the 12th amino acid position was confirmed in the IGV.

Copy number variation inference

To distinguish LLC tumors, we used the chromosomal gene expression pattern to infer CNVs. We established a reference group comprising all cell types except for LLC tumors and then created an inferCNV object (https://github.com/broadinstitute/inferCNV). We predicted copy number alteration (CNA) with the following settings: cutoff = 0.1, HMM_type = ‘i6’, tumor_subcluster_partition_method = “leiden”.

Cell–cell interaction analysis

CellPhoneDB (version 4.1.0) was used to infer cell-cell interactions (23). The UMI count matrix assigned to all cell types was loaded as an input onto the CellPhoneDB. Mouse gene symbols were converted to human ensemble gene IDs using the biomaRt package (version 2.46.3). Cells expressing specific ligands and receptors were set to more than 25%.

Immunohistochemistry staining

Lung tissues were fixed in 10% neutral buffered formalin, embedded in paraffin, and cut at 4-μm thickness. IHC staining was performed with Anti-CD74 (ab289885 from Abcam, Cambridge, UK, 1:2,000) and Anti-c-Myc antibodies (ab32072, Abcam, 1:50), followed by ImmPRESS HRP Universal Antibody (Horse Anti-Mouse/Rabbit IgG) Polymer Detection Kit (Vector Laboratories). Signals were detected with diaminobenzidine (DAB) and counterstained with hematoxylin. Stained slides were observed using an automatic digital slide scanner (Pannoramic SCAN II, 3DHISTECH) and analyzed with CaseViewer Software.

Quantitative RT-PCR

Total RNA was extracted using RNeasy Mini Kit (QIAGEN). RNA concentration was determined using a NanoDrop 2000/2000c Spectrophotometer (Thermo Scientific). Quantitative RT-PCR (qRT-PCR) was performed using the SuperScript First-Strand Synthesis System (Invitrogen) and IQ Sybr Green SuperMix (BioRad). Primer Sequences: Sftpc forward 5’-GTCCTCGTTGTCGTGGTGATTG-3’, reverse 5’-AAGGTAGCGATGGTGTCTGCTC-3’; Myc forward 5’-TCGCTGCTGTCCTCCGAGTCC-3’, reverse 5’-GGTTTGCCTCTTCTCCACAGAC-3’; Cxcl1 forward 5’-TCCAGAGCTTGAAGGTGTTGCC-3’, reverse 5’-AACCAAGGGAGCTTCAGGGTCA-3’.

Data availability

Raw and processed scRNA-seq data have been deposited at the Gene Expression Omnibus (GEO) with accession number GSE246710.

This work was supported by the National Research Foundation of Korea (NRF) grant funded by the Korea government (MSIT) (2022R1A2C1091451, RS-2023-00220840).

The authors have no conflicting interests.