Immunoglobulin superfamily genes are novel prognostic biomarkers for breast cancer

Breast cancer progression is associated with dysregulated expression of the immunoglobulin superfamily (IgSF) genes that are involved in cell-cell recognition, binding and adhesion. Despite widespread evidence that many IgSF genes could serve as effective biomarkers, this potential has not been realized because the studies have focused mostly on individual genes and not the entire network. To gain a global perspective of the IgSF-related biomarkers, we constructed an IgSF-directed neighbor network (IDNN) and an IgSF-directed driver network (IDDN) by integrating multiple levels of data, including IgSF genes, breast cancer driver genes, protein-protein interaction (PPI) networks and gene expression profiling data. Our study shows that IgSF genes in the PPI network have important topological features related to cancer. Most IgSF genes are either cancer driver genes themselves or associated with them. We also identified a 21-gene IgSF network module with enriched mutations that are associated with overall survival based on 450 breast cancer patient samples extracted from The Cancer Genome Atlas (TCGA) and multiple independent microarray validation datasets. These results highlight the potential of IgSF genes as novel diagnostic, prognostic and therapeutic targets for breast cancer.

INTRODUCTION
Breast cancer is the leading cause of cancer death among women worldwide. In Chinese women, breast cancer is the most prevalent form of cancer with more than 1.6 million people diagnosed and 1.2 million people dying every year. The most common type of breast cancer is invasive ductal carcinoma (IDC) that can spread from the ducts or the lobules to the surrounding tissue. Prognostic biomarkers are useful to choose the appropriate treatment for IDC, and they significantly affect the process of cancer therapy [1][2][3]. Studies have shown that the genetic diversity in breast cancer impacts response to treatment and patient outcomes. This is exemplified by the estrogen receptor negative (ER−) and positive (ER+) subtypes that have different prognostic gene signatures and responses to treatment [4]. Therefore, there is scope to identify novel signatures that can enhance predicting the prognostic and clinical behavior.
Gene expression of many IgSF members is altered in breast cancer, and hence, they are promising candidates as prognostic biomarkers. ALCAM (CD166) is a potential breast cancer biomarker and a therapeutic target due to its role in induction of programmed cell death, apoptosis and autophagy in breast cancer [5]. Down regulation of CXCR4 inhibits cell migration in breast cancer cells [6]. The expression of MUC18 (CD146) promotes the progression of human breast cancer cells by increasing their motility, invasiveness and tumorigenesis [7]. L1CAM is potentially an early diagnostic biomarker in breast cancer progression as it promotes cell adhesion and migration in vitro [8]. Although these findings demonstrate the important role of IgSF members in breast cancer progression and metastasis, these studies focus Research Paper on one or a few IgSF members analyzed in either cell lines or in limited patient samples and therefore do not present a global perspective of the entire immunoglobulin superfamily.
Current advances in cancer biology and genomic methods have generated large-scale gene expression profiling datasets (such as TCGA) and other OMICs and provide an opportunity to study the entire network of IgSF genes as well as identify novel biomarkers for breast cancer. Previously, Li and others used cancer gene microarray and network data to develop a network-based method for cancer prognostic biomarker identification [9]. Similarly, Chuang and others used a protein-networkbased approach using data from protein interaction databases to identify markers as sub networks [10]. These studies demonstrate that integrating gene expression and protein-protein interaction data can improve prediction performance in network biomarker identification.
Our aim was to study the role of the IgSF network in breast cancer and test its diagnostic, prognostic and therapeutic potential. Towards this goal we constructed an IgSF-directed neighbor network (IDNN) and an IgSFdirected driver network (IDDN) to address the role of the IgSF network in breast cancer. We identified IgSF genes in the PPI network with hub topological features connected to breast cancer. We also identified a 21-gene module from the IDDN network that was associated with the overall survival of breast cancer patients. This module included several key IgSF and breast cancer driver genes with enriched mutations, demonstrating the functional significance of the IgSF genes to predict breast cancer. Our findings highlight the novel role of the IgSF-directed network in breast cancer. It also highlights their potential for biomarker-guided development of preclinical and clinical therapeutic modalities.

IgSF genes play a crucial role in breast cancer
We constructed an IgSF-directed neighbor network (IDNN) that had 1050 nodes including 283 IgSF genes and their 767 neighbors from the PPI network ( Figure 1A, Supplementary Table S1). We found 33 IgSF genes (KIT, PDGFRB, KDR, FGFR1, CD28, NTRK1, etc.) that were themselves breast cancer driver genes as well as 250 IgSF genes (including PTPN11, LCK, GRB2, ABL1, APP, etc. Supplementary Table S2) that associated with the breast cancer driver genes. The top 10 IgSF genes (KIT, PDGFRB, KDR, FGFR3, FGFR1, CD4, NTRK1, IGSF21, HSPG2, and PDGFRA) had a direct connection in the network, indicating that IgSF genes play hub roles in the sub-networks ( Figure 1B). We also found a greater degree of genes (mean = 16.76) that were both IgSF and breast cancer driver genes in the IDNN suggesting a intricate link between the IgSF genes and breast cancer ( Figure 1C, Supplementary Figure S1). Also, enrichment of IgSF cancer driver genes in specific modules suggests that the IgSF driver genes may play consequential roles in special biological modules (Supplementary Figure S2,  Supplementary Table S3).
Notably, six of the top 10 IgSF genes in the PPI network were well-known breast cancer driver genes, including KIT, whose high expression occurs infrequently in breast cancer [11]. In the triple negative breast cancer, immunohistochemical expression of C-kit and mutations of PDGFRA are frequent suggesting that they are good candidates for molecular targeted therapy [12]. Our findings also suggest that CD4 + follicular helper T cells (Tfh) may be prognostic indicators as they are found in the breast tumors [13]. Furthermore, down-regulation of KDR expression induces apoptosis in breast cancer cells [14]. Also, FGFR1 activity is required for the survival of an FGFR1-amplified breast cancer cell line [15]. Among the IgSF neighbor genes, ABL1 and PDGFR are well known breast cancer driver genes that promote acquired resistance to aromatase inhibitor (AI) therapy in ER+ breast cancers [16]. The activation of the PDGFR and ABL1 pathways is associated with long-term estrogen deprivation in MCF7 breast cancer cells and decreased anti-proliferative response to AI treatment in primary ER-positive breast carcinomas.

IgSF genes directly interact with breast cancer driver genes
To further explore the relationship between IgSF and breast cancer driver genes, we constructed a network called IgSF-directed driver network (IDDN). This included IgSF genes and breast cancer driver genes that were extracted from the IDNN (Figure 2A, Supplementary  Table S4). The IDDN contained 253 genes, of which 103 were IgSF and 121 were breast cancer driver genes. Among these, 29 of the IgSF genes were also breast cancer genes. Compared to the nodes in IDNN, the nodes in the IDDN had higher degrees, betweenness centrality and closeness centrality (avg.  Figure S3C; Wilcoxon rank sum test). This suggested that the IDDN obtained from IDNN was closer in structure and played a crucial role in the biogenesis of breast cancer.
To explore the biological functions of these genes, we conducted pathway analysis using the Subpathway Miner [17]. The T cell receptor signaling pathway emerged as the most significant pathway from this analysis with several of the IgSF genes located in key positions, such as CBL, LCK, ZAP70, FYN, CD28, PD-1, CTLA4, ICOS, and CD4/8 ( Figure 2B). Since the expression of CBL can inhibit LCK and ZAP70 gene expression, CBL gene may play a crucial role in the pathway [18]. www.impactjournals.com/oncotarget We then performed a GO enrichment analysis using IgSF and breast cancer driver genes. Some common GO terms that we found associated with breast cancer included cell surface receptor-linked signal transduction, immune response, regulation of cell proliferation, regulation of cell activation, regulation of T cell activation, regulation of T cell differentiation, and negative regulation of the immune system process. Significance of this analysis is exemplified by the fact that mutations of BTLA, CD28, CD4, and CD8A genes in the GO term cell surface receptor-linked signal transduction contribute to sporadic breast cancer risk. More importantly, we found overlapping GO terms between the IgSF genes and the breast cancer driver genes, including phosphate metabolic process, protein amino acid phosphorylation, enzyme linked receptor protein signaling pathway, transmembrane receptor protein tyrosine kinase signaling pathway, cell surface receptor-linked signal transduction, and regulation of cell proliferation ( Figure 2C). This indicated a close connection between IgSF genes and breast cancer. Additionally, the degree of the common genes of IgSF and the breast cancer driver genes derived from the IDDN was also greatest ( Figure 2D).

IgSF genes associate with driver genes in functional modules
To understand the communication between the IgSF and the breast cancer driver genes, we performed a module analysis in the IDDN. We found six significant IgSF-associated modules in which IgSF and breast cancer driver genes were closely associated ( Figure 3). We conducted a GO terms enrichment analysis of the module genes using DAVID [19] and identified the biological process terms that enriched the module genes using a cutoff of FDR < 0.05 (Supplementary Table S5). The top-ranked GO terms included immune response, cell adhesion, biological adhesion, regulation of alpha-beta T cell proliferation, transmembrane receptor protein tyrosine kinase signaling pathway, and T cell activation. The GO term cell adhesion enriched by the second and sixth module genes is a very important biological process of the IgSF genes (VCAM1, ICAM1, F11R, ITGAL,  EZR, ICAM4, ICAM5, ICAM2, ICAM3, ITGB5, ITGB2, and CD226) that participate in the immune response. Among these, the expression of members of the ICAM family regulates tumorigenesis and are potential diagnostic biomarkers and therapeutic targets for breast cancer [20]. Moreover, both the first and fourth modules enriched the common GO term phosphate metabolic process, which is differentially regulated in breast cancer [21].

IgSF-related modules are enriched with cancer mutations
Since genetic mutations can cause cancer, we analyzed if IgSF genes in our modules are mutated. Previously, the IgSF genes in the module that includes well-known breast cancer genes have been shown to be mutated. Some SNPs in six immunological genes, BTLA, ITGAL, CTLA4, ICOS, PDCD1, and VTCN1 were reported as breast cancer risk mutations in previous studies [22][23][24][25][26][27] (Supplementary Table S6). We performed gene mutation enrichment analysis and evaluated if these module genes were enriched in the top 15% of mutated genes in breast cancer. We identified three module genes that overlapped with the top 15% mutated genes in breast cancer (p = 0.001 for module 1, Supplementary Figure S4A; p = 0.001 for module 4, Supplementary Figure S4B; p = 1.28e-07 for module 5, Supplementary Figure S4C; Hypergeometric Test).
Next, we identified the somatic mutations of the genes in the IDDN using TCGA breast cancer somatic mutation data and determined if they overlapped with the top 15% mutated genes in the IDDN using hypergeometric test. We found that the second module with 21 genes was enriched in mutated genes ( Figure 4A, P < 0.001). Four of the six common genes (PTPN11, TRIM2, FGR and ZAP70) were IgSF genes ( Figure 4B). PTPN11, a HER2inhibition up-regulated PTP (protein tyrosine and dualspecificity phosphatase) transduces positive signals and is an oncogene [28]. TRIM2 is a diagnostically significant and conserved element of the SOX10 signature in BBC (breast basal-like carcinomas) cell lines [29]. FGR and ZAP70 have diagnostic and therapeutic potential due to their relationship with breast cancer development and progression [30,31].  these 20 genes to generate OS (overall survival) prediction (see the Material and Methods section). The Cox regression coefficients of the total samples, ER+ samples and ER-samples are as listed in Supplementary Table S7. Using the median risk score of the test series as the cutoff point, we calculated the risk scores for the 20 genes for each patient and then ranked the patients according to their risk score. The patients grouped into a high-risk (N = 225) or a low-risk (N = 225) category by using the median risk score of the test series as the cutoff point. Patients in the high-risk group had significantly shorter median OS than those in the low-risk group ( Figure 5A, HR = 3.82, P = 2.81e-05). In addition, patients with high risk in the ER+ series (HR = 2.78, P = 0.00507) as well as the ERseries (HR = 8.77, P = 0.0124) had significantly shorter median OS than those in the low-risk group ( Figure 5A). Based on the distribution of gene risk score, the survival status and the gene expression signature of the breast cancer patients, patients with high risk scores expressed higher levels in the ten risk genes, whereas patients with low risk scores expressed higher levels in the remaining ten protective genes ( Figure 5B-5D).

Validating the prognostic potential of the IgSF module in independent patients
To confirm the prognostic value of the IgSF module, we validated the 21 gene signatures in four independent microarray datasets (Table 1). Using the same risk score formula, we classified patients in GSE4922 into high-risk and low-risk groups using the median score of the test series as the cutoff point. Consistent with our previous findings, patients in the high-risk group had significantly shorter median OS than those in the low-risk group ( Figure 6A, all samples P = 0.076; Figure 6B, ER+ samples P = 0.012; Figure 6C, ER-samples P = 0.021). Similarly, we classified the patients in GSE7390 into a high-and a low-risk group ( Figure 6D, all patients, P = 0.091; Figure 6E, ER+ patients, P = 0.018 and Figure  6F, ER-patients, P = 0.054) and obtained similar results.
After further adjusting for other markers, the univariate analysis indicated that the IgSF module, as an independent risk factor, was significantly associated with the overall survival of breast cancer patients from TCGA (HR = 2.71, 95% CI: 2.02-3.64, P = 2.75e-11), GSE4922 (HR = 1.69, 95% CI: 0.94-3.05, P = 0.079), and GSE7390 (HR = 1.58, 95% CI: 0.92-2.73, P = 0.094) datasets (Table 2). Further, when multivariate analysis was performed to investigate the independence of the module to other clinical factors, the high-and low-risk groups remained independent of other clinical factors in TCGA patients (HR = 2.60, 95% CI: 1.84-3.67, P = 4.32e-8), GSE4922 (HR = 1.86, 95% CI: 1.01-3.43, P = 0.045), and GSE7390 (HR = 1.70, 95% CI: 0.97-2.96, p = 0.060) datasets (Table 2). Data stratification analysis on TCGA patients also indicated that the module was independent of PR, age and tumor stage. It performed similarly in PR category (log-rank test p = 0.005 for PR+ group and log-rank test p = 0.002 for PR-group) as well as the age category older or younger than 60 (log-rank test p = 0.01 for the older patients and log-rank test P = 9.64e-4 for the younger patients) and tumor stage III/IV category (logrank test P = 0.012) ( Figure 7A-7E). Similar trend was observed in GSE4922 patients with G3 grade (log-rank   test P = 6.48e-07; Figure 7F). Based on all these data, we conclude that the IgSF-related module is a strong prognostic indicator for breast cancer.

DISCUSSION
Elucidating the functional significance and molecular mechanism of IgSF members in breast cancer may provide new opportunities for the early detection and treatment of cancer. A variety of molecular biomarkers in breast cancer have been previously identified. Van de Vijver and others categorized breast cancer based on the gene expression profile of 70 genes and found that their classification was better than standard systems that were based on clinical and histological criteria [32] . Pawitan and others developed a 64-gene signature to predict the response of breast cancer patients to therapy [33]. Venet and others showed that many random gene expression signatures are significantly associated with breast cancer survival, although the underlying mechanism of the biomarkers was unclear implying the need for more effective biomarkers [34]. We performed a global analysis of the IgSF genes in breast cancer and developed a network-based strategy to identify the IgSF-related PPI network and modules. Our study shows that the IgSF network genes closely interact with the cancer driver genes in the mutations enriched module and can serve as a survival-associated biomarker for breast cancer. Also, the IgSF genes participate in the dysregulation of breast cancer driver genes. These results provide a novel understanding of the role of IgSF network in breast cancer progression. Our study also highlighted the importance of immune-related IgSF genes. Immune genes play a key role in cell-cell communication, and dysfunctional immune response cause various diseases in humans including cancers [9]. In breast cancer, the immune genes may modulate the communication between malignant cells and normal cells. We had postulated that the IgSF-related modules mediate breast cancer progression by regulating cancer metastasis. However, functional enrichment analysis showed that the IgSF-related modules were involved in a variety of biological events, including immune response, cell adhesion, biological adhesion and regulation of alphabeta T cell proliferation. This suggests that in the early stage of breast cancer, IgSF-related modules restrict cell proliferation by regulating signal transduction and the immune response.
In summary, our data shows that the IgSF-related module that is of prognostic importance. Our data also reveal novel functional insights for dysregulated IgSFrelated members and modules in the etiology of breast cancer.

Breast cancer driver gene dataset
We downloaded cancer associated genes from the Catalogue of Somatic Mutations in Cancer (COSMIC v70; Aug 2014), which is the world's largest and most comprehensive resource for exploring the impact of somatic mutations in human cancer. Firstly, we obtained the breast cancer associated genes from the Cancer Gene Census by searching the keyword "breast" in the tumor types. In addition, we searched cosmic mutation data with the keyword "breast" in the primary site, "y" if the entire genome/exome is sequenced and "CANCER" if the mutations affected the tumor generation in pathology. In summary, we obtained 1,307 unique genes acting as the breast cancer driver genes for subsequent analyses.

Human protein-protein interaction data
The protein-protein interaction (PPI) data was downloaded from the Human Protein Reference Database (HPRD Release9, http://www.hprd.org/) [35]. It contained more than 42,000 manually curated interactions between 9,826 human genes.

Analysis of breast cancer gene expression and somatic mutations in clinical datasets
In this study, we focused only on invasive ductal carcinoma (IDC) patients with breast cancer. The breast cancer gene expression, somatic mutation data and the corresponding clinical data was downloaded from The Cancer Genome Atlas (TCGA, http://cancergenome. nih.gov/). The gene expression profiling from 591 breast cancer patients was generated by the UNC__ AgilentG4502A_07_3 microarray platform. The samples that contained only expression profiling without clinical data were removed from the analysis. We generated 450 IDC samples that had both expression profiling and clinical data, involving 17,814 genes. Also, the patients were further stratified into the estrogen receptor negative (ER−) and positive (ER+) groups.

Validation of microarray datasets
The four independent microarray breast cancer datasets [36][37][38][39] that were used in this study were obtained from the Gene Expression Omnibus (GEO) database. All the datasets were produced by the HGU133A or HGU133-PLUS2 platform. The datasets were chosen based on the criteria of no less than 100 samples and the availability of clinical outcome data. Raw microarray datasets were normalized using Robust Multichip Average [40]. The three main steps were as follows: background correction, quantile-normalization and log2-transformation. All the probes were mapped based on their EntrezGeneID. When multiple probes were mapped to the same gene ID, the mean value was used to represent the expression value of the single gene. To account for differences in systematic measurement between different datasets, each dataset was standardized independently by transforming the expression of each gene into a mean of 0 and a standard deviation of 1.

Network construction and analysis of topological features
We constructed a human PPI network based on the HPRD data. Then, we distinguished all the IgSF genes in this network using HGNC IgSF family data. Finally, using IgSF genes and their direct interacting genes in the network (referred to as IgSF neighbor genes), we constructed a sub-network of the human PPI network named the IgSF-directed neighbor network (IDNN). Furthermore, we extracted all the IgSF genes and the breast cancer driver genes to construct an IgSF-directed driver network (IDDN) for subsequent analysis. The Cytoscape software was used for the construction of networks. The topological properties of the IgSF genes were analyzed in both IDNN and IDDN.

Identifying functional modules in the IgSF networks
Based on the IDDN constructed by the IgSF and breast cancer driver genes, we identified all the network modules using GraphWeb [41], a web server for identifying the network-based biomarkers that best represent the property of the network. The GraphWeb web server was made of three component processes: (i) Network datasets (to input human protein-protein interaction pairs; IDDN in this study), (ii) Network algorithm (we used the betweenness centrality clustering method and the default values were set) and (iii) Network settings (including default edge settings, node setting and module settings with less than 3 nodes and insignificant modules hidden).

Statistical analysis
Hypergeometric test was used to explore the overlap between the genes in the IgSF-related modules and the top mutated genes in the sub-network for breast cancer. We also studied the enrichment of breast cancer driver genes in IgSF by using this method. Univariate and multivariate analyses were performed using Cox proportional hazards regression model to determine whether the IgSF-related prognostic module was independent of other clinical variables, and adjusted for ER, PR, HER2, age, stage and grade. Hazard ratio (HR) and 95% confidence intervals (CI) were estimated by Cox proportional hazards regression model. To verify if the modules we identified are associated with patient survival, we determined the regression coefficient of every gene in the module related to patient survival using the gene expression data. The classifier was built as a linear combination of the gene expression values of select immune-related genes with the standardized Cox regression coefficient as the weight. A risk score formula for each patient was established by including the expression values of each selected gene, weighed by their estimated regression coefficients in the multivariate Cox regression analysis [42]. Finally, the patients were divided into high-and low-risk groups using the median of the risk score as the threshold. The patients with high-risk scores were classified as poor outcomes. Kaplan-Meier survival plots and log-rank tests by R package "survival" were used to assess the differences in overall survival (OS) time between the high-and low-risk patients. Bioinformatic analysis was performed with R 3.0.0 statistical software.

Functional enrichment analysis
Functional enrichment analysis at the GO and KEGG levels were performed using DAVID bioinformatics resources (http://david.abcc.ncifcrf.gov/, version 6.7) [43]. The DAVID enrichment analysis was limited to KEGG pathways and GO-FAT biological process (BP) terms with the whole human genome as background. Functional categories with a p-value of < 0.05 were considered statistically significant.

CONFLICTS OF INTEREST
No potential conflicts of interest were disclosed.


s2fieldsofstudy:
Biology
Agricultural And Food Sciences

text:
Partial Sequencing of 16S RRNA Gene of Selected Staphylococcus Aureus Isolates and Its Antibiotic Resistance

The choice of primer used in 16S rRNA sequencing for identification of Staphylococcus species found in food is important. This study aimed to characterize Staphylococcus aureus isolates by partial sequencing based on 16S rRNA gene employing primers 16sF, 63F or 1387R. The isolates were isolated from milk, egg dishes and chicken dishes and selected based on the presence of sea gene that responsible for formation of enterotoxin-A. Antibiotic susceptibility of the isolates towards six antibiotics was also tested. The use of 16sF resulted generally in higher identity percentage and query coverage compared to the sequencing by 63F or 1387R. BLAST results of all isolates, sequenced by 16sF, showed 99% homology to complete genome of four S. aureus strains, with different characteristics on enterotoxin production and antibiotic resistance. Considering that all isolates were carrying sea gene, indicated by the occurence of 120 bp amplicon after PCR amplification using primer SEA1/SEA2, the isolates were most in agreeing to S. aureus subsp. aureus ST288. This study indicated that 4 out of 8 selected isolates were resistant towards streptomycin. The 16S rRNA gene sequencing using 16sF is useful for identification of S. aureus. However, additional analysis such as PCR employing specific gene target, should give a valuable supplementary information, when specific characteristic is expected.

INTRODUCTION
Staphylococcal food poisoning is occurred when people consume food contaminated by Staphylococcus aureus that produce enterotoxins. Various staphylococcal enterotoxins (SE) incriminated in staphylococcal food poisoning have been reported. They were SEA, SEB, SEC, SED and SEE, and recently new serological types of SEs (SEG, SEH, SEI, SEJ, SEK, SEL, SEM, SEN, SEO, SEP, SEQ, SER and SEU) were also identified (Argudin et al., 2010;Xie et al. 2011, Roussel et al., 2015. Staphylococcal enterotoxins-A (SEA), which is found on most food poisoning by S. aureus, is expressed in the mid-exponential phase, and its gene appears to be transferred by temperate bacteriophage (Argudin et al., 2010). The prevalence and genetic diversity of S. aureus has been investigated in raw and pasteurized milk (Rall et al., 2008), ready-to-eat foods (Huong et al., 2010), milk and food product (Salasia et al., 2011), and street-vend foods (Rohinishree & Negi, 2011). For epidemiological purposes, the accuracy of identification of Staphylococcus species isolated from various food products is critical. In this regard, molecular characterization is reported to be more accurate (Becker et al., 2004) than phenotypic identification (Rohinishree & Negi, 2011). The 16S rRNA gene is extensively used as taxonomic marker molecules during molecular characterization (Janda & Abott, 2007), particularly during the sequence analysis for differentiating species and sub species of bacteria (Rohinishree & Negi, 2011).
Generally, during sequencing analysis, the DNA target is firstly amplified using a pair primers followed by the sequencing using a single primer (SenGupta & Cookson, 2010). Identification based on highly conserved genes such as 16S rRNA usually uses long sequences primers (≥ 500 bp to about 1500 bp) (Janda & Abbott, 2007), although species specific shorter sequences can also be applied. For this purposes, universal primers for amplification of 16S rRNA genes are widely available, such as primers 63F and 1387R (1350bp) (Marchesi et al., 1998) and 27f and 1492r (Frank et al., 2008). Relatively shorter primer, i.e. 16sF and 16sR3 were also used (Lee et al., 2007). Since the choice of primers will affect the diversity of bacterial species that will be detected (Fredriksson et al., 2013), in this study three different primers targeting the 16S rRNA gene, i.e.
16sF, 63F and 1387R, was used separately to sequence DNA of selected S. aureus isolates from milk, egg dishes and chicken dishes. Milk, egg dishes and chicken dishes were reported contaminated by S. aureus in previous study (Handayani et al., 2014).
Furthermore, as information on the spreading of antibiotic resistance strains among S. aureus in food is also important, the antibiotic susceptibility testing of the isolates against antibiotics was also conducted in this study. The antibiotics tested were antibiotics that usually used to control human infections as well as to control and treat infections on farms. Many S. aureus strains that demonstrated resistance to different antibiotics were isolated from hospitals (Schmitz et al., 1999;Brown & Ngeno, 2007;Xie et al., 2011), hospital waste waters (Thompson et al., 2012), as well as from animal based food product, such as raw milk and dairy products (Jamali et al., 2015), poultry retail meat (Teramoto et al., 2016), and goat milk powder (Xing et al., 2016).

Bacterial Strain and Isolates
The wild-type of S. aureus from raw milk (S1, S4 and S10), egg dishes (TB1), sautéed chicken cuts (UA2 and UA13) and chicken cuts satay (SJ1 and SJ4) were isolated in previous study (Handayani et al., 2014). S. aureus ATCC 25923 was used as a reference bacterium. All bacteria were grown in a tryptic soy broth or agar (Difco Laboratory, Detroit, MI, U.S.A.), and incubated at 37 o C for 18h-24h. For confirmation, the isolates were spread onto Baird-Parker Agar (Oxoid Ltd., Hampshire, UK) supplemented with egg-yolk tellurite. Plates were incubated at 37°C for 18-24 h, thereafter the colonies were picked and streaked on Mannitol Salt Agar (Oxoid Ltd., Hampshire, UK). Typical colonies were then tested for production of catalase using Staphylase test kit (Oxoid Ltd., Hampshire, UK) and biochemical identification using API Staph (bioMérieux Inc., North Carolina, USA) according to manufacturer's instructions.

DNA Extraction
Genomic DNA was isolated as described previously by Mason et al. (2001) with slight modification, as reported by Handayani et al. (2014), i.e. lysostaphin (10 mg/mL) was substituted by lysozyme (Bio Basic Canada Inc., Ontario, Canada) solution (10 mg/ml). The concentration of genomic DNA was determined by the Spectrophotometer UV -1800 (Shimadzu, Japan) at 260 nm while the quality was assessed based on the ration of the reading at 260/280 nm. The integrity of the DNA was checked by running in 1.5% agarose gel at 75V for 40 min electrophoresis (Bio-Rad Laboratories Pte. Ltd, Singapore).

Detection of 16S rRNA and Sea Gene by PCR
The amplification of the gene encoding 16S rRNA and sea was performed using primers listed in Table 1 at Thermal Cycler 2720 (Applied Biosystems, California, USA). PCR master mix consisted of 12.5 µL of DreamTaq Green master mix (Thermo Fisher Scientific, Massachusetts, USA), 1 µL of each primer (10 µM), 2 µL of DNA template, and 8.5 µL nuclease free water (Thermo Fisher Scientific, Massachusetts, USA).

Sequence Analysis of 16S rRNA Gene
The genotypic characterization of bacteria isolates was made through partial sequence analysis of 16S rRNA gene, using single primer 16sF, 63F, or 1387R. The PCR products that were amplified with primers 16sF/16sR3 were sequenced using primer 16sF. The PCR products amplified with primers 63F/1387R were sequenced using primer 63F and 1387R separately. The process of DNA sequencing was conducted by BigDye Applied Biosystem sequencer engine model 3730 at Macrogen inc., Singapore. Partial sequence data obtained, in FASTA format, was then submitted to the BLAST process at NCBI (National Center for Biotechnology Information) database for the identification of isolates (http://blast.ncbi.nlm.nih.gov/ Blast.cgi).
The BLAST process was conducted using nucleotide collection searching setting (nr/nt) with Staphylococcus aureus subsp. aureus (taxid:46170) as organism of choice. The sequencing results were then compared to the reported 16S rRNA gene sequences of Staphylococcus species available in the GenBank database. The isolates were identified based on the highest homology percentages to the reported sequences.

Genotypic Characteristic of Isolates
Ten isolates of 78 presumptive S. aureus isolates were positive for S. aureus (Handayani et al., 2014). Eight of these isolates were reconfirmed in this study carrying sea gene that responsible for formation of staphylococcal enterotoxin-A (SEA). The presence of sea gene was indicated by the occurrence of 120 bp amplicon, after PCR amplification using primer SEA1/SEA2 ( Table 2). The reference strain, S. aureus ATCC 25923, did not show this gene. In addition to the work of Handayani et al. (2014), all isolates also demonstrated 1350 bp amplicon after amplification by 63F/1387R primer (universal primer). The amplified PCR products of some isolates are shown in Figure 1.
The BLAST results are listed in Table 3. Sequencing by primer 16sF resulted in high identity percentages (almost all achieved 99%) towards the existing genome of S. aureus strains found in database of NCBI GenBank, in comparison to the sequencing by 63F and 1387R. Sequencing by 63F showed similarity to S. aureus strains in a range between 94% and 99%, while by 1387R in a range of 90% to 99%. In addition, using 16sF also resulted in low E-value. The lower the E-value, or the closer it  is to zero, the match is more significant (Pearson, et al., 2014). The expect value (E) is a parameter that describes the number of hits one can "expect" to see by chance when searching a database of a particular size. As shown in  Table 4. Since all isolates were detected carrying sea gene (Table 2), the isolates were in agreeing to S. aureus subsp. aureus ST288 that able to form enterotoxin-A. The present study has been indicating that sequencing with 16S rRNA as gene target has been successfully identifying the isolates to specific strains. Additional PCR analysis employing SEA1/SEA2 primers increased accuracy of characterization of the isolates.

Antibiotic Resistance Among the Isolates
Four isolates showed resistance to streptomycin (Table 5). All isolates, however, were susceptible to gentamycin and oxytetracycline. Interestingly, the resistant strains to streptomycin were isolated from different food sources, i.e. from milk (S10), egg dishes (TB1), sautéed chicken cuts (UA13) and chicken cuts satay (SJ1). Next to resistance to streptomycin, isolate TB1 showed also intermediate resistance to kanamycin, chloramphenicol and tetracycline.

DISCUSSION
Among the SEs, SEA is reported as the most common cause of staphylococcal food poisoning worldwide, but the involvement of other SEs has been also found. SEA is considered as the main cause of staphylococcal food poisoning, probably due to its high resistance to proteolytic enzymes (Argudin et al., 2010). In the present study, all (eight) isolates were identified carrying sea gene, suggesting the potential risk of these strains due to production of SEA. Another study conducted in Indonesia by Salasia et al. (2011), however, did not find this gene in 11 food isolates (i.e. fermented milk product, sausage, meat ball, cakes and cheese), but found seb, sec, see and seg genes. They found, however, sea gene in 5 isolates of 19 milk isolates.
Regarding the sequencing results, as predicted, the use of 3 different primers resulted in various identity percentages. Sequencing by primer 16sF (5'-CCGCCTGGGGAGTACG-3') resulted in higher identity percentages (almost all achieved 99%) towards the existing genome of S. aureus strains in comparison to the sequencing by 63F (5'-CAGGCCTAACACATGCAAGTC-3') and 1387R (5'-GGGCGGWGTGTACAA GGC-3'). As was presented in Table 3, sequencing by 16sF also resulted in higher query coverage than that achieved by 63F and 1387R, except for isolate SJ4. These results highlighted that 16S rRNA sequencing by short sequence could provide sufficient identification amongst S. aureus strain. The 16sF has been used before to detect 16S rRNA genes of S. aureus isolates from food sample (Lee et al., 2007;Lee & Park 2016).
This study found that all isolates demonstrated 99% homology to the sequence genome of 4 strains of S. aureus by sequencing using primer 16sF. Based on the highest total score of the BLAST results, all isolates showed similarity to the genome of S. aureus strain MSSA476 and S. aureus subsp. aureus ST772-MRSA-V strain DAR4145. Strain MSSA476 is an invasive community acquired methicillin-susceptible S. aureus (Holden et al. 2004). On the other hand strain DAR4145 is a multidrug resistant strain of ST772-MRSA-V (Steinig et  al., 2015). These results suggested that the isolates in this study showed equal likelihood whether they were corresponding to the MRSA or the MSSA. This finding is possible since these strains also occupied many people in Indonesia. Recent epidemiological study showed MRSA carriage rate of 4.3%, and MSSA carriage rate of 1.5% among 1,502 patients in hospitals in Java and Bali (Santosaningsih et al., 2014).
As discused above, when sequenced by 16sF all isolates showed good homology (99%) to some strains found in the NCBI GenBank data base with different characteristics on its antibiotic resistance. Next to this characteristic, the coresponding strains also showed differences on capability to form enterotoxin. Strain MSSA476 and DAR4145, as well as strain DSM 20231 did not produce enterotoxin-A (NCBI, 2016). On the other hand, S. aureus subsp. aureus ST288 isolate 10338,10497,15532,16035,18341,18412 and 18583 are known as enterotoxin-A producers. Therefore, additional information to confirm the presence of a gene that responsible for formation of enterotoxin-A is considerably important. Since all selected isolates were confirmed carrying sea gene, indicated by the presence of 120 bp amplicon after PCR amplification using primer SEA1 and SEA2 as presented in Table 2, all isolates were most in agreeing with S. aureus subsp. aureus ST288. The sea gene is 771 bp in size encoding an enterotoxin A precursor of 257 amino acid residues (Huang et al., 1987). Specific primers SEA1 and SEA2 were frequently used in PCR analysis to detect the presence of sea gene in S. aureus isolates from food, such as in raw and pasteurized milk (Rall et al., 2008); Kérouanton et al., 2007) and ready-to-eat Kimbap (Lee et al., 2007). Specific PCR primers have commonly been employed to confirm the presence or absence of specific characteristics associated with target microorganisms such as virulence factors.
Next to genotypic characteristic, information on antibiotic resistance amongst S. aureus strains found in food is also important for surveillance and epidemiology study. This study found that four isolates from different food sources (from milk, egg dishes, sautéed chicken cuts and chicken cuts satay) were resistant to streptomycin. Streptomycin resistance among S. aureus isolates was also reported in other study. Jamali et al. (2015) found that amongst S. aureus isolates from raw milk and dairy products (n=328), 5.8% demonstrated resistance to streptomycin, 4% to kanamycin, 3.7% to chloramphenicol, and 2.1% to gentamicin. Most isolates were resistant to tetracycline (56.1%) followed by to penicillin (47.3%). The high percentage of resistant isolates to these last two antibiotics could be due to the widespread use of these antibiotics to control and treat infections on dairy farms (Jamali et al., 2015). Moreover, the fact that streptomycin resistant strains were found in milk, egg dishes and chicken dishes indicated possible  occurrence of cross contamination from human or food vendor. Schmitz et al. (1999) found that 21% of the S. aureus isolates (n=699) collected from different hospitals in Europe were resistant to streptomycin. They also found, however, 23% of the S. aureus isolates were resistant to gentamycin, 29% to tobramycin, and 31% to kanamycin. In more recent study, Onwubiko & Sadiq (2011) also found that 55.8% of S. aureus from clinical isolates in a tertiary health institution in North-western Nigeria (n= 129) showed resistance to streptomycin, 68.8% to tetracycline, 38.1% to chloramphenicol, and 7.6% to gentamycin. The spreading of S. aureus strains that resistant to antibiotics has become a global concern. Continued surveillance of S. aureus producing enterotoxin-A in milk, egg and poultry food products at genotypic levels is necessary to understand and limit further increases of staphylococcal food poisoning incidences.

CONCLUSION
This study has demonstrated that in order to increase the accuracy of the identification results, next to the sequencing of S. aureus targeting 16S rRNA gene, PCR analysis using specific primer is considerably important. All eight isolates were carrying sea gene, detected by PCR analysis, indicating that they can produce staphylococcal enterotoxin-A. Genotypic characterization of the selected strains by sequencing using 16sF, showed agreeing to the sequence genome of S. aureus subsp. aureus ST288 that also produce enterotoxin-A. This study also found that 4 of 8 selected isolates were resistant to streptomycin.


s2fieldsofstudy: Biology Medicine

text:
Peripheral blood Th9 cells are a possible pharmacodynamic biomarker of nivolumab treatment efficacy in metastatic melanoma patients

ABSTRACT Although nivolumab is associated with a significant improvement in overall survival and progression-free survival, only 20 to 40% of patients experience long-term benefit. It is therefore of great interest to identify a predictive marker of clinical benefit for nivolumab. To address this issue, the frequencies of CD4+ T cell subsets (Treg, Th1, Th2, Th9, Th17 and Th22), CD8+ T cells, and serum cytokine levels (IFNγ, IL-4, IL-9, IL-10, TGF-β) were assessed in 46 patients with melanoma. Eighteen patients responded to nivolumab, and the other 28 patients did not. An early increase in Th9 cell counts during the treatment with nivolumab was associated with an improved clinical response. Before the first nivolumab infusion, the responders displayed elevated serum concentrations of TGF-β compared to non-responders. Th9 induction by IL-4 and TGF-β was enhanced by PD-1/PD-L1 blockade in vitro. The role of IL-9 in disease progression was further assessed using a murine melanoma model. In vivo IL-9 blockade promoted melanoma progression in mice using an autochthonous mouse melanoma model, and the cytotoxic ability of murine melanoma-specific CD8+ T cells was enhanced in the presence of IL-9 in vitro. These findings suggest that Th9 cells, which produce IL-9, play an important role in the successful treatment of melanoma patients with nivolumab. Th9 cells therefore represent a valid biomarker to be further developed in the setting of anti-PD-1 therapy.

Introduction
Effective immune checkpoint blockades have improved the overall survival of patients with metastatic melanoma. The monoclonal antibody nivolumab blocks programmed cell death 1 (PD-1), an inhibitory immune checkpoint receptor expressed on activated T cells. 1 Nivolumab is associated with a significant improvement in overall survival and progression-free survival, and 20 to 40% of patients experience long-term benefit. [1][2][3] Although the expression of programmed cell death 1 ligand 1 (PD-L1) in tumor cells has been associated with responsiveness to the blockade of this immune checkpoint, 4 the objective measurement of PD-L1 protein levels reveals heterogeneity within tumors and prominent interassay variability or discordance. 5 The pharmacodynamic biomarkers of nivolumab, however, remain unknown to date. It is therefore of great interest to find a predictive marker of clinical benefit for nivolumab and a parameter that can be validated as a surrogate marker of response or survival benefit.
Recently, interleukin (IL)-9-producing CD4 C T helper cells (Th9) have been identified as a new subset of CD4 C T helper cells mediating both proinflammatory events and the induction of tolerance. 6,7 Although Th1 cells, which produce interferon (IFN)-g, have a clear role in cancer immune surveillance and in promoting antitumor responses, 8 reports on the role of Th9 cells in tumor development remain contradictory. 9 For instance, IL-9 is known to promote the proliferation, migration, and adhesion of human lung cancer cells, 10 but IL-9 seems to have the opposite effect on melanoma proliferation and migration in the B16 melanoma murine model. 11 Herein, we analyzed the immunological profile of peripheral blood in patients receiving nivolumab treatment in order to identify clinically useful biomarkers. We discovered that Th9 cells but not Th1 cells are increased in melanoma patients who were successfully treated with nivolumab. Using the autochthonous mouse melanoma model, we found that IL-9 treatment suppressed melanoma progression and increased granzyme B and perforin in CD8 C T cells. We therefore propose that elevated levels of Th9 cells may represent a valid biomarker in the setting of anti-PD-1 therapy. Our data support the idea that boosting IL-9 in itself might be a therapeutic avenue.

Th9 cell frequency is increased in responders to nivolumab treatment
Forty-six melanoma patients who received nivolumab were prospectively included in this study ( Table 1). The group contained 18 males and 28 females. The median age of the patients was 66 y (ranging from 34 to 89 y). Patients were divided into two groups, responders (SD, PR, and CR) and non-responders (PD) to nivolumab treatment. Eighteen patients responded to the treatment, and the other 28 patients did not. We found no difference in the total number of whole blood cells and lymphocytes, and in serum lactate dehydrogenase levels between responders and nonresponders (Table S1).
To investigate the candidates for biomarkers related to treatment response, we first compared the frequencies of CD8 C T cells, CD4 C T cells, and Tregs in the peripheral blood between responders and non-responders. We found no significant differences among these cell populations before and after treatment (Fig. 1A). We further compared CD4 C T cell subsets in responders and non-responders. In addition to Th1 and Th2 subsets, we investigated Th9, Th17, and Th22 subsets since recent studies showed that these Th subsets might also play some roles in tumor immunity. [12][13][14] Although there was no significant difference in Th1, Th2, Th17, and Th22 cells, Th9 cells were significantly increased in the responder group (Fig. 1B). IFNgproducing CD8 C T cell numbers before and after nivolumab treatment were comparable between responders and nonresponders (Fig. 1B). These results suggest that Th9 cells may play some role in nivolumab-induced antitumor immunity.

PD-1/PD-L1 blockade promotes Th9 differentiation
Our data showed that Th9 cells and serum TGF-b levels are increased in responders to nivolumab treatment. We therefore hypothesized that nivolumab enhanced tumor immunity by promoting Th9 differentiation. To this end, we evaluated the effect of anti-PD-1 antibody on Th9 induction in vitro. Human PBMCs were stimulated with recombinant IL-4 and recombinant TGF-b in the presence or absence of anti-PD-1 blocking antibody. After 48 h of stimulation, the frequency of Th9 cells was evaluated by flow cytometry. IL-4 and TGF-b-induced Th9 differentiation was enhanced by anti-PD-1 antibody in a dose-dependent manner (Fig. 2C). In addition, anti-PD-L1 antibody also enhanced Th9 differentiation (Fig. 2C). These results suggest that PD-1 signaling blockade by nivolumab may promote IL-4 and TGFb-dependent Th9 differentiation: Anti-IL-9 neutralizing antibody downregulates granzyme B and perforin expression in CD8 C T cells in vitro Since our data suggest that Th9 cells may promote antitumor immunity, we next investigated the effects of IL-9 on immune cells. A previous report demonstrated that Th9 cells promote the recruitment of dendritic cells (DCs) to the tumor tissues. 11 We hypothesized that Th9 cells also promoted the recruitment of CD8 C T cells to melanoma tissues. We evaluated the effect of IL-9 on the expression levels of chemokine receptors, which are responsible for CD8 C T cells. Because CXCR3 and CCR5 are reportedly important chemokine receptors for the recruitment of T cells to melanoma, 16,17 we investigated the expression levels of these receptors on CD8 C T cells cultured with or without anti IL-9 neutralizing antibody. However, we found that CXCR3 and CCR5 expression on CD8 C T cells was not affected by anti-IL-9 neutralizing antibody (Fig. 2D).
Granzyme B and perforin are cytolytic molecules produced by cytotoxic CD8 C T cells, and they have activity against a variety of tumors. 18 We therefore evaluated next the effect of IL-9 on granzyme B and perforin expression in CD8 C T cells in vitro. Human PBMCs were cultured with or without anti-IL-9 neutralizing antibody, and their expression levels were evaluated analyzed via flow cytometry. The expression levels of granzyme B and perforin in CD8 C T cells were reduced in the presence of anti-IL-9 neutralizing antibody (Fig. 2E). These results suggest that IL-9 promotes the expression of granzyme B and perforin in CD8 C T cells.
In vivo IL-9 blockade augments tumor progression in melanoma-bearing mice To investigate the effect of IL-9 on anti-melanoma immunity in vivo, we next evaluated melanoma progression through the administration of anti-IL-9 neutralizing antibody to melanoma-bearing mice. We found that anti-IL-9 neutralizing antibody administration promoted tumor progression in a melanoma cell line, the B16 melanoma cell line injection model (Fig. 3A). The subcutaneous inoculation of B16 cells mirrors human disease development poorly because tumor cells are an artificially inoculated and so it has low immunogenicity. 19 We next used the Baf/Pten autochthonous mouse melanoma model, in which melanoma develops de novo within the murine skin. 20 Consistent with the B16 injection model, the administration of anti-IL-9 neutralizing antibody also promoted tumor progression in the Braf/Pten model (Fig. 3B, Table S2). To exclude the possibility that IL-9 directly inhibits melanoma progression but does not modulate tumor immunity, we next performed a tumor proliferation assay. B16 melanoma cells were cultured with or without recombinant murine IL-9.
We found that there was no significant difference in the proliferation of B16 cells between the two groups (Fig. S1). These results support the notion that IL-9 suppresses melanoma progression via immune modulation.

IL-9 blockade leads to the downregulation of granzyme B and perforin in CD8 C T cells but not in NK cells in mice
To further investigate the mechanism of IL-9, we used the Braf/Pten melanoma model and analyzed the immune cells infiltrating into the tumor in mice treated with or without anti-IL-9 neutralizing antibody. First, the expression of granzyme B and perforin in the whole melanoma tissues was investigated by means of real-time polymerase chain reaction (RT-PCR). We found that granzyme B and perforin expression were reduced in mice treated with anti-IL-9 neutralizing antibody (Figs. 3C and D), suggesting that IL-9 promotes the expression of granzyme B and perforin in the melanoma tissues. Since both CD8 C T cells and NK cells produce granzyme B and perforin, we next analyzed the effect of IL-9 on granzyme B and perforin expression in CD8 C T cells and NK cells by flow cytometry. We already demonstrated that the expression levels of chemokine receptor responsible for tissue infiltration were not changed by anti-IL-9 treatment with human samples (Fig. 2D), suggesting that lymphocytes infiltration into the skin is not regulated by IL-9. Consistent with the above findings, there was no significant difference in the frequency of CD8 C T cells (Fig. 3E, left) or NK cells (Fig. 3F, left) infiltrating into murine melanoma tissues treated with or without anti-IL-9 antibody. Next, we evaluated the expression levels of granzyme B and perforin in CD8 C T cells and NK cells in melanoma tissues. The mean fluorescence intensity (MFI) levels of granzyme B and perforin in CD8 C T cells were significantly lower after IL-9 blockade (Fig. 3E, right), whereas MFI levels of granzyme B and perforin in NK cells were unaltered (Fig. 3F, right). These results suggest that IL-9 causes an increase in granzyme B and perforin in tumor-infiltrating CD8 C T cells.
IL-9 enhances cytotoxicity of tumor-specific mouse CD8 C T cells We evaluated the effect of IL-9 on the cytotoxic ability of tumor specific CD8 C T cells in vitro. First, we co-cultured B16 murine melanoma cells stably transduced with OVA (termed MO4 cells) and whole lymph node cells from OT-I mice (which have OVA-specific CD8 C T cells), in the presence or absence of recombinant IL-9. The MFI levels of granzyme B and perforin in CD8 C T cells increased significantly after recombinant IL-9 exposure (Fig. 4A). In addition, the cytotoxic assay showed that immune cells from inguinal lymph nodes of MO4 tumor-bearing OT-I mice killed melanoma cells more effectively when cultured with recombinant IL-9 (Fig. 4B). Furthermore, we prepared purified CD8 C T cells from inguinal lymph nodes of MO4 tumor-bearing OT-I mice, and confirmed that the tumorspecific cytotoxic ability of CD8 C T cells increased by recombinant IL-9 (Fig. 4C), while IL-9 did not directly affect the proliferation of B16 (Fig. S1). We next evaluated that the effect of IL-9 on the ability of proliferation and antigenicity of human melanoma cell lines using two human melanoma cell lines, A375 and SK-MEL-28. We found that IL-9 did not affect the proliferation of these cell lines (Fig. S2). The expression levels of HLA-ABC and HLA-DR were quantified by means of flow cytometry to evaluate the antigenicity. The expression of PD-L1 and IL-9 receptor (IL-9R) was also assessed to evaluate the effect of IL-9 on these cell lines. The results showed that IL-9 did not affect these expression levels (Fig. S2).

IL-9 is highly expressed in human melanoma lesions
Finally, we assessed 10 melanoma samples by immunohistochemistry to evaluate the localization of IL-9 C cells and CD8 C T within the tumor before nivolumab treatment. We analyzed sequential sections for these two stains using 10 samples selected from both responders and non-responders. All samples showed that high IL-9 expression and CD8 C T cell infiltration were observed in the peritumoral lesion (Fig. 4D,  Fig. S3). These findings suggest that IL-9 may be related to some extent to antitumor immunity by CD8 C T cells in the lesional area of human melanoma.

Discussion
In this study, we demonstrated that Th9 cells in peripheral blood were significantly increased in the responders to nivolumab treatment. In addition, the serum level of TGF-b, which contributes to the development of Th9, was significantly higher in responders compared to non-responders before nivolumab treatment. Moreover, anti-PD-1 antibody enhanced Th9 differentiation in vitro. Using melanoma-bearing mice, we showed that anti-IL-9 antibody decreased granzyme B and perforin in CD8 C T cells. Furthermore, in vitro experiments revealed that IL-9 enhanced the cytotoxic ability of tumor-specific CD8 C  T cells in mice. Finally, we demonstrated that IL-9-positive cells existed near CD8 C T cells in human melanoma tissues. These results suggest that Th9 cells may play an essential role in antimelanoma immunity and that anti-PD-1 antibody may elicit an antitumor effect through upregulating Th9 differentiation in melanoma patients successfully treated with nivolumab. This is the first report showing that Th9 cells in peripheral blood can be a pharmacodynamic biomarker for nivolumab efficacy in melanoma patients. The literature suggests the diverse effects of anti-PD-1 antibody treatment on antitumor immunity. Anti-PD-1 antibody caused an increase in the frequency of CD8 C T cells in the melanoma lesion of patients who responded to the treatment. 21,22 In addition, anti-PD-1 antibody mediates antitumor effects through augmented T cell proliferation, increased IFNg, and IFNg inducible chemokine production at the tumor site. 23 We also demonstrated that nivolumab increases Th9 cells both in vivo and in vitro. It has been reported that the PD-1 blockade in itself enhanced T cell migration into the tumor lesion, 23 and our study proposes the possibility that nivolumab may promote antitumor immunity in the melanoma lesion via CD8 C T cell activation by Th9 cells. This is supported by our observation that IL-9 C cells and CD8 C T cells co-localized in melanoma lesion. We are currently working to evaluate whether IL-9 C cells are increased in the melanoma lesion in responders using available samples, such as in-transit metastases.
Recent studies have shown that genetic markers might be useful as biomarkers for the effectiveness of immune checkpoint blockades in the treatment of melanoma. [24][25][26] However, no serum cytokines have been reported as useful and effective biomarkers. We demonstrated that serum TGF-b levels were significantly higher in responders compared to non-responders before nivolumab treatment. Our study therefore suggests that TGF-b may also be a good biomarker for nivolumab treatment. Several human cancer cells express high levels of TGF-b, which influences the microenvironment and promotes tumor growth, invasiveness, and metastases. 27 Increased expression and secretion of TGF-b in melanoma cell lines have been reported. 28,29 Interestingly, the expression levels of TGF-b in melanoma are different in among melanoma cell lines and patient disease stages. 29,30 Although the mechanisms that caused high serum TGF-b levels in responders remain unclear in our study, one hypothesis is that TGF-b expression may be related to mutations in the melanoma cells. In human hepatocellular carcinoma, Kras mutation deregulates the TGF-b signaling pathway. 31 Above all, the differences in TGF-b expression in melanoma may have a genetic basis. Further studies are needed to investigate the relationship between TGF-b expression and mutation burden.
The role of Th9 cells has remained controversial in tumor immunity. It has been reported that IL-9 acts directly to drive tumor growth and contributes to the establishment of an immunosuppressive environment. 9 For example, IL-9 promotes the proliferation of human lymphoid tumors, such as Hodgkin's lymphoma, diffuse large B-cell lymphoma, and NK T-cell lymphoma. [32][33][34] Furthermore, IL-9 is known to inhibit adaptive immunity and promote tumor progression using colon carcinoma and breast cancer in mice. 9 Other studies have suggested the beneficial roles of IL-9 in preventing tumor progression through a multivariate effector response. 11,12 In this study, we observed an antitumor effect of IL-9 against melanoma in two different murine melanoma models, the B16 injection model and the Braf/Pten model. This is the first report demonstrating the beneficial roles of IL-9 for tumor immunity using the autochthonous mouse melanoma Braf/Pten model. Consistent with our study, some recent reports refer to the importance of IL-9 and Th9 cells in suppressing melanoma progression. For example, one report showed that exogenous rIL-9 inhibited melanoma growth in Rag1 ¡/¡ mice but not in mast-cell-deficient mice, suggesting that mast cells are essential for IL-9-mediated antitumor immunity. 12 Another report demonstrated that Th9 cells elicited strong cytotoxic T cell responses by promoting the recruitment of DCs to the tumor tissues. 11 As for the origin of Th9 cells, one study showed that IL-1b induced Th9 cells, which in turn exerted potent anticancer functions in an interferon regulatory factor 1 (IRF1)-and IL-21-dependent manner. 35 Our study proposes a new mechanism whereby IL-9 directly enhances tumorspecific cytotoxic activity of CD8 C T cells by increasing their levels of granzyme B and perforin.
In addition, we observed that PD-1/PD-L1 blockade promotes Th9 differentiation in the present study. Several studies report that PD-1 modulates the metabolic program of T cells. 36 For example, PD-1 ligation is known to prevent T cell development by altering metabolic reprogramming of cells. 37 Since glycolytic activation is required for the differentiation of Th9 cells, 38 the PD-1/PD-L1 blockade might promote Th9 differentiation by modulating the cell metabolism. Further investigation is required on this subject.
In conclusion, we present data to show that the frequency of Th9 cells can serve as a pharmacodynamic biomarker for anti-PD-1 therapy. Of note, we found an increase in Th9 cells after the third infusion of nivolumab and provided evidence that Th9 cells promote anti-melanoma immunity. Therefore, we propose that Th9 cells may represent a biomarker to be further developed in the setting of anti-PD-1 therapy.

Patients, treatment, and clinical evaluation
This observational immunomonitoring study included 46 metastatic melanoma patients receiving nivolumab (Ono Pharmaceutical) at Kyoto University Hospital and other collaborating hospitals in Japan. This study was approved by the ethic committee of the Kyoto University Graduate School of Medicine (R0251). Patients were included if they (i) had a confirmed diagnosis of stage IV melanoma according to the 2009 American Joint Committee on Cancer (AJCC) melanoma staging and classification, (ii) were alive 12 weeks after the first nivolumab perfusion, and (iii) were receiving at least four courses of nivolumab over 90 min at a dose of 2 mg/kg of body weight every 3 weeks. 39 Other inclusion criteria were: at least 20 y of age and no specific melanoma therapy during the previous 28 d. All histological types of melanoma, including mucosal and uveal melanoma, were eligible for inclusion. Exclusion criteria were the presence of an autoimmune disease, HIV, hepatitis B or C, pregnancy, or concomitant systemic therapy or any history of prior immunotherapy for melanoma. Treatment efficacy was assessed using contrast-enhanced computed tomography (CT), magnetic resonance imaging (MRI), or positron emission tomography-CT (PET-CT) after the third nivolumab infusion and clinical response was defined based on response evaluation criteria in solid tumors (Response Evaluation Criteria in Solid Tumors, Version 1.1 (RECIST, v1.1)). A clinical response was defined as complete response (CR), partial response (PR), or stable disease (SD).

Collection of human samples and analysis of serum
Peripheral blood was taken one to 7 d before the first nivolumab infusion (pre) and within one to 3 weeks after the third infusion (post). In most cases, nivolumab was administered every 3 weeks in Japan in conformity to the national health insurance and thus the PBMCs after treatment were obtained from 10 to 12 weeks after the first infusion. Peripheral blood mononuclear cells (PBMCs) were obtained from venous blood anticoagulated with EDTA by density gradient centrifugation using lymphocyte separation solution (Nacalai Tesque). Isolated cells were cryopreserved in Bambanker (Nippon Genetics) at ¡80 C or washed with FACS buffer (PBS containing BSA and sodium azide) for flow cytometry as previously described. 40,41 Serum levels of IFNg, IL-4, IL-9, IL-10, IL-17, and transforming growth factor (TGF)-b were measured by enzyme-linked immunosorbent assay (ELISA) kits for quantitative detection of each cytokine (eBioscience, San Diego). All measurements were made in duplicate and mean values were obtained.

Analysis of peripheral blood samples
The following fluorescent-labeled monoclonal antibodies were used for surface or intracellular staining: TCRab-PerCP/Cy5. To detect intracellular protein, cells were permeabilized, fixed, and stained according to the manufacturer's instructions using the Cytofix/Cytoperm kit (BD Biosciences) for cytoplasmic targets or the Foxp3 Staining Buffer Set (eBioscience, San Diego, CA) for the nuclear target FoxP3. Tregs were defined as FOXP3 C CD45RO C C D4 C TCRab C T cells.
Acquisition was performed by eight-color flow cytometry using FACS Fortessa with FACS Diva software (both from BD Biosciences). The compensation control was performed with BD CompBeads (BD Biosciences). FlowJo software (Tree Star) was used for analysis. Data were expressed as dot plots.
B16 melanoma cells were injected into C57/B6N mice subcutaneously at each lateral region by 10 £ 5 cells in 100 mL of PBS per place. The tumor volume was measured every 3 d. Braf/Pten mice were treated topically with 10 mL of 1.9 mg/mL (5 mM) 4-Hydroxytamoxifen (4-HT, 70% Z-isomer, Sigma) in acetone at 6 to 8 weeks of age on both ears, the abdomen, and back. The tumor size (length £ width, in mm 2 ) was first measured at week 2 of tamoxifen treatment and every subsequent 5 d. Change in tumor size was expressed as percentage change compared to the first measurement. Anti-IL-9 neutralizing antibody (9C1, Bio X Cell, NH, USA) administration was started one day before the application of tamoxifen. It was injected intraperitoneally at 200 mg per mouse every 3 d for 2 weeks. Braf/Pten mice receiving control IgG or a-IL-9 every 2 d were analyzed on day 15 followed by a previous study using IL-9-neutralized melanoma mice. 11
Human IL-9 blocking assay and mouse cytotoxicity assays with recombinant IL-9 PBMCs from healthy donors were cultured with anti-CD3 antibody (0.5 mg/mL, BioLegend) for 48 h with or without anti-IL-9 neutralizing antibody (10 ng/mL, MH9D1, BioLegend). The expression levels of CXCR3, CCR5, Granzyme B, and Perforin of CD8 C T cells were evaluated by flow cytometry. Mouse T cell cytotoxicity with or without recombinant IL-9 was assayed using a Cytotox 96 nonradioactive kit (Promega) following the instructions provided.
OT-I mice, MHC class I-restricted TCR transgenic line specific for ovalbumin (OVA), 44 were inoculated with 1 £ 10 6 MO4 tumor cells into the skin of both lateral regions. On day 10 after tumor inoculation, draining lymph nodes (inguinal) were excised, mashed, and washed with RPMI in order to obtain single-cell suspensions. Purified CD8 C T cells were obtained from the lymph node cells using CD8a C T cell isolation kit (Miltenyi Biotec KK, Tokyo, Japan). Whole lymph node cells or purified CD8 C T cells were plated as the effector cells in 96-well plates at the effector/target ratio of 10/1 using 5 £ 10 3 MO4 target cells per well in RPMI lacking phenol red for 4 h at 37 C. Lactate dehydrogenase release was subsequently assessed by incubation of the supernatants with the provided substrate for 30 min and the absorbance was read at 490 nm using a Thermomax plate reader (Molecular Devices). Percentage cytotoxicity was calculated as follows: (experimental effector spontaneous ¡ target spontaneous/target maximum ¡ target spontaneous) £ 100. All cytotoxicity assays were reproducible in at least three separate assays.

Immunohistochemistry
Ten paraffinized human primary melanoma samples from 10 melanoma patients were cut into 5-mm-thick sections. Their clinical information is shown in Table S3. Antigens were retrieved by boiling in citrate buffer, pH 6.0, using a microwave. Non-specific binding of immunoglobulin G was blocked by normal goat serum (Vector Laboratories, Burlingame, CA). The sections were incubated with rabbit anti-IL-9 antibody (polyclonal, Abcam, Tokyo, Japan) overnight at 4 C. Then, they were incubated with biotinylated goat-anti-rabbit secondary antibody (Vector Laboratories, Burlingame, CA). Secondary antibodies were visualized using the Vectastain ABC-AP kit (Vector Laboratories, Burlingame, CA).

Statistical analysis
Unless otherwise indicated, data are presented as the means § standard deviation (SD) and are a representative of three independent experiments. p-values were calculated with the twotailed Student's t-test or one-way analysis of variance (ANOVA) followed by the Dunnett multiple comparison test. p-values less than 0.05 were considered to be statistically significant and are denoted by asterisks ( Ã ) in the figures.

Disclosure of potential conflicts of interest
No potential conflicts of interest were disclosed.