ScRNA-seq revealed targeting regulator of G protein signaling 1 to mediate regulatory T cells in Hepatocellular carcinoma

Abstract

BACKGROUND:

Regulatory T cells (Tregs) are central to determine immune response, thus targeting Tregs for immunotherapy is a promising strategy against tumor development and metastasis.

OBJECTIVES:

The objective of this study was to identify genes for targeting Tregs to improve the outcome of HCC.

METHODS:

We integrated expression data from different samples to remove batch effects and further applied embedding function in Scanpy to conduct sub-clustering of CD4 $+$ T cells in HCC for each of two independent scRNA-seq data. The activity of transcription factors (TFs) was inferred by DoRothEA. Gene expression network analysis was performed in WGCNA R package. We finally used R packages (survminer and survival) to conduct survival analysis. Multiplex immunofluorescence analysis was performed to validate the result from bioinformatic analyses.

RESULTS:

We found that regulator of G protein signaling 1 (RGS1) expression was significantly elevated in Tregs compared to other CD4 $+$ T cells in two independent public scRNA-seq datasets, and increased RGS1 predicted inferior clinical outcome of HCC patients. Multiplex immunofluorescence analysis supported that the higher expression of RGS1 in HCC Tregs in tumor tissue compared to it in adjacent tissue. Moreover, RGS1 expression in Tregs was positively correlated with the expression of marker genes of Tregs, C-X-C chemokine receptor 4 (CXCR4), and three CXCR4-dependent genes in both scRNA-seq and bulk RNA-seq data. We further identified that these three genes were selectively expressed in Tregs as compared to other CD4 $+$ T cells. The activities of two transcription factors, recombination signal binding protein for immunoglobulin kappa J region (RBPJ) and yin yang 1 (YY1), were significantly different in HCC Tregs with RGS1 high and RGS1 low.

CONCLUSIONS:

Our findings suggested that RGS1 may regulate Treg function possibly through CXCR4 signaling and RGS1 could be a potential target to improve responses for immunotherapy in HCC.

Keywords

Tregs RGS1 immunotherapy tumor microenvironment ScRNA-seq

1. Introduction

Hepatocellular carcinoma (HCC), accounting for $\sim$ 90% of all primary liver cancer, is one of most frequent cause of cancer-related death worldwide [1]. HCC arises from transformed hepatocytes and is often diagnosed late. With traditional treatment for HCC, the 5-year survival rate of patients is less than 20% [2], moreover, the recurrent rate of HCC can reach up to 70–80% after 5 years [3, 4]. Immunotherapy provides a promising treatment to improve the clinical outcome of patients with HCC. Recent studies proved the safety and efficacy of CTLA-4 and PD-1-targeted immunotherapy [5], and the median overall survival after treatment with either immune checkpoint inhibitors (ICIs) or combination of those two inhibitors has been improved compared to traditional treatments [5, 6, 7].

Even though targeting the immune checkpoints for HCC treatment seems promising, only $\sim$ 15–20% HCC patients benefited from ICI therapy that inhibits either cytotoxic T-lymphocyte (CTLA-4) or programmed cell death-1 (PD-1) [8]. One of the main problems leading to the low response of those immunotherapies is the immune microenvironment dominated by immunosuppressive cells, such as regulator T cells (Tregs) and tumor-associated macrophages (M2) [9, 10]. These cells, normally correlate with poor prognosis, cooperate with dysfunctional CD8 $+$ T cells with highly expressed PD-1 to determine the response and resistance to immunotherapies. Single cell RNA sequencing (scRNA-seq) technology is extremely useful for elucidating the immune components. Previous analyses at single-cell levels have confirmed that immunosuppressive cells, especially regulatory T cells (Tregs), were infiltrated in the HCC immune microenvironment [11, 12].

Considering the evidence of Tregs’ infiltration in HCC microenvironment and the critical role of Tregs in suppressing the beneficial immune response in tumor microenvironment, studies on how to inhibit the suppressive function of infiltrated Tregs in tumor microenvironment became increasingly important. Recently, PTEN has been reported to affect Treg function through regulation of PI3K-AKT activity [13]. In non-small-cell lung cancer, IRF4 determines differentiation of effector Treg and immune suppressive function [14]. Moreover, enzymatic activity of HPGD, catabolized PGE ${}_{2}$ into 15-keto-PGE ${}_{2}$ , would mediate Treg suppression with a tissue- and context-dependent suppressive mechanism in human and mouse [15]. These previous studies all supported that targeting Tregs by key molecules for their functions has largely potential application for cancer immunotherapy and the efficiency of these targeting genes in Tregs needs to be improved. Taking advantage of sequencing technology at single cell level, some groups found a few Tregs feature genes in cancers, such as CARD16 and RGS1 [16, 17]. But these studies have been unclear in understanding mechanism of Tregs’ suppressive function in tumor microenvironment (TME) and how to better target Tregs in cancer immunotherapy.

To investigate the candidate genes for targeting Tregs in HCC as well as the potential mechanisms, we first utilized two independent public scRNA-seq datasets to profile CD4 $+$ T cells and identified Tregs cluster in HCC patients based on its marker genes. Surprisingly, we found that RGS1 expression was dramatically higher in Tregs compared to other CD4 $+$ T subpopulations. In addition, expression of RGS1 was positively correlated with the expression of Tregs marker genes, and elevated RGS1 was associated with short survival time in HCC using bulk RNA-seq data. Functional bioinformatic analyses suggested that four Treg function-related genes (CARD16, CXCR4, DUSP4, and PMAIP1) and two transcription factors (RBPJ and YY1) were all differentially expressed in RGS1 high Tregs compared to RGS1 low Tregs. Our results indicate that RGS1 potentially regulates Tregs in HCC and targeting RGS1 for Tregs is a promising strategy for HCC treatment.

2. Materials and methods

2.1 Data preprocessing

Two independent raw count matrixes for single-cell RNA-seq data were retrieved from the NCBI Gene Expression Omnibus database (accession code: GSE156337 and GSE151530) [18, 19]. We applied the same quality-control filtering used in original study to reproduce similar results. The annotation of cells was directly retrieved from original results. Additionally, we conducted a re-clustering analysis to validate annotations for CD4 $+$ and CD8 $+$ cells. We included quality-control filtering cells (containing all cell types) to conduct clustering analysis, and further used canonical markers (CD3D, CD3E, IL7R, and CD4 for CD4 $+$ T cell; CD3D, CD3E, CD8A, and CD8B for CD8 $+$ T cell) to annotate CD4 $+$ and CD8 $+$ T cells, respectively. We then compared the number of annotated CD4 $+$ and CD8 $+$ T cells between original study and our re-clustering analysis. The annotations from our analyses were consistent with original studies. To integrate expression data of CD4 $+$ T cells from patients with different status, we only included samples with number of CD4 $+$ T cells $>=$ 200. Finally, all CD4 $+$ T cells after above filtering were used for further analyses.

For two independent bulk RNA-seq data, the raw count files generated by HTSeq and clinical information for The Cancer Genome Atlas Liver Hepatocellular Carcinoma (TCGA-LIHC) cohort were obtained from GDC Data Portal (https://portal.gdc.cancer.gov/) [20]. We then applied TMM (trimmed mean of M values) normalization to raw count matrix for the following analysis. Another quantitated gene expression matrix came from RSME software with upper quartile normalization were downloaded from NODE (accession code: OEP000321) (https://www.biosino.org/node/) [21]. The log2-transformed was further conducted for matrix of Gao cohort.

2.2 Sub-clustering CD4 $+$ T cell

To investigate CD4 $+$ T cell sub-populations in HCC, we integrated expression data from different samples and removed batch effects by BBKNN (batch balanced k nearest neighbors) method. We further applied embedding function in Scanpy (V1.8.2) [22] to reduce the dimensionality of the scRNA-seq data and assigned cells to clusters with the Leiden graph-clustering method. Above analysis was performed to each of two public scRNA-seq data, respectively, to obtain CD4 $+$ T sub-population information for each dataset.

2.3 Differential gene expression analysis using scRNA-seq data

To compare gene expression signatures among different sub-clusters or groups, we first computed a ranking for the highly differential genes in each CD4 $+$ T sub-cluster using a Wilcoxon rank-sum test in Scanpy (V1.8.2). The differential expression of gene for cluster was defined as log2fc $\geqslant$ 0.25 and adjusted $p$ -value $\leqslant$ 0.05 between cluster and the rest of all clusters. Finally, we inferred the activity of transcription factors (TFs) for special group using DoRothEA [23]. The regulons of DoRothEA are gene regulatory networks including the collection of a transcription factor and its transcriptional targeting genes. Coupled with statistical methods, these regulons have been used to infer TF activities by transcriptomic data.

2.4 Validation of differential gene expression based on bulk RNA-seq

We applied the k-Means clustering to classify group based on expression matrix of selected genes from bulk RNA-seq data using scikit-learn module in python [24]. The difference of gene expression between groups from bulk RNA-seq was calculated by Wilcoxon rank-sum test in R project (https://www.r-project.org/). The heatmap was plotted by pheatmap R package (https://CRAN.R-project.org/package=pheatmap).

2.5 Weighted gene correlation network analysis (WGCNA)

To identify co-expression genes on HCC bulk RNA-seq data, we performed a gene expression network analysis with the WGCNA R package [25]. Top 10,000 genes with high variance among bulk RNA-seq samples were included into the further analyses for each dataset. The soft-thresholding power and co-expression adjacency were determined by recommended method from WGCNA. We then detected modules and calculated the correlation of genes for each module with default parameters.

2.6 Survival analysis for RGS1

To determine the optimal cut-point for the expression of RGS1, we applied surv_cutpoint function in survminer R package (https://CRAN.R-project.org/package=survminer) to calculate a value of RGS1 expression that correspond to the most significant relation with outcome of HCC patients in Gao and TCGA cohort, respectively. The Kaplan-Meier method with a log-rank test was used to generate survival curves (https://CRAN.R-project.org/package=survival). Hazard ration (HR) between two groups was estimated by the Cox regression model.

2.7 Multiplex immunofluorescence analysis

To label the expression spatial distribution and intensity of RGS1, CD4 and FOXP3 in HCC tumor tissues and adjacent tissues, the tissues were fixed with 4% paraformaldehyde within 0.5 h after being excised, and then dehydrated and embedded in paraffin using routine methods. The paraffin blocks of HCC tumor tissues or adjacent tissues were cutting into 4 $\mu$ m thick sections. The FFPE tissue sections melted and dehydrated at 60 ${}^{\circ}$ C for 20 min, and then subjected to deparaffinization and rehydration using xylene and alcohol, respectively. The paraffin sections were placed in EDTA buffer or citrate buffer in a microwave oven for antigen retrieval. Then a blocking buffer was used to block the sections for 30 min. The sections were incubated with the RGS1 primary antibody, HRP-conjugated secondary antibodies, and tyramine signal amplification (TSA). Following TSA deposition, antibody stripping and antigen retrieval were performed to label CD4 and FOXP3. Finally, 4’,6-diamidino-2- phenylindole (DAPI) was used to stain DNA. A Vectra 3 pathology imaging system was used to obtain images, and these images were then analyzed using the CaseViewer software.

3. Results

3.1 Landscape of CD4 $+$ T cells in HCC

To characterize the comprehensive spectrum of CD4 $+$ T cells in HCC, we analyzed scRNA-seq data (GSE156337) consisting of 14 patients with multiple time points (total 57 samples). After data preprocessing, mainly including quality control filtering, 17 HCC samples with more than 200 CD4 $+$ T cells were included, and BBKNN method was applied to integrate gene expression matrix and remove batch effects. Our database finally comprised 14,614 high-quality singe-cell transcriptomes. Uniform manifold approximation and projection (UMAP) of CD4 $+$ T cells by Scanpy (V1.8.2) identified eight distinct clusters in HCC (Fig. 1A).

Figure 1.

The profile of CD4 ${}^{+}$ T cell populations from hepatocellular carcinoma (HCC) patients. (A). Plot of eight CD4 ${}^{+}$ T cell subpopulations with 14,614 cells from 17 HCC samples (GSE156337). The dashed line highlights Tregs clusters. (B). Stacked-violin plot of eight CD4 ${}^{+}$ T cell subpopulations with signature genes. (C). Plot of eight CD4 ${}^{+}$ T cell subpopulations with 4,712 cells from 14 HCC samples (GSE151530). The dashed line highlights Tregs clusters. (D). Dot plot of eight CD4 $+$ T cell subpopulations with signature genes. (E) The fractions of four CD4 ${}^{+}$ main clusters in each sample in GSE156337. (F). The fractions of four CD4 $+$ main clusters in 14 samples in GSE151530.

To annotate CD4 $+$ T cell sub-populations, we applied Scanpy (V1.8.2) to conduct wilcoxon rank-sum test to detect differentially expressed marker genes in each of the eight clusters. Based on canonical markers, CD4 $+$ T cells with eight subpopulations represented major immune cell types, including naïve, cytotoxic, and regulatory states, mainly defined by CCR7, GZMA, TIGIT, CTLA4 and FOXP3 gene expression (Fig. 1B), respectively. Notably, CD4 $+$ cytotoxic T cells have been separated into three subpopulations based on their characterization of transcriptomes, respectively. In addition, there were three clusters without significant high expression of those marker genes, however, IL7R, CD69, and FOS genes were highly expressed in those three clusters.

To validate the classification of CD4 $+$ T cells, we analyzed another scRNA-seq data (GSE151530) including 4,712 high quality CD4 $+$ T cells from 14 HCC samples using the same method and identified eight CD4 $+$ T cell subpopulations. In addition to naïve, cytotoxic (3 subpopulations) and regulatory CD4 $+$ T cells, we also identified 3 CD4 $+$ T cell subpopulations without high expression of CCR7, GZMA, and FOXP3 marker genes in an independent scRNA-seq dataset (Fig. 1D). Moreover, the gene expression signature in each subpopulation of CD4 $+$ T cell was similar across the two independent scRNA-seq datasets. We further investigated the proportions of CD4 $+$ T cell subpopulations and these results showed the diverse composition of CD4 $+$ T cells in HCC patients from two datasets (Fig. 1E and F). Taken together, our analyses on two scRNA-seq datasets revealed the robust clustering of CD4 $+$ T cell subpopulations and heterogeneity of immune cell types in HCC.

3.2 High expression of RGS1 in HCC Tregs is correlated with poor clinical outcome

As one subpopulation of functionally distinct CD4 $+$ T cells, Tregs are essential for the regulation of inflammation and immune responses to determine the development of cancer. We performed differential gene expression analysis between Treg cells and the rest of CD4 $+$ T cells on single-cell transcriptomics to elucidate the gene expression signature of Tregs in HCC. As expected, the expression of Tregs marker genes (CTLA4, FOXP3, IL2RA and TIGIT) were all significantly increased in Treg cells compared to other CD4 $+$ T cells in both scRNA-seq datasets (Fig. 2A), while naïve (CCR7) and cytotoxic (GZMA and GZMK) marker genes were minimally expressed in Treg cells. Notably, RGS1 was among the top highly differentially expressed genes in HCC Tregs. Further investigated analysis showed that RGS1 were highly expressed in all Tregs subpopulations (Fig. 2B).

Figure 2.

RGS1 highly expressed in Tregs compared to other CD4 ${}^{+}$ T subpopulations in HCC. (A). Volcano plot of differentially expressed genes (DEGs) between Tregs and other CD4 ${}^{+}$ T cell clusters in GSE156337 and GSE151530. (B). Violin plots comparing the expression of RGS1 across eight CD4 $+$ T cell subpopulations in GSE156337 and GSE151530. (C). RGS1 expression positively correlated to the average expression of three Tregs marker genes (TIGIT, CTLA4, and FOXP3) in HCC from Gao cohort and TCGA cohort. (D) Increased expression of RGS1 related to poor overall survival (OS) in Gao cohort and TCGA cohort. Log-rank test is performed for survival analysis and curves are plotted by Kaplan-Meier method, and the hazard ratio (HR) is measured at 95% confidence interval. (E). Representative image showing multiplex immunofluorescence staining of RGS1, CD4, and FOXP3 in HCC tumor tissues and adjacent tissues. Each marker was represented by a different color, as indicated in the panel. White squares indicated represent image based in the RGS1, CD4, and FOXP3 marker expression. Scale bars: 50 $\mu$ m (Upper panel), 20 $\mu$ m (Lower panel).

To investigate the correlation between RGS1 and Tregs marker genes, we then analyzed bulk RNA-seq data from two independent HCC cohorts. After normalization of gene expressions, we first calculated the mean expression of three Tregs marker genes (CTLA4, FOXP3, and TIGIT) and then conducted correlation analysis of gene expression between RGS1 and the mean of the three Tregs marker genes. Results showed that RGS1 expression is positively associated with Tregs marker genes expression in both Gao cohort ( $R=$ 0.51, $P<$ 2.2e-16; Fig. 2C) and TCGA cohort ( $R=$ 0.51, $P<$ 2.2e-16; Fig. 2C), respectively.

Given the significant high expression of RGS1 in Tregs and the role of Tregs in facilitating oncogenesis, we next asked whether RGS1 expression predicts overall survival of HCC patients. Using the optimal cutpoint determined by survminer package, HCC patients in cohort were classified into low and high RGS1 expression groups. Survival analysis result indicated that high expression of RGS1 was significantly associated with poor outcome in both Gao cohort ( $P=$ 0.0002, HR $=$ 0.38; Fig. 2D) and TCGA cohort ( $P=$ 0.0066, HR $=$ 0.52; Fig. 2D). Multiplex immunofluorescence imaging further confirmed that the expression level of RGS1 in Tregs is higher in the HCC tumor tissues than it in adjacent tissues (Fig. 2E).

3.3 RGS1 may regulate Treg function through CXCR4 signaling

To investigate the role of RGS1 in Tregs, we classified Treg cells into RGS1 low and RGS1 high groups based on the 25 ${}^{\text{th}}$ percentile and 75 ${}^{\text{th}}$ percentile of RGS1 expression, respectively. Because of the high dropouts for scRNA-seq data, here we excluded Treg cells with RGS1 expression $=$ 0. We further identified differential expression genes between those two groups by Scanpy (V1.8.2). Interestingly, at least eight genes were significantly higher expressed (log2fc $\geqslant$ 0.25, $p$ -value $<$ 0.05) in RGS1 high Tregs as compared to RGS1 low Tregs in both scRNA-seq datasets (Fig. 3A and B). Here, we excluded ribosomal genes and 13 mitochondrial protein-coding genes.

Figure 3.

The expression of RGS1 significantly associated genes involved in Tregs function by CXCR4 signaling among scRNA-seq data. (A) Highly different expression genes between RGS1 high and low among Tregs in HCC (GSE156337). (B) Validation of DEGs in another dataset (GSE151530). (C) Stacked violin plot showing selectively expression of four genes involved in HCC Tregs function compared to other conventional CD4 $+$ cells in GSR156337. (D) Validation of selective expression of four key genes in GSE151530. (E) The expression of RGS1 positively correlated to expression of Tregs marker gene FOXP3 in GSE156337 ( $P=$ 1.2e-07). (F). Validation of positive correlation between RGS1 and FOXP3 in another dataset ( $P=$ 0.0015) (GSE151530).

As chemokines play a critical role in maintenance of Tregs homeostasis, we asked whether chemokine-related genes are differently expressed between RGS1 low Tregs and RGS1 high Tregs. Notably, the expression of chemokine receptor gene (CXCR4) was significantly enriched in RGS1 high Tregs in both scRNA-seq datasets (Fig. 3A and B). Our analysis suggested that RGS1 and CXCR4 may cooperate to affect the function of Tregs in HCC TME.

To further assess the role of RGS1 in Tregs function, we identified candidate genes selectively expressed in Tregs rather than in other conventional CD4 $+$ T cells from the eight differentially expressed genes. Although CXCR4 is ubiquitously expressed in different CD4 $+$ T subpopulations due to its key functions, half (4 of 8) of the differential expression genes mainly expressed in Tregs compared to other CD4 $+$ T cells (Fig. 3C and D). Among these four genes, at least two genes were shown to regulate Tregs function in TME either directly or indirectly. DUSP4 knockdown resulted in decreasing of FOXP3 expression and thus inhibited the suppressive function of activate Tregs [26]. NOXA, encoded by PMAIP1, is one of the pro-apoptotic BH3- only proteins, and affects the infiltration of Tregs in TME by BCL-XL dependent manner [27]. BTG3 and CARD16 have also been reported to be the signature genes for Tregs in multiple cancers [16, 17, 28]. Most importantly, the expressions of these four genes are either directly regulated by CXCR4 signaling or significantly correlated with the expression of CXCR4 [29, 30].

Co-inhibitory receptor TIGIT is predominantly expressed on the surface of Tregs and is shown to play a vital role in the suppressive function of Tregs [31, 32]. Our differentially expressed genes (DEGs) analysis identified that TIGIT expression was significantly increased in Tregs with RGS1 high compared to RGS1 low (Fig. 3A and B), suggesting an important role of RGS1 in regulation of Tregs function. Interestingly, similar positive correlation was not found on FOXP3, another Treg marker gene. To figure out why, we investigated the expression of FOXP3 for each individual Tregs and found that the FOXP3 expression values in more than half of Tregs were zero due to high drop-out rate in the scRNA-seq data we used here. However, it should be noted that RGS1 expression was positively related to FOXP3 expression (GSE156337, $P$ -value $=$ 1.2e-07; GSE151533, $P$ -value $=$ 0.0015) in Tregs with expressed both genes (Fig. 3E and F). Our analysis revealed that expression of RGS1 was closely related to two functional marker genes of Tregs in HCC. Taken together, our analysis found expressions of Tregs signature marker gene FOXP3 and TIGIT were positively correlated to RGS1 expression, and CXCR4 and related genes were significantly increased in HCC Tregs with RGS1 high compared to its low.

3.4 Validated correlation using bulk RNA-seq

To validate the above analysis results found in scRNA-seq data, we conducted WGCNA using two bulk RNA-seq datasets. Although bulk RNA-seq data was mixed by diverse cell types and the proportion of Tregs in bulk RNA-seq was comparatively small, the results from bulk RNA-seq data were consistent with our analyses from scRNA-seq, revealing that the correlation network of RGS1 not only included Treg marker genes (CTLA4, IL2RA and TIGIT), but also part of key genes with CXCR4 signaling (CARD16, DUSP4, and PMAIP1) (Fig. 4A and B). The Pearson correlation analysis also suggested that RGS1 was highly related to CXCR4 and three genes involved in CXCR4 signaling (Fig. 4C and D).

Figure 4.

Validation of association between RGS1 and Tregs signature genes by two independent bulk RNA-seq datasets. (A) Weighted gene co-expression network analysis (WGCNA) showing the association of RGS1 and Tregs signature genes in bulk RNA-seq data from Gao cohort. Red: top 100 genes associated with RGS1; yellow: top 200 associated genes; blue: top 500 associated genes; green: top 600 associated genes. (B) Consistent association found in TCGA cohort. (C) Significantly positive association between RGS1 and Tregs signature genes (CARD16, CXCR4, DUSP4 and PMAIP1) in Gao cohort. The correlation coefficients were measured by Pearson correlation. (D) Validation of positive association by TCGA cohort. (E) Heatmap showing diverse expression of genes used for k-means clustering. Three (low, middle, and high) groups are determined by k-means clustering method based on the expression of RGS1, CXCR4, CARD16, DUSP4 and PMAIP1. Boxplot showing significantly different expression of Tregs marker genes between RGS1 low and high group from Gao cohort. (F) Validation of significantly different expression of Tregs marker genes in TCGA cohort.

Given the high correlation between RGS1 and the genes related to CXCR4 signaling, we asked whether the expressions of Treg marker genes can be predicted by RGS1 and CXCR4 signaling genes. We first applied k-means clustering to classify HCC patients into three groups based on expression matrix of RGS1 and CXCR4 signaling genes. As shown in Fig. 4E and F, there were three distinct groups with significantly different expression of their genes. We further compared the expression of Treg marker genes (CTLA4, FOXP3, IL2RA, and TIGIT) in low and high groups defined by RGS1 and CXCR4 signaling genes. Our data suggested that these four Treg marker genes were all significantly increased in RGS1 high group compared to its low group ( $P$ values from Wilcoxon rank-sum test for each comparison all $<$ 0.05). All results supported that the expression of RGS1 and CXCR4 signaling genes were positively correlated with Tregs in HCC.

3.5 Potential transcription factors upstream of RGS1 in HCC Tregs

To identify the key transcription factors downstream of RGS1, we conducted functional TF analysis between RGS1 high and RGS1 low cells in Tregs using DoRothEA. This approach predicts the activity of transcription factors by the expression of their targets. Here, we excluded the low confidence levels D and E to only maintain high quality regulons in DoRothEA for further analysis. The prediction results revealed that the TF activity of RBPJ (Fig. 5A and B) was significantly increased in Tregs with RGS1 high in both scRNA-seq datasets (mean change $=$ 0.40, adjusted $P$ -value $=$ 6.29e-23 in GSE156337 and mean change $=$ 0.48, adjusted $P$ -value $=$ 7.76e-23 in GSE151530, respectively), whereas YY1 activity (Fig. 5C and D) was suppressed in RGS1 high Tregs compared to RGS1 low Tregs (mean change $=-$ 0.27, adjusted $P$ -value $=$ 0.0077 in GSE156337 and mean change $=-$ 0.49, adjusted $P$ -value $=$ 2.31e-42 in GSE151530, respectively). By digging into public available ChIP-seq data, we found that RGS1 is a direct downstream target of both RBPJ and YY1 [33]. Additionally, transcription factor YY1 is showed to regulate CXCR4 expression in a calcineurin-dependent manner in human T lymphocytes [34]. Thus, we predicted that RBPJ and YY1 might be the key upstream transcription factors of RGS1 in regulating Treg function in TME via CXCR4 signaling.

Figure 5.

The inference of activation of transcription factors for RGS1 low and high expressed groups. (A) RBPJ activity significantly elevated in RGS1 high expressed group compared to it low group in GSE156337. (B) Validation of elevated RBPJ activity with RGS1 high expression in another dataset (GSE151530). (C) Decreased YY1 activity in Tregs RGS1 high expressed group compared to it low group in GSE156337. (D) Validation of decreased activity of YY1 in another dataset (GSE151530). (E) Elevated RBPJ activity in group with high expression of RGS1 and five related genes in Gao cohort ( $P=$ 0.0005, Wilcoxon rank-sum test). Three (low, middle, and high) groups are determined by k-means clustering method. (F) Validation of elevated RBPJ activity in TCGA cohort ( $P=$ 2.8e-5, Wilcoxon rank-sum test). (G) No significant difference of YY1 activity in Gao cohort ( $P>$ 0.05, Wilcoxon rank-sum test). (H) YY1 activity decreases in patient group with high expression of RGS1 and five related genes in TCGA cohort ( $P=$ 2.1e-5, Wilcoxon rank-sum test).

To validate the candidate TFs identified at single cell level, we directly compared the expression level of RBPJ and YY1 between low and high groups defined by RGS1 and its five related genes from bulk RNA-seq data with k-means clustering. We identified the significantly elevated RBPJ in HCC patients with high expression of RGS1 and its related genes in both cohorts (Fig. 5E and F; $P=$ 0.0005 in Gao cohort, $P=$ 2.8e-5 in TCGA cohort, Wilcoxon rank-sum test). Even though the expressions of YY1 were not significantly different in two groups in Gao cohort (Fig. 5G; $P>$ 0.05, Wilcoxon rank-sum test), YY1 in TCGA cohort was significantly downregulated in HCC RGS1 high patients (Fig. 5H; $P=$ 2.1e-5, Wilcoxon rank-sum test). In total, our TF functional analysis suggested that RBPJ may positively regulate the expression of RGS1 in Tregs, while negatively regulate YY1 transcription activity in HCC.

4. Discussion

Immunotherapy for HCC mainly including ICIs has revolutionized the treatment paradigm and improved overall survival of patients. However, one major challenge of immunotherapy is the low response rate, primarily caused by the accumulation of immunosuppressive subset of cells such as Tregs in tumor microenvironment [9, 35, 36]. These immune suppressive cells can hamper the antitumor response in TME and eventually lead to poor outcome of HCC patients. In this study, our comprehensive bioinformatic analyses suggested that high expression of RGS1 may determine the function of infiltrated Tregs in HCC through CXCR4 signaling and unveiled a new role of RGS1 in regulating Tregs in HCC.

Our study appears to expand the understanding the role of RGS1 in regulation of immune cells and provide a new way to improve the efficacy of targeting immunotherapy in cancers. In HCC, high expression of RGS1 was positively associated with bad outcomes, supported by previous studies in few cancers, such as gastric cancer [37] and multiple myeloma [38]. Moreover, emerging evidence suggested that RGS1 regulating inflammation and accumulation of immune cells was a potential target for immunotherapy [39]. Recent studies described a novel mechanism where the trafficking of antitumor T cells, such as CD8 $+$ CTLs and CD4 $+$ T ${}_{\text{H}}$ 1, is be inhibited by elevated RGS1 expression [40, 41]. Here, we found that RGS1 also promotes the suppressive function of infiltrated Treg in tumor TME. Thus, inhibiting RGS1 could theoretically not only facilitate the infiltration of antitumor T cells to tumor sites, but also suppresses the function of immunosuppressive Tregs in TME.

Previous studies demonstrated that inhibition of CXCR4 affecting Treg function in multiple cancers, such as renal cancer and pancreatic cancer [42, 43]. More importantly, RGS1 is shown to bound to CXCR4 and determine homeostasis of CTLs and T follicular helper cells [44]. Consistently, our analysis suggested that RGS1 and CXCR4 may cooperate to affect the function of Tregs in HCC TME. Although we conducted a systematic bioinformatic analysis and revealed a potential role of RGS1 in regulating Tregs functions, it is should be mentioned that the detailed function of RGS1 in regulating Tregs and the mechanism behind it still need large experiment evidence to better guide clinical application of RGS1. In future, validation from functional experiments could help us to determine the role of CXCR4 signaling pathway in RGS1 regulating Treg function.

Metabolic reprogramming in Tregs is closely related to their suppressive function and targeting alternation of metabolic genes may provide a promising therapeutic strategy by Tregs for treatment of cancers. Interestingly, our study identifies that CARD16 and PMAIP1, two metabolism related genes, were significantly high expressed in Tregs with RGS1 high compared to RGS1 low in all HCC scRNA-seq and bulk RNA-seq data. CARD16 and PMAIP1 should be high confidently candidate genes for further experiment study to valid their clinical application. Due to these two metabolic genes related to Tregs function through metabolic changes, targeting them may largely improve efficiency of HCC treatment.

In addition, we identified potential upstream transcription factors of RGS1 in regulating Treg function in the TME of HCC. Based on RNA-seq data at single cell and bulk levels, two TFs (RBPJ and YY1) are inferred to be the candidate transcription factors upstream of RGS1 in regulating Treg function. Consistent with previously experimental studies [45, 46, 47], our study revealed RBPJ and YY1 both involved in the regulation of Tregs function, these two TFs are also candidate molecules to conduct validated study.

In conclusion, we characterized gene expression signature of CD4 $+$ T subpopulations in HCC and observed dramatically elevated RGS1 expression in Tregs. Functional bioinformatic analysis predicted a novel function of RGS1 in regulating Tregs suppressive function through CXCR4 signaling. Our study provides a potential therapeutic target to boost efficacy of antitumor response and improve the clinical outcome of HCC patients.

Author contributions

Conception: CX Zhao.

Interpretation or analysis of data: LH Zou, KH Liu, YZ Shi, GW Li.

Preparation of the manuscript: KH Liu, LH Zou, CX Zhao.

Revision for important intellectual content: LH Zou, HY Li, CX Zhao.

Supervision: CX Zhao.

All authors contributed to the article and approved the submitted version.

Funding

This study was funded by the Key Project of Hunan Provincial Science and Technology Innovation (No. 2020SK1010).

Footnotes

Acknowledgments

We would like to thank the Dr. Linh T. Bui for critical reading of the manuscript and discussions about the project.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

References

Llovet

J.M.

et al., Hepatocellular carcinoma, Nature Reviews Disease Primers 7 (2021), 1–28.

Villanueva

, Hepatocellular carcinoma, New England Journal of Medicine 380 (2019), 1450–1462.

Sherman

, Recurrence of hepatocellular carcinoma, New England Journal of Medicine 359 (2008), 2045–2047.

Saito

et al., Prediction of early recurrence of hepatocellular carcinoma after resection using digital pathology images assessed by machine learning, Modern Pathology 34 (2021), 417–425.

Llovet

J.M.

et al., Immunotherapies for hepatocellular carcinoma, Nature Reviews Clinical Oncology (2021), 1–22.

Pinter

Scheiner

and Peck-Radosavljevic

, Immunotherapy for advanced hepatocellular carcinoma: A focus on special subgroups, Gut 70 (2021), 204–214.

Sangro

et al., Advances in immunotherapy for hepatocellular carcinoma, Nature Reviews Gastroenterology & Hepatology 18 (2021), 525–543.

Cheng

A.-L.

et al., Challenges of combination therapy with immune checkpoint inhibitors for hepatocellular carcinoma, Journal of Hepatology 72 (2020), 307–319.

Zhu

et al., Nanomedicines modulating tumor immunosuppressive cells to enhance cancer immunotherapy, Acta Pharmaceutica Sinica B 10 (2020), 2054–2074.

10.

et al., Multifunctional nanoparticles boost cancer immunotherapy based on modulating the immunosuppressive tumor microenvironment, ACS Applied Materials & Interfaces 12 (2020), 50734–50747.

11.

et al., Current perspectives on the immunosuppressive tumor microenvironment in hepatocellular carcinoma: Challenges and opportunities, Molecular Cancer 18 (2019), 130.

12.

D.W.-H.

et al., Single-cell RNA sequencing shows the immunosuppressive landscape and tumor heterogeneity of HBV-associated hepatocellular carcinoma, Nature Communications 12 (2021), 3684.

13.

Lam

A.J.

et al., PTEN is required for human Treg suppression of costimulation in vitro , European Journal of Immunology 52 (2022), 1482–1497

14.

Alvisi

et al., IRF4 instructs effector Treg differentiation and immune suppression in human cancer, Journal of Clinical Investigation 130 (2020), 3137–3150.

15.

Schmidleithner

et al., Enzymatic activity of HPGD in Treg cells suppresses Tconv cells to maintain adipose tissue homeostasis and prevent metabolic dysfunction, Immunity 50 (2019), 1232–1248.e14.

16.

Luo

et al., Single-cell transcriptomic analysis reveals disparate effector differentiation pathways in human Treg compartment, Nature Communications 12 (2021), 3913.

17.

Venema

W.T.U.

et al., Single-cell RNA sequencing of blood and ileal T cells from patients with Crohn’s disease reveals tissue-specific characteristics and drug targets, Gastroenterology 156 (2019), 812–815e22..

18.

Sharma

et al., Onco-fetal reprogramming of endothelial cells drives immunosuppressive macrophages in Hepatocellular carcinoma, Cell 183 (2020), 377–394.e21.

19.

et al., Single-cell atlas of tumor cell evolution in response to therapy in hepatocellular carcinoma and intrahepatic cholangiocarcinoma, Journal of Hepatology 75 (2021), 1397–1408.

20.

Cancer Genome Atlas Research Network, Comprehensive and integrative genomic characterization of Hepatocellular carcinoma, Cell 169 (2017), 1327–1341.e23.

21.

Gao

et al., Integrated proteogenomic characterization of HBV-related Hepatocellular carcinoma, Cell 179 (2019), 1240.

22.

Wolf

F.A.

Angerer

and Theis

F.J.

, SCANPY: Large-scale single-cell gene expression data analysis, Genome Biology 19 (2018), 15.

23.

Garcia-Alonso

et al., Benchmark and integration of resources for the estimation of human transcription factor activities, Genome Research 29 (2019), 1363–1375.

24.

Pedregosa

et al., Scikit-learn: Machine learning in Python, Journal of Machine Learning Research 12 (2011), 2825–2830.

25.

Langfelder

and Horvath

, WGCNA: An R package for weighted correlation network analysis, BMC Bioinformatics 9 (2008), 559.

26.

Yan

et al., Imbalanced signal transduction in regulatory T cells expressing the transcription factor FoxP3, Proceedings of the National Academy of Sciences 112 (2015), 14942–14947.

27.

Kolb

et al., Proteolysis-targeting chimera against BCL-XL destroys tumor-infiltrating regulatory T cells, Nature Communications 12 (2021), 1281.

28.

Dias

et al., Effector regulatory T cell differentiation and immune homeostasis depend on the transcription factor Myb, Immunity 46 (2017), 78–91.

29.

Kremer

K.N.

et al., CXCR4 chemokine receptor signaling induces apoptosis in acute myeloid leukemia cells via regulation of the Bcl-2 family members Bcl-XL, Noxa, and Bak, Journal of Biological Chemistry 288 (2013), 22899–22914.

30.

Zhang

et al., Tissue Treg secretomes and transcription factors shared with stem cells contribute to a Treg niche to maintain Treg-ness with 80% innate immune pathways, and functions of immunosuppression and tissue repair, Frontiers in Immunology 11 (2021), 632239.

31.

Joller

et al., Treg cells expressing the co-inhibitory molecule TIGIT selectively inhibit pro-inflammatory Th1 and Th17 cell responses, Immunity 40 (2014), 569–581.

32.

Kurtulus

et al., TIGIT predominantly regulates the immune response via regulatory T cells, Journal of Clinical Investigation 125 (2015), 4053–4062.

33.

Rouillard

A.D.

et al., The harmonizome: A collection of processed datasets gathered to serve and mine knowledge about genes and proteins, Database (Oxford) 2016 (2016), baw100.

34.

Cristillo

A.D.

and Bierer

B.E.

, Regulation of CXCR4 expression in human T lymphocytes by calcium and calcineurin, Molecular Immunology 40 (2003), 539–553.

35.

Waldman

A.D.

Fritz

J.M.

and Lenardo

M.J.

, A guide to cancer immunotherapy: From T cell basic science to clinical practice, Nature Reviews Immunology 20 (2020), 651–668.

36.

Tang

et al., Advantages of targeting the tumor immune microenvironment over blocking immune checkpoint in cancer immunotherapy, Signal Transduction and Target Therapy 6 (2021), 1–13.

37.

et al., High expression of regulator of G-protein signalling 1 is associated with the poor differentiation and prognosis of gastric cancer, Oncology Letters 21 (2021), 322.

38.

Roh

et al., RGS1 expression is associated with poor prognosis in multiple myeloma, J Clin Pathol 70 (2017), 202–207.

39.

Bai

et al., Single-cell transcriptome analysis reveals RGS1 as a new marker and promoting factor for T-cell exhaustion in multiple cancers, Frontiers in Immunology 12 (2021), 767070.

40.

Fercoq

and Carlin

L.M.

, “Mind the GAP”: RGS1 hinders antitumor lymphocytes, Nature Immunology 22 (2021), 802–804.

41.

Huang

et al., Targeting regulator of G protein signaling 1 in tumor-specific T cells enhances their trafficking to breast cancer, Nature Immunology 22 (2021), 865–879.

42.

Santagata

et al., CXCR4 and CXCR7 signaling pathways: A focus on the cross-talk between cancer cells and tumor microenvironment, Frontiers in Oncology 11 (2021), 591386.

43.

Santagata

et al., Targeting CXCR4 reverts the suppressive activity of T-regulatory cells in renal cancer, Oncotarget 8 (2017), 77110–77120.

44.

Caballero-Franco

and Kissler

, The autoimmunity-associated gene RGS1 affects the frequency of T follicular helper cells, Genes & Immunity 17 (2016), 228–238.

45.

Delacher

et al., Rbpj expression in regulatory T cells is critical for restraining TH2 responses, Nature Communications 10 (2019), 1621.

46.

Balkhi

M.Y.

et al., YY1 upregulates checkpoint receptors and downregulates type I cytokines in exhausted, chronically stimulated Human T cells, iScience 2 (2018), 105–122.

47.

Kosasih

F.R.

and Bonavida

, Involvement of Yin Yang 1 (YY1) expression in T-cell subsets differentiation and their functions: implications in T cell-mediated diseases, Critical Reviews in Immunology 39 (2019), 491–510.

ScRNA-seq revealed targeting regulator of G protein signaling 1 to mediate regulatory T cells in Hepatocellular carcinoma

Abstract

BACKGROUND:

OBJECTIVES:

METHODS:

RESULTS:

CONCLUSIONS:

Keywords

1. Introduction

2. Materials and methods

2.1 Data preprocessing

2.2 Sub-clustering CD4 + T cell

2.3 Differential gene expression analysis using scRNA-seq data

2.4 Validation of differential gene expression based on bulk RNA-seq

2.5 Weighted gene correlation network analysis (WGCNA)

2.6 Survival analysis for RGS1

2.7 Multiplex immunofluorescence analysis

3. Results

3.1 Landscape of CD4 + T cells in HCC

Author contributions

Funding

Footnotes

Acknowledgments

Conflict of interest

References

2.2 Sub-clustering CD4 $+$ T cell

3.1 Landscape of CD4 $+$ T cells in HCC