Deciphering the potential pathogeny of rare tracheal adenoid cystic carcinoma by single-cell RNA-sequencing

Introduction

Primary tracheal carcinoma is a rare sub-type of cancer, with an annual incidence of 0.1 to 0.26 cases per 100,000 persons (1). Squamous cell carcinoma (SCC) and adenoid cystic carcinoma (ACC) account for approximately 75% and 15% of primary tracheal carcinoma cases, respectively (1). Traditionally, tracheal adenoid cystic carcinoma (TACC) has been classified as a sub-type of salivary gland-type tumors, and was thought to originate from submucosal glands of the trachea (2,3). It usually arises between the ages of 40–60 years, and its gender distribution is equal (2,3). Due to its indolent behaviors of slow growth and late distant metastasis, TACC patients have a 5-year survival rate of 52% and a 10-year survival rate of 29%, both of which are better than those of other types of lung cancer (1). However, its symptoms, including dyspnea, cough, wheezing, and hemoptysis, not only seriously affect the quality of life of patients, but may also lead to misdiagnoses of asthma or bronchitis, resulting in treatment delays of months and even years (4).

At present, surgery is the first primary treatment recommendation for TACC. Nevertheless, only 42–57% of patients achieve complete resection (1-4), as the occult invasion of submucosal and perineural spaces is very common, making it difficult to determine the surgical margin properly. Radiotherapy offers an alternative treatment option for patients with unresectable tumors or positive surgical margins, but the findings of previous studies are controversial, and the benefits of radiotherapy remain uncertain (5). Systemic treatments like chemotherapy and target therapy are not as widely used in the treatment of TACC as in the treatment of other lung cancers. Despite the effect of some agents had been demonstrated in sporadic reports, the evidence from randomized and prospective clinical trials is still lacking. Moreover, no standard chemotherapy regimen has been established (6). Further, to date, no data on immunotherapy in TACC are available. This lack of evidence creates a dilemma in the treatment of TACC.

The molecular profile of TACC is not yet known, hindering the investigation of novel targets and therapies. Research has been conducted on the well-known driver genes of lung cancer in TACC by immunohistochemistry (IHC), fluorescence in situ hybridization, and polymerase chain reaction (PCR), but no “hot” mutations have been identified (7,8). Using high-throughput next generation sequencing technology, research has confirmed that “hot” mutations in well-known driver genes are absent in TACC; however, these reports lack overlapping evidence (9-12). This may be because data originating from bulk sequencing are the average result of all tested cells, and thus cannot reflect the exact molecular characteristics of TACC.

Single-cell RNA-sequencing (RNA-seq) is a novel and emerging tool with great power to recognize the heterogeneity of tumors. Its ability to reveal interactions between tumor cells and the immune microenvironment has been successfully demonstrated in previous studies (13,14). Therefore, single-cell RNA-seq was employed to investigate TACC in the present study. We sought to determine the pathogeny of TACC by deciphering the molecular characteristics of this rarely seen cancer in typical young patients, and to verify the findings in patients by fluorescent multiplex immunohistochemistry (mIHC). We present this article in accordance with the MDAR reporting checklist (available at https://tlcr.amegroups.com/article/view/10.21037/tlcr-2025-832/rc).

Methods

Patients

In July 2021, a 32-year-old male, who had been suffering from continuous hemoptysis for about one month, was admitted to The First Medical Center of Chinese PLA General Hospital. He had suffered from the same but less serious symptoms one year before, and had recovered quickly after receiving treatments such as anti-inflammatory and hemostatic agents. Other than a smoking habit of 20 cigarettes each day before his second hemoptysis, he reported no other significant past medical history. A chest computed tomography (CT) scan indicated a space-occupying lesion in the middle lobe of the right lung. Bronchoscopy was subsequently performed, revealing a neoplasm blocking the bronchial lumen, with obvious congestion and edema in the adjacent mucous layer. The pathological diagnosis of the biopsy was ACC.

The patient then underwent thoracoscopic surgery, and a tumor (40 mm × 25 mm × 15 mm) was successfully removed. A small part of the tumor tissue and its adjacent normal tissue were sent for single-cell sequencing. The other part of the surgical tissue was sent to the pathological department, and the diagnosis was confirmed to be consistent with the biopsy results. Specifically, the whole layer of the bronchial wall and nerve tissue had been invaded by the tumor, but the surgical margin and adjacent lymph nodes were all negative (Figure 1).

Figure 1 The results of the bronchoscopy and pathological diagnosis. The bronchoscopy before surgery (A,C) revealed a neoplasm blocked in the bronchial lumen, with obvious congestion and edema in the adjacent mucous layer. A bronchoscopy in the same position after surgery (B,D) showed that the neoplasm had been successfully removed. The pathological diagnosis after surgery (E,F) confirmed that the neoplasm was an adenoid cystic carcinoma invading the whole layer of the bronchial wall and nerve tissue, but the surgical margin and adjacent lymph nodes were all negative [H&E staining; ×40 (E) and ×100 (F)]. H&E, hematoxylin and eosin.

The study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments. The study was approved by the Ethics and Scientific Committee of the Chinese PLA General Hospital (No. S2021-216-01), and informed consent was obtained from all the patient (the study workflow is shown in Figure 2A).

Figure 2 Workflow of the present study and the integrated RNA-seq results. As the workflow (A) shows, the tumor tissue, adjacent tissue, and blood of one patient were obtained for RNA-seq. The RNA-seq results were analyzed using the bioinformatics method and confirmed by fluorescent multiplex immunohistochemistry. All the sequenced cells could be distinctly divided into immune, epithelial, and stromal clusters (B) based on the RNA expression of well-recognized markers (E). The distribution and percentage of these cells in different samples is shown in (C) and (D), respectively. RNA-seq, RNA-sequencing; UMAP, Uniform Manifold Approximation and Projection.

Single-cell isolation

Fresh tumor tissue (15 mm × 5 mm × 5 mm) and adjacent normal tissue (6 mm × 8 mm × 5 mm) were removed from the surgical tissue and immediately immersed in MACS^®Tissue Storage Solution (Cat# 130-100-008, Miltenyi Biotec, Bergisch Gladbach, Germany). Additionally, 5 mL of venous blood was collected in an ethylene diamine tetraacetic acid (EDTA) anticoagulant vacuum tube before surgery. The tissue and blood samples were transported to laboratory at proper temperature and processed within 1 hour.

The tissue samples were dissociated using a Tumor Dissociation Kit (Cat# 130-095-929, Miltenyi Biotec) in accordance with the manufacturer’s instructions, and then transferred to a MACS SmartStrainer (70 µm, Cat# 30-098-462, Miltenyi Biotec) and placed in a 50-mL tube. The SmartStrainer was washed with 20 mL of Roswell Park Memorial Institute (RPMI) 1640 (Cat# 11875-093, Gibco, New York, USA), and the eluate was centrifuged at 300 ×g for 7 minutes. The re-suspended samples and the blood mononuclear cells that had been isolated by ficoll density gradient centrifugation (Cat# P9011, Solarbio, Beijing, China) were both treated with a 10-fold volume of 1×Red Blood Cell Lysis Solution (Cat# 30-098-462, Miltenyi Biotec). The viability, concentration, and clumping rate of the cells was determined using Acridine Orange/Propidium Iodide staining agent (Cat# F23001, Logos Biosystems, Gyeonggi-do, South Korea) on a LUNA-FL^TM Automated Fluorescence Cell Counter (Cat# L20001, Logos Biosystems). The dead cells were removed using the Dead Cell Removal Kit (Cat# 130-090-101, Miltenyi Biotec) as necessary in accordance with the manufacturer’s instructions. Before the next procedure was performed, the samples were centrifuged and re-suspended with an adequate volume of 1× phosphate-buffered saline, 0.04% bovine serum albumin, and 1 U/µL of ribonuclease inhibitor to obtain a final concentration of 700–1,200 viable cells/µL.

Single-cell capture and cDNA library construction

The samples were prepared using the Single-Cell 3’ Kit v3.1 and Chip G on Chromium Single-Cell Controller (10x Genomics, California, USA) and ETC811 Thermal cycler (Eastwin, Beijing, China) in accordance with the manufacturers’ instructions. In brief, appropriate volumes of nuclease-free water, Master Mix, and single-cell suspension for targeting 10,000 single cells were mixed and added to the Chip G to generate the gel bead-in-emulsions (GEMs) with gel beads and partitioning oil on the Chromium Single-Cell Controller. The GEMs were dispensed into the tube strip on ice and sent for reverse transcription within 1 hour. Recovery Agent was then added to the products to break the single-cell droplets, and the first-strand complementary DNA (cDNA) was isolated and cleaned with Dynabeads Cleanup Mix. Afterwards, the resulting samples underwent the steps of cDNA amplification, fragmentation, adaptor ligation, and sample index PCR to generate the DNA libraries, with a cleanup procedure using the SPRI select Reagent Kit performed at the end of each step. Finally, the libraries were assessed using the LabChip GX Touch Nucleic Acid Analyzer (Caliper Lifesciences, California, USA) and KIT iQ SYBR GRN (Bio-Rad, California, USA) to ensure the sample quality met the requirements for sequencing.

Single-cell sequencing and bioinformatics analysis

The libraries were sequenced on the Illumina NovaSeq6000 platform (Illumina, San Diego, CA, USA) using the following parameters: paired-end 150, dual indexing, and a depth of no less than 23,000 read pairs per cell. Adopting the method used in previous reports, the cellular barcodes were demultiplexed by Cell Ranger software pipeline (10x Genomics, version 7.0.1), and the reads were mapped to the genome and transcriptome using the STAR aligner. The down-sample reads were used to generate normalized aggregate data across samples, producing a matrix of gene counts versus cells.

The single-cell RNA data was filtered, processed, and analyzed using the Seurat software (version 4.3.0.1) with R (version 4.2.2). First, cells with fewer than 200 or more than 100,000 genes, or with fewer than 1,000 unique molecular identifiers were filtered out, thus removing the low-quality cells and likely multiplet captures. Meanwhile, the low-quality cells for which >20% of the counts belonged to mitochondrial genes were also discarded. Second, library size normalization was performed using the NormalizeData function to obtain the normalized count. Specifically, the gene expression measurements for each cell by total expression were normalized by “LogNormalize”, the results were log transformed after been multiplied by a scaling factor (10,000 by default). Third, the top variable genes across single cells were identified using the method described by Macosko et al. (15). In brief, the most variable genes were selected using the FindVariableGenes function (mean.function = FastExpMean, and dispersion.function = FastLogVMR); A principal component analysis (PCA) was conducted to reduce dimensionality using the RunPCA function. Graph-based clustering was applied to group cells based on their gene expression profiles using the FindClusters function.

Cells were visualized by a two-dimensional Uniform Manifold Approximation and Projection (UMAP) algorithm using the RunUMAP function. Marker genes for each cluster were identified using the FindAllMarkers function (test.use = presto). Differentially expressed genes (DEGs) were detected using the FindMarkers function (test.use = presto). The thresholds for significant differential expression were set as a P value <0.05 and a |log₂fold change| >0.25. Both Gene Ontology enrichment and Kyoto Encyclopedia of Genes and Genomes pathway enrichment analyses of the DEGs were performed using R based on the hypergeometric distribution.

To infer the genomic copy number structure, InferCNV (version 1.8.1) was used according to the developer’s suggestions with standard parameters (https://github.com/broadinstitute/inferCNV). For the minimum average read counts per gene among reference cells, a cut-off value of 0.1 was used. Meanwhile, alveolar epithelial type 1 (AT1) and alveolar epithelial type 2 (AT2) cells were used as internal reference controls. Further, the Monocle2 (16) package (version 2.8.0) was used to analyze single-cell trajectories to discover the cell state transitions with the following parameters: average expression R0.125, num_cells_expressed R10, and qval <0.01 (differentialGeneTest function). The trajectory was visualized as a two-dimensional t-distributed stochastic neighbor embedding (tSNE) graph, and dynamical expression heatmaps were generated using the plot_pseudotime_heatmap function.

Subsequently, a protein-protein interaction (PPI) network was constructed using the STRING online database (https://string-db.org/). The PPI pairs with selected larger scores were used to construct the PPI network. Finally, CellChat (17) (version 0.5) was employed to conduct the cell-cell communication analysis. In brief, cell annotation labels and normalized gene expression levels generated through the Seurat workflow were integrated as inputs for CellChat. The expression levels of ligand and receptor genes in each cell population were projected onto the PPI network. Permutation testing of randomized network connections was used to identify significant source-target network connections.

Result verification by fluorescent mIHC

In addition to the patient that underwent single-cell sequencing, the formalin-fixed and paraffin-embedded tissues of three other TACC patients (2 males, 1 female, aged 64, 73 and 42 years, respectively) were examined as follows: The sections were cut into 5 µm-thick slides, which were deparaffinized in xylene for 30 minutes and sequentially rehydrated in absolute ethyl alcohol for 5 minutes (twice), 95% ethyl alcohol for 5 minutes, and 75% ethyl alcohol for 2 minutes. The slides were then washed three times with distilled water. Heat-induced epitope retrieval was performed using a microwave oven, during which the slides were immersed in boiling EDTA buffer (ZLI-9079, ZSBio, Beijing, China) for 15 minutes. Blocking was conducted using Antibody Diluent/Block from Alpha X Bio (Beijing, China). The mIHC experiments were performed with the AlphaXPainter^® X30 (Alpha X Bio) and analyzed according to one panel, in which the following primary antibodies were used: panel 1: CD8 (ZA0508, ZSGB-Bio), CD16a (ab227665, Abcam, UK), CD56 (ZM0057, ZSGB-Bio), and PANCK (ZM0069, ZSGB-Bio); Panel 2: CD68 (ZM0060, ZSGB-Bio), Vimentin (HA721174, Huabio, Hangzhou, China), MIF (87501S, CST, Massachusetts, USA), CD74 (ET1702-51, Huabio), CXCR4 (ab181020, Abcam, Cambridge, UK), and PANCK (ZM0069, ZSGB-Bio). All the primary antibodies were incubated for 1 hour at 37 ℃. Subsequently, the slides were incubated with Alpha X Polymer HRP Ms+Rb (Alpha X Bio) for 10 minutes at 37 ℃.

Visualization was performed using the Alpha X 7-Color IHC Kit (catalog no. AXT37100031; Alpha X Bio). Following each staining cycle, heat-induced epitope retrieval was carried out to remove all antibodies, including both primary and secondary antibodies. The slides were counterstained with diamidino-phenyl-indole (DAPI) for 5 minutes and mounted in Antifade Mounting Medium (catalog no. I0052; NobleRyder, Beijing, China). The images were scanned using a Zeiss Axioscan 7 (Zeiss, Oberkochen, Germany) and analyzed with HALO software (version 3.6; Indica Labs, New Mexico, USA). Based on the hematoxylin and eosin staining results, pathologists reviewed the slides and delineated the tumor and its adjacent area. Using HALO, the positivity of the cytoplasm/membrane and/or nucleus of each cell was measured separately, and the number of positive cells and average cell intensity for each marker was obtained. The number of positive cells for a specific marker per square millimeter (i.e., the density) was calculated as the number of positive cells for a specific marker in the relevant area divided by the area of that region. Further, the average distances (ADs) between the Vimentin⁺PANCK⁺ cells and the other cells were calculated.

Statistical analysis

The P values for the cell density and ADs were calculated using the rank-sum test. The analyses were carried out using SPSS Statistics 24.0 (IBM Corporation, USA), and P values <0.05 were considered statistically significant.

Results

The integrated analysis of the RNA-seq data

In total, 33,337 cells met the quality control requirements, including 13,720 cells from the blood samples, 5,123 cells from the adjacent tissue samples, and 14,494 cells from the tumor tissue samples. As the violin plot in Figure S1 shows, the overall quality of the data from the three samples was satisfactory. As Figure 2 and Figure S1 show, these cells could be distinctly divided into immune, epithelial, and stromal clusters based on the RNA expression of well-recognized markers like PTPRC, EPCAM, KRT18, KRT19, PECAM1, and DCN. The distribution of these cells in different samples is shown in Figure 2C. For each type of sample, the constitution of the cell types is shown in Figure 2D. Notably, the blood samples comprised only immune cells, while the adjacent tissue comprised 19.7% epithelial cells, 72.5% immune cells, and 7.8% stromal cells; and the tumor tissue comprised 42.1% epithelial cells, 46.6% immune cells, and 11.3% stromal cells.

The construction of epithelial cells in the tumor and adjacent tissues

Based on the RNA expression of well-recognized markers (as shown in Figure S2), the epithelial cells in the tumor tissue and adjacent tissue were divided into 12 subgroups (Figure 3A). Among them, the AT1 and AT2 cells were the main cell types in the adjacent tissue (accounting for more than 60% of the cells) and were not found in the tumor tissue. Conversely, the main cell types in the tumor tissue were Basal_II, Basal_III, and Basal_V (accounting for more than 60% of the cells), but these cell types were not found in the adjacent tissue (Figure 3A).

Figure 3 The RNA-seq results for the epithelial cells. The epithelial cells in the tumor and adjacent tissues were divided into 12 subgroups (A). The CNV analysis revealed the genomes of the Basal_I, Basal_II, Basal_III, Basal_IV, Basal_V, ciliated, and ciliated-like cells were unstable (B). A further analysis of these subgroups in the tumor tissue revealed different degrees of activation in the tumor-related signaling pathways, indicating that these cells possessed a malignant phenotype (C). The pseudo-time analysis showed that the ciliated cells were at the beginning of the curve (D). A further analysis revealed eight genes with abnormal expression in the ciliated cells had close interactions with most of the cancer driver genes, and thus warrant further investigation (E). CNV, copy number variation; RNA-seq, RNA-sequencing; UMAP, Uniform Manifold Approximation and Projection.

The copy number variation (CNV) analysis took the AT1 and AT2 cells as the reference cells and revealed that the genomes of the Basal_I, Basal_II, Basal_III, Basal_IV, Basal_V, ciliated, and ciliated-like cells were unstable (Figure 3B and Figure S2D). A further analysis of these subgroups in the tumor tissue observed a different degree of activation in the tumor-related signaling pathways, indicating that these cells possessed a malignant phenotype (Figure 3C).

Subsequently, the pseudo-time analysis of the above subgroups demonstrated that the ciliated cells were at the beginning of the curve, followed by the Basal_I cells (Figure 3D). The other cells were divided into two distinct branches, one of which mainly comprised Basal_V and ciliated-like cells, and the other of which mainly comprised Basal_II, Basal_III, and Basal_IV cells (Figure 3D). The expression of the commonly seen stem-cell markers also confirmed the primitive characteristics of the ciliated cells (Figure S2E), indicating that they might be linked to the origin of the disease. Further, the gene expression of the ciliated cells in the tumor tissue was compared with the same type of cells in the adjacent tissue to identify the potential biomarkers playing crucial roles in this process. As Figure S2F shows, the comparison identified 20 and seven significantly up- and down-regulated genes, respectively. A further analysis revealed that eight of these genes (MUC16, SERPINB3, KTR17, KTR19, KRT14, IGFBP2, SAA2, and MB) interacted closely with most of the cancer driver genes and thus may warrant further investigation (Figure 3E).

The construction of immunocytes in the blood, and tumor and adjacent tissues

As Figure S3 shows, the immunocytes were divided into seven subgroups according to the RNA expression of well-recognized markers; no significant difference in their proportions in the blood, tumor tissue, and adjacent tissue were observed. In each type of sample, T cells and natural killer (NK) cells, which are responsible for target cell destruction, represented the majority of the cells. Based on the RNA expression of well-recognized markers (Figure S3C), these T cells and NK cells were further subdivided into 12 subgroups as shown in Figure 4A. Among them, the proportion of CD8 effector memory cells was significantly higher in the tumor tissue than the blood and adjacent tissue. Additionally, the proportion of CD8 experience cells was similar in the tumor tissue and adjacent tissue, but was higher in the tumor tissue than the blood. These data indicate that the local tissue was undergoing immune activation in response to the tumor. However, the proportion of cytotoxic NK cells and recently activated CD8 effector memory cells, which reflect the immediate immune response status, was significantly lower in the tumor tissue than the blood and adjacent tissue, indicating their inability to prevent tumor occurrence. Notably, the expression of the cytotoxic-related genes was lower in the NK cytotoxic cells from the tumor tissue than those from the blood, with FCGR3A (the encoded protein also known as CD16a) showing the most significant difference (Figure 4B). Meanwhile, the pseudo-time analysis of four subgroups of CD8 T cells (Figure 4C) revealed that these cells were mostly concentrated at the end of the curve in the tumor, suggesting they had finished or were undergoing an immunological reaction.

Figure 4 The RNA-seq results for the T cells and NK cells. The T cells and NK cells that made up the majority of the immunocytes were subdivided into 12 subgroups (A). Among them, the proportion of cytotoxic NK cells and recently activated CD8 effector memory cells was significantly lower in the tumor tissue than the blood and adjacent tissue. Notably, the expression of the cytotoxicity-related genes was lower in the cytotoxic NK cells of the tumor tissue than the blood, with FCGR3A showing the most significant difference (B). The pseudo-time analysis of the four subgroups of CD8 T cells revealed that these cells were mostly concentrated at the end of the curve in the tumor tissue (C). NK, natural killer; RNA-seq, RNA-sequencing; UMAP, Uniform Manifold Approximation and Projection.

The construction of stromal cells in the tumor and adjacent tissues

The stromal cells in the tumor and adjacent tissues were divided into five subgroups (Figure 5) according to the RNA expression of well-recognized markers as shown in Figure S4. Among them, lymphatic endothelial cells were rarely seen in the tumor tissue compared with the adjacent tissue, indicating the tumor cells may have less chance to metastasize through the lymph nodes. Conversely, the proportion of endothelial cells was slightly higher in the tumor tissue than the adjacent tissue, suggesting that the tumor was undergoing active angiogenesis. Notably, a subgroup of mesothelial cells that should belong to epithelial cells were identified as stromal cells in the tumor tissue. The CNV analysis in which endothelial and lymphatic endothelial cells served as the references revealed that the genomes of these mesothelial cells were unstable (Figure 5B). A further analysis confirmed these mesothelial cells were distributed along the pseudo-time curve of the tumor cells (Figure S4B) and expressed a higher level of genes that were related to invasion and metastasis (Figure S4C). Thus, it was speculated that these cells were tumor cells that had undergone epithelial-mesenchymal transition (EMT), and might play significant roles in the recurrence and metastasis of TACC.

Figure 5 The RNA-seq results for the stromal cells. The stromal cells were divided into five subgroups, with lymphatic endothelial cells nearly absent and mesothelial cells exclusively present in the tumor tissue (A). The copy number variation analysis revealed that the genomes of these mesothelial cells were unstable (B). The crosstalk analysis between different cell types demonstrated that these mesothelial cells had significant relationships with myeloid cells, followed by endothelial cells (C). The gene expression analysis revealed that MIF-(CD74 + CXCR4) and MIF-(CD74 + CD44) played significant roles in the crosstalk between the mesothelial and myeloid cells (D). NK, natural killer; RNA-seq, RNA-sequencing; UMAP, Uniform Manifold Approximation and Projection.

Meanwhile, the crosstalk analysis between different cell types demonstrated that these mesothelial cells had significant relationships with myeloid cells, followed by endothelial cells (Figure 5C). More specifically, these mesothelial cells mainly crosstalked with macrophage 1 cells and NK-like monocytes (Figure S4D). The gene expression analysis revealed MIF-(CD74 + CXCR4) and MIF-(CD74 + CD44) played significant roles in the crosstalk between the mesothelial and myeloid cells (Figure 5D). Similar results were obtained when the tumor cells were analyzed, further confirming the speculation that these mesothelial cells were tumor cells had undergone EMT (Figure S4E-S4G).

No matter which mesothelial cells or tumor cells were analyzed with the sub-types of myeloid cells, the importance of MIF-(CD74 + CXCR4) and MIF-(CD74 + CD44) was repeatedly emphasized (Figure S4H-S4I). For the mesothelial cells and endothelial cells, the interaction between ACKA1 and chemokines (e.g., CCL2, CXCL2, CXCL3, and CXCL8) should be noted, as the results suggested that the endothelial cells were recruited by mesothelial cells to participate in angiogenesis (Figure 5D). Further analyses confirmed the above markers were highly expressed in the corresponding cells (Figure S4J-S4K).

The mIHC results for preliminary verification

Due to the availability of the antibodies and the limited number of TACC patients, the mIHC verification was only performed on some critical biomarkers in four samples (including those of the patient who undergone single-cell sequencing, patient #3 in Figure 6). Figure 6A,6B show a merged picture of biomarkers in panels 1 and 2, respectively. Although the statistical analysis was challenging, the results from the four patients were relatively consistent. The most noticeable result was the density of the CD56⁺CD16a⁺ cells (representing cytotoxic NK cells), CD68⁺ cells (representing macrophages), and CD68⁺CD74⁺ cells. For the CD56⁺CD16a⁺ cells, the density was lower in the tumor tissue than the adjacent tissue (Figure 6E), which is consistent with the results in Figure 4A,4B. For the CD68⁺ cells and CD68⁺CD74⁺ cells, the density of both was significantly higher in the tumor tissue than the adjacent tissue (Figure 6F,6G).

Figure 6 The mIHC result for preliminary verification. The representative mIHC results of biomarkers in panels 1 and 2 are shown in (A) and (B), respectively. Taking the CD56⁺CD16a⁺ cells as an example, the density was calculated as the number of positive cells in the relevant area (C). Additionally, the AD between the Vimentin⁺PANCK⁺ cells and other cells was calculated using the method shown in (D). The biomarkers with a significant difference between the tumor and adjacent tissue are shown in the box plots, including the density of CD56⁺CD16a⁺ cells (E), CD68⁺ cells (F), and CD68⁺CD74⁺ cells (G), as well as the AD between the Vimentin⁺PANCK⁺ cells and CD68⁺CD74⁺ cells (H). AD, average distance; mIHC, multiplex immunohistochemistry; UMAP, Uniform Manifold Approximation and Projection.

More importantly, the AD of the Vimentin⁺PANCK⁺ cells, which were considered tumor cells that had undergone EMT, to the CD68⁺CD74⁺ cells was shorter in the tumor tissue than the adjacent tissue (Figure 6H). The above data confirmed the findings in Figure 5 and Figure S4 that CD74 might play significant roles in the crosstalk between macrophages and mesothelial cells. In relation to the other biomarkers, the density of the CD8⁺, CD68⁺CXCR4⁺, PANCK⁺MIF⁺, and Vimentin⁺PANCK⁺MIF⁺ cells, as well as the AD from the Vimentin⁺PANCK⁺ cells to CD68⁺ and CD68⁺CXCR4⁺ cells, were also explored, but no significant results were found (data not provided).

Discussion

TCAA is a rare type of cancer that is often misdiagnosed as asthma or bronchitis, leading to treatment delays of months or even years (1-4). Due to the limited understanding of its pathogenesis and molecular characteristics, a highly effective therapy for this disease has yet to be established (5,6). Thus, the present study employed single-cell RNA-seq to examine the molecular characteristics of TACC. Notably, our work revealed three main findings.

First, the pseudo-time analysis revealed TACC originates from ciliated cells that display stem-cell properties. This finding contradicts the traditional view that TACC originates from the mucous gland (2,3), but it remains plausible. It is well recognized that normal ciliated cells in the trachea are located on the luminal surface, with dense cilia growing from their top side. Through the rapid and rhythmic swing of the cilia, the mucus and the foreign matters in the trachea (e.g., dust, bacteria, virus, and cigarette tar) can be transported to the pharynx and coughed out (18). As the first defense line of the lung, the ciliated cells come into contact with these harmful substances directly, enhancing the possibility of cancerization.

Given that the patient in present study used to smoke frequently (about 20 cigarettes per day), his pathogeny may involve smoke-induced changes in the gene expression of ciliated cells. The comparison of the ciliated cells in the tumor and adjacent tissues showed that the expression of the eight genes requires further investigation (Figure S2F and Figure 3E). Among them, KTR14, KTR17, and KRT19 are known to be involved in the formation and stabilization of the cytoskeleton, the abnormality of which may indicate that the morphology or rigidity of these cells is changing (19). MUC16-encoded cancer antigen-125 is a biomarker that is widely used in tumor screening, diagnosis, and evaluation (20). SERPINB3-encoded human SCC antigen 1 is known to regulate the immune reaction to tumor cells by acting as the inhibitor of papain-like cysteine proteases (21). The protein encoded by IGFBP2 is a member of the insulin-like growth factor (IGF)-binding protein family, and it can extend the half-life of IGF, and exhibits dual functions by inhibiting or stimulating the effect of IGF (22). SAA2-encoded serum amyloid A2 has been reported to plays significant roles in the regulation of immunity and inflammation (23). MB-encoded myoglobin is a protein in charge of oxygen storage and distribution, and may have potential implications in cancer (24). Interestingly, the aforementioned genes do not overlap with the genetic events of TACC reported in previous studies, including MYB-NFIB fusion, NOTCH1 mutations, and EGFR/c-KIT signaling (7-12). This discrepancy may be attributed to our use of single-cell RNA-seq specifically on the ciliated cells (from which TACC might originate); previous reports have used bulk sequencing to analyze all cell populations. However, the specific mechanisms of the aforementioned genes in the development of TACC need to be further elucidated. The integration of findings from previous studies will help us better understand the pathogenesis of TACC, and provide an important basis for the effective prevention and precision treatment of this disease.

Second, the suppression of immune cells in tumor tissue might lead to the expansion of TACC. The immune cells capable of directly killing tumor cells mainly include NK cells and cytotoxic T cells (25,26). We found that the functions of these two types of cells in TACC were both weakened to some extent. The quantity and the expression of the cytotoxicity-related genes in the NK cells were both lower in the tumor tissue than the peripheral blood and adjacent tissue. Conversely, most of the cytotoxic T cells were located at the end of the differentiation curve. Notably, the number of recently activated CD8 effector memory cells was significantly lower in the tumor tissue than the blood and adjacent tissue. These findings indicate that the immune response was activated in the tumor tissue but did not exert its appropriate functions to eliminate the tumor cells. According to our results, immune cell therapy like NK cells or chimeric antigen receptor-T (CAR-T) cell transfusions (27) may benefit patients with TACC, but this needs to be verified in future clinical trials.

Third, a distinct subgroup of mesothelial cells disguised as stromal cells was identified in the present study, which has not been reported in previous research. These cells had an unstable genome, were distributed along the pseudo-time curve of tumor cells, and expressed higher levels of genes related to invasion and metastasis. Therefore, it was speculated that they were tumor cells that had undergone EMT, and might play significant roles in the recurrence and metastasis of TACC. Further analysis revealed they had close interactions with macrophage-1 and endothelial cells, indicating that they had participated in immune reaction and angiogenesis, probably via the MIF-(CD74 + CXCR 4 or CD44) and CXCL2-ACKR1 pathways, respectively.

Notably, the above findings were confirmed when using epithelial tumor cells in the analysis and verified by the mIHC method, suggesting that these molecular alterations may serve as drug targets for TACC treatment. Fortunately, such drugs are already in development or approaching clinical application. For example, some small molecule inhibitors [e.g., ISO-1 for MIF (28) and AMD3100 for CXCR4 (29)], monoclonal antibodies [e.g., RG7356 for CD44 (30), and PF-06747143 for CXCR4 (31)], and antibody-drug conjugates [e.g., STRO-001 for CD74 (32)] have exerted promising effects in many kinds of tumors. Among them, research should focus on CD74, as the mIHC results showed that the CD68⁺CD74⁺ cells had a significantly higher density and a closer distance to the Vimentin⁺PANCK⁺ cells in the tumor tissue than the adjacent tissue (Figure 6G,6H), indicating that it could serve as a target or biomarker in treatment (33). If these drugs are validated in future clinical trials, they could address the current therapeutic gap for TACC by offering targeted treatment options for patients.

Collectively, this study puts forward a number of targeted prevention and treatment recommendations for patients with TACC. First, long-term exposure to harmful gases containing substances such as cigarette smoke, lampblack, and dust should be avoided to minimize adverse stimulation to tracheal ciliated cells—a key measure for reducing TACC susceptibility. On the onset of warning symptoms (e.g., chest tightness and hemoptysis), prompt bronchoscopy is strongly recommended to rule out TACC, addressing the critical challenge of delayed or misdiagnosis for this rare disease. Second, for patients with confirmed TACC, NK or CAR-T cell therapy, either as monotherapy or in combination with conventional treatments (e.g., surgery and radiotherapy), is advised. This strategy aims to enhance tumor clearance efficacy and mitigate the risk of tumor recurrence and distant metastasis, overcoming the limitations of traditional TACC therapies. Finally, to expand TACC’s therapeutic arsenal, preclinical studies using validated TACC models and large-scale multi-center clinical trials in TACC patients should be conducted. These investigations should focus on evaluating the real-world efficacy and safety of drugs targeting key molecular pathways (MIF, CD74, CD44, and CXCR4) closely linked to TACC cell proliferation, invasion, and immune evasion.

While the clinical significance of these recommendations awaits further verification through rigorous prospective studies, we anticipate they will contribute to improving TACC patient care. In essence, these recommendations are not merely clinical suggestions but a roadmap to redefine TACC care by reducing disease risk, and innovating treatments and advancing research. Their successful validation could not only improve TACC patient outcomes but could also serve as a model for addressing unmet needs in other rare thoracic malignancies.

The present study had two main limitations. First, single-cell sequencing was performed on only one patient, and validation experiments were conducted on a total of four patients (including the sequenced patient). This constraint primarily stems from the relative rarity of TACC, coupled with the strict requirement for an adequate volume of fresh samples to ensure reliable sequencing results. Consequently, identifying eligible patients for inclusion has proven challenging in clinical practice. To our knowledge, aside from the present study, there is only one other case report documenting the single-cell sequencing of TACC (34). However, the sampling site in that report was the throat, and its analysis focused exclusively on intratumoral T cell receptors, making it less comprehensive than the multi-dimensional analysis conducted herein. Given this limited sample size, the conclusions of our study should be interpreted with caution and require further validation in larger, well-powered cohorts.

The second limitation is the unavailability of a suitable TACC cell line, which precluded the completion of relevant functional mechanism research. Notably, organoid culture technology has advanced significantly in recent years (35); thus, the establishment of a qualified TACC organoid bank could provide valuable support for both mechanism exploration and drug development. This avenue will be a key direction for future research by our team.

Conclusions

This study demonstrated that TACC originates from ciliated epithelial cells. The abnormal gene expression triggered by external factors, along with the dysfunction of NK and T cells, are pivotal factors in the genesis and progression of TACC. Notably, the MIF-(CD74 + CXCR4) and MIF-(CD74 + CD44) signaling pathways are likely to be of critical importance in the recurrence and metastasis of TACC. Therefore, the development of drugs targeting these pathways holds great promise for achieving breakthroughs in the treatment of this type of tumor.

Acknowledgments

We sincerely acknowledge the technical assistance from Xiaochen Pan, Li Zhou, and Huchao Xu (EASYRESEARCH, Beijing, China) for single cell sequencing library and bioinformatics analysis, as well as Zuofu Peng, Huanrong Shang and Yun Jia (Alpha X Bio, Beijing, China) for mIHC verification.

Footnote

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

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

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

Funding: The study was supported by National Natural Science Foundation of China (No. 82241230).

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

Ethical Statement: The authors are accountable for all aspects of the work in ensuring that questions related to the accuracy or integrity of any part of the work are appropriately investigated and resolved. The study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments. The study was approved by the Ethics and Scientific Committee of the Chinese PLA General Hospital (No. S2021-216-01) and informed consent was taken from all the patients.

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

References

1. Ran J, Qu G, Chen X, et al. Clinical features, treatment and outcomes in patients with tracheal adenoid cystic carcinoma: a systematic literature review. Radiat Oncol 2021;16:38. [Crossref] [PubMed]
2. Maziak DE. Biology of Adenoid Cystic Carcinoma of the Tracheobronchial Tree and Principles of Management. Thorac Surg Clin 2018;28:145-8. [Crossref] [PubMed]
3. Falk N, Weissferdt A, Kalhor N, et al. Primary Pulmonary Salivary Gland-type Tumors: A Review and Update. Adv Anat Pathol 2016;23:13-23. [Crossref] [PubMed]
4. Cantù G. Adenoid cystic carcinoma. An indolent but aggressive tumour. Part B: treatment and prognosis. Acta Otorhinolaryngol Ital 2021;41:296-307. [Crossref] [PubMed]
5. Yang Y, Ran J, Wang Y, et al. Intensity modulated radiation therapy may improve survival for tracheal-bronchial adenoid cystic carcinoma: A retrospective study of 133 cases. Lung Cancer 2021;157:116-23. [Crossref] [PubMed]
6. Chae YK, Chung SY, Davis AA, et al. Adenoid cystic carcinoma: current therapy and potential therapeutic advances based on genomic profiling. Oncotarget 2015;6:37117-34. [Crossref] [PubMed]
7. Macarenco RS, Uphoff TS, Gilmer HF, et al. Salivary gland-type lung carcinomas: an EGFR immunohistochemical, molecular genetic, and mutational analysis study. Mod Pathol 2008;21:1168-75. [Crossref] [PubMed]
8. Huo Z, Wu H, Li S, et al. Molecular genetic studies on EGFR, KRAS, BRAF, ALK, PIK3CA, PDGFRA, and DDR2 in primary pulmonary adenoid cystic carcinoma. Diagn Pathol 2015;10:161. [Crossref] [PubMed]
9. Li M, Zhao BR, Liu SQ, et al. Mutational landscape and clonal diversity of pulmonary adenoid cystic carcinoma. Cancer Biol Ther 2018;19:898-903. [Crossref] [PubMed]
10. Pei J, Flieder DB, Patchefsky A, et al. Detecting MYB and MYBL1 fusion genes in tracheobronchial adenoid cystic carcinoma by targeted RNA-sequencing. Mod Pathol 2019;32:1416-20. [Crossref] [PubMed]
11. Ho AS, Ochoa A, Jayakumaran G, et al. Genetic hallmarks of recurrent/metastatic adenoid cystic carcinoma. J Clin Invest 2019;129:4276-89. [Crossref] [PubMed]
12. Wang F, Xie X, Song M, et al. Tumor immune microenvironment and mutational analysis of tracheal adenoid cystic carcinoma. Ann Transl Med 2020;8:750. [Crossref] [PubMed]
13. Zhang Y, Wang D, Peng M, et al. Single-cell RNA sequencing in cancer research. J Exp Clin Cancer Res 2021;40:81. [Crossref] [PubMed]
14. Tirosh I, Suva ML. Cancer cell states: Lessons from ten years of single-cell RNA-sequencing of human tumors. Cancer Cell 2024;42:1497-506. [Crossref] [PubMed]
15. Macosko EZ, Basu A, Satija R, et al. Highly Parallel Genome-wide Expression Profiling of Individual Cells Using Nanoliter Droplets. Cell 2015;161:1202-14. [Crossref] [PubMed]
16. Qiu X, Mao Q, Tang Y, et al. Reversed graph embedding resolves complex single-cell trajectories. Nat Methods 2017;14:979-82. [Crossref] [PubMed]
17. Jin S, Guerrero-Juarez CF, Zhang L, et al. Inference and analysis of cell-cell communication using CellChat. Nat Commun 2021;12:1088. [Crossref] [PubMed]
18. Fujisawa T, Tanaka Y, Ikegami K. Analysis of motility and mucociliary function of tracheal epithelial cilia. Methods Cell Biol 2023;176:159-80. [Crossref] [PubMed]
19. Dmello C, Srivastava SS, Tiwari R, et al. Multifaceted role of keratins in epithelial cell differentiation and transformation. J Biosci 2019;44:33.
20. Song Y, Yuan M, Wang G. Update value and clinical application of MUC16 (cancer antigen 125). Expert Opin Ther Targets 2023;27:745-56. [Crossref] [PubMed]
21. Cagnin S, Pontisso P, Martini A. SerpinB3: A Multifaceted Player in Health and Disease-Review and Future Perspectives. Cancers (Basel) 2024;16:2579. [Crossref] [PubMed]
22. Zhang B, Hong CQ, Luo YH, et al. Prognostic value of IGFBP2 in various cancers: a systematic review and meta-analysis. Cancer Med 2022;11:3035-47. [Crossref] [PubMed]
23. Sack GH Jr. Serum Amyloid A (SAA) Proteins. Subcell Biochem 2020;94:421-36. [Crossref] [PubMed]
24. Elkholi IE, Elsherbiny ME, Emara M. Myoglobin: From physiological roles to potential implications in cancer. Biochim Biophys Acta Rev Cancer 2022;1877:188706. [Crossref] [PubMed]
25. Chiossone L, Dumas PY, Vienne M, et al. Natural killer cells and other innate lymphoid cells in cancer. Nat Rev Immunol 2018;18:671-88. [Crossref] [PubMed]
26. Iwahori K. Cytotoxic CD8(+) Lymphocytes in the Tumor Microenvironment. Adv Exp Med Biol 2020;1224:53-62. [Crossref] [PubMed]
27. Han J, Zhang B, Zheng S, et al. The Progress and Prospects of Immune Cell Therapy for the Treatment of Cancer. Cell Transplant 2024;33:9636897241231892. [Crossref] [PubMed]
28. Cheng B, Wang Q, Song Y, et al. MIF inhibitor, ISO-1, attenuates human pancreatic cancer cell proliferation, migration and invasion in vitro, and suppresses xenograft tumour growth in vivo. Sci Rep 2020;10:6741. [Crossref] [PubMed]
29. Xie XF, Wu NQ, Wu JF, et al. CXCR4 inhibitor, AMD3100, down-regulates PARP1 expression and Synergizes with olaparib causing severe DNA damage in BRCA-proficient triple-negative breast cancer. Cancer Lett 2022;551:215944. [Crossref] [PubMed]
30. Menke-van der Houven van Oordt CW. First-in-human phase I clinical trial of RG7356, an anti-CD44 humanized antibody, in patients with advanced, CD44-expressing solid tumors. Oncotarget 2016;7:80046-58. [Crossref] [PubMed]
31. Liu SH, Gu Y, Pascual B, et al. A novel CXCR4 antagonist IgG1 antibody (PF-06747143) for the treatment of hematologic malignancies. Blood Adv 2017;1:1088-100. [Crossref] [PubMed]
32. Abrahams CL, Li X, Embry M, et al. Targeting CD74 in multiple myeloma with the novel, site-specific antibody-drug conjugate STRO-001. Oncotarget 2018;9:37700-14. [Crossref] [PubMed]
33. Fey RM, Nichols RA, Tran TT, et al. MIF and CD74 as Emerging Biomarkers for Immune Checkpoint Blockade Therapy. Cancers (Basel) 2024;16:1773. [Crossref] [PubMed]
34. Ye W, Clark EA, Sheng Q, et al. Primary tracheal adenoid cystic carcinoma: A case report and analysis of the tumor immune microenvironment using single cell RNA sequencing. Head Neck 2024;46:E91-8. [Crossref] [PubMed]
35. Yan HHN, Chan AS, Lai FP, et al. Organoid cultures for cancer modeling. Cell Stem Cell 2023;30:917-37. [Crossref] [PubMed]

(English Language Editor: L. Huleatt)