GARP Correlates With Tumor-Infiltrating T-Cells and Predicts the Outcome of Gastric Cancer

Accepting the crucial role of the immune microenvironment (TME) in tumor progression enables us to identify immunotherapeutic targets and develop new therapies. Glycoprotein A repetitions predominant (GARP) plays a vital part in maintaining regulatory T cell (Treg)-mediated immune tolerance. The impact of GARP in TME of gastric cancer is still worth exploring. We investigated public genomic datasets from The Cancer Genome Atlas and Gene Expression Omnibus to analyze the possible role of GARP and its relationship with TME of gastric cancer. Fluorescence-based multiplex immunohistochemistry and immunohistochemistry for T-cell immune signatures in a series of tissue microarrays were used to validate the value of GARP in the TME. We initially found that GARP expression was upregulated in gastric carcinoma cells, and diverse levels o3f immune cell infiltration and immune checkpoint expression were detected. Gene expression profiling revealed that GARP expression was related to the TME of gastric cancer. GARP upregulation was usually accompanied by increased FOXP3+ Treg and CD4+ T cell infiltration. In addition, GARP expression had positive relationships with CTLA-4 and PD-L1 expression in gastric cancer. Cox regression analysis and a nomogram highlighted that the probability of poor overall survival was predicted well by GARP or GARP+CD4+ T cell. Taken together, this research underlines the potential effect of GARP in regulating survival and tumor-infiltrating T-cells. In addition, the function of CD4+ T cell immune signatures in the prognosis can be clinically meaningful, thereby providing a new idea for the immunotherapeutic approach.


INTRODUCTION
Gastric cancer (GC) patients which received the conventional treatment at the same stage usually showed heterogeneous clinical prognosis (1,2). Therefore, we need a prognostic signature that is different from the previous staging system to accurately predict the outcome of patients and better guide adjuvant therapy (2)(3)(4). Tumor-associated immune cells in the tumor microenvironment have been demonstrated to play a vital part in tumor development and affect the clinical outcomes of patients (5,6). Although remarkable progress has been made in cancer treatment through the blockade of CTLA-4 or PD-1 signaling using monoclonal antibodies (mAbs), most patients do not respond to immunotherapy because of primary or acquired drug resistance (7,8). Therefore, a better understanding of the markers associated with T cells in the TME is meaningful for deciphering the mechanisms of immunotherapy and identifying new therapeutic targets (9,10).
The transforming growth factor-b (TGF-b) superfamily is an important family of regulatory cytokines with multiple functions in development, immunity, and cancer (11). GARP (commonly known as leucine-rich repeat-containing 32) is a cell surface docking receptor for latent TGF-b, and also has been studied as a non-signal receptor on the surface of Tregs, platelets, and certain cancer cells (12)(13)(14). GARP forms a complex with integrin and releases active TGF-b from the cell surface, thereby enhancing the inhibitory phenotype of Tregs (15)(16)(17). It has been reported that GARP is overexpressed in colon, lung, and breast cancers, and patients with high GARP expression tend to have a poor prognosis (12,18). Therefore, the roles of GARP in the immune microenvironment of gastric cancer and prognosis are worthy of further exploration. In the present study, we combined experiment and bioinformatic technique to further characterize the potential impact of GARP in regulating survival and the TME of gastric cancer, thereby finding TME-associated prognostic signature.

Bioinformatic Analysis
Evidence From the Public Database TCGA clinical and RNA-Seq data for GC patients, including 375 tumor samples, 27 paracancerous samples, and 32 normal samples, were download from Genome Data Commons (https://portal.gdc. cancer.Gov/). We excluded data missing key information, such as overall survival (seven cases), age (four cases), and lymph node metastasis (two cases). Our research meets the publishing requirements provided by TCGA. We also obtained an additional GEO dataset, GSE84437, which contained 434 GC patients with survival information (https://www.ncbi.nlm.nih.gov/geo/).

TISIDB Analysis
TISIDB is a website for comprehensive research on the immune microenvironment that integrates tumor immunology with multiple types of data resources (http://cis.hku.hk/TISIDB/) (19). In TISIDB, we can use literature mining and high-throughput data analysis to clarify the roles of genes of interest in tumor-immune interactions. We analyzed the effect of GARP expression on the prognosis of patients with gastric cancer and its connections with the clinicopathological parameters and immune subtypes of gastric cancer.

TIMER Database Analysis
TIMER (v.2.0.) used a deconvolution statistical method to infer the prevalence of tumor-infiltrating immune cells (TIICs) based on the gene expression profile (https://cistrome.shinyapps.io/ timer/), The database used TCGA data from 10897 samples of 32 cancers to approximate the abundance of TIICs. We performed a gene module to assess the association between GARP expression in gastric cancer and TIICs, including B cells, CD4 + T cells, CD8 + T cells, neutrophils, macrophages, and dendritic cells.

CIBERSORT Estimation
CIBERSORT is an analytical tool for deconvolution of the expression matrix of immune cell subtypes based on the principle of linear support vector regression (https://cibersort. stanford.edu/index.php) (20). We used the CIBERSORT database to explore the infiltration levels of 22 immune cells in gastric cancer. Standard annotation files were utilized to generate gene expression datasets. CIBERSORT approximates the p-value via Monte Carlo sampling and deconvolution to determine the credibility of the results. We grouped the data downloaded from the TCGA database according to the immune subtypes obtained via single-sample Gene Set Enrichment Analysis (ssGSEA) to evaluate the infiltration of immune cells in different immune subtypes.

Gene Set Enrichment Analysis (GSEA) and Unsupervised Clustering
GSEA is a tool for analyzing genome-wide expression profiling data (21). The basic idea is to use a predefined set of genes, sort the genes according to the degree of differential expression between two sample types, and then test whether the predefined set of genes is enriched at the top or bottom of the sorted table. The samples were first grouped according to phenotypes, and then the differential gene sets were selected according to the group. GSEA determined which group the gene sets assembly chose. In this case, the gene sets were associated with the phenotypic grouping. We downloaded RNA-Seq data for gastric cancer from the TCGA database. Then, we performed GSEA using R (v.3.5.3) to identify signaling pathways that were differentially activated in gastric cancer. The threshold was determined using the following parameters: false-discovery rate (FDR) < 0.05 and P < 0.05.
The infiltration levels of the different immune cell populations were determined via single-sample GSEA (ssGSEA) using the R Bioconductor package Gene Set Variation Analysis with the default parameters. The ssGSEA algorithm is a rank-based method that defines a score representing the degree of absolute enrichment of a particular gene set in each sample. GSEA was performed on each sample using transcriptome data and clinical data downloaded from the TCGA database. We obtained the immune cells, immune-related gene sets, and immune-related pathways of each sample, thereby permitting the immune activity of each sample to be evaluated using the CIBERSORT and ESTIMATE algorithms. Unsupervised clustering classifies the samples into distinct subtypes according to the immune cell infiltration pattern of each sample. The unsupervised clustering "Pam" method in accordance with Euclidean and Ward's linkage was used in our analysis, executed by using the "ConsensuClusterPlus" R package, and repeated 1,000 times to ensure the classification stability.

Protein-Protein Interactions
STRING is a database of known and predicted protein-protein interactions (https://string-db.org/) (22). The interactions include direct (physical) and indirect (functional) associations. They stem from computational prediction, knowledge transfer between organisms, and interactions aggregated from other (primary) databases. We used the STRING database to build a protein network of interactions between GARP and related immune signatures.

Human Tissue Samples and Patient Clinical Information
The tissue microarray (TMA) (176 gastric cancer tissues, 52 normal gastric mucosa tissues) used in this research was prepared by the Department of Clinical Biobank of the Affiliated Hospital of Nantong University. A core on the TMA represents a sample with a diameter of 2 millimeters. We averaged the results of multiple samples from the same patient. This research retrospectively analyzed the clinicopathological features and prognoses of the patients. We collected clinicopathological information from the patients' medical records. The patients had not received radiotherapy, chemotherapy, or biological immunotherapy before surgery. This research protocol was approved by the Human Research Ethics Committee of the Affiliated Hospital of Nantong University (Jiangsu, China).

Immunohistochemistry (IHC)
Formalin-fixed, paraffin-embedded TMA sections were deparaffinized and rehydrated using alcohol and xylene. TMA sections were heated using a microwave in sodium citrate buffer (0.01 M, pH 6.0) to repair antigen. The sections were incubated with 5%BSA to quench endogenous peroxidase activity and then with rabbit anti-PD-L1 (13684S, Cell Signaling Technology) and mouse anti-CTLA-4 antibodies (NB10064849, NOVUS). An EliVision Plus DAB Kit (Kit-0015; Maxim Biotechnologies, Fuzhou, China) was used to analyze the result of antibody binding. The results of TMA staining were assessed using the semiquantitative H-score method by a pathologist who was blinded to the clinical information of the patients. The staining intensity score was multiplied by the percentage of positively stained cells to calculate the total score, which ranged from 0 to 300.

Fluorescence-Based Multiplex Immunohistochemistry (mIHC)
TMA sections were heated using a microwave in AR6 buffer (AR600, AKOYA) to repair antigen. MIHC staining was performed after the secondary antibody was added, and then the antigen was repaired via heat induction and cooling. The nucleus was stained with DAPI and sealed. The slides were scanned using the Vectra 3.0 automated quantitative pathology imaging system to detect and measure the positive rate of biomarkers. The cores containing both tumor and stroma were captured with a ×20 Olympus lens objective. Using inForm ® Cell Analysis software, we train machine-learning algorithms to segment the images into tissue areas of cancerous cells and stromal cells, to segment individual cells by DAPI counterstaining, and to accurately identify and quantify the phenotypes of those cells in all high-power fields within the entire tissue section.

Statistical Analysis
Student's ?-test was used to compare GARP protein expression between tumor and non-tumor tissue samples. Pearson's c 2 test was performed to determine the correlation between GARP expression and clinicopathologic parameters. Cox regression models were used to identify prognostic factors. We used the "rms" R package to formulate nomograms, which can predict the probability of 1-year, 3-year, and 5-year overall survival for GC patients. R software (v.3.6.0), SPSS (v.17.0), GraphPad Prism (v.5.0), and Strawberry Perl (v.5.30.1) were used in the early data processing of this study. For all tests, P < 0.05 was considered statistically significant.

Immune Microenvironment Grouping of Patients With Gastric Cancer
TME is mainly composed of tumor-infiltrating immune cells, extracellular matrix, and secreted factors that are highly related to overall survival and treatment response (23). We used the CIBERSORT and ESTIMATE algorithms to evaluate each TCGA sample by scoring immune cells, stromal cells, and tumor purity. Besides, we divided TCGA samples into high and low immunity groups via the unsupervised clustering "Pam" method. Compared with that in the high immunity group, the immune score was significantly lower in the low immunity group, whereas the tumor score was higher in the low immunity group (P < 0.05) ( Figures 1A, B). In addition, GARP expression was significantly increased in the high immunity group (P < 0.001) ( Figure 1C). Then, we investigated the link between GARP expression and the markers of CD4+ T cells, CD8+ T cells, and Tregs in different immune groups. We found that GARP expression was related to CD4, CD8A, and FOXP3 expression, and the correlation was stronger in the high immunity group than in the low immunity group (P < 0.05) ( Figures 1D-F). To determine whether GARP was involved in the activation of Treg in gastric cancer, we subdivided the GARP high/low group into a TGF-b1 low and high group, and then made a Kaplan-Meier curve with overall survival (Supplementary Figure 1). We also divided the GARP high/low group into a FOXP3 low and high group. However, the P-value of Kaplan-Meier curve with overall survival was no statistical significance. In TIMER database, there was a significant link between GARP expression and the levels of CD4+ T cells (Pearson correlation = 0.450, P < 0.05), CD8+ T cells (Pearson correlation = 0.290, P < 0.05), macrophages (Pearson correlation = 0.617, P < 0.05), neutrophils (Pearson correlation = 0.322, P < 0.05), and dendritic cells (Pearson correlation = 0.506, P < 0.05) ( Figure 1G). These findings illustrated the relationship of GARP with the immune microenvironment in gastric cancer.

High Expression of GARP Is Associated With Clinicopathological Features
Although histological classification or clinical staging can well help to predict the prognosis of GC patients, other markers are needed to detect tumor progression. GARP expression was correlated with the tumor grade (P < 0.05) and stage (P < 0.05) (Figures 2A, B). In the TISIDB database, the samples were classified into distinct subtypes according to the median GARP mRNA levels. High GARP expression in patients with gastric cancer was associated with a worse prognosis (P < 0.05) ( Figure 2C). A great amount of transcriptomic data may not be translated into proteins. We next performed fluorescence-based mIHC using TMA and determined that GARP protein levels significantly differed between tumor and normal tissues (P < 0.05) ( Figure 2D). MIHC staining was combined with multispectral image analysis to estimate the positive rate of GARP in a cohort of GC patients. Cytokeratin (CK) was used to identify epithelial cells in tumor samples and to define tumor and stroma, and DAPI was used to stain nuclei. Machinelearning algorithms were trained to distinguish between different tissues (tumor tissue, stroma, and no tissue) and cell phenotypes (tumor cell, and immune cell). Then, we analyzed whether GARP expression levels were associated with clinicopathological features, including gender, age, tumor size (T), lymph node metastasis (N), distant metastasis (M), TNM stage, tumor differentiation, preoperative serum carcinoembryonic antigen (CEA) levels, and preoperative serum carbohydrate antigen 19-9 (CA19-9) levels. 176 GC patients were divided into the GARP-high group (88 cases) and GARP-low group (88 cases) based on the median GARP expression. From our analysis, we observed marked correlations of GARP expression with tumor size (P < 0.05), distant metastasis (P < 0.05), and TNM stage (P < 0.001) ( Table 1). We performed Bonferroni adjustment for multiple comparisons of clinicopathologic characteristics in Supplementary Table 1.

The Relationship Between GARP and TIICs
We performed computational imaging techniques to evaluate multiple lymphocyte markers at the same time, allowing spatial analysis of different T cell populations in the same sample tissue section ( Figures 3E, F). TMA sections were developed to visualize CD3, CD4, CD8, FOXP3, GARP, and CK simultaneously on a cohort of gastric cancer samples. These markers were indicated as signatures for T cells (CD3, CD4, CD8, and FOXP3). Our results demonstrated that nearly all samples had varying degrees of immune cell infiltration. Then, we analyzed whether TIIC counts differed between patients with gastric cancer according to GARP expression. Our samples were divided into two groups based on the median GARP expression level. CD3+ T cell infiltration was significantly suppressed in the high GARP expression group compared with that in the low GARP expression group (P < 0.05) ( Figure 3A). In addition, CD4+ T cell (P < 0.05) and FOXP3+ Treg infiltration (P < 0.05) were obviously enhanced in the high expression group, whereas CD8+ T cell infiltration did not differ significantly between the two groups (P = 0.728) ( Figures 3B-D). In addition, our analysis illustrated that the levels of immune cell infiltration were vastly correlated with tumor size (T) (Figures 3G, H). According to GARP high/low expression, we subdivided our data to reveal the effect of GARP within T groups on levels of immune cell infiltration. CD4+ T cell and FOXP3+ Treg infiltration were slightly higher in GARP high group (Supplementary Figure 2). From our exploration, we can conclude that GARP, as a surface molecule of Tregs, is associated with the infiltration of CD4+ T cell and FOXP3+ Treg but not that of CD8 + T cell.

GARP Upregulation or GARP+CD4+ T Cell Is an Independent Prognostic Factor for Poor Overall Survival
We utilized a cohort of 434 GC patients (GSE84437) to further comprehend the survival mechanism associated with the relationship between GARP and T-cell immune signatures.
Cox regression analysis showed only GARP can be used as an independent factor affecting the prognosis of gastric cancer compared with other immune molecules (P < 0.001) ( Figures  4A, B). In our research cohort, GARP upregulation in gastric cancer was a prognostic factor for poor overall survival (P < 0.001) (Figures 4C, D). We evaluated GARP expression levels in stroma and tumor cells, respectively. Immunofluorescence results showed the positive staining of GARP in CD4+ T cells ( Figures 4F, G). Additionally, high proportions of GARP+CD4+ T cells from all T cells translated to the inferior outcome (P < 0.05) ( Figure 4E). On the contrary, CD3+ T cell, CD4+ T cell, CD8+ T cell, and FOXP3+ Treg were not associated with the survival of gastric cancer.

Construction and Evaluation of a Nomogram for Overall Survival
Cox regression analyses were conducted to exhibited that GARP or GARP+CD4+ T cell could serve as an independent predictor of patients' overall survival after adjusted by TME-associated signatures in multiple GC cohorts (Figure 4). Based on logistic regression, we generated a nomogram that integrated GARP, GARP+CD4+ T cell, and other clinicopathological features, including tumor size, lymph node metastasis, distant metastasis, TNM stage, tumor differentiation to predict the probability of 1-year, 3-year, and 5-year overall survival for GC patients with the GSE84437 and the experimental cohort ( Figures 5A, C). The calibration plots revealed the probability of a 3-year survival rate is well predicted in the GSE84437 cohort and the experimental cohort ( Figures 5B, D).

Exploration of the Molecular Mechanism of GARP
We investigated whether the prognostic effect of GARP is related to immune checkpoints in gastric cancer. IHC was performed to explore the expression of PD-L1 and CTLA-4 in gastric cancer and then analyzed their correlations with GARP expression. As presented in Figures 6A-F, GARP expression was associated with CTLA-4 (P < 0.05) and PD-L1 expression (P < 0.05). Meanwhile, a positive relationship between CTLA-4 and PD-L1 expression was noted in gastric cancer (P < 0.05).
We divided 176 GC patients into the GARP-high group and GARP-low group based on the median GARP expression. GSEA analysis screened the differential genes according to the sample groups, and then enriched the genes. The results showed that gene sets were enriched with GARP upregulation in GC samples. We selected 10 KEGG pathways with significant differences according to the normalized enrichment score (FDR < 0.25, NOM P < 0.05). Specifically, the following pathways were significantly enriched in the high expression phenotype: ECMreceptor interaction, GAP junction, leukocyte transendothelial migration, the JAK-STAT signaling pathway, the MAPK signaling pathway, pathways in cancer, the TGF-b signaling pathway, the Toll-like receptor signaling pathway, the intestinal immune network for IgA production, and natural killer cell-mediated cytotoxicity ( Figure 6K).

DISCUSSION
The genome resources of the public database provide a unique platform for us to further explore the molecular characteristics of different cancers (24). Our research explored the relationship between GARP and the immune microenvironment of gastric cancer based on TCGA and GEO data. The success of cancer immunotherapy has revealed that immune cells, especially T cells, can be helpful in eliminating tumor cells (26). Wang Yu Cai et al. revealed that evaluating the TME components can predict survival time and provide a new idea for the immunotherapeutic approach of gastric cancer (GC) (27). Compared with a low density of T cells, a higher density of T cells in the TME can better predict the prognosis of gastric cancer (28,29). Salem et al. demonstrated that GARP restrains antitumor immunity by adjusting the function of Tregs in colorectal cancer (30). our analysis illustrated that the significant connection between GARP and TME was still worth exploring. T helper cells, cytotoxic T cells, and Tregs are associated with T cell-mediated immune responses in the TME (31). Bioinformatic analysis illustrated that GARP expression was relevant to the immune groups of gastric cancer. In the colon cancer model, the loss of GARP in Treg leads to spontaneous inflammation and enteritis with high activation of CD4+ and CD8+ T cell, which has an important impact on immune surveillance (30). Our research illustrated GARP expression correlated with FOXP3+ Treg and CD4+ T cell infiltration and CTLA-4, and PD-L1 expression, whereas GARP expression had no remarkable connection with CD8+ T cell infiltration. Lucas et al. reported that the enhancement of this combination therapy did not depend on increasing the number of CD8+ T cells (7). Our study undoubtedly provided evidence for this result and also raised a question, specifically whether the antineoplastic effect of GARP affects the infiltration of CD4 + T cells. Despite the high expression of T cell markers indicate the improvement of progress in many cancers, no statistically significant correlation was observed in gastric cancer (24). A recent report has shown that T cells expressing immune checkpoints such as LAG3, PD1 may represent exhausted T-cell (32). In our analysis, the levels of CD4+ T cell and FOXP3+ Treg infiltration were significantly higher in patients with T3 gastric cancer than patients with T1 and T2 gastric cancer. We further found the expression of GARP in CD4+ T cells and analyzed that GARP+CD4+ T cells play a significant role in the prognosis of gastric cancer. The nomogram also showed that GARP+CD4+ T cell can predict the survival rate of GC patients together with other clinicopathological parameters. CD4 positive cells are likely to annotate as Tregs, dendritic cells, macrophages, or NK cells (27,33). GARP may identify a special subgroup of T cells related to inferior prognosis. The variations in the phenotype of tumor-infiltrating immune cells and their relationship with prognosis highlight the clinical significance of the crosstalk between tumor cells and TME (34). However, more efforts are needed to determine whether T cell dysfunction is associated with GARP expression and poor prognosis of GC patients.
By analyzing the STRING database and our study, we also found an interaction between GARP, PD-L1 and, CTLA-4 expression. Researches have also found a significant correlation between PD-1+PD-L1+ T cells and Tregs. In animal experiments, combined therapy targeting GARP and PD-1 has achieved positive results (7,35). A recent report suggests that combining checkpoint inhibitors with chimeric antigen receptor T-cells may also be of great significance in the treatment of cancer (27). GARP upregulation was related to the TGF-b signaling pathway. Evidence of TGF-b signaling in T cells has been found in melanoma specimens infiltrated by GARP-expressing T cells, indicating that inhibiting the activity of Treg-derived TGF-b1 using anti-GARP: TGF-b1 mAbs may effectively enhance CD8+ T cell-mediated antitumor immunity (7). Our study illustrated that the TGF-b signaling pathway was differentially enriched in the GARP high expression phenotype. Whether GARP affects the immune microenvironment of gastric cancer through the TGF-b signaling pathway is worthy of further exploration.
This study had several limitations. For example, the number of samples was limited, thus limiting the strength of our conclusions. In addition, we have not yet identified the regulatory pathway connecting GARP expression and CD4+ T cell, nor have we confirmed whether there are interactions among GARP, CTLA-4, and PD-L1. The combination of anti-GARP: TGF-b mAbs and PD-1 inhibitors can significantly enhance the effector ability of tumor T cells (36)(37)(38). Although remarkable progress has been made in cancer treatment by blocking CTLA-4 or PD-1 pathway with mAbs, most patients do not respond to therapy because of T cell-mediated primitive or acquired immune resistance to anti-tumor drugs (8,39,40). Combining checkpoint inhibitors with chimeric antigen receptor T-cell may also be of great significance in the treatment of cancer.
Taken together, our study proved that GARP is an independent influencing factor that is significantly upregulated in gastric cancer. We underlined the relationship between GARP and tumor-infiltrating T-cell. The phenotype of tumorinfiltrating immune cells was clinically meaningful, thereby providing a new idea for the immunotherapeutic approach of gastric cancer (30,41).

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.