The tumor ecosystem of ACC was examined using scRNA-seq of digested living cells derived from ACC paracarcinoma and carcinoma tissues using a 10x Genomics-based platform ( Figure 1A ). A total of 49,948 cells from the parcarcinoma and carcinoma tissues of three ACC patients were captured using the sequencer and 42,714 cells were retained for downstream analysis after quality control filtering ( Supplementary Table 1 ). The results of the staining pathologic sections of three patients are shown in Supplementary Figure 1 . The clinicopathological characteristics of three patients are listed in Supplementary Table 2 .

Unsupervised clustering of the cells identified 13 cell types, including fibroblasts, epithelial cells, endothelial cells, basal cells, T or natural killer cells, smooth muscle cells, myeloid cells, muscle satellite cells, mast cells, skeletal muscle cells, lymphatic endothelial cells, Schwann cells, and an unknown cell cluster, which was a mixture of plasma cells, B cells, and fibroblasts ( Figure 1B and Supplementary Figure 2B ).

Uniform manifold approximation and projection (UMAP) visualization of the hallmark genes expressed in each cell subtype was performed ( Figure 1C ). The dot plots show the well-expressed marker genes in corresponding cell types ( Supplementary Figures 2A, C ). The complex tumor ecosystem comprised 40.1% epithelial cells, 27.9% fibroblasts, 11.3% T/NK cells, 10.8% endothelial cells, and a small number of other cells ( Figure 1D ). Notably, basal cells only existed in sample 0222 paracarcinoma tissue, whereas smooth muscle cells, muscle satellite cells, unknown cells, skeletal muscle cells, and Schwann cells only existed in sample 7420 paracarcinoma and carcinoma tissues ( Figure 1E ).

The proportion of cell types in the paracarcinoma and carcinoma tissues were examined. The proportion of epithelial cells was 9.3% in paracarcinoma, while in carcinoma was 40.1%. In contrast, the proportion of endothelial cells in paracarcinoma and carcinoma were 24.0% and 10.8%, respectively ( Figures 1D, E ). Considering ACC is a malignant tumor that originates from epithelial cells, the transcriptome characteristics of epithelial cells were further evaluated in the study.

The 9,685 epithelial cells present were categorized as three subtypes ( Figure 2A ): myoepithelial-like cells (MECs, ACTA2⁺/MYH11⁺/CNN1⁺), intercalated duct-like cells (Inter-Duct 1–7, KRT19+/AQP5+/KIT+), and duct-like cells (Duct 1 and Duct 2, KRT19+/AQP5-/KIT-). MUC5B, which is a mucinous acinar marker, was expressed at lower levels in a few subtypes, such as Inter-Duct 5. Moreover, MUC7, which is a serious acinar marker, was barely expressed in any subtypes ( Figure 2B ). In the paracarcinoma tissues, Inter-Duct 1, Inter-Duct 3, Duct 1, Inter-Duct 4, and MECs accounted for 34.7%, 21.9%, 15.6%, 12.2%, and 8.0% of the epithelial cells, respectively. In the carcinoma tissues, MECs, Inter-Duct 5, and Duct 2 accounted for 18.7%, 7.0%, and 3.4% of the epithelial cells, respectively ( Figure 2C ).

Each cluster was analyzed to identify differentially expressed genes. HES1 and ID4 were highly expressed in Inter–Duct 3 cells. CENPF, TOP2A, and MK167 were highly expressed in Inter–Duct 4 cells. ACTA2 and TP63 were highly expressed in MECs. SCGB3A1, MMP7, and ZG16B were highly expressed in Inter–Duct 5 cells. SRGN, CCL5, and CREM were highly expressed in Inter–Duct6 cells. PLVAP and PECAM1 were highly expressed in Inter–Duct 7 cells ( Figure 2D ).

Gene set variation analysis (GSVA) revealed the function of the different epithelial clusters. Notch signaling pathway, MYC targets, DNA replication process, oxidative phosphorylation, Wnt β-catenin, and Pi3k-Akt-mTOR signaling pathway were upregulated in Inter-Duct 3–4 cells. Duct 2 and MECs were involved in epithelial–mesenchymal transformation, angiogenesis, myogenesis, and apical junction. Inter-Duct 5–7 cells showed high expression of immune-related genes and were enriched for pathways including interferon α/γ response, allograft rejection, and IL6 JAK STAT3 signaling ( Figure 2E ).

These results suggested that each epithelial cluster was enriched for different signaling pathways. Inter-Duct 3–4 cells were enriched for activation of tumor progression pathways and Inter-Duct 5–7 cells may be involved in the immune response.

Differences in gene expression programs of three main epithelial cell clusters were performed by non-negative matrix decomposition. Hierarchical clustering identified five expression programs that varied within the Inter-Duct cells, including cell cycle, extracellular matrix organization, epidermal development, stress response, and regulation of neuron death. Three gene expression programs were found in MECs, including epithelial cell migration and proliferation, neurotransmitter transport, and muscle system process ( Figure 2F ). We identified the dominant gene sets in the different expression programs of three main epithelial cell clusters and performed GO (Gene Ontology) enrichment analysis on specific genesets to support above findings ( Supplementary Table 3 , Supplementary Figures 3, 4 ).

These results suggested that Inter-Duct cells may be directly involved in tumor development and myoepithelial cells play a contributory role in this process.

The potential molecular basis driving the distinct epithelial clusters was examined using SCENIC (Single-cell Regulatory Network Inference and Clustering) to identify the underlying regulatory network in epithelial cells.

Each cluster was driven by different transcription factors (TFs). Notably, coexpression of MYB and EN1 was observed in the upstream TFs of Inter-Duct 3–4 cells ( Figure 3A ). In addition, genes involved in DNA repair and cell cycle transition, such as E2F2, TFDP1, E2F3, BRCA1, POLE3, E2F1, EZH2, and RB1, showed aggregation ( Figure 3B ). Prediction of genes downstream of MYB and EN1 identified 10 coregulated downstream target genes, including COLEC12, ELAVL2, FRS2, IL17RD, ITGA6, LAMB1, LRIG1, NAV2, NCALD, and HOMER3 ( Figure 3C ), which are predominantly associated with the nervous system (15–21). Genemania (22) was used to predict possible functions coregulated by MYB and EN1, which are involved in cell–substrate junctions and neuronal guidance. ( Figure 3D )

Gene Ontology (GO) enrichment analysis of the TFs drove different expression programs. Duct 1 and Duct 2 cells were enriched for G1/S transition of the mitotic cell cycle, protein–DNA complex subunit organization, and intracellular receptor signaling pathway. Inter-Duct 3 cells were enriched for regulation of neuron death, histone deacetylation, and eyelid development in camera-type. Inter-Duct 4 were enriched for negative regulation of G0/G1 transition, cell cycle regulation, and transcription regulation involved in the G1/S transition of the mitotic cell cycle. Inter-Duct 5 cells were enriched for cellular response to type I interferon and transcription initiation from RNA polymerase II promoter. MECs were enriched for maintenance of the somatic cell population and steroid hormone-mediated signaling pathway ( Figure 3E ). These results indicate that Inter-Duct 3 cells were regulated by the upstream TFs, MYB and EN1.

The large-scale chromosome CNV status of all the cells from the paracarcinoma and carcinoma tissues was examined using infercnv ( Figure 4A ). Comparing with the endothelial cells and fibroblasts from the paracarcinoma, all epithelial cells showed complex CNV changes, suggesting a malignant tendency in most cells ( Supplementary Figure 5A ). Moreover, different clusters of epithelial cells from the paracarcinoma tissue underwent the massive copy number amplification in 6q, 8q, 12q, and 17p, and chromosomal deletions in 14q ( Figure 4C ). We assume that the majority of tumor tissue derived from epithelial cells with similar CNV pattern are malignant, while epithelial cells that do not conform to that pattern, which were mainly identified in the paracarcinoma tissue are considered to be premalignant (Pre-M) cells. ( Figure 4B ).

Examination of the paracarcinoma and carcinoma tissues revealed that 47.6% of epithelial cells in paracarcinoma tissues were pre-M cells ( Figure 4D ). Moreover, among the pre-M cells, 26.2% were Duct 1 cells, 40.3% were Inter-Duct 1 cells, 5.4% were Inter-Duct 3 cells, 1.4% were Inter-Duct 5 cells and 9.7% were MECs. In comparison, the proportion of Inter-Duct 3 cells in malignant cells was 19.9%, which was significantly higher than the pre-M cells 5.4%. In addition, nearly 17.3% malignant cells were MECs ( Figure 4E ). The percentage of the pre-M cells in the Duct cells, Inter-Duct cells, and MECs were 32.8%, 12.8%, and 8.9%, respectively. On the contrary, the proportion of malignant cells in each cell clusters were 67.2%, 87.2%, and 91.1%, respectively. ( Figure 4F ).

MYB and EN1 were identified as the upstream TFs and MYB was identified as the hallmark ACC gene. MYB/MYBL1/EN1 were only highly expressed in epithelial cells ( Supplementary Figure 5B ). MYB and EN1 were highly expressed in malignant Inter-Duct 3–7 cells, and MYBL1 was highly expressed in Inter-Duct 1 cells, Inter-Duct 2 cells and MECs. ( Supplementary Figure 5C ). Moreover, among three main epithelial cell types, MYB was highly expressed in malignant Inter-Duct cells, MYBL1 was highly expressed in malignant MECs, and EN1 was highly expressed in all of them ( Supplementary Figure 5E ). These findings were consistent with previous studies that reported that the expression of MYB and MYBL1 were mutually exclusive (11) ( Supplementary Figures 5C, D ).

Importantly, the stem cell gene score of pre-M cells was higher than that of malignant cells, reflecting the accuracy of our method ( Supplementary Figure 5F ).

Differentiation trajectory analysis was used to examine the evolution of epithelial cells. The nodes of pre-M cells with a higher stemness score were set as the root node of pseudo-time. Most of the pre-M cells were gathered initially and evolved into malignant cells by the end of the pseudo-time. Notably, MECs, which were used as an independent differentiation group, were not involved in these processes ( Figure 5A ).

The pseudo-time axis was divided into four stages according to its lower quartile, median, and upper quartile. Graph-autocorrelation analysis revealed coregulated genes in the modules within the four stages. MYB and MYBL1 were highly expressed in modules 5 and 13 ( Figure 5B , Supplementary Table 4 ), originated from the same root node, and showed increased expression along the pseudo-time trajectory, evolving toward two completely different branches ( Figure 5C ). Consistent with the above-mentioned results, during tumor evolution, gene MYB and MYBL1 showed opposite trends in expression ( Figure 5D ).

In addition, GO enrichment analysis of these two gene modules showed that the genes in Module 5 were involved in the rhythmic process, regulation of cellular component size, and dopaminergic neuron differentiation ( Figure 5E ). Furthermore, the genes in Module 13 were associated with diverse and complex biological processes, such as cellular response to external stimulus ( Figure 5F ). Taken together, these findings indicate that genes coregulated with MYB or MYBL1 may drive the two different evolutionary directions of ACC.

The tumor microenvironment is essential for the proliferation, migration, survival, anti-immune killing of malignant cells. Thus, we applied CellPhoneDB to infer cell-cell communication from combined expression of multi-subunit ligand-receptor complexes from scRNA-seq data. We found that malignant_inter_duct, Basal cells, endothelial cells, fibroblasts, myeloid cells, malignant_MEC were the dominant communication hubs ( Figure 6A ).

We classified cell communication pairs into 12 categories according to malignancy grade, predominant epithelial subset, and molecular origin. A significant upregulation in the expression of a number of receptors and ligands has been observed in malignant Inter-Duct and malignant MECs. Moreover, compared with preM cells, MIF-TNFRSF14, CD74-MIF and CD74-APP were generally highly up-regulated in the communications between malignant Inter-Duct cells or MECs and other cells.

MIF (Macrophage Migration Inhibitory Factor) plays an import dual role of pro-inflammatory and pro-oncogenic (23). TNFRSF14 is a member of tumor necrosis factor (TNF) receptor superfamily, functioning in activating inflammatory and inhibitory T-cell immune response (24). The up-regulation of MIF-CD74 were found in many cancers, such as cervical squamous cell carcinoma (25), hepatocellular carcinoma (26) and prostatic cancer (27). Studies have shown that MIF-CD74 may enhance the proliferation and inhibit the apoptosis of tumors by promoting angiogenesis of tumor microenvironment. In addition, it has been reported that MIF may be a novel prognostic marker for human oral squamous cell carcinoma (28). This suggested that inflammation and immunosuppression in tumor microenvironment may provide favorable conditions for ACC proliferation and invasion ( Figures 6B, C ).

The upregulated receptor–ligand genes in different groups were implicated in different Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways strongly associated with cancerization ( Figure 6D ). Overall, the enhanced communications between malignant and stromal cells in the tumor microenvironment further indicate the complexity of tumor behavior. The tumor microenvironment is essential for the proliferation, migration, survival, and anti-immune killing of malignant cells and the combined expression of multisubunit ligand–receptor complexes identified in the scRNA-seq data may infer cell–cell communication.

The dynamic changes in the homologous subpopulations of the Inter-Duct 1 and 7 cells from paracarcinoma and carcinoma tissues may be related to the different biological functions of each cell subpopulation in tumor progression.

GSVA showed that each Inter-Duct cell clusters were enriched for different signaling pathway, including Notch signaling pathway, MYC, and immune response-related pathways (31, 32). The results indicated that different homologous cell clusters had different biological functions.

We used SCENIC to identify the TFs for each epithelial cell cluster. GO enrichment analysis of the dominant TFs showed that they controlled different expression programs.

MYB and EN1, as upstream TFs, regulated downstream genes related to domain neural activity together (33), histone deacetylation, and other tumor processes (34, 35) in Inter-Duct 3–4 cells. This is the first time we identified EN1 as a TF in ACC, may regulate the activity of different cell types in scRNA-seq at the transcription level.

The clustering of TFs, POLE3 (DNA Polymerase Epsilon 3), TFDP1, RB1, and E2F transcription factor families (E2F3, E2F3, and EF7) were associated with cell cycle transition (36), suggesting possible disruption of the normal epithelial cell division cycle ( Figure 3B ). Duct 1 and Duct 2 cells shared many biological processes, which were mostly related to transcription. Inter-Duct 1 cells were enriched for stress response-related pathways, such as response to peptide and laminar shear stress. Inter-Duct 2 cells were associated with stem cell proliferation and differentiation. Inter-Duct 5 and Inter-Duct 6 cells were associated with immune processes, such as T cell activation and interferon response. MECs were enriched for maintenance of somatic cell populations and steroid hormone-mediated signaling pathways ( Figure 3E ).

Our findings contribute to further understanding of the heterogeneity of biological functions by different TFs. Furthermore, these can explain the neural invasion from a scRNA-seq perspective in ACC.

Large-scale CNV analysis validated an intermediate state between precancer and cancer in the transcriptional profile of samples from exocrine glandular malignancies, such as breast cancer, prostate cancer, and their paracarcinoma tissues (37–39).

In our study, the complex CNVs were found in the epithelial cells from the paracarcinoma tissues. We speculated that there was an abundant population of specialized cells with abnormal transcription. They will transform into malignant cells in the paracarcinoma tissues in the future but cannot be accurately determined at the cellular or protein levels. This may be relevant to the frequently positive incisal edge in ACC. We tried to identify and define this unique cell population as pre-M cells.

Sample 0222 and 0329 showed a long arm amplification at chromosomes 6 (40) and 8 (41), respectively. As previously reported (42), MYB and MYBL1 were located at 6q23 and 8q13.1, respectively. This suggested that the expression of these two genes may affect copy number changes at the chromosomal level. In addition, translocation fusion of MYB family genes with NFIB is the gold standard for diagnosis of ACC, in which the expression rate is about 65%–85% (43). Previous study reported that EN1 was a potential biomarker for worse prognosis in ACC (44). This may explain the association of EN1 with poor prognosis from the single cell transcriptome level.

The distribution of the pre-M and malignant cells in each cell cluster was observed.

We attempted to calculate the stemness gene score of pre-M and malignant cells as it may be useful to determine the potential for transforming cells from pre-malignant to malignant states. Assuming that about 95% of the epithelial cells in carcinoma tissues are malignant, we observed a relatively equal proportion of pre-M and malignant cells in paracarcinoma tissues. This may suggest that pre-M cells may transform into malignant cells. Comparing with the different cluster, more malignant cells were found in Inter-Duct 3, which activated tumor-associated signaling pathways by upstream TFs, MYB and EN1.

Subsequently, MYB family genes and EN1 were expressed in both pre-M and malignant Inter-Duct cells. MYB was previously shown to be highly expressed in only duct cells by immunohistochemistry and fluorescence in situ hybridization (45). Notably, the stemness gene score of pre-M was higher, suggesting its potential cancer transforming characteristics. These findings fully validate the pre-M cells in the paracarcinoma tissues. It can be a more accurate explanation for the false-negative incision edge and high recurrence rate from a single-cell transcriptome perspective in ACC.

The abnormal manifestation of complex CNVs verified the presence of specific pre-M cells in the paracarcinoma tissues and the expression of MYB verified its accuracy. We used scRNA-seq pseudo-time analysis to infer the differentiation trajectory of the epithelial cells and evolution of the cell subtypes to explain how the pre-M cells eventually become malignant throughout tumor progression.

The pseudo-time was divided into four stages. At the beginning of the time, a considerable number of cells were pre-M, which was mainly enriched in the Duct 1 cluster. As the timeline progressed, they gradually transformed into malignant cells and overlapped with Inter-Duct 3/4/6. Inter-Duct 3–4 cells were regulated by the TFs, MYB and EN1, and underwent activation of tumor-associated pathways. Myoepithelial cells were separated from the rest of the epithelial cell population and similarly transformed over time into a malignant cell population. This finding is fully consistent with the pathological features of ACC.

Graph-autocorrelation analysis revealed that MYB gene expression was progressively enhanced during the transition from stage 3 to 4, whereas MYBL1 showed high gene expression during the transition from stage1 to 2. UMAP (Uniform Manifold Approximation and Projection) analysis of the pseudo-time showed that MYB and MYBL1 showed exclusive expression at the end of the proposed time, which is consistent with the findings of previous studies (41, 46). GO enrichment analysis of these two gene modules showed that these two modules may drive two different evolutionary directions of ACC. This could explain the further development of ACC under the regulation of MYB homologous genes from a single-cell transcriptional perspective.

To sum up, the study is the first to report the use of scRNA-seq to examine ACC at the transcriptome level. We identified a special population of Inter-Duct cells. Inter-Duct 3–4 cells were coregulated by the upstream TFs, MYB, and EN1, which were enriched for tumor progression pathways. In addition, pre-M cells were found in paracarcinoma tissues, which was verified by the presence of MYB gene, and were highly expressed mainly in malignant Inter-Duct cells, including Inter-Duct 3–7. In addition, MYBL1 was highly expressed in malignant Inter-Duct and myoepithelial cells. Finally, pseudo-time analysis revealed that different cell clusters eventually transformed from a pre-M to malignant state in the ACC progression. Two modules containing MYB and MYBL1 drove different trajectories. Our findings further explain the high recurrence rate and unique biological characteristics of ACC.

The present study was approved by the Ethics Committee of Chinese People's Liberation Army (PLA) General Hospital and was performed according to the guidelines of the Declaration of Helsinki. Written informed consent was obtained from each patient prior to sample collection.

Preoperatively, the CT data of patients were imported into the Robotic-assisted navigation system. The surgeon could locate the tumor precisely by the robot. We removed the paracarcinoma and carcinoma tissues by enlargement. Samples were obtained from the paracarcinoma and carcinoma tissues of three ACC patients who underwent surgery at the Chinese PLA General Hospital Stomatology department. The primary foci were on the minor salivary gland. The sample 7420 and 0329 were both cribriform types, in addition sample 0222 was cribriform type in 2019 and recurrence to be solid type in 2021 in pathology.

The tissues from sample 0222 for scRNA-seq was incised for recurrence in 2021. And then we compared pathological sections of 0222 samples between 2019 and 2021 for further study. Each sample was carefully reviewed by two experienced pathologists to confirm the pathology. Single-cell data information of the six samples in shown in Supplementary Table 1 . The clinical data are summarized in Supplementary Table 2 .

Differentially expressed genes were identified using the FindAllMarkers function of the Seurat package for each subcluster. Specific genes were selected to serve as a basis for artificially defining cell populations as follows (52–58): fibroblasts (COL3A1, DCN, COL1A1, LUM, COLA2, COL6A2, FBN1), epithelial cells (EPCAM, KRT19, CLDN3 and KRT8), endothelial cells (PECAM1, ENG, VWF), basal cells (KRT5, KRT14, TP63), T/NK cells (GZMA, HCST), smooth muscle cells (ACTA2, MYH11, ACTG2, CNN1, CALD1, MCAM, TAGLN, PDGFRB, MYL9), myeloid cells (AIF1, CD163, LYZ), muscle satellite cells (PAX7, CD82, NCAM1, MYF5), mast cells (MS4A2, TPSB2, GATA2), skeletal muscle cells (ACTA1, NEB, MYL2), lymphatic endothelial cells (PDPN, PROX1, LYVE1), Schwann cells (NGFR, SOX10, GAP43, CDH19), plasma cells (JCHAIN, IGKC, IGHA1, IGHA2), T cells (CD3G, CD3D, CD3E), and B cells (MS4A1, BANK1, CD37).

Large-scale chromosomal copy number alterations within cells were detected using the R package inferCNV (v1.8.1). Immune and stromal cells from paracarcinoma tissues were used as presumptive “normal” cells as a reference, and their CNV scores were set as baseline. A previously described method was used to distinguish malignant cells from all epithelial cells (59). In brief, 1,600 fibroblasts and endothelial cells were stochastically picked from paracarcinoma samples, among which 1,000 were regarded as a reference and the remaining cells and epithelial cells were considered as an observation. CNV analysis was performed using inferCNV. When most of the spiked “normal” stromal cells were gathered in a specific cluster on a dendrogram, the other cells in this cluster were considered as “normal epithelial cells”. Correspondingly, cells not belonging to this cluster were identified as “malignant epithelial cells.” Since only 29 epithelial cells were identified as normal epithelial cells, accounting for only 0.003% of all epithelial cells, we concluded that all epithelial cells show malignant characteristics. Epithelial cells with malignant features were further subclassified. The median CNV estimates of all fibroblasts in the paracancerous samples were defined as the baseline, and the Euclidean distance between the CNV estimates of all observation cells and the above baseline was calculated. We determined the range of Euclidean distances (median ± 2 SD) by assuming that most epithelial cells (95%) in the carcinoma samples were “malignant” in the true sense, whereas cells whose distance from the baseline was not within this range were considered pre-M cells. The signature genes used to distinguish subtypes of epithelial cells or to calculate the score of stemness of different class of epithelial cells are as follows (60–62): intercalated duct cells (KIT or AQP5), myoepithelial cell genes (ACTA2, MYH11, CNN1), and stemness genes (ALDH1A1, CD44, PROM1, NANOG, KIT, NES, KLF4, CD55, ALCAM, NOTCH4, WNT7A, PDPN). Each cell was scored for stemness activity using the function AddModuleScore in Seurat with default settings. A re-clustered heatmap of CNV scores of all epithelial cells was plotted using pheatmap (v1.0.12).

The expression matrix and meta information of epithelial cells migrated from Seurat object were subsequently imported into R and analyzed using monocle3 (v1.0.0) (63) using the default parameters and a standard pipeline. After graph-autocorrelation analysis, the top 3,000 genes that varied with the trajectory were selected and collected into the gene module divided by pseudo-time. The graph of the change of key gene expression over the pseudo-time was drawn using the plot_genes_in_pseudo-time function.

The specific transcriptional drivers in each epithelial cell subtype were analyzed by using the R package, SCENIC (v1.2.4) (64), under the guidance of a tutorial from https://github.com/aertslab/SCENIC. The binary score matrix of intracellular regulator activities was processed using limma (v3.48.3) (65) to find TFs that were differentially expressed in different cell subpopulations (logFC > 0 & p.adjust < 0.05), and the results were displayed using pheatmap (v1.0.12). The interaction network construction of high confidence dominating transcription factor and their target genes in epithelial cells and final visualization were implemented using Cytoscape (v3.9.1) (66).

In order to predict the possible functions of genes co-regulated by MYB and EN1, we used an online tool called Genemania (http://genemania.org/), which can find other genes that are related to a set of input genes by using a large set of genome and proteome association data. The result of GO enrichment analysis of above genes can be download from the website directly.

Intercellular communication within the tumor microenvironment was predicted using the python package, CellPhoneDB (v3.0.0) (67), based on the expression of interacting ligand and receptor genes. The number of iterations for the statistical analysis was set as 1,000. Significant (p < 0.05) mean results were further processed as follows. Cell pairs containing Duct/InterDuct/MEC were selected and divided into “source” and “target” groups according to the sequence of which in cell pairs. The expression scores of receptor–ligand pairs were considered as “counts” and constructed as a Seurat object and the differentially expressed receptor–ligand pairs in each group were revealed using the FindAllMarkers function with default parameters. Receptor–ligand pairs were selected for plot as those with p.adjust value <0.05, average log2FC >1. Heatmaps and dot plots were plotted using CellPhoneDB.

All epithelial cells were extracted from six samples and divided into three subgroups and the expression matrix was then normalized and decomposed using python package cNMF (v1.3) (68). The high-quality expression programs in each sample were manually selected, the Pearson’s correlation coefficients between them were calculated, and those with coefficients >0.18 and clustered together were considered as characteristic programs. In addition, correlation coefficients >0.6 were corrected to 0.6. The parameters, methods, and code used for data processing were adapted from a GitHub tutorial (https://github.com/dylkot/cNMF).

GSVA was executed on hallmark gene sets, which were obtained from MSigDB (The Molecular Signatures Database) (69) and contained 50 well-defined pathways, using the R package, GSVA (v1.40.1) (70), with default settings. Differences between the pathway enrichment scores of the different groups was calculated by R package Limma (v3.48.3).

GO and KEGG term enrichment analysis were performed using R packages clusterProfiler (v4.0.5) and org.Hs.eg.db (v3.13.0) with default settings (71).

The software, methods, and thresholds used for the statistical analyses are detailed in the Materials and Methods section. Wilcoxon test was used to reveal the statistical differences in Figures 1E , 2C and 4B . Two-tailed t test was used to reveal the statistical difference in Supplementary Figures 5E, F . A p.adjust-value or q-value <0.05 was considered significant.