Biomarkers Journal Open Access

Keywords

Personalized immunotherapy; Tumor-infiltrating lymphocytes (TILs); CD8+cells; Quantitative pathway; Network; Gene expression signature (GES); Biomarkers; Heterogeneous responses

Introduction

Quiescent T-cell has been defined by small cell size, lack of spontaneous proliferation and low metabolic rate in some special environment such as in tumor tissue [1]. In order to study quiescent T-cell in tumor tissue, we have used single-cell technique to isolate the T-cells in tumor tissue which have already contacted tumor cells and then performed single cell genomics analysis to study their genomic profiles related quiescent T-cell network [2]. Our studies have reached the following conclusions [3]: (1) TIL CD8+cells are an actively maintained quiescent status in tumour microenvironment as compared to the default status in the absence of the stimulating signals; (2) CD8+cells maintain the quiescent status by activating several regulatory pathways; and (3) these activated signalling pathways finally inhibit cell cycling to maintain the quiescent state of TIL. In molecular biology. some scientists supported the genomic profiles related network. For example, lung-Krüpple-like factor (LKLF, a zinc finger-containing transcription factor) has been discovered to maintain T cell quiescence [4]; Tob (nuclear protein) was uncovered to have anti-proliferative activities by blocking T cell receptor (TCR) engagement in the presence of either CD28 co-stimulation or IL-2 [5].

Because some or most of CD8+ cells are quiescent in tumor microenvironment, the conventional TIL immunotherapy includes harvesting T-cell from patients’ tumor tissues, stimulating T-cell by IL-2 under ex vivo culture and then administering to patients in vivo for immunotherapy. However, current evidences have indicated that responses of CD8 cells are “heterogeneous”, that is, T-cells will produce different genomic profiles in different patients under GWAS analysis [6,7]. In this study, we performed in silico analyses of special biomarkers [also called as gene expression signatures (GES) as genomic terms] related quiescent networks which they were obtained from genomic data of quiescent CD8+ cells harvested by single-cell technique from a pair of tumor specimens with similar liver cancer patients. The manual purpose is to study how to select specific therapeutic targeting molecules or stimulating agents to perform individual “heterogeneous” immune responses for further personalized immunotherapy.

Methods

mRNA expression validation

We have reported genomic results of quiescent CD8+ cells or CD3+ cells from single TIL of two liver cancers. Specimen-1 had primary liver cancer from male, 43 years old with 7 years HBV positive carrying and specimen-2 was primary liver cancer from male, 37 years without HBV and HCV carrying history [8]. All quiescent CD8+cell expression profile was confirmed by quantitative real time PCR. Briefly, as our previous report [9], the real time PCR assay was performed in triplicates for each gene. The PCR reaction mixture is totally 50 ul containing 25 ul 2x SYBR Green (BioRad), 500 nM each primer, RNA extract and 1 ul iScript reverse transcriptaseTM. According to the primer conditions and manufacture’s recommendations, one step realtime PCR was performed by incubating at 10 min at 50°C and 5 min at 95°C, followed by 45 cycles of denaturation for 10 s at 95°C and annealing/extension for 30 s at 55°C. The SYBR fluorescent signal intensities were recorded and analyzed during PCR in an ABI Prism 7700 sequence detector system (Applied Biosystems) using the SDS (Ver. 1.91) software.

Network construction and layout

Definition

To characterize the network and layout [10], we first define the terms as below:

Seed-proteins

Proteins (nodes) whose genes were differentially expressed in single cell mRNA differential display and were involved in quiescent status or active status.

Non-seed-proteins (neighbours)

Proteins (nodes) selected as a consequence of in silico experimentally validated interactions starting from seed-proteins.

Network

Network that includes seed-proteins and non-seed-proteins (neighbours) and their interactions. Only seed-proteins linked to neighbouring proteins were included in the network analysis.

Network construction

After gene expression signatures of quiescent or active CD8 cells were profiled using single cell mRNA differential display, the differentially expressed genes were used in pathway analysis and network construction. Briefly, seed proteins (seed nodes) were input into Pathway Studio 5.0 software. Starting from seed proteins involved in either quiescent status or active status, we obtained a network through the interaction of these proteins with their direct neighbours. A general scheme of our approach included (1) experiments were used as Expand methods (both edge in and edge out); (2) node filter (entity filter) was regulatory connectivity (protein-protein reaction) and pathway; (3) parameters used to generate “network” in the “String database” were one step method (only direct neighbours); (4): parameters used to generate the relation filter were proteinprotein interaction (PPI). This configuration implies that only the experimental evidence with direct interactions was extracted from the database as valid links for each quiescent or active CD8 network. A detailed description of each parameter can be found in Pathway Studio 5.0 [11]. We did not consider either the indirect interaction of proteins or indirect self-interactions. All datasets with the regulatory connectivity and calculated betweenness of each protein node within these networks were observed from the Pathway Studio 5.0 software.

Topology analysis for genomic signatures

Definition

Connectivity

Node connectivity (node degree) is defined as the total number of edges that connect to a given node.

Regulatory connectivity (or Degree Centrality, DC)

Regulatory connectivity is defined as the total number of protein edges that connect to a given protein node.

Betweenness (or Betweenness Centrality, BC)

The network diameter is defined as the shortest optimal path between any pair of nodes divided by average pathway length in the network.

The network parameters and calculation

DC

DC is protein-protein reaction in the analytic network regarding pathway and signal transduction as previously defined. Briefly, networks are generated by loading the proteins which are related with quiescent CD8 or active CD8 cells into the Pathway Studio 5.0 software with the following parameters: “Expand Methods” (both edge in and edge out), parameters, “One Step Method” (only direct neighbors) and “regulatory connectivity” (protein-protein reaction). The number of all paths in the regulatory networks was recorded for each gene/ node which was differentially expressed in quiescent or active CD8 cells. A detailed description of the regulatory connectivity analysis is available as the Pathway studio software 5.0 manual.

BC

To calculate node BC within networks, we used the traditional formulation defined as below [12].

BC has been established as an important quantity to characterize the communications between each pair of nodes. The communication paths: σst (v) is between a pair of nodes (s, t), which are the shortest pathway through v point and σst is all number of such pathways from s points to t points. Betweenness (BC) means σst (v) denoted by σ (s, t). Betweenness centrality is, in some sense, a measure of the influence that a node has over the flow of information through the network. In theory, high betweenness nodes lie on a large number of non-redundant shortest paths between other nodes. In fact, they can thus be thought of as “choke point”, which is very important parameter to control the pathway regulation or network.

Biomarker analysis by GES

After topology analysis by using BC and DC, the set of differentially expressed biomarkers is identified as “topologically significant”, i.e., “gene expression signatures, GES”. Therapeutic identification was performed using GES according to the GeneGo software manual. Briefly, the candidate genes were input into the GeneGo software. GESs were then used to construct networks that connected corresponding nodes in the global database of drugs, small molecules and other agents.

Software availability

GES were processed by the Pathway Studio 5.0. Topology study is performed by the public software Baysialab in the Pathway studio 5.0. The GeneGo software [13] was used for therapeutic identification and drug targeting.

Experimental confirmation

Harvest and culturing

TIL of tumor tissues were isolated and harvested as our previous publication [14]. Briefly, after tissue was minced into ~ 1 mm³ pieces and collagenase IV enzyme digestion medium for 4°C overnight, filtered cells were separated by Ficoll gradient centrifugation (600 g with 30 min. at 25°C) and the leukocyteenriched layer collected for TIL growth. The cell suspensions were placed at 5 × 10⁵ total cells per ml CM (RPMI-1640 containing 0.1 ug/ml of human recombinant IL-2 and 10% human serum with 100 U/ml penicillin and 100 μg/ml streptomycin) in a 37°C, 5% CO₂ humidified incubator for the first four days and then added inhibiting compounds relied on GES for further 72 hours inhibition.

Cell and molecule characteristics of heterogeneous responses

After 72 hours of inhibitor addition, cultured cells were analyzed by cell viability and cell growth inhibition, quiescent gene expression inhibition assay and TIL cytotoxic activity. mRNA expression was studied by Q-rtPCR assay and measurement were shown as previously published [15]. Briefly, after T-cells were harvested and RNA was extracted by Trizol solution, Q-rtPCR analysis was used to study targeting gene to study mRNA expression level.

To study TIL cytotoxic activity in heterogeneous responses of quiescent status, cytotoxic T lymphocyte (CTL) assays with cocultured with 5 × 10³ target cells were defined by MTT (MTT Cell Growth Assay Kit, cat, CT02, Millipore Inc, USA) with a calculation formula as our previous reports [16]. Briefly, after TIL were cultured at day-7, 5 × 10³ autologous tumor cell was loaded in each well. After cultured for 24 hours for tumor cells, 1 × 10⁵ TIL (at a 50:1 effector: target ratio against target cells) were loaded into each wells and then incubated for additional 24 hour. Mechanism of TIL cytotoxic activity for heterogeneous responses in quiescent status also was measured by Q-rtPCR and Western blot using Anti-Granzyme B antibody (ab53097, Abcam Inc, USA), Anti-Perforin antibody (ab97305, Abcam Inc, USA) to confirm the protein expression level of the heterogeneous responses in quiescent status.

Results

Subtle difference of gene expression in a pair of specimens

In order to study individual “heterogeneous” immune responses, we selected 8 candidate genes to study gene expression level from a pair of two liver cancer patients which were already confirmed CD8+ cell quiescent genes from single-cell libraries of TIL CD8+ cell (312 and 288 colonies in two specimens, respectively). After measuring mRNA level change by Q-rtPCR obtained from quiescent CD8+ cells of the two liver tumor specimens, eight candidate genes were subtly different as summarized in Table 1 and Figure 1A and 1B. Five genes (Myc, Tob, ERF, Ets-2 repressor factor, REST and TGF-β) had a high level of gene expression in specimen-1 and seven genes (KLF2, Ski, Sno-A, Myc, Tob, ERF and REST) were confirmed a high level of gene expression in Specimen-2 after 1.5 fold cut off.

Table 1: Genomic expression in both specimens.

Figure 1: mRNA expression level from both quiescent CD8 cells..

Network construction

Although mRNA level from quiescent CD8+cells obtained from the two liver tumor specimens are subtle change, in order to study individual “heterogeneous” immune responses in the system biology, three network were constructed by Pathway Studio 5.0 software, control-group (all of eight quiescent genes), specimen-1 (five quiescent genes) and specimen-2 (seven quiescent genes). In control-group, after all eight (8) seed proteins were loaded into the Pathway Studio 5.0 software, totally 196 nodes were discovered in the network profiles of quiescent CD8 cells, including 8 seed-proteins and 188 neighbours which are derived (Figure 2). As Figure 2, one of eight seed protein, Sno-A, has not any links resulting in non-involvement in the network analysis (i.e., there was no in silico evidence of interactions). As Figure 2A and 2B, TGFB1, one of seed protein, is growth factor. Tob, one of seed is membrane expression protein and other four of them (c-MYC, REST, ERF and KLF2) are transcription factors located at nucleus obtained from the networks. As Figure 2C and 2D, network is constructed symmetrical platform and force power platform for further topology analysis.

Figure 2: Network and layout of common pathways from quiescent CD8 cells- A and B are networks constructed at cellular structure level (A is plain cellular platform and B is circular cellular platform); C is the network building from symmetrical platform and D is the network construction using force power platform..

In order to study heterogeneous responses in the networks for individualized therapeutic targeting, we further analysed networks from specimens-1 and 2, respectively. After five seed proteins from specimen-1 were input into the platform, the resulting network contains 5 seed proteins and 155 neighbours derived from direct interaction (Figure 3). As Figure 3A and 3B, TGFB1, one of seed protein, is growth factor. Tob, one of seed is membrane expression protein and other three of them (c-MYC, REST, ERF) are transcription factors located at nucleus obtained from the networks. Similarly to specimen-1, we studied networks from specimen-2, which contain 7 seed proteins. Eventfully, 107 neighbours of totally 114 proteins were constructed as Figure 3C and 3D. Again, Sno-A, a seed protein, had no links in this network analysis as shown in Figure 2. Although only four proteins are different between these two samples, i.e., TGFB1 in specimen-1 and Ski, Sno and KLF2 in specimen-2, the resulting networks (Figure 3A and 3B vs. Figure 3C and 3D) have shown some differences. The results revealed that only a few of gene change maybe result in network change in different individuals.

Figure 3: Network and layout of individual pathways in Specimen 1- A and B are networks constructed at cellular structure level (A is plain cellular platform and B is circular cellular platform) from specimen-1; C and D are networks constructed at cellular structure level (C is plain cellular platform and D is circular cellular platform).

According to our previous report, after activation of IL2, all four genes (pre-TCR, TRAIL, Perforin and TNF receptor) have a high expression in activated CD8 cells. In order to study heterogeneous responses in an activated network, as compared to the quiescent CD8 cells from both specimens, we also constructed network from four genes obtained from activated CD8+ cells. In this analysis, the network contains 4 seed proteins, 81 neighbours and totally 85 proteins in the whole networks as Supplement figure 1A-1D.

Topology analysis

We first studied eight quiescent genes listed in the common network from both specimens with 196 nodes, including 8 seed proteins and 188 non-seed proteins. The results of topology regarding regulatory connectivity are summarized in Table 2 and displayed in Figure 4. The higher regulatory connectivity of a protein means more direct interactions with other proteins, such as c-Myc (85) and TGFB1 (63). The lower regulatory connectivity of a protein means fewer direct interactions with other proteins, such as Tob1 (4), KLF2 (2), REST (13) and Ski (18).

Table 2: Regulatory connectivity and betweenness centrality.

Figure 4: Regulatory connectivity of eight seed proteins- Sno has not any linking so that it is omitted from the figure..

Betweenness centrality (BC) is very important measurement for downstream targeting hub. Here we used a formula to study BC, or a pair of node (s, t), the shortest pathway through v point is denoted by σ (s, t). Although regulatory connectivity is similar to shared proteins in both specimens as results shown above, the values for BC are greatly different in some hubs such as c-Myc, ERF, REST and Tob1 between specimen-1 and specimen-2 as shown in Table 2 and Figure 5A-5C.

Figure 5: Betweenness from Specimen 1 and Specimen 2- All A, B and C are rebuilt by force power platform. A is control networks from common pathways; B is individual networks of Specimen 1 and C is individual networks from Specimen 2..

GES analysis

Recently, therapeutic targeting such as drug targeting, small molecule targeting, Ab targeting and RNA-interfering therapy has focused on GES based on topology analysis. Most publications focused on two parameters in topology, DC and BC as a measurement. Theoretically, both BC and DC all play an important role in cell function while DC is also likely to be toxic due to their system-wide influence, thus we search GES with higher BC and low DC. A high BC value indicates a significant targeting node because it is very important hub in the cell functional network and low DC means very few branches without their system-wide influence to cause whole cell functional toxicity. This study has similar DC so that we analyzed higher betweenness as key nodes. After GES were defined by topology analyses, the nodes were loaded into the GeneGo softwere to mine drugs, small molecule and other molecular therapy agents. The resulting drug candidates were combined with GES and loaded into the Pathway Studio software to predict therapeutic targeting. For example, we predicted that drugs like etoposide and taurine could inhibit quiescent CD8 cells in specimen-1 as Figure 6A and loratadine and taurine could directly inhibit quiescent CD8 cells in specimen-2 as Figure 6B.

Figure 6: Combination results from mining therapeutic targeting and mining gene expression signature- A is a diagram of combination results of therapeutic targeting mining with GES mining from individual networks of Specimen-1; B is for Specimen-2; C is after drug treatment GES change for Specimen 1 and D is after drug treatment GES change for Specimen-2..

Experimental confirmation

According to previous evidence for TIL CD8+cells quiescent status, we studied two set of compounds at culture day-4 with 72 hour induction. After culture at day-7, inhibition of targeting gene expression were confirmed that etoposide could, respectively, inhibit 65% and 87% REST and ERF gene expression but without SKI gene inhibition; taurine can inhibit 88% and 87% REST and ERF gene expression but without SKI gene inhibition while loratadine have not any obvious inhibition for REST, ERF and SKI in specimen-1as Figure 6C. Taurine inhibited ERF and SKI expression and loratadine only inhibit SKI in specimen-2 as Figure 6D.

Following inhibiting analysis of CD8+cells, cytotoxic T lymphocyte (CTL) assays were used to study cytotoxicity of TIL activity. Cytotoxic T lymphocyte (CTL) was measured for effector:target ratio by 1:50 using autogenously tumor cell as Figure 7A. The results of killing activity showed 61 + 5.6%, 45 + 9% and 21.2 + 3% inhibited by etoposide, taurine and loratadine, respectively, in specimen-1 and 19.2 + 3.2%, 48 + 6% and 45 + 4%, respectively, in specimen-2. Mechanism of heterogeneous responses of TIL was confirmed by a Western blot using anti-granzyme-B antibody and anti-perforin antibody. Results of Q-rtPCR and Western blot demonstrated granzyme-B and perforin were highly expressed by etoposide and taurine to inhibit REST and ERF at specimen-1 while granzyme-B and perforin were highly expressed by inhibition of loratadine and taurine to SKI and ERF in specimen-2 as Figure 7B and 7C. As most of inhibiting compounds, we also see that drugs interfere TIL growth curve to compare normal TIL cells as Figure 7D, and viability as Supplement Figure 2, especially by Loratadine inhibition.

Figure 7: Experimental confirmation - A is a diagram of CTL combination results of therapeutic targeting mining with GES mining from individual networks of Specimen 1 and 2. B is mRNA expression change after inducing quiescent pathway for Specimen 1-2. C is Western blot after inducing quiescent pathway for Specimen 1-2. D is cell growth after drug treatment GES change for Specimen 1 and 2..

Discussion

In order to study heterogeneous responses in network, we tool to study the complexity in heterogeneous responses of CD8 cell quiescent and active status. Many genomics publications have extensively reported network construction. However, as we have demonstrated that a tiny change (including noise) of key nodes can be magnified in network analysis. This drawback has demanded scientists and physicians to develop more accurate genomics analyses. In order to overcome this kind of challenge, we applied for data generated from single cell mRNA differential display to analyse network. We selected the seed proteins obtained from differentially expressed genes in single cell level and analysed networks based on these seed proteins and their neighbour proteins.

Up to now, genomics profiles have been successfully used in analysis for clinical application [17]. New generation medicine such as “precision medicine” has encouraged physicians to develop the accurate genomics diagnosis with network construction to apply for the new fields. Here, we first used the data generated from “pure” genomic data to study a network. This research aim is subject to study subtle individual immune response to further apply for personalized immunotherapy.

Network studies including BC and DC play an important role in therapeutic targeting [18]. In these data from two specimens, we obtained some similar nodes (including four genes: c-Myc, Tob-1, ERF and REST) from two specimens indicating a shared response in quiescent network. After network analysis to study the complexity of individual responses of CD8+ cell quiescent, once one or more proteins are added into network, such as TGFB1 increase from specimen-1, it could result in individual pathway shifts in network. Although DC and BC all play an important role in therapeutic targets, here a high BC value is defined as significant targeting nodes because they belong to very important hubs in cell functional network without more toxicity [19]. Following system modeling results, we defined REST and ERF as GES for therapeutic targets in specimen-1 and Ski, Tob1 and ERF in specimen-2. Compounds inhibiting GES further support our network analysis. For example, etoposide and taurine could inhibit REST and ERF in specimen-1 and loratadine and taurine could directly inhibit Ski, Tob1 and ERF in specimen-2. Moreover, etoposide and taurine can significantly increase CTL in specimen-1 while taurine and loratadine increase CTL in specimen-2. Mechanism analysis of individual responses of TIL demonstrated granzyme-B and perforin were highly expressed by etoposide and taurine at specimen-1 while granzyme-B and perforin were highly expressed by inhibition of loratadine at specimen-2 as Figure 8. Our results demonstrated that an in silico network analysis can be employed for personalized immunotherapy [20].

Figure 8: Diagram of Gene expression signature for heterogeneous responses in specimen 1 and 2..

Acknowledgment

Under the support of former Dr. H. D. Preisler, we set up methods to analyse CD3, CD4 and CD8 from TIL. The work is supported by National Cancer Institute IRG-91-022-09 for Dr. Biaoru Li as P.I. This work was supported in part by USDA CRIS Project No. 1265- 31000-098-00D (to GEL). Mention of trade names or commercial products in this article is solely for the purpose of providing specific information and does not imply recommendation or endorsement by the U.S. Department of Agriculture.

Author Contributions

BL conceived, designed and performed the experiments. BL, HH and GEL performed genomic analyses. HH, JQ and BL performed bench work. BL, GEL, JZ and AT wrote the paper.

Competing Interest Statement

No competing financial interest is declared.

References

1. KilincMO, Rowswell-Turner RB, Gu T, Virtuoso LP, Egilmez NK (2009) Activated CD8+ T-effector/memory cells eliminate CD4+CD25+Foxp3+T-suppressor cells from tumors via FasL mediated apoptosis.JImmunol
2. Xu Y, Hu H, ZhengJ, Li B (2014) Feasibility of whole RNA sequencing from single-cell mRNA amplification. Genet Res In 2013: 1-8.
3. Zhang W, Ding J, Qu Y, Hu H, Lin M, et al. (2009) Genomic expression analysis by single-cell mRNA differential display of quiescent CD8 T cells from tumour-infiltrating lymphocytes obtained from in vivo liver tumours. Immunology 127:83-90.
4. HaalandRE, Yu W, Rice AP (2005) Identification of LKLF-regulated genes in quiescent CD4+ T lymphocytes. MolImmunol 42:627-641.
5. Tzachanis D, Li L, Lafuente EM, Berezovskaya A, Freeman GJ, et al. (2007) Twisted gastrulation (Tsg) is regulated by Tob and enhances TGF-beta signaling in activated T lymphocytes. Blood 109:2944-2952.
6. Chan WK, Rujkijyanont P, Neale G, Yang J, Bari R, et al. (2013) Multiplex and genome-wide analyses reveal distinctive properties of KIR+ and CD56+ T cells in human blood. J Immunol, 191:1625-1636.
7. Obara W, Ohsawa R, Kanehira M, Takata R, Tsunoda T, et al. (2012) Cancer peptide vaccine therapy developed from oncoantigens identified through genome-wide expression profile analysis for bladder cancer. J ClinOncol 42:591-600.
8. Li B, Chang T, Larson A, Ding J (2000) Identification of mRNAs expressed in tumor-infiltrating lymphocytes by a strategy for rapid and high throughput screening. Gene 255:273-279.
9. Li B, Perabekam S, Liu G, Larson A (2002) Experimental and bioinformatics comparison of gene expression between T cells from TIL of liver cancer and T cells from UniGene.J Gastroenterol 37:275-282.
10. Daraselia N, Wang Y, Budoff A, Lituev A (2012) Molecular signature and pathway analysis of human primary squamous and adenocarcinoma lung cancers. Am J Cancer Res 20:93–103.
11. Lim D, Kim NK, Park HS, Lee SH (2011) Identification of Candidate Genes related to Bovine Marbling using Protein-Protein Interaction Networks. Int JBiolSci 7:992-1002.
12. Brandes U (2001) A faster algorithm for betweenness centrality. Journal of Mathematical Sociology25:163–177.
13. Bessarabova M, Ishkin A, JeBailey L, Nikolsky T (2012) Knowledge-based analysis of proteomics data. BMC Bioinformatics13:13.
14. Li B, Tong SQ, Zhang XH (1994) A new experimental and clinical approach of combining usage of highly active tumor-infiltrating lymphocytes and highly sensitive antitumor drugs for the advanced malignant tumor. Chin Med J (Engl) 107:803-807.
15. Perambakam S, Li B, Preisler H (2001) Quantitation of interferon regulatory factor transcripts in patients with acute myeloid leukemia. Cancer Invest 19: 346-351.
16. Wang JH, Tong SQ, Li B (2000) Immunologic character of tumor infiltrating lymphocytes in ovarian carcinoma. Chinese Journal of Cancer Research 12: 99-104.
17. MortazaviCA (2014) The common ground of genomics and systems biology. BMC SystBiol 8:1-2.
18. Bonaz C (2013)Inflammatory bowel diseases: a dysfunction of brain-gut interactions. Minerva GastroenterolDietol. 59:241-259.
19. Chassey L, Meyniel-Schicklin A, Aublin-Gex P, Andre V Lotteau (2012) New horizons for antiviral drug discovery from virus–host protein interaction networks. Current Opinion in Virology2:606-613.
20. Rosenberg SA,Restifo NP (2015) Adoptive cell transfer as personalized immunotherapy for human cancer. Science348:62-68.