Alteration of specific cytokine expression patterns in patients with breast cancer

Systemic inflammation has been associated with aggressive tumor growth, invasion, and metastasis. Here we performed a comprehensive analysis of 26 kinds of inflammatory cytokine expression patterns among 185 patients with breast cancer and 54 healthy volunteers followed by chemometric analysis. We identified the specific cytokine expression patterns of breast cancer patients compared to healthy volunteers with (1) VEGF, IL-9, GM-CSF, IL-13, IL-4, and IFNγ, (2) IL-8, IL-10, IL-12, IL-5, IL-7, IL-1α, GCSF, IL-1β, and TNFα and (3) IL-2, Eotaxin, MIP1β, MIP1α, IL-17, and bFGF. Among the patients with breast cancer, we identified the specific cytokine signature of metastatic patients compared to non-metastatic patients. We also established a mathematical model for distinguishing patients with breast cancer from healthy volunteers and metastatic patients from non-metastatic patients. This cytokine network analysis could provide new insights into early intervention and effective therapeutic strategy for patients with breast cancer.

Expression pattern of cytokines in BC patients was distinct from healthy volunteers. PCA of 185 BC patients and 54 healthy volunteers showed that the former were clearly distinguishable from the latter (Fig. 1a). Almost all healthy volunteers were located in Area 1; there was only one healthy volunteer in Area 4. In contrast, BC patients were located in broader zones. A small overlap area between healthy volunteers and BC patients can be observed. Figure 1b shows the contribution of each cytokine to Primary Component 1 (PC1) and PC2. Cytokines, which existed in Area 2 (see Fig. 1b), contributed to Areas 1 and 2 (see Fig. 1a), such as in the form of IFN-γ and IL-13, and also existed in Area 4 (see Fig. 1b) and contributed to Areas 2 and 4 (see Fig. 1a), such as in the form of IL-17 and bFGF.
Expression pattern of cytokines in cancer progression and phenotypes. Next, we investigated the expression pattern in BC patients by subgroup plotting. Figure 2a depicts the PCA of the plot according to tumor staging. The patients with Stages 1, 2 and 3 (non-metastatic patients) were primarily located in Areas 3 and 4. In contrast, the patients with Stage 4 were primarily located in Areas 1 and 2. Figure 2b depicts the PCA of the plot of the BC phenotype. There was no signature for the luminal type, whereas the HER2 type tended to exist in the upper area (Areas 1 and 2) and the TNBC type tended to exist in the lower area (Areas 3 and 4). There was no significant pattern analyzed according to Ki67 and tumor grade (Fig. 2c,d).
Correlation mapping of cytokines in BC patients and healthy volunteers. Next Prediction of BC prognosis and therapeutic effect by PCA and IL-17 expression pattern. We investigated the correlation between PCA expression patterns and BC prognosis. We followed up Stage 1-3 patients, with distant metastasis occurring in 14 out of 139 patients (Fig. 4a). Interestingly, these metastatic patients have distributed in the PC1 < 0 areas. Kaplan Mayer analysis showed that disease-free survival (DFS) was significantly shorter in patients in PC1 < 0 (Areas 2 and 4) compared to PC1 > 0 (Areas 1 and 3) (hazard ratio: 0.21; 95% confidence interval [CI], 0.075 to 0.681; P = 0.009 with the use of the log-rank test) (Fig. 4b). www.nature.com/scientificreports www.nature.com/scientificreports/ In subgroup analysis, the only grade was involved in DFS. Higher stage and HER2 and TNBC subtypes had also trend, but there was not significant due to small sample size ( Supplementary Fig. S2). We also investigated the correlation between IL-17 and BC prognosis due to IL-17 being the strongest effector about PC1. There was significantly poor prognosis in the high IL-17 group, which was defined by median levels of IL-17 (hazard ratio: 0.23; 95% CI, 0.078 to 0.662; P = 0.006 with the use of the log-rank test) (Fig. 4c). Next, we investigated the correlation between PCA and IL-17 expression patterns and therapeutic effects among BC patients who received neoadjuvant chemotherapy (Tables 2 and 3). A pathological complete response (pCR) was achieved in 36.6% of all cases (30/82), 25.0% of PC1 < 0 cases (11/44), 50.0% of PC1 > 0 cases (19/38) (Fig. 4d), 24.4% of high IL-17 cases (10/41) and 48.8% of low IL-17 cases (20/41) (Fig. 4e). We checked alteration between cytokine profile at BC diagnosis and after neoadjuvant therapy using eight patient's samples ( Supplementary Fig. S3). Our findings indicate that there were significant alteration patterns in non-pCR with recurrence case. These results suggest that PCA based on serum cytokine levels could predict BC prognosis and therapeutic effects.

Modeling of orthogonal projections to latent structures (OPLS) based on cytokine array provides information on patients with BC.
Our data based on PCA and correlation analysis suggested the existence of three different types of cytokine patterns in the patients with BC. Based on the cytokine pattern, we performed OPLS-discriminant analysis (OPLS-DA) to predict the patients' status and classification. In this analysis, healthy volunteers, patients with BC, operable patients with BC (Stages 1-3), metastatic patients with BC (Stage 4), TNBC and other types (luminal and HER2) are used as classification variants.
As shown in Fig. 5a-c, we showed the scatterplot for each variant according to the trained model based on all sample data. The red and blue plot indicates each sample and the ellipse line represents the variance in each category. Intuitively, the smaller overlapped area of variance ellipses indicates the more precise prediction of classification. Our data show that the OPLS-DA model could distinguish BC patients from healthy volunteers www.nature.com/scientificreports www.nature.com/scientificreports/ (Fig. 5a). In the analysis on operable and metastatic BC, even the overlap of variance ellipses is slightly larger than in the previous plot, while the data also suggest that our OPLS-DA model may predict both operable and metastatic BC (Fig. 5b).
In contrast, in the analysis on TNBC and other types, a large part of each ellipse line overlap, indicating that our OPLDS-DA model makes it difficult to distinguish TNBC from other types (Fig. 5c). To validate our OPLS-DA models, we performed k-fold cross-validation and predicted classification performance, here we used k = 5. The analysis of TNBC and other types could not be performed because of the insufficient sample size. Precision was calculated by using the following formula: precision = (true positive)/(true positive + false positive). The recall was also calculated using the following formula: recall = (true positive)/(true positive + false negative). As shown in Fig. 5d, the OPLDS-DA model for predicting BC patients and healthy volunteers shows satisfactory precision and recall values.
Moreover, the OPLDS-DA model for predicting both operable and metastatic BC shows potent precision and recall values. These results indicate that OPLS-DA analysis based on serum cytokine levels can precisely predict www.nature.com/scientificreports www.nature.com/scientificreports/ patients with BC. More importantly, our data showed that OPLS-DA analysis based on serum cytokine levels could predict operable or metastatic BC with a 90% probability.
Correlation between serum cytokine expression and immune profile on BC tissues. The previous study revealed that immune profile, such as tumor infiltrating lymphocytes (TILs) and FOXP3, on BC tissues has both predictive and prognostic values in BC 11,12 . However, the correlation between serum cytokine expression www.nature.com/scientificreports www.nature.com/scientificreports/ and immune profile in BC tissues is still unclear. We analyzed this correlation by evaluating TILs and FOXP3 in breast cancer tissues to uncover the immune response in breast cancer tumor microenvironment. We evaluate 83 paired samples in this study and characteristic of patients is in Supplementary Table S1. Interestingly, the majority of cytokines have a negative correlation with TILs and have a positive correlation with FOXP3 ( Fig. 6a-c). These results suggest that high TILs correlates with low pro-inflammatory cytokine expression status and high FOXP3 correlates with high pro-inflammatory cytokine expression status in BC.  www.nature.com/scientificreports www.nature.com/scientificreports/

Discussion
Tumor tissue is the gold standard for clinical and investigational sequencing. However, major barriers exist regarding acquisition and utility, such as inconvenience from a scheduling perspective, the cost of patient care, and invasive procedures for patients. The primary limitation of tissue biopsy is heterogeneity, which characterizes most advanced cancers 13 , with heterogeneity existing between metastases within the same patient. To overcome the limitations of tissue biopsy, less invasive techniques capable of capturing tumor heterogeneity and the molecular changes that cancer cells undergo when they are exposed to therapy are needed 14 . The blood test, which is the appropriate and most straightforward measure for mass screening, is widely used in the testing of several diseases including metabolic syndrome.
In this study, we presented a novel approach of cytokine networks based on the comprehensive analysis of 26 kinds of expression inflammatory cytokine expression patterns, followed by chemometric analysis. We revealed the specific cytokine expression patterns between patients with BC and healthy volunteers and the specific cytokine signature in metastatic patients with BC.
Recent studies have suggested the existence of complex crosstalk between the tumor and host, which plays a crucial role in tumor initiation, growth, metastasis, and drug resistance. The inflammatory cytokine is the main component of this crosstalk. For instance, several inflammatory cytokines are crucial in the pre-metastatic niche www.nature.com/scientificreports www.nature.com/scientificreports/ in several animal models. There were several reports, which described cytokines and BC progression. Wang et al. reported that serum IL-6, IL-8, IL-10, SCC-Ag, and CYFRA 21-1 may serve as potential markers for predicting metastasis and prognosis of BC 15 . Circulating tumor cells (CTCs), which is increased with BC progression and metastasis, involved in cytokine expression in BC. Vilsmaier et al. reported that IL-1α is a candidate marker for the release of tumor cells and IL-1β, IL-12, sFlt1 and PlGF appears to be related to CTCs release in patients with BC 16,17 . Another group showed that pro-inflammatory cytokines TNFα, IL-6, CSF-1, and IFNγ correlated with disseminated tumor cells in BC 18 . We demonstrate that cytokine expression pattern correlated with both immune profiles on BC tissues (Fig. 6a). Inflammation in a tumor microenvironment comprises infiltrating immune cells that secrete cytokines, chemokines, and growth factors to which the tumor responds 1,19 . Previously, we revealed that BC progression alters immune profile in peripheral blood mononuclear cells (PBMCs) 9 . FOXP3 expression in PBMCs is one of the genes correlated with BC progression, and we discovered that FOXP3 in BC tissues also is correlated with serum cytokine expression pattern. Overall, our results may uncover significant immune response in BC tumor microenvironment.
Consistent with these findings, clinical studies have shown the upregulation of several inflammation cytokines including G-CSF, IL-6, and IL-17 in patients with BC 5 . Meanwhile, recent reports have shown that IL-17-producing cells increased in tumor tissues and peripheral blood from different cancer patients 20 {Numasaki, 2004 #707 . In contrast, a blood test to measure IL-17 levels is the more feasible way to evaluate the significance of IL-17 in tumor microenvironments, compared with measuring IL-17-producing cells. In non-small-cell lung cancer and hepatocellular carcinoma, serum IL-17 levels were correlated with poor prognosis 21,22 . Here, we are the first to show that serum IL-17 is up-regulated in BC patients and that high levels of IL-17 are also correlated with poor clinical outcomes.
Recently, preclinical models have shown that IL-17 mediates cancer progress via various mechanisms 23-28 . Coffelt et al. reported that increasing systemic IL-17 levels leads to the up-regulation of G-CSF, which subsequently causes neutrophil expansion and alteration of the neutrophil phenotype. Another study reported the effects of IL-17 and bFGF on mesenchymal stem cell (MSC) growth and elucidated the signaling pathways involved. These results indicated that IL-17 increased the frequency of colony-forming unit fibroblast (CFU-F), as well as MSC proliferation in a dose-dependent manner, while bFGF supplementation induced an increase in cell proliferation 27 . Consistent with these findings, our data indicated that the interaction of serum IL-17 levels with an increased level of both bFGF and G-CSF are correlated with poor DFS among BC patients. While IL-17and IL-17-producing cells are significantly insensitive towards chemotherapy in experimental settings 29 , little www.nature.com/scientificreports www.nature.com/scientificreports/ is known about the correlation between systemic IL-17 levels and sensitivity towards chemotherapy in clinical settings, especially concerning BC. Although there are limitations to our study, given that patient characteristics include all subtypes and involve various chemotherapy regimens, our data indicate that systemic IL-17 levels are correlated with pathological complete response rates in BC patients. To overcome such limitations, we need to further investigate the correlation between serum IL-17 levels and clinical outcomes in larger-volume clinical trials.
We acknowledge that this study has several other limitations. Notably, this study is retrospective, meaning that it is difficult to evaluate correlations between cytokine diversity and clinical outcome due to treatment and phenotype bias. We should validate using individual cohort to enhance our conclusion. To overcome this limitation, we are planning to validate using two ongoing clinical trial data set. The first one is to investigate the efficacy and safety of dual-HER2 blockage therapy in HER2-positive breast cancer (Clinical trial information: UMIN000007576). The other one is to evaluate the safety, antitumor effect and improved prognosis of Nivolumab with radiation therapy in patients with HER-2-negative metastatic breast cancer have bone metastasis that possible to irradiate (Clinical trial information: UMIN000026046). These clinical trials will cover all phenotypes. Additionally, we will investigate not only cytokine expression pattern but also tumor immune microenvironment status using CyTOF analysis and RNA-seq analysis of PBMCs. Then, we will confirm the correlation between cytokine expression pattern and tumor microenvironment status.
In this study, at least in the evaluations of BC development and progression, our analysis could be potentially relevant as we revealed them by OPLS-DA (Fig. 5). To further strengthen our model, we performed k-cross-validation as it is used for the development of biomarker classifiers from high dimensional data set 30 . Interestingly, results of k-fold cross-validation showed high accuracy to distinguish BC from HV and stage IV patients from stage I-III patients. Additionally, the standard deviation value is less than 5% both recall and precision rate. Therefore, we think the risk of over-learning is very low. Considering the application, our data are crucial because our results suggest that serum analysis is a useful approach to provide information, for example, on the subtype of cancer, inflammation status and prognosis by a 12.5-µl serum sample. Regarding accessibility, for example, in terms of low invasion, established method and low cost, the comprehensive analysis of inflammatory cytokines in serum could be useful for monitoring inflammation status in the use of antagonizing/neutralizing antibodies to address inflammatory cytokine signaling.

Methods
Serum samples. Serum from healthy volunteers, early-stage BC patients, and metastatic BC patients was collected by the Department of Breast Surgery, Kyoto University Hospital, and the Resource Center for Health Science, Kyoto, Japan. All the BC patient's samples were collected at the BC diagnosis between April 2011 and March 2015. Eight patients had two-point samples, both at diagnosis and after neoadjuvant therapy. Written informed consent was given by all patients before collection. All study protocols were approved by the Ethics Committee for Clinical Research, Kyoto University Hospital (Authorization Number G424) and performed according to the Declaration of Helsinki and the Ethical Guidelines for Epidemiological Research of the Ministry of Education, Culture, Sports, Science, and Technology and the Ministry of Health, Labour, and Welfare of Japan. Samples were frozen in liquid nitrogen and stored at −80 °C within 2 hours from sample collection. We checked if sample storage period effects on cytokine expression in this study ( Supplementary Fig. S4). All of the value of correlation coefficient were within 0.2, and there was no significant distribution in the characteristic of patients.
Analysis of the immune profile on BC tissues. The paraffin sections were rehydrated and immersed in 3% hydrogen peroxide for 10 min to quench endogenous peroxidase activity. Epitope retrieval was performed by heating the sections at 95 °C in Target Retrieval Solution (pH 9, S2375; Dako) for 30 min in an electric pressure cooker (SR-P37; Panasonic, Tokyo, Japan). All incubations after the blocking step were performed at room temperature. The wash buffer was 0.05 M Tris-buffered saline with 0.05% Tween 20 (pH 7.6). Following blocking with 5% normal goat serum (Abcam) in PBS for 10 min, the sections were incubated with anti-FOXP3 antibodies (236 A/E7; Abcam) for 1 h.
Histopathologic assessment of the percentage of TILs and FOXP3 was performed on one representative immunohistochemical section of a tumor using methods recommended by the International TILs Working Group 2014 18 . % TILs was defined as the percentage of lymphocyte area per tumor area in vivo samples. Areas of in situ carcinomas, normal lobules, necrosis, hyalinization, and crush artifacts were not included. Histopathologic evaluation of TILs was performed by two breast pathologists who scored each case independently in a blind manner. The mean values of both observers were obtained as the final scores for each case.
Statistics and data analysis. Multivariate statistical analysis (PCA and OPLS-DA) was performed in order to evaluate the measured cytokine data using R version 3.3.1 (R Foundation for Statistical Computing, 2016) and the R ropls package version 1.8.0. PCA can project higher-dimension data into lower-dimensional space in order to visualize the distribution of complex data. Typically two-dimension is used as the lower-dimensional space. Usually, PCA results are represented by a score plot and a loading plot. In a score plot, samples expressing www.nature.com/scientificreports www.nature.com/scientificreports/ similar cytokine values are plotted in the near position, while the spatial groups of similar samples are recognized as a pattern. In a loading plot, the positions of each point corresponding to the cytokine variables describe the correlation with the PCs, such as PC1 and PC2. The loading plot implies which cytokines mainly contribute to forming the pattern in the score plot. In our PCA calculation, a few outlier samples, which were more than 10-100 times larger than the average of the distribution, were omitted because substantial outlier values collapse the pattern. OPLS-DA is a potent statistical tool, which provides insights into distinguishing between groups of high-dimensional spectral data 31 . The score axis is the predictive component, which is the variation correlated to the factor of interest. The orthogonal score axis is the first component of the uncorrelated variations. A χ 2 test (two-sided, α = 0.05) was used to compare the populations between groups in the score plot, using SPSS 17.0 (Chicago, IL, USA).
k-fold-cross-validation: The dataset was randomly partitioned into k = 5 disjoint subsets to be used as a training and testing sets to validate the Delphi exercise-approved diagnostic criteria. Four subsets were used as the training set to build a diagnostic model, and the remaining subset was used as the testing set to test the model. This process is repeated five times for each dataset, each round with a different training set and complementary testing set. Model performance was estimated by calculating precision and recall which definitions were described at the results section. The averages and standard deviations over the 5 datasets were calculated.

Data Availability
The datasets generated during and analyzed during the current study are not publicly available due to protecting participant confidentiality but are available from the corresponding author on reasonable request.