Single-cell data revealed the regulatory mechanism of TNK cell heterogeneity in liver metastasis from gastric cancer.

Introduction
Gastric cancer (GC) is one of the most frequent malignancies worldwide and the third major cause of cancer-related mortality [1]. Cancers spread through the lymphatic system, blood or directly to other organs [2], and liver metastasis is one of the most common forms of GC metastasis [3]. The presence of liver metastasis often indicates a quick progression of GC, with a 3-year overall survival of lower than 10% [4, 5]. The treatment of GC with liver metastasis is often highly challenging because GC is biologically heterogeneous with a variety of genetic mutations [6, 7]. At present, treatment and prognosis of GC has not been greatly improved, therefore revealing the molecular mechanisms of liver metastasis from GC is crucial for the prevention of the cancer and its treatment.
Tumor microenvironment (TME) is a complex ecological niche with a dynamic network of extracellular matrix, tumor-associated fibroblasts, immune cells, etc [8]. Understanding the interactions between tumor cells and the components of tumor immune microenvironment (TIME) will facilitate the exploration of the mechanisms and therapeutic targets related to GC with liver metastasis. For example, the most abundant cells in the TME are cancer-associated fibroblasts (CAFs), which are capable of secreting various biological mediators including vascular endothelial growth factor (VEGF), chemokines, TGF-β, and stromal cell-derived factors [9]. Tumor cells nutritionally supported by CAFs exhibit greater invasiveness compared to those lacking nutrition, moreover, tumor cells can also be reversely transformed or induced to increase CAFs, thereby generating a positive feedback loop of bi-directional signaling between the tumor cells and CAFs [10]. Neutrophils are also abundant in the TME of GC, and many patients with advanced GC exhibit neutrophilia, while patients with liver metastases from GC and a low neutrophil/lymphocyte ratio normally have a relatively longer survival [11, 12]. Moreover, neutrophils can secrete various pro-inflammatory factors to activate ERK pathway in GC cells, inducing epithelial-mesenchymal transition (EMT) and enhancing GC cell invasion and metastasis [13]. Recent study employed the analysis of scRNA-seq data to reveal a highly proliferative cell subpopulation related to cell cycle pathways in GC, providing guidance for personalized treatment of GC [14]. Similarly, based on the relevant scRNA-seq data, another study revealed the critical role of tumor endothelial cells in GC progression and developed a prognostic model to help predict the immune response of GC patients [15]. Hence, analysis on TIME of GC through single-cell transcriptomics could help reveal the mechanisms of GC with liver metastasis, contributing to targeted therapeutic treatment of the cancer.
In this study, we classified the cellular subpopulations that mediated the liver metastasis from GC at the single-cell level and determined the downstream pathways through which these subpopulation cells regulated GC progression. Subsequently the interactions between key cell subpopulations that mediated GC liver metastasis and other immune cells were analyzed. The present discoveries provided guidance for understanding the molecular mechanisms of liver metastasis from GC and improving clinical treatment strategies.

Materials and methods

Collection and preprocessing of scRNA-seq data for GC
The GSE163558 dataset containing three GC samples and two GC samples with liver metastasis was downloaded from the GEO (https://www.ncbi.nlm.nih.gov/geo/) database. The scRNA-seq data were preprocessed by reading the data with the Read10X function using the Seurat package to obtain cells with less than 10% of mitochondrial genes and the number of expressing genes between 200 and 6000 [16]. The SCTransform function was used for normalization. After PCA downscaling, the batch effect between samples (max.iter.harmony = 20, lambda = 0.5) was removed applying the harmony package [17]. Next, we performed UMAP downscaling according to the first 20 principal components and clustered the cells using the FindNeighbors and FindClusters functions. The CellMarker2.0 database (http://bio-bigdata.hrbmu.edu.cn/CellMarker) provided marker genes to annotate each cell type [18].

Functional enrichment analysis
The FindMarkers was sed to calculate differentially expressed genes (DEGs) between subgroups with p < 0.05, |logFC| > 1 as the screening criteria. Subsequently, differentially enriched biological process (BP) and KEGG pathways between the two groups were visualized by used the “gseGO()” and “gseKEGG()” functions in the clusterprofiler package, respectively. Finally, the GseaVis package was applied to show the genes in the differentially enriched pathways [19].

SCENIC analysis
Cellular heterogeneity in tissues is a result of the differences in cellular transcriptional states, while the specificity of transcriptional states is in turn determined and maintained by transcription factors (TFs)-dominated gene regulatory networks (GRNs). Thus, analysis on single-cell GRNs could help delve into the biological significance behind cellular heterogeneity. SCENIC [20] is a highly reliable algorithm developed specifically for single-cell analysis to identify TF-dominated GRNs. Following the official tutorial [5], the GENIE3 method was used to calculate the potential target genes of each TF and the top10perTarget method was used to construct the TF regulatory network. The AUCell function was used to calculate the degree of regulon activity in each cell after identifying highly significant TF-target gene relationship pairs. A regulon is a gene set of TFs and their directly regulated target genes. Based on the expression value of the gene, the activity of each regulon in individual cells was scored by the SCENIC package, with a higher score indicating greater activation of the gene set.

Cellular communication analysis
Cell-to-cell interactions between GC samples and GC samples with liver metastasis were analyzed applying the CellChat software package [21]. The createCellChat function was used to construct objects, and identifyOverExpressedGenes and identifyOverExpressedInteractions were used to identify ligand–receptor pairs overexpressed in the cellular subpopulations. The expression values of the ligand–receptor pairs were mapped to protein interaction networks using the projectData function, and the probability of ligand–receptor pair interactions between cellular subpopulations was inferred with the computeCommunProb function. The probability of the pathways locating in each ligand–receptor pair was calculated by the computeCommunProbPathway. Finally, the netVisual_bubble function was used to visualize the results.

Cell culture and transfection
Dulbecco’s Modified Eagle Medium (Gibco, 11965-092) added with 1% antibiotics (Gibco, 15070-063) and 10% fetal bovine serum (Gibco, 26140-095) were utilized for culturing human gastric mucosa cell line GES-1 purchased from SUNNCELL (Wuhan, China) and GC cell line AGS (BNCC338141) purchased from BNCC (Beijing, China) at 37 °C and 5% CO2. Negative control (NC) and JUNB siRNA (5ʹ-CTTCTACCACGACGACTCATAC-3ʹ) (Segun, China) were transfected with Lipofectamine 2000 (Invitrogen, USA).

QRT-PCR
cDNA was synthesized from RNA using the Qiagen One-Step RT-PCR kit (Qiagen Gmbh, Germany) to perform qRT-PCR. An ABI 7500 system (Thermo Fisher Scientific, USA) and SYBR Green were used for amplification. The experimental procedure was referred to previous studies reported [22]. The primer sequence of JUNB was Forward sequence: 5′-CGATCTGCACAAGATGAACCACG-3′, Reverse sequence: 5ʹ-CTGCTGAGGTTGGTGTAAACGG-3′, the Forward sequence of FOS was 5ʹ-GCCTCTCTTACTACCACTCACC-3ʹ, Reverse sequence: 5′-AGATGGCAGTGACCGTGGGAAT-3′ and the primer sequence of GAPDH was 5ʹ-GTCTCCTCTGACTTCAACAGCG-3ʹ, Reverse sequence: 5′-ACCACCCTGTTGCTGTAGCCAA-3′.

Wound healing assay
The GC cell line were seeded into a 6-well plate to grow until they covered the entire bottom. Scratches were created vertically with a 200 µL pipette tip. The cells were washed by phosphate-buffered saline twice and imaged with an inverted microscope at 0 h and 48 h after scratching. The rate of wound healing was calculated as [(0 h width − 48 h width)/0 h width] × 100%. The experiment was conducted in triplicate.

Transwell analysis
Cells were seeded into the upper Transwell chamber (Corning, USA) containing serum-free medium, while the lower chamber was filled with DMEM medium containing 10% FBS. After incubation at 37 °C for 1 day, cells remaining on the upper chamber were eliminated and the migrated cells were fixed by 4% paraformaldehyde and dyed with crystal violet for 30 min. Finally, a microscope was utilized to observe and count the cells.

Statistical analysis
Differences in continuous variables between two groups were compared by the Student’s t-test. All the calculations were performed in R language (version 4.3.1) and GrpahPad Prism 8.0. P < 0.05 was considered as statistically significant.

Results

Cell landscape of GC and GC samples with liver metastasis
A total of 3627 cells from GC samples were screened by the Seurat package after data pre-processing and we classified six major cell subpopulations, namely, B cells, Myeloid cells, Fibroblasts, TNK cells, Mast cells, and Epithelial cells (Fig. 1A, B). The distribution of genes, unique molecular identifiers, and the proportion of mitochondrial genes in each cell subpopulation was shown in Fig. 1C–E. The expressions of marker genes in each cell subpopulation were specifically as follows: Mast cells high-expressed MS4A2 and TPSAB1; Epithelial cells high-expressed KRT18 and KRT8; Myeloid cells high-expressed S100A8 and IL1B; B cells high-expressed MS4A1 and CD79A; Fibroblasts high-expressed LUM, DCN; TNK cells high-expressed CCL5 and NKG7 (Figs. 1F and 2A). Finally, we quantified the proportion of each cell type in GC and GC with liver metastasis groups, respectively, and observed that the GC liver metastasis group had more TNK cells than that of the GC group (Fig. 2B).

Differences in TNK cell-associated pathways between GC and liver-metastatic GC samples
DEG and enrichment analyses were performed to investigate differentially enriched biological process (BP) of TNK cells in GC cancer and GC liver metastasis groups. It was found that TNK cells were enriched in the pathways such as IL-17 signaling pathway, Pi3k–Akt Signaling Pathway, etc. in the GC liver metastasis group, while in the GC group TNK cells were enriched in Protein Export and some other pathways (Fig. 3A–D). This revealed that these pathways may be associated with tumor cell inflammatory responses, cancer cell survival and migration. Meanwhile, the present study analyzed the expression of genes related to Pi3k–Akt and IL-17 signaling pathways in the two GC groups and found that these genes (such as MAPK12, IL-1β, and MMP9) were significantly upregulated in GC liver metastasis group (Fig. 3E, F). Overall, these results demonstrated that the development of liver metastases from GC may be associated with intracellular signal transduction and inflammatory responses.

Landscape of TNK cell subpopulations in the GC and GC liver metastasis groups
To explore the heterogeneity of TNK cells, the cells were further clustered into CD8 T cells, Exhausted T cells, NKT cells, NK cells, and Treg cells (Fig. 4A). Specifically, CD8 T cells high-expressed CD8A; Exhausted T cells high-expressed CXCL13; NK cells high-expressed GNLY; NKT cells high-expressed WDR74; Treg cells high-expressed FOXP3 (Fig. 4B, C). Differential analysis of cellular infiltration between the samples showed that the infiltration ratio of both CD8 T cells and NKT cells in the GC liver metastasis group was significantly higher than that in the GC group (Fig. 4D, E).

Identification of TNK cell-related TFs associated with GC liver metastasis
The SCENIC package was used to identify TFs in TNK cells and the activity of each TF within a cellular subpopulation was calculated by the AUCell algorithm so as to reveal differentially expressed TNK-associated TFs in the GC liver metastasis group in comparison to the GC group (Fig. 5A). Among them, JUNB and FOS, which are associated with EMT of tumor cell lines in previous studies, were important regulators of the metastatic phenotype of cancer cells [23, 24]. Subsequent enrichment analysis on the two genes showed that FOS was mainly enriched in BP such as nuclear-transcribed mRNA catabolic, RNA catabolic, translational initiation, MAPK, Ribosome, Spliceosome, Apoptosis and some other pathways (Fig. 5B, C), and that JUNB was mainly enriched in nuclear-transcribed mRNA catabolic, mRNA catabolic process, RNA catabolic process, and some other BP (Fig. 5D). KEGG enrichment analysis on JUNB also revealed its involvement in ribosome, RNA transcription and other pathways (Fig. 5E). These results supported the role of the two TFs as potential therapeutic targets and provided a theoretical basis for in-depth study of their specific functions in tumors.

TNK cell-related TFs in regulating the metastasis of GC cell line
For the two key TFs mined in the above study, this study further explored their regulatory effects on GC cell lines. Molecular assay results showed that JUNB and FOS were remarkably upregulated in GC cell line, and that JUNB was a significant factor during GC progression (Fig. 6A). Wound healing assay revealed that JUNB significantly promoted the migration of GC cells (Fig. 6B). Transwell results showed that the invasion of JUNB-silenced GC cell lines was significantly inhibited (Fig. 6C). These findings demonstrated the regulatory role of JUNB in the invasion and migration of GC.

TNK cell-to-cell communication in GC with liver metastasis
To investigate the differences in intercellular communication of TNK cells in GC and GC liver metastasis groups, this study compared the differences in the enrichment of signaling pathways related to cell communication between the two groups. It was found that the signaling pathways differed significantly between the GC and GC metastasis groups, and that the number of these factors and cellular communication intensity were relatively increased in the liver metastasis group of GC (Fig. 7A, B). However, the interactions between different types of immune cells were remarkably reduced in the GC liver metastasis group (Fig. 7C, D). In the GC group, the cellular interaction sites mainly involved BAFF, CDH5, PARs, NRG, VWF, etc. (Fig. 7E), while in the GC liver metastasis group, the intercellular interaction sites mainly involved IL1, RESISTIN, L1CAM, FASLG, SELL, etc. (Fig. 7F).

Discussion
GC is the fifth most commonly diagnosed cancer and the fourth major cause of cancer-related mortality all over the world. The liver is the primary site for distant metastasis from GC but the specific mechanisms remains unclear. Using a series of in vivo screening techniques and transcriptome analysis together with qRT-PCR and tissue array assessments, Li et al. discovered that downregulated mitogen-activated protein kinase 4 (MAPK4) in the GC tissues was strongly related to liver metastasis and an unfavorable prognosis [25]. The present study performed single-cell transcriptome analysis to reveal the cell subpopulations that affected GC liver metastasis, elucidated their functions in GC progression, and also determined the key downstream pathways involved in regulating GC liver metastasis.
In this study, single-cell clustering of GC samples and GC liver metastasis samples classified six major cell subpopulations, of which the TNK cell subpopulation showed the highest infiltration in the GC liver metastasis group. Further clustering of TNK cell subpopulations subdivided the cells into NK/T cells and CD8+ T cells, both of which were significantly high-expressed in GC liver metastasis group. The regulatory role of immune cells in the TME favoring cancer progression has been previously reported [26]. Different from T and B cells that express the antigen recognition receptors TCR and BCR, NK cells are lymphocytes with the ability to lyse tumor cells and do not require antigenic pre-sensitization to destroy target cells [27]. At the same time, activated NK cells secrete a range of cytokines and chemokines to fulfill immunomodulatory functions. In addition to an anti-tumor effect, NK cells also secrete cytokines such as IFN-γ and TNF-α to promote inflammation, which can lead to oncogenic transformation of cells [28, 29]. Furthermore, NK cells can exhibit a pro-tumorigenic phenotype when the TME attempts to evade surveillance from NK cell, which explained the positive correlation between high infiltration of NK cells and the metastatic phenotype of GCs observed in this study [30]. CD8+ T cells are the main effector cells of antitumor immune response, and primary CD8+ T cells are activated to differentiate into cytotoxic effector T cells through MHC-I upon specific recognition of antigenic peptides presented by APC [31]. Similar to NK cells, CD8+ T cells, which could be divided into multiple subpopulations (TC1, TC2, TC17), are also pro-tumorigenic but not all the subpopulations are cytotoxic [32]. TC17 cells play dual roles in different cancers through the secretion of IL-17 A, IL-17 F, and IL-22. For instance, TC17 cells shows anti-tumor and pro-tumor activity in esophageal squamous carcinoma and head and neck squamous cell carcinoma, respectively [33, 34]. The important regulatory role of these cell types in the malignant phenotype of cancer cells, especially the metastatic phenotype of cancer cells, also further supports the conclusion of the present study that multiple immune cells were significantly upregulated in GC patients with liver metastasis.
Our enrichment analysis indicated that IL-17 and Pi3k–Akt signaling pathways were differentially enriched between the GC liver metastasis group and GC group. Studies showed that aberrant activation of IL-17 not only contributes to the development of autoimmune diseases, but is also involved in tumorigenesis and progression. In diseases such as GC, IL-17 activation involved in the regulation of PD-L1 expression mediates cancer immune escape to facilitate cancer progression [35]. In addition, high-expressed IL-17 positively promotes the JAK2/STAT3 signaling pathway, which is often associated with malignant phenotypes such as cancer cell proliferation, migration, and invasion [36]. Through the regulation of different downstream effectors, PI3K/AKT can be combined to form different cellular pathways to jointly participate in the regulation of different cellular processes such as cell cycle, cell survival, inflammation, metabolism and apoptosis [37]. PI3K/AKT pathway together with different growth factor receptors could effectively regulate cancer cell survival, proliferation and migration, especially in GC, where it associates with the inactivation of miR34a to enhance the metastatic phenotype of GC cells and cancer progression [38]. These discoveries could all support the mechanism of the regulatory role of TNK cells in GC liver metastasis discovered by the present research.
In this study, we identified a variety of TFs in TNK cells, in particular, FOS and JUNB showed significantly different expressions between the GC and GC liver metastasis groups, suggesting that these two genes might be important regulators in mediating the metastatic phenotypes of GC cells. FOS is a structural gene that encodes a tyrosine protein kinase involved in the promotion of cellular differentiation and also promotes the proliferative and migratory phenotype of cancer cells through the mechanisms including activating cellular inflammatory response [39]. FOS enhances angiogenesis in GC by elevating VEGFD expression and is significantly associated with a poor prognosis of patients with GC [40]. JUNB is a proto-oncogene that promotes cancer progression through the modulation of immunosuppressive phenotype of immune cells such as T-lymphocytes, macrophages, etc [41]. High expression of JUNB promotes the HGF-mediated migration phenotype of GC cells through cascade regulation with factors such as NF-kappaB and MMP-9 [23]. Thus, FOS and JUNB mediated the activity of TNK cells to affect the heterogeneity of the TME in GC and regulate the metastatic phenotype of GC cells.

Conclusion
In conclusion, we analyzed the single-cell transcriptomic data in GC samples and GC liver metastasis samples to reveal the regulatory role of TNK cell subpopulations in GC liver metastasis. The present results demonstrated that TNK cells affected GC progression mainly through Pi3k–Akt and IL-17 signaling pathways, and that CD8+ T cells, NKT cells, and their TFs (FOS and JUNB) promoted the metastatic phenotype of GC. Our findings are expected to provide guidance for the improvement of clinical treatment strategies for GC.

Limitation
However, this study also had certain limitations to be noted, in particular, the data in this study were mainly derived from public databases, which required tissue experiments to further validate the regulatory mechanisms of the cancer-related pathways and TFs during liver metastasis from GC.