Genomic and transcriptional Profiling of tumor infiltrated CD8⁺ T cells revealed functional heterogeneity of antitumor immunity in hepatocellular carcinoma

ABSTRACT

As key players in HCC antitumor response, the functions of tumor infiltrated CD8⁺ T cells are significantly affected by surrounding microenvironment. A detailed profiling of their genomic and transcriptional changes could provide valuable insights for both future immunotherapy development and prognosis evaluation. We performed whole exome and transcriptome sequencing on tumor infiltrated CD8⁺ T cells and CD8⁺ T cells isolated from other tissue origins (peritumor tissues and corresponding PBMCs) in eight treatment-naive HCC patients. The results demonstrated that transcriptional changes, rather than genomic alterations were the main contributors to the functional alterations of CD8⁺ T cells in the process of tumor progression. The origins of CD8⁺ T cells defined their transcriptional landscape, while the tumor infiltrated CD8⁺ T cells shared more similarity with peritumor-derived CD8⁺ T cells compared with those CD8⁺ T cells in blood. In addition, tumor infiltrated CD8⁺ T cells also showed larger transcriptional heterogeneity among individuals, which was modulated by clinical features such as HBV levels, preoperative anti-viral treatment and the degree of T cell infiltration. We also identified multiple inter-connected pathways involved in the activation and exhaustion of tumor infiltrated CD8⁺ T cells, among which IL-12 mediated pathway could dynamically reflect the functional status of CD8⁺ TILs and activation of this pathway indicated a better prognosis. Our results presented an overview picture of CD8⁺ TILs’ genomic and transcriptional landscape and features, as well as how the functional status of CD8⁺ TILs correlated with patients’ clinical course.

KEYWORDS:

• CD8⁺ T cells
• Hepatocellular carcinoma (HCC)
• transcriptome analysis
• IL-12 signaling pathway
• tumor infiltrated lymphocytes (TILs)

Introduction

Hepatocellular carcinoma (HCC), which is considered as the most common primary liver malignancy in China, ranks the third among the causes of cancer mortality.¹ At present, treatment of HCC, still faced great challenges.² Cancer immunotherapies are under rapid development, providing a promising strategy to deal with various cancers by boosting the body’s natural immune response.³ CD8⁺ T cells are the essential effector cells in antitumor immunity⁴ and generation of CD8⁺ cytotoxic T lymphocytes (CTLs) is one of the major immunity response for encountering tumor progression and possible hepatitis virus infection. When migrated into tumor, these CD8⁺ Tumor-infiltrating lymphocytes (TILs) become the first line of defense, with their abundance closely correlating to patients’ prognosis and clinical outcomes.⁵

To resist adaptive immunity, HCC developed multiple mechanisms to hamper the induction of functional TILs. One of the mechanisms is cancer immunoediting, which reshapes tumor immunogenicity via evolutionary selection of mutations with weak immunogenicity,⁶ leading to tumor escape from immune surveillance.⁷ Besides, when entering into the tumor microenvironment, under the pressure provided by tumor cells as well as inflammatory background, the tumor-infiltrating effector T cells often become ‘exhausted’.⁸ Exhaustive T cells only retained impaired effector functions, which often lead to dysfunction in tumor cytotoxicity. Studies have been conducted to investigate how T cell exhaustion come to place, and several intrinsic mechanisms have been discovered. The lack of proper expression of co-stimulatory molecules (such as CD80 and CD86) can prevent T cells’ normal activation, while the existence of high level inhibitory receptors further worsens the scenario.⁹ Several different types of co-inhibitory signals, such as PD-1and CTLA-4, have already been identified in past decades, and the blockade of these inhibitory receptors can restore or even further enhance T cell function in cancers.¹⁰ At present, immunotherapies to block the checkpoints of T cells, have achieved effective immune response in several cancer types, but there are limited clinical response in HCC, which may be resulted from the chronic HBV-mediated immunosuppression.¹¹

To improve the immunotherapy based on tumor-infiltrated CD8⁺ T cells, the functional diversity of CD8⁺ TILs in HCC and the underlying mechanisms should be fully investigated. Recent study from Chunhong Zheng et al. discovered the functional composition of T cells in HCC by transcriptome analysis, demonstrating a comprehensive picture of T cell subpopulations. However, crosstalk between pathways involved in functional alteration of CD8⁺ TILs and their clinical relevance have not been thoroughly discussed. Here, we present an extensive analysis of CD8⁺ TILs’ multi-dimensional profile including both genomic and transcriptional levels to reveal the functional heterogeneity, clinical relevance and transcriptional features of CD8⁺ TILs. Tumor, peritumor tissues and match PBMCs were collected from 8 HCC patients with clinical variety to investigate the dynamic changes of CD8⁺ TILs’ cell function, which could provide valuable knowledges for future prognosis evaluation and immunotherapy development in HCC treatment.

Results

Genomic and transcriptional landscape of CD8⁺ T cells in HCC patients

To investigate the multi-dimensional landscape of CD8⁺ T cells in HCC patients, eight treatment-naive patients diagnosed with HCC were enrolled. Detailed information regarding to patients’ clinical information is available in Table 1. Tumor tissues, peritumor tissues and matched PBMCs were collected and the CD8⁺ T cells from these origins were isolated respectively (Figure 1(a)). The presence of infiltrating CD8⁺ T cells in HCC tissues as well as peritumor tissues was confirmed using immunohistochemistry (IHC) with CD3 and CD8 antibodies (Figure 1(b) and Figure S1). Isolated CD8⁺ T cells with high purity (≥ 90%) were subjected to whole exome and transcriptome sequencing (Figure 1(c)).

Table 1. Clinical characteristics of HCC patients in present study.

Patient ID	Gender	Age	Tumor location	Tumor thrombus	TNM stage	HBV-DNA (U/ml)	Cirrhosis	Hepatitis B	Antiviral therapy
1377	Female	55	V,VII	Positive	IV	Undetectable	Positive	Negative	No
1967	Female	57	II,III	Positive	II	6.26 × 10^4	Positive	Positive	No
2129	Male	43	VI,VII	Positive	II	Undetectable	Negative	Positive	No
2834	Male	43	VII	Positive	II	1.12 × 10^5	Positive	Positive	No
6423	Female	64	IV,VI,VII,VIII	Positive	II	4.22 × 10^6	Positive	Positive	No
7898	Male	56	VIII	Positive	II	3.75 × 10^5	Positive	Positive	No
7971	Female	65	V	Positive	II	Undetectable	Positive	Positive	Yes
8306	Male	61	V,VII	Positive	II	Undetectable	Positive	Positive	No

Figure 1. Isolation and purification of CD8⁺ T cells from HCC patients. (a) Flow chart of study design. (b) Representative immunohistochemistry staining images of CD3 and CD8 for HCC and peritumoral biopsies from patient 2129. (c) Representative flow cytometry plots of CD8⁺ T cells collected from different tissue origins after MACs purification.

Somatic alteration identification revealed extremely low burden of mutations in both CD8⁺ TILs and peritumor-derived CD8⁺ T cells (Figures 2(a) and 2(b)). Detailed information of somatic SNPs and CNVs was listed in online supplementary Table S1. The low burden and relatively low frequency (88% of the somatic mutations occur at mutation frequencies of less than 0.2, Table S2) of the somatic mutations in CD8⁺ TILs and peritumor-derived CD8⁺ T cells indicated that the genome of CD8⁺ T cells did not change significantly compared with corresponding PBMCs and thus might have limited contribution to functional variation.

Figure 2. Genomic and transcriptomic landscapes of CD8⁺ T cells from different origins in HCC patients. (a and b) Circos plots of somatic SNVs and CNVs. Somatic SNVs and CNVs identified in CD8⁺T cells from tumor and peritumor tissues of all the eight patients were shown. The outermost ring shows the chromosome location. Black spots indicate somatic mutations, blue/red lines represent somatic deletions/duplications respectively. (c) The X axis represents the fold changes of gene expression between peritumor-derived CD8⁺ T cells versus PBMC-derived CD8⁺ T cells while the Y axis represents the fold changes between CD8⁺ TILs versus peritumor-derived CD8⁺ T cells. Genes significantly differentially expressed (adjusted P < 0.05) in both comparisons were colored with red, and genes only showed significant different expression in either CD8⁺ TILs vs peritumor-derived CD8⁺ T cells or peritumor-derived CD8⁺ T cells vs PBMC-derived CD8⁺ T cells were colored with blue and golden respectively. Other genes were colored with grey. (d) PID pathway analysis of DE genes. The outline color of each node represents the normalized significance for comparison between peritumor-derived CD8⁺ T cells and PBMC-derived CD8⁺ T cells while the inner color represents that of comparison between CD8⁺ TILs and peritumor-derived CD8⁺ T cells. Nodes are connected if the genes within two pathways were significantly overlapped (Fisher’s exact test P < 0.05). (e) Heatmap of normalized FPKM of DE genes within IL-12 mediated signal pathway.

To explore whether changes on transcriptional level affect CD8⁺ T cells’ function in different microenvironments, differentially expressed (DE) genes among CD8⁺ T cells from different tissue origins were extracted. The existence of large amount of DE genes between CD8⁺ TILs/peritumor-derived CD8⁺ T cells and corresponding PBMC-derived CD8⁺ T cells clearly demonstrated significant differences on transcriptome level (2200 and 3919 DE genes, respectively, Figure 2(c)), while CD8⁺ TILs and peritumor-derived CD8⁺ T cells exhibited more similarity on expression patterns (162 DE genes, Figure 2(c)). Furthermore, for DE genes identified from comparisons of either CD8⁺ TILs vs PBMC-derived CD8⁺ T cells or peritumor-derived CD8⁺ T cells vs PBMC-derived CD8⁺ T cells, pathway analysis was conducted to understand how the expression profile changes impacted different biological pathways. Significant enrichment of multiple inter-connected pathways focusing on “downstream signaling of naive CD8⁺ T cells” was discovered (Figure 2(d)), indicating the functional changes of CD8⁺ TILs and peritumor-derived CD8⁺ T cells were triggered by naive CD8⁺ T cell related signal alterations. Interestingly, “IL-12 mediated signaling events” was the only pathway enriched in DE genes identified in both comparisons, suggesting the importance of this pathway during the entire process of adaptive immunity. Most DE genes involved in IL-12 mediated signaling events were up-regulated in peritumor-derived CD8⁺ T cells while down-regulated in CD8⁺ TILs (Figure 2(e)). Intriguingly, expression pattern of this pathway was in well consistence with T cell functional status since peritumoral immune cells exhibited activated phenotypes for battling against tumor cells,¹² while T cells in HCC tumor tissue were more likely to display exhausted phenotypes.⁸

Principal component analysis revealed how microenvironment affect CD8⁺ TILs’ transcriptome profiles

Given the complicated interactions among the biological pathways identified in T cells from different tissue origins, principal component analysis (PCA) of transcriptome data was carried out to provide a more in-depth description of the CD8⁺ T cells’ function. As demonstrated in Figure 3(a), CD8⁺ TILs, peritumor-derived CD8⁺ T cells and PBMC-derived CD8⁺ T cells could be well separated by the first three principal components (PCs), indicating the inherent transcriptional similarities of CD8⁺ T cells within the same tissue origin. PC1, which accounted for 23.9% of the total variance, could clearly separate PBMC-derived CD8⁺ T cells from the CD8⁺ T cells resided in tumor or peritumor tissue (P < 0.001). Although PC2 accounted for 16.8% of the total variance, no significant difference could be observed among different groups (P = 0.49). Furthermore, PC3 (accounted for 11.5% of the total variance) could nicely separate CD8⁺ TILs from peritumor-derived CD8⁺ T cells, with all CD8⁺ TILs from different patients have higher PC3 scores (P < 0.001). This excellent separation of CD8⁺ T cells from different tissue origins indicated that PC1 and PC3 might reflect the transcriptome changes affected by surrounding microenvironment.

Figure 3. PCA analysis on transcription profiles of CD8⁺ T cells. (a) The left subplot showed the distribution of CD8⁺ T cells from different tissue origins according to the first three PCs. Blue dots represent PBMC-derived CD8⁺ T cells, green dots represent peritumor-derived CD8⁺ T cells and red dots represent CD8⁺ TILs. The right subplot showed the relationship between tissue origin of CD8⁺ T cells and the three PCs, the significance was examined by one-way ANOVA. (b) Heatmap of normalized FPKM of genes that most correlated to PC1 and PC3. For better visualization, only top ten most correlated genes were shown in left subplot. Additionally, histone genes, HSPs and cell proliferation regulator genes positively correlated to PC3 were also shown in the right subplot. (c) The left subplot showed PC3 score of peritumor-derived CD8⁺ T cells (green) and PBMC-derived CD8⁺ T cells (red) from each patient. The right subplot displayed the correlation between PC1/PC3 scores and expression of known signatures of CD8⁺ T cell subtypes from HCC patients.¹³ Spearman’s rank correlation was performed and significant correlations (P < 0.05) were colored according to their correlation coefficient.

To further explore these two PCs, their most correlated genes were extracted. No surprising, genes significantly positively correlated with PC1 were overexpressed in PBMC-derived CD8⁺ T cells. These genes included two categories of T cell markers: naive T cell markers such as CCR7, SELL and LEF1, and effector T cell markers such as FGFBP2 and CX3CR1 (Figure 3(b)).^13–15 As expected, PC1 scores were positively correlated with naive and effector CD8⁺ T cell signatures (Figure 3(c)). These results supported the idea that naive and effector T cells were two main components of CD8⁺ T cells from PBMCs, which is in consistence with previous report.¹³

PC3 correlated genes could be divided into two groups with distinct expression tendency in CD8⁺ TILs and peritumor-derived CD8⁺ T cells. Genes negatively correlated to PC3 were significantly up-regulated in peritumor-derived CD8⁺ T cells, which included genes acting as inhibitors of inflammation or immune system, such as C1QA, C1QC, CEBPD,^16,17 together with some pro-inflammatory factors such as EGR1, CD161.^18,19 These results indicated that this part of PC3 could reflect the balance of pro- and anti-inflammatory process, which has been reported to play a crucial role in immunity response involved in protecting human body from inflammatory disease.^20,21 Interestingly, PC3 was also negatively correlated with HCC mucosal-associated invariant T (MAIT) cell signature (Figure 3(c)), indicating an enrichment of MAIT cells within CD8⁺ T cells from peritumor. On the other hand, genes positively correlated with PC3 showed significantly higher expression in CD8⁺ TILs, which consisted of genes related to both T cell exhaustion and proliferation. In this cluster of genes, several T cell exhaustion associated genes including HAVCR2 (also known as TIM-3), CTLA4 and TIGHT were identified.²² Furthermore, PC3 also showed positive correlation with exhaustive CD8⁺ T cell signature (Figure 3(c)). Moreover, other genes in this cluster were closely related to DNA synthesis and T cell proliferation, including histones, heat shock proteins (HSP), and many other key regulators involved in cell proliferation such as TOP2A and MKI67 (also known as Ki-67).²³

Interestingly, three patients with undetectable levels of serum HBV-DNA have the highest PC3 scores (sample 8306, 1377, 2129), with relatively higher expression levels of genes representing T cell exhaustion and proliferation. The higher expression levels of proliferation associated genes in CD8⁺ TILs may reflect the existence of more HBV-specific CD8 ⁺ T cells, which is required to sustain low levels of HBV. Likewise, hepatitis B patients with lower levels of HBV have also been reported to possess higher fraction of HBV-specific CD8 ⁺ T cells.²⁴ Moreover, co-expression of proliferation and exhaustion markers indicate that these two important biological processes of CD8⁺ T cells might be somehow linked. Additionally, the CD8⁺ TILs from the patient who received antiviral therapy (entecavir) before surgery (sample 7971) have the lowest PC3 scores, with the lowest expression of exhaustion genes and the highest expression of effector T cell markers. This result is in consistent with previous reports that T cells could be activated by antiviral drugs such as abacavir.²⁵ Collectively, our analysis revealed that HBV burden and antiviral therapy might affect T cell functional status in HCC and further studies with larger sample size may be valuable.

Transcriptional similarity and heterogeneity of CD8⁺ T cells in HCC patients

As demonstrated above, T cell function alterations showed great variety among different HCC patients. To intuitively visualize the transcriptional similarities and differences of each individuals, SOM analysis was conducted by translating transcriptional data into 1600 metagenes representing sets of co-expressed genes. The landscape of expression portraits for each CD8⁺ T cells was shown in Figure 4(a). Clustering samples according to these 1600 metagenes demonstrated that CD8⁺ T cells from the same tissue origin were usually transcriptionally similar (Figures 4(a) and 4(b)), while CD8⁺ TILs exhibited significant transcriptional heterogeneity, which is related to their clinical features.

Figure 4. Transcriptional similarity and heterogeneity of CD8⁺ T cells. (a) SOM analysis results showing the expression pattern of each sample. Grey rectangles indicate significantly upregulated spots for each sample. (b) Sample clustering according to all the metagenes. (c) Distribution of all identified significantly regulated spots. (d) Sample level enrichment scores of each spot.

In general, seven statistically differentially expressed spots (labeled A-G) were identified (Figure 4(c)), and each of them were named according to the gene set enrichment analysis of the protein coding genes within the spots, namely: A) naive T cell, B) Effector T cell C) memory T cell, D) MAIT cell, E) T cell proliferation, F) drug metabolism, G) infiltration (Detailed annotation information of the spots were provided in Figure 4(c) and Table S3). Sample level enrichment of each spot was shown in Figure 4(d).

Generally, all plasma CD8⁺ T cells displayed similar expression patterns, with all of them consistently overexpressing spot A and most (7/8) of them also overexpressing spot B. As to spot A, enrichment analysis showed significant overlap with several known naive CD8⁺ T cell associated signatures (Table S3). Spot B shows significant overlap with HCC effector CD8⁺ T cell signature, accompanied with significant enrichment in multiple integrin related pathways, suggesting that integrins are crucial for effector T cell function. This observation is in consistent with previous report that several members of integrin family take part into effector T cells’ migration.²⁶ Furthermore, the specially presence of spot A and spot B in plasma CD8⁺ T cells confirmed the idea that naive and effector T cells were prevalent in PBMCs.

On the other hand, peritumor-derived CD8⁺ T cells were characterized by increased expression of spot C and spot D. Genes from spot C were significantly overlapped with several memory CD8⁺ T cell signatures. The higher density of memory T cells in peritumor has been reported in previous studies.²⁷ Furthermore, spot D, which was specially expressed in peritumor-derived CD8⁺ T cells, represented the functional status of MAIT cells. MAIT cells are highly enriched in liver tissue and was considered as the first line of defense against viral infections.²⁸ These results suggested that MAIT cells were enriched in peritumor tissue and diminished in HCC tissue, indicating an impaired immune microenvironment in HCC tissue. Interestingly, peritumor-derived CD8⁺ T cells from patients with undetectable levels of serum HBV-DNA (patient 8306, 1377, 2129) and the patient received antiviral therapy (patient 7971) have higher expression levels of spot D, further proving the importance of MAIT cells in viral infection control.

Large heterogeneity was observed for CD8⁺ TILs, showing four spots of differentially expressed genes (C, E, F and G). For 3 of 4 patients with high level of serum HBV (patient 6423, 1967 and 2834), spot C (memory T cell spot) was reserved but expressed at a lower level in CD8⁺ TILs compared to corresponding peritumor-derived CD8⁺ T cells, indicating high level of serum HBV might benefit the retainment of memory T cells in tumor microenvironment

Meanwhile, CD8⁺ TILs from patients with undetectable level of serum HBV all have increased expression level of spot E, which was characterized by enriched cell proliferation related genes, including numerous histones, HSPs genes, both are critical for DNA synthesis, and some classical proliferation markers such as TOP2A and MKI67.²³ Spot E was largely resembled to PC3 positive correlated gene cluster, and both of them were upregulated in patients with undetectable serum HBV level. In addition, spot E was significant overlapped with two CD8⁺ effector T cell associated gene clusters,²⁹ suggesting a more cytotoxic phenotype of CD8⁺ TILs from patients with undetectable serum level of HBV.

Interestingly, patient 7898, the only patient with large amounts of T cell infiltration in tumor tissue (Figure S1 and Figure 5(a)), showed dramatic upregulation of spot F. No CD8⁺ T cell signature was significantly enriched for spot F, while PID pathway enrichment analysis revealed that this spot contained genes involved in T cell infiltration associated pathways including “CXCR3-mediated signaling events” and “IL6-mediated signaling events”.^30,31

Figure 5. Downregulation of IL12-mediated signal events is associated with worse clinical outcome. (a) T cell signature enrichment of peritumor-derived CD8⁺ T cells (green) and CD8⁺ TILs cells (red) from each patient using expression profiles of bulk cells. (b) Correlation of IL-12 signature enrichment between CD8⁺ T cells and corresponding bulk cells. P value and coefficient of correlation were calculated using spearman’s rank correlation. (c) Scatter plot of T cell signature enrichment of TCGA samples included in survival analysis (N = 124). (d) Expression levels of IL-12 signature in HCC tissues and paired peritumor tissues from TCGA data set (N = 50). P values were evaluated by Wilcoxon matched-pairs singed rank test. (e) Kaplan–Meier plot of HCC patients stratified by IL-12 signature expression levels (N = 124). Patients expressing higher levels (higher than 75th percentile) of IL-12 signature had better disease-free survival than other patients (log-rank P = 0.021).

The only patient who have received antiviral therapy before surgery (patients 7971) was characterized by overexpression of spot G. Though no known CD8⁺ T cell signature was significantly overlapped with spot G, GO analysis indicated that this spot was closely related to “drug metabolic process”. Interestingly, spot G contained multiple Cytochrome P450s (CYP) genes such as CYP2D6, CYP2E1 and CYP3A4, all of which play critical roles in drug metabolism.³² Their existence may reflect the influence of the antiviral treatment.

Intriguingly, CD8⁺ TILs from patient 6423 showed overexpression of three spots (C, E and F) simultaneously. Since sample 6423 had largest tumor and his tumor had occupied in multiple regions of liver including liver segments IV, VI, VII and VIII), we speculated that his clinical conditions could lead to higher heterogeneity of TILs.

The above results demonstrated that the CD8⁺ TILs’ functional phenotype was greatly influenced by patients’ clinical course, which in turn could result in different clinical outcomes.

Clinical correlation between IL12-mediated pathway and HCC prognosis

Since IL-12 mediated pathway has been identified in two of these spots (spot C and spot D) that closely related to CD8⁺ T cell phenotype, and its vital function has been also demonstrated in DE gene analysis, we further investigated how this pathway correlated with patients’ prognosis. Genes involved in IL-12 mediated pathway within spot C and spot D were extracted and defined as IL-12 signature. Since most of IL-12 signature genes were interleukins, chemokines and TNFs, which are mostly exclusively expressed in T cells, we hypothesized that bulk sequencing of corresponding tissue samples could reflect their expression in CD8⁺ T cells. This hypothesis was confirmed by the significant correlation of the expression patterns in CD8⁺ T cells and tissue samples except patient 7898 (P = 0.005, correlation coefficient = 0.701, Figure 5(b)). However, extremely high pan-T cell infiltration in HCC tissue of this patient was observed by immunohistochemistry (IHC) staining (Figure S1), which could explain the disruption of correlation. A pan-T cell signature consisted of 13 genes was also adopted to evaluate the T cell infiltration degree via bulk sequencing data,³³ which provided consistent results to show high pan-T cell infiltration level in HCC tissue of patient 7898 (Figure 5(a)). Further analysis with TCGA dataset also suggested the existence of a very small proportion of HCC patients with high level of T cell signature related to T cell infiltration (Figure 5(c)). We then compared the expression levels of IL-12 pathway with patients’ prognosis using TCGA data. Interestingly, we find significantly decreased expression levels of IL-12 signature in tumor tissue compared with peritumor tissue (P = 1.2 × 10^–6,Figure 5(d)) and lower expression levels of IL−12 signature indicated a worse disease-free survival (log rank P = 0.021). The expression of IL−12 signature could reflect the functional activation of CD8+ T cells, which could be further utilized for both immunotherapy development and patients’ prognosis evaluation.

Discussion

Serving as an essential component of adaptive immunity, CD8⁺ T cells play important roles in antitumor immune responses.⁴ Activation and dysfunction of T cells during the process of tumor progression have become the major focus of related research fields, particularly after breakthrough of immunotherapy in cancer treatment. Our analysis presented a detailed depiction for the dynamic change of CD8⁺ T cells from different origins on both genomic and transcriptional levels, showing that the functional alterations of TILs could greatly vary among patients.

After infiltrating into the inflammatory or cancer tissues, CD8⁺ T cells face the influence of the changed microenvironment.⁸ The interaction between T cells and their surroundings can significantly affect their functions on multiple levels. Except the commonly known transcriptional alterations, the genome of CD8⁺ T cell is also under the pressure of mutational processes which is similar to somatic mutation accumulation in cancer cells. Our analysis did prove that small number of somatic mutations, including both SNVs and CNVs, can indeed take place on CD8⁺ T cell’s genome from both tumor tissues and peritumor tissues, however, the functional alterations caused by these mutations seem to be weak compared with transcriptional changes. On transcriptional level, similar transcriptome changes were observed when CD8⁺ T cells entered into tumor or peritumor tissues, indicating similar microenvironment surrounding T cells could shape the CD8⁺ T cells into similar transcriptomes or phenotypes. Our results identified multiple connected pathways involved in the functional regulation of T cells, among which IL-12 signaling pathway was one of the most important ones. The expression pattern of IL-12 pathway was in consistent with the expected T cell activity, since CD8⁺ T cells in peritumor were considered to be activated to against tumor invasion while become exhausted in HCC microenvironment.^12,34 This correlation suggested that this pathway might dynamically reflect T cell functional status and could be utilized for T cell functional evaluation or immunotherapy development. Noteworthily, IL-12 also plays critical roles in natural killer (NK) cell activation,^35,36 which might provide a paralleled/joint effect in tumor immunology. However, whether or how these two immunity mechanisms interact with each other still needs further investigations.

Substantial functional difference among TILs from different patients has been reported,³⁷ while our analysis further revealed that the difference of TIL function was greatly affected by patients’ clinical conditions including HBV-DNA level, anti-viral treatment and degree of infiltration. One intriguing fact is that HBV infection level might somehow reflect antitumor immune activity, since patients with low level of HBV tended to show higher expression of T cell proliferation signature. One possible explanation is that proliferation of CD8⁺ T cells might contribute to the low HBV-DNA microenvironment while excessive T cell proliferation lead to exhaustion. Consistent with our hypothesis, CD8⁺ T cells have been reported to play roles in viral control during hepatitis B virus infection, and higher frequency of HBV-specific CD8⁺ T cells have been detected in hepatitis B patients with a lower level of HBV replication.³⁸ Moreover, antiviral treatment also contribute to T cell functional alteration, since TILs from the patient who received antiviral treatment before surgery showed increased expression of effector T cell markers and reduced expression of T cell exhaustion markers. Previous reports have suggested that T cell would be activated by anti-HBV drugs such as abacavir,²⁵ and entecavir which has also been reported to improve the function of T cell immunity.³⁹ Additionally, several epidemiological investigations suggested that antiviral treatment could improve the survival and decrease the recurrence of HBV-related HCC.⁴⁰ Taken together, serum HBV level might be utilized for patients’ immunity response evaluation, while combine therapy with antiviral drugs in antitumor treatment could greatly benefit for T cell activation in HBV-related HCC patients.

In conclusion, our analysis revealed the multi-dimensional profiles of tumor infiltrated CD8⁺ T cells in HCC patients and demonstrated that how their functional change was related to patients’ clinical course. Our study also could further improve the fundamental understanding of CD8⁺ T cell functional alteration in HCC by discovering multiple inter-connected pathways in T cells activation. Selectively targeting the common hub genes in these pathways may provide a more effective strategy for future immunotherapy development.

Materials and methods

Sample collection

Eight HCC patients who received surgical operation in Mengchao Hepatobiliary Hospital of Fujian Medical University were recruited. None of these patients received any chemotherapy and radiation before surgical resection. Primary tumor tissue, paired peritumoral tissue (at least 2 cm away from the primary HCC lesion) and corresponding peripheral blood were collected after hepatectomy. A written informed consent was obtained from each participant. All of the study protocols were approved by the Institution Review Board of Mengchao Hepatobiliary Hospital of Fujian Medical University and were conducted according to the principles of the Declaration of Helsinki.

Single cell suspension preparation and CD8⁺ T cell purification

Fresh tumor and peritumor excisions were minced to small pieces in 35mm dishes, and then dispersed in 10 ml RPMI 1640 containing 0.1% type I collagenase (Solarbio, China, C8140), 0.01% hyaluronidase (Solarbio, China, D8030) and 0.002% DNA enzyme (Solarbio, China, D8070) at 37 ◦C for 30–40 min. Further filtering was conducted with 40μm BD Falcon Cell Strainer to remove undigested tissues. The cell suspensions were then washed and re-suspended with PBS. Isolation of mononuclear cells were conducted with discontinuous percoll density gradient centrifugation (density gradient: 30% and 60%) in 50ml tubes. The tubes were centrifuged at 800g for 20 min. Enriched cells between 60% percoll and 30% percoll solution were gathered, and washed with PBS. Matched peripheral blood mononuclear cells (PBMCs) were isolated from peripheral blood collected before surgery by Ficoll-Paque PLUS kit (GE Healthcare).

Afterwards, single-cell suspensions were incubated with anti-CD8 magnetic beads (Milteny, Germany, 130–045-201), according to the manufacturer’s protocol. The cells were then washed and separated using charged MACS columns (Milteny, Germany, 130–042-201). The columns were washed and the flow-through was discarded. The columns were then removed from the magnet and the CD8⁺ T cells were eluted from the column. Finally, the cells were stained with CD3⁺ and CD8⁺ antibody and subjected to FACS analysis for cell purity evaluation and CD8⁺ T cells with high purity (≥ 90%) were subjected to subsequent sequencing.

Exome and transcriptome sequencing

Total genomic DNA and RNA from acquired CD8⁺ T cells and corresponding tissue samples were extracted using QIAamp DNA mini Kit (Qiagen) and RNeasy Mini Kit (Qiagen) respectively. Exome regions of qualified DNA were enriched using Agilent SureSelect Human All Exon V6 Kit. Genomic DNA were subjected to DNA library preparation, and exome regions of qualified DNA were enriched using Agilent SureSelect Human All Exon V6 Kit. RNAseq libraries were prepared from total RNA with removal of ribosomal RNA (rRNA) according to manufacturer’s instruction. Both enriched exome and transcriptome libraries were sequenced by Ruibo (Beijing, China) on Illumina HiSeq3000 platform (paired end, 150 bp).

WES data analysis

All sequencing reads were first evaluated for sequencing quality to remove reads with low quality based on the following criteria: containing adaptor sequence; with > 10% ambiguous bases (N) of total read length; containing > 20% low quality base (base quality ≤ 20) of total read length. WES Sequencing read pairs were aligned to UCSC human reference genome hg19 (GRCh37) using BWA (v0.7.15) with default parameters. Possible PCR duplicates were marked and removed using Picard (v2.9.3).⁴¹ Both local insertion/deletion (indel) realignment and base quality score recalibration were performed using Broad’s Genome Analysis Toolkit (GATK, v3.6). Somatic mutations of T cells were detected with Mutect2 in GATK (v3.6), using matched T cells from PBMC as normal reference.⁴² Mutations were then further filtered with following criteria: ≥ 20X depth in both samples of comparison; variant reads ≥ 5 and variant allele frequency ≥ 10%. CopywriteR (v2.10.0) and cn.MOPS (v1.24.0) were used to infer somatic CNVs.^43,44 For CopywriteR, bin size of 20 kb was used and the depth of coverage were corrected for both GC-content and mappability. The output of CopywriteR was further analyzed by circular binary segmentation (CBS) algorithm to determine segment mean value of each region. The cut-off value for duplications was set as segment mean value > 0.25, and deletions was set as segment mean value < −0.25. Meanwhile, cn.MOPS was run with default parameters and regions < 1kb were excluded. Finally, CNVs identified by CopywriteR with more than half segment overlapping with cn.MOPS results were included for further analysis.

Transcriptome data analysis

The quality control for raw sequencing reads were performed using the same criteria as WES data. Qualified reads were aligned to ribosomal rRNA sequences download from RNAcentral database.⁴⁵ Then unmapped reads were aligned to human genome reference (GRCh37) with GENCODE gene annotation using STAR.⁴⁶ LncRNAs were further defined as those with transcript biotypes of “antisense”, “3prime_overlapping_ncRNA”, “lincRNA”, “misc_RNA”, “processed_transcript”, “sense_intronic” and “sense_overlapping”. Transcripts shorter than 200 bp or showed low expression level (average reads count < 10 within all samples) were excluded. Protein coding genes were defined as those with gene biotype of “protein_coding” and genes with no mapped reads across all samples were removed. To quantify the expression levels for each gene, fragments per kilobase of exon per million mapped fragments (FPKM) were calculated. Differentially expressed protein coding genes and lncRNAs between different group of samples were identified using DESeq2,⁴⁷ and limma.⁴⁸ Significantly changed genes with adjusted p value < 0.05 (Benjamini-Hochberg correction) were selected for Pathway Interaction Database (PID) pathway enrichment using the web tool ConsensusPathDB.⁴⁹

Principal component analysis (PCA)

PCA was performed on log₂-transformed expression matrix consisting of all the qualified protein coding genes and lncRNAs. The relationships between tissue origin of CD8⁺ T cells and the first three PCs were examined by one-way ANOVA. Top genes positively correlated with PC1/PC3 and negatively correlated with PC3 were extracted according to their coefficients of eigenvectors. To better characterize different PCs, the correlation between overall expression of reported signatures representing different subtypes of CD8⁺ T cells and the PC scores was evaluated by spearman’s rank correlation. Single-cell RNA-seq data for different subtype of HCC lymphocytes was download from Gene Expression Omnibus database (accession number: GSE98638).¹³ Expression data of naive, effector, mucosal-associated invariant T (MAIT) cells and exhausted CD8⁺ T cells were subjected to subtype specific signature gene identification using prediction analysis for microarrays (PAM).⁵⁰ Top 50 genes with highest PAM scores within each subtype were then defined as signature genes for corresponding subtype of CD8⁺ T cells.

Self-organizing map (SOM) analysis

SOM analysis was conducted using R package oposSOM⁵¹ with default parameters. A total of 20,000 RNAs including top 10,000 protein coding genes and 10,000 lncRNAs showing highest variance across all samples were selected for dimensionality reduction and clustering. Input transcriptional data were translated into 1600 metagenes, each of which represented a set of co-expressed genes and can be visualized by a two-dimensional 40 × 40 mosaic grid. Metagenes of similar expression patterns were clustered and localized in neighboring regions, forming so-called spots of co-expression modules. In addition, all CD8⁺ T cell samples were clustered by pairwise correlation clustering according to the expression of all the 1600 metegenes. Genes from each spot were extracted (FDR< 0.05) for PID pathway enrichment analysis using ConsensusPathDB.⁴⁹ Meanwhile, gene sets from aforementioned signatures of HCC CD8⁺ T cell subtypes,¹³ CD8⁺ T cells related C7-immunologic signatures and C5-gene ontology (GO) gene sets taken from the Molecular Signatures Database (MSigDB) were also included for gene set enrichment by Fisher’s exact test. Gene sets with P < 10^–10 were considered indicative. In addition, sample level enrichment score of each spot was calculated by single sample gene set enrichment (ssGSEA) using R package GSVA.

Prognostic assessment using TCGA data

RNA-seq and clinical data of HCC patients in The Cancer Genome Atlas (TCGA) was download from TCGA website (http://www.cancergenome.nih.gov). Patients with available transcriptome data for both tumor and adjacent normal tissues were extracted (N = 50) to examine the difference of IL-12 signature expression levels by Wilcoxon matched-pairs singed rank test.

In addition, patients with more than two years’ follow-up were extracted from TCGA dataset for survival analysis (N = 124). IL-12 signature expression level of each individual sample was calculated as the mean z-score of all signature genes. Survival curves of TCGA patients with high (> 75th percentile) or low (< 75th percentile) expression levels of IL-12 signature were estimated using the Kaplan-Meyer method, and the log-rank test was performed to examine the significance.

Author contributions

XL Liu, JF Liu, and ZX Cai - study concept and design; ZL Li, G Chen, and XQ Dong - patient selection and sample collection; HP Xu, YY Zeng, G Chen, LM Qiu and XQ Dong - clinical information collection; ZX Cai and XQ Dong - sequencing data acquisition; ZL Li and G Chen - analysis and interpretation of data; XL Liu and JF Liu - technical support, obtained funding and study supervision; all authors have read and edited the manuscript.

Declaration of interest statement

All authors declare no competitive interests related to this study.

Ethics approval

All of the study protocols were approved by the Institution Review Board of Mengchao Hepatobiliary Hospital of Fujian Medical University.

Supplemental Material

Supplemental data for this article can be accessed here.

Additional information

Funding

This work was supported by the Science and Technology major project of Fuzhou (Grant No. RongKe 16 in 2016); the Scientific Foundation of Fujian Health and family planning Department (Grant No. 2018-1-92); the Scientific Foundation of Fuzhou City (Grant No. 2017-S-131-3); the Startup Fund for scientific research of Fujian Medical University (Grant No. 2016QH085, Grant No. 2017XQ1165); the Scientific Foundation of Fuzhou Health Department (Grant No. 2017-S-wq38); the Science and Technology development project of central government guiding local government (Grant No. 2017L3017, Grant No. 2018L3016);Startup Fund for scientific research, Fujian Medical University [2016QH085];Startup Fund for scientific research, Fujian Medical University [2017XQ1165];