Key biomarkers within the colorectal cancer related inflammatory microenvironment

Therapeutic approaches focused on the inflammatory microenvironment are currently gaining more support, as biomolecules involved in the inflammatory colorectal cancer (CRC) tumor microenvironment are being explored. We analyzed tumor and paired normal tissue samples from CRC patients (n = 22) whom underwent tumor resection surgery. We assessed 39 inflammation-involved biomolecules (multiplex magnetic bead-based immunoassay), CEA and CA19-9 (ELISA assay) and the tissue expression levels of occludin and also pErk, STAT1 and STAT3 transcriptional factors (western blot). Tumor staging has been established by histopathological evaluation of HE stained tumor tissue sections. We report 32 biomarkers displaying statistically significant differences in tumor vs. control. Additionally, positive statistical biomarker correlations were found between MMP2–IL8 and BAFF–IL8 (Pearson correlation coefficients > 0.751), while APRIL–MMP2, APRIL–BAFF and APRIL–IL8 were negatively correlated (correlation coefficients < − 0.650). While APRIL, BAFF, IL8 and MMP2 did not modulate with tumor stage, they were inversely related to the immune infiltrate level and CD163 tissue expression. We conclude that the significantly decreased APRIL and increased BAFF, IL8 and MMP2 expression were tumor-specific and deserve consideration in the development of new treatments. Also, the positive correlation between Chitinase 3-like 1 and IL8 (0.57) or MMP2 (0.50) suggest a role in tumor growth and metastasis pathways.

www.nature.com/scientificreports/ therapeutic targets in CRC 7 . Several studies employing extensive multiplex assays for inflammatory profiling led to a wide diversity of encouraging results, for instance IL-6 [8][9][10][11] and IL-8 [10][11][12] , expressions are higher in serum patients with CRC, and are currently being considered as candidates in immunotherapies. Other biomarkers with diagnostic and prognostic value in CRC reported by several studies are IL-8 programmed death ligand 1 (PD-L1) 10 , CEA 11,13 , CA19-9 12 , MMP-9 13 etc. Despite such a multitude of studies, to date there is no unique panel of biomarkers capable to distinguish between different types of cancers and more studies are needed to find a reliable panel of biomarkers in terms of sensitivity and specificity for predictive and diagnostic purposes. Moreover, a combination of these markers could improve the understanding of the complex chain of molecular phenomena leading to this phenotype.
In this work, we describe through statistical correlations some possible molecular pathways involving tumorrelated inflammation markers and their association to CRC progression. Molecules that are simultaneously involved in CRC progression through common pathways might be either eligible targets for new immunotherapy drugs (aiming to improve life quality and lifespan of CRC patients) or markers in a screening panel that will monitor CRC patients evolution and improve of post-surgical care.
All the remaining raw data sets detected beyond LOD (see Supplementary Table S1), have failed the initial Shapiro-Wilks test. After log-transformation, the following biomarkers displayed a normal distribution and statistically significant differences in tumor vs. control (in alphabetical order): APRIL, BAFF, CA19-9, pERK, IL-1β, IL-1RA, IFN-β, IL-26, MMP-2, and TWEAK ( Fig. 1). Showing similar distribution as in the previous series, TSLP exhibited downregulation when comparing its log transformed levels in tumor vs. control (p = 0.001).
All these markers except APRIL, IFN-β and TWEAK, were upregulated in tumor samples compared to controls. The log-transformed levels which were normally distributed but exhibited no difference between tumor tissue and control were those of CD30 (significance value p = 0.115).
The log-transformed levels which do not exhibit a normal distribution but exhibit statistically significant tumor vs. control differences via Wilcoxon tests were (in alphabetical order): CEA, Chitinase 3-like 1 (CHI3L1), IFN-γ, IL-8, IL-22, IL-35, LIGHT/TNFSF14, OCLN, MMP-1, MMP-3, sIL-6Rα, STAT1, STAT3, sTNFR1 and sTNFR2 (Fig. 2). All these markers except sIL-6Rα, OCLN, STAT1 and STAT3, are upregulated in tumor samples compared to controls. Additionally to the biomarkers presented in Fig. 2 Histopathology and immunohistochemistry. After assessing hematoxylin and eosin slides in order to evaluate the immune infiltrate levels, our samples fell into three qualitative categories: reduced (54.5%), moderate (41%) and large (4.5%) immune infiltrate. Moreover, we have noticed that CD163 protein expression ratios (tumor/control) were sample-wise positively correlated to immune infiltrate level in hematoxilin and eosin (HE) slides i.e. reduced immune infiltrates were correlated with low CD163 ratios (0.0-0.5), moderate immune infiltrate with intermediate CD163 ratios (0.5-4.0) and large immune infiltrate with CD163 ratios (4.0-6.0). We did not have sufficient data to verify if such a correlation can be generalized to other conditions (e.g. other types of solid tumors or other CRC patient cohorts). Samples have been also analyzed by immunohistochemistry and tested for microsatellite instability. All the tumors have been found microsatellite stable.

Correlation matrix.
A new set of variables has been defined as the ratio between the biomarker levels in tumors and in control samples. Correlations between biomarker levels and ratios were gathered into a combined multi-type variables correlation matrix (see Fig. 3). This matrix has been compiled by reporting the most significant correlation coefficients as follows: log-transformed biomarker levels ratios which were normally distributed are represented by Pearson correlation coefficients, the untransformed (linear scale) biomarker level ratios are represented by Spearman coefficients, the log-transformed biomarker tumor levels which are normally distributed are represented by Pearson coefficients, and untransformed (linear scale) biomarker tumor levels are represented by Spearman coefficients.
Also, we note the strongest five negative correlations (less than − 0.  Fig. 3) have been grouped in clusters represented as triangles and quadrangles ( Fig. 4) with the sides, representing the actual two-by-two www.nature.com/scientificreports/ correlations, drawn using different color-coded types of lines corresponding to the type of variables (biomarker levels ratio or biomarker tumor levels, in linear or log scale). The four biomarkers levels ratios (tumor/control) involved in the strongest above-mentioned correlations i.e. APRIL, BAFF, IL8 and MMP2 were patient-wise negatively correlated to the CD163 levels ratios (Fig. 5).
For the normally distributed biomarkers (presented in Fig. 1) demonstrating reliable strong correlations, we have plotted their absolute levels and grouped the samples according to the tumor AJCC stages (I, II, III, IV) with their corresponding paired normal tissue controls (see Supplementary Figure S6). An Independent-Samples Kruskal-Wallis Test has been applied to the biomarker tumor levels in order to find any possible differences across staging groups. The test has been applied to all these biomarkers and also to IL-8 which was involved in several strong correlations highlighted in Fig. 4. With the exception of IL 26 (p = 0.044), there was no statistically significant difference in tumor biomarker expression across the groups containing samples of the same stage (Supplementary Table S2).

Discussion
Underexpressed and low level biomarkers. An important number of cytokine levels were below the detection limit of the method (LOD and/or LLOQ) such as IL-2, IL-10, IL-11, IL-12 (p40), IL-12 (p70), IL-19, IL-27 (p28), IL-28A/IFN-λ2, IL-29/IFN-λ1, IFN-α2 and osteocalcin in both tumor and normal tissue control samples. Most of these cytokines were previously reported as participants in the regulation of the anti-tumor immune response. These findings might suggest the existence of a pro-malignant phenotype in these CRC patients. Another potential explanation for such downregulation is that these cytokines are active mainly during very early stages of tumor formation. This group of cytokines deserves more attention in the future because it might reveal potent biomarkers for early stage tumor detection. For example, IL-27, whose concentrations were below our method detection limit, was previously shown to posess a potent antitumor activity, related not only to the induction of tumor-specific Th1 and cytotoxic T lymphocyte responses but also to direct inhibitory effects on tumor cell proliferation, survival, invasiveness, and angiogenic potential 14 . In addition, IL-27 together with IL-12 (another undetectable marker in our samples) mediates the activation of both STAT1 and STAT3 www.nature.com/scientificreports/ and also it can enhance CD4+T cell proliferation, Th1 cell differentiation, and IFN-γ production 14 . Moreover, IFN-γ production was low in both normal and tumor tissue (Supplementary Table S1) and STAT1 and STAT3 were downregulated in tumor versus normal tissue (Fig. 2n,o). Hence, IL-27 and his co-enhancing partner IL-12 may not be active in the stages were neoplastic tissue is well differentiated into a tumor, as it is in the case of our cohort formed by patients with tumors beyond the IIA stage. Both IFN-α and IFN-β, together with IL-28 have been reported to have potent antitumor activity and are currently used in the clinical treatments of several malignancies 15 . Our results show that IFN-α and IL-28 are poorly expressed and IFN-β is significantly downregulated in tumors compared to controls (Fig. 1c).
Finally, we emphasize the low IL-2 levels in our samples. The antitumor effect of high-dose IL2 therapy was demonstrated more than three decades ago 16 . Now, IL2 is used in cancer immunotherapy for the expansion of immune cells such as NK cells, T-cells, NKT-cells, cytokine-induced killer cells and is also used as an adjuvant in the treatment of patients with melanoma, advanced colorectal cancer or ovarian cancer with autologous dendritic cells stimulated by autologous tumor lysate 6 . The efficacy of these therapies is easy to understand as the suppressed IL-2 levels are a barrier in cellular immune response. CHI3L1-related pathways. CHI3L1 has been proposed both as a biomarker and a potential therapeutic target in gastric and colorectal cancer, being overexpressed in serum and tumor tissue [17][18][19][20] being involved in promoting cancer cell proliferation, invasion and metastasis 21 . As shown in Fig. 2a, the protein expression levels in cancer tissues were significantly higher than in adjacent normal tissues of the same patients (p < 0.001). Our results are in accordance with the CHI3L1 plasma levels observed in CCR patients vs. control and mRNA expression levels detected by Kawada et al. in colorectal cancer tissue versus normal tissues 17 . In addition, their results revealed that a strong CHI3L1 expression was significantly correlated with T stage, lymphatic invasion, vascular invasion and lymph node metastasis and the in vitro studies on SW480 cells indicated that CHI3L1 overexpression induces IL-8 secretion and up-regulation of Erk and JNK pathway 17 . In our cohort, immunoblot analysis of pErk1/2 ( Fig. 1j) revealed that this pathway is active, the pErk1/2 levels were significantly increased in tumor   4). The CHI3L1-Erk correlation is in agreement with the previously proposed mechanisms involved in tumor biology. Thus, CHI3L1 binding to CD44v3 induces the activation of Erk, Akt, and β-catenin signaling pathways enhancing cancer metastasis 22 . In cell lines, CHI3L1 stimulation results in the phosphorylation of Erk1/2 23 and a recombinant CHI3L1 was reported to significantly enhance the proliferation of SW480 cells, through the activation of MAPK/Erk signaling pathway 24 . Overall, CHI3L1-induced MMP2 overexpression plays a crucial role in ECM regulation promoting cancer cell growth, proliferation, invasion, and metastasis 24 . The CHI3L1, MMP2, and IL8 form a triad (Fig. 4) based on the correlation matrix ( Fig. 3), which seems to play a central role in tumor local and distal development. Another important triad correlation (Fig. 4) revealed the link between the upregulation of pErk (Fig. 1a) and CHI3L1 (Fig. 2a) and the downregulation of OCLN (Fig. 2m). The possible involvement of the Ras-Raf-MEK-ERK signaling module in regulating OCLN expression was previousley examined in Pa-4 cells transiently transfected with either an oncogenic K-RAS, an active mutant of Raf-1, Raf BXB or with constitutively active MEK1, pFC-MEK1. In all these transfected cells, northern blots revealed decreased OCLN mRNA levels. In addition, downregulation of OCLN expression induced by pFC-MEK1 was blocked by PD98059, a selective inhibitor of Erk activation. Furthermore in A549 cell model, with a high Ras-Raf signaling activity due to an oncogenic K-ras mutation, an upregulation of OCLN was demonstrated by transfecting a dominant negative Raf-1 construct or by treating the cells with PD98059 25 . These results are in agreement with our study showing that elevated pErk levels are inversely correlated with OCLN protein expression (Pearson correlation − 0.54; p < 0.01), reinforcing that the decrease in OCLN is due to the activation of the Erk signaling pathway. Moreover, 45% of patients were diagnosed with K-RAS mutation (data not shown). In addition, a Pearson correlation coefficient (− 0.49) was determined between CHI3L1 and OCLN (Fig. 3).
Selected biomarkers pairs with positive Pearson correlation coefficients (ranging between 0.47 and 0.83) and negative Pearson correlation coefficients (ranging from − 0.76 and − 0.48) (as shown in Fig. 3) have been grouped in clusters represented as triangles and quadrangles ( Fig. 4) with the sides representing the actual twoby-two correlations. Corresponding to the type of variables (biomarker levels ratio and biomarker tumor levels in linear and log scales) for which the correlation has been calculated, the sides have been drawn using different color-coded types of lines (as explained in the legend).

IL-8 an important cross-link point in cytokines constellation.
We emphasize the fact that in our tumor samples IL-8 cytokine had the highest levels vs. controls (fourfold increase, Fig. 2c, Supplementary  Table S1). Our study showed that this highly expressed biomarker is a member of a highlighted correlation triad formed by IL-8, IL-β1 and MMP-2 with Pearson correlation coefficients of 0.56, 0.48 and 0.88, respectively www.nature.com/scientificreports/ (Figs. 3, 4). In agreement with the IL-8 and MMP2 positive correlation revealed by our data, Pengjun et al. (2013) showed that compared with the healthy controls, the colorectal adenoma patients exhibited a concomitant increase of IL-8 and MMP-2 26 . Consecutively, the importance of MMP-2 was demonstrated by Groblewska et al. 27 who showed that positive tissue expression of MMP-2 was a significant prognostic factor for CRC patients survival being involved in the invasion and metastasis of CRC. An important implication in EMT was also revealed by IL-22 in relation with IL-1β 28 , which was shown to be critically involved in CRC cell metastasis 29 . We showed an important correlation between IL-8, IL-1β, MMP-2 and IL-22 (Figs. 3, 4). IL-22 has tumor-promoting properties, enhances tumor-cell proliferation, protects against apoptosis, and mediates the attraction of immunosuppressive immune cells and the release of other pro-and anti-inflammatory cytokines 30 . In addition, cancer cells induce IL-22 production by memory CD4+T cells via IL-1β in order to promote tumor growth 31 .
A study of the immune condition of CRC patients 1 has suggested a direct relation between IL-8 and CA19-9 on one hand and IL-8 and CEA on the other hand. Our study suggests that such relation might develop through intermediary correlations such as IL8-IL22-CA19-9 and IL8-IL1 β-CEA, respectively (see Fig. 4).
As discussed above, β-catenin pathway activation through IL-8 induces EMT and the increases of MMP-2 expression and activity 32 . We found a downregulation of IFN-β in tumor tissue compared to control (Fig. 1c) and a negative Pearson correlation with MMP-2 (− 0.501). Moreover, in CRC it was demonstrated that β-catenin could inhibit the expression of IFN-β and interferon-stimulated gene 56 (ISG56) by interacting with the central transcription factor, interferon-regulatory factor 3, blocking its nuclear translocation, responsible for the induction of IFN-β and hence being essential for the activation of interferon responses 33 . This result might be relevant for immunotherapy development because it has recently been shown that IFN-β can sensitize CRC cells to 5-fluorouracyl treatment with a potent effect on the reduction of tumor mass, suggesting a novel strategy to selectively target CRC 34 .

APRIL, BAFF, IL-8 and MMP-2 cluster as potential therapeutic targets.
Only scarce data exist about APRIL and BAFF roles in CRC tumor biology. The breast cancer is one of the few solid tumors described to display a differential role of BAFF and APRIL, which are produced in important quantities either by stromal cells or infiltrating neutrophils 35 . BAFF is constantly present in both normal and tumor tissue, while APRIL is preferentially expressed by the non-cancerous breast epithelial tissue, while its expression was shown to be decreased in breast tumor cells. In addition, APRIL expression was inversely correlated with the tumor stage of breast cancers 35,36 .
Recently the possible role of APRIL-BAFF and of their receptors in solid tumors has been studied by using Oncomine resources to investigate The Cancer Genome Atlas (TCGA) in order to compare tumor vs. normal mRNA expression in the whole spectrum of the samples collection. All the tumor samples analyzed exhibited APRIL and BAFF expression, together with those of their receptors. APRIL and BAFF were usually downregulated in tumors as compared to normal control. CRC was the only exception: BAFF is upregulated in tumor samples as compared to their non-transformed counterparts 36  The levels ratios (tumor vs. control) of APRIL, BAFF, IL-8 and MMP-2 were patient-wise negatively correlated in log scale with CD163 protein expression ratios (tumor/control) (as shown in Fig. 5). CD163 is a potential inflammation biomarker, which can be found either free (sCD163) or bound to the plasma membrane of macrophage cells 37 , and it is not known to be expressed in the normal colon tissue 38 . The CD163 protein tissue level can thus be a good indicator of the presence of active inflammatory foci. Moreover, CD163 has also been reported to bind TWEAK 37 that has a downregulated expression in tumor tissues vs. control (Fig. 1). The immune infiltrate levels (reduced, moderate or large) in HE slides relate in a patient-wise manner to the CD163 expression levels ratio: reduced immune infiltrate is correlated with low values of the CD163 ratio while large immune infiltrate is related with high values of the CD163 ratio. This suggests that the four marker cluster might be expressed mostly by the CRC tumor cells and less by the immune infiltrate cells.

The STAT 1 and STAT3 downregulation plays a role in tumor growth. The STAT1 and STAT3
transcription factors have been identified as major players in tumor genesis and they are being considered promising targets for cancer therapy. They are thought to play opposite roles, namely, STAT1 is involved in activating immunosurveillence and it is considered a tumor suppressor 39 , while STAT 3 is considered an oncogene, and its persistent signaling contributes to stimulate cell proliferation and prevent apoptosis 2 . STAT1 inhibition was reported in diverse tumor types, along with overexpressed STAT3 40 . However, the role of STAT3 in CRC development remains controversial, as there are reports that show elevated STAT3 had tumor promoting activity and also tumor inhibitory effects 41 . Interestingly, more recently, concomitant absence of STAT1 and STAT3 was reported in CRC tumor tissue, which was found to be significantly correlated with shorter overall survival of CRC patients 41 . In vitro experiments on various colon cancer cell lines revealed that STAT3 activity is subjected to particular changes in the inflammatory tumor microenvironment 42 , most notably IL6 high levels stimulate STAT3.
In vitro functional studies revealed that IL6 utilizes the GP130/IL6 hexameric signaling complex, which includes IL6Rα 43 , a receptor chain component which is typically upregulated in cancer cell lines. Here we report, both diminished levels of STAT1 and STAT3 in CRC as compared to non-transformed colon tissue, but also reduced IL6Rα, which may contribute to explain this phenomenon. Interestingly, recent developments of www.nature.com/scientificreports/ therapeutic approaches are considering the pharmacological blockage of STAT3 signaling 43 , however, given our findings of already reduced STAT1 and 3 in CRC tumors and also the association of low STAT1 concomitant with STAT3 with reduced survival of CRC patients 41 , great caution should be taken when designing such therapies. As a major final by-product of multiple oxidation pathways that occur in the cell, protein carbonylation is an appropriate marker of oxidative stress. Modification of the carbonylated protein levels is associated with the pathogenesis of several diseases, including colorectal cancer (CRC) 44 . The current understanding researchers in this field share is that activated oncogenes (specifically K-Ras G12D , B-Raf V619E and Myc ERT2 ) contribute to enhance the transcription Nrf2, leading to a more intense antioxidant program, that in turn contributes to decrease tumor cell reactive oxygen species (ROS) levels 45 . In support of this idea, our findings indicate that protein carbonylation levels significantly decreased in CRC tumor tissue. This may also contribute to explain why such perturbations of redox homeostasis could favor worse outcomes in CRC. Interestingly, we found positive correlations between CP and both STAT1 and STAT3 (Figs. 3, 4). Although there is tentative data suggesting that in some tumor cells both STAT1 and STAT3 post-transcriptional regulation are stimulated by ROS 46 , and specifically STAT3 phosphorylation can be activated 47 , there is insufficient data regarding colon cancer STAT1/3 ROS regulation, and this may be a good direction for further study as it could contribute to elucidate the role of ROS in CRC, and potentially clarify the controversial role of antioxidants in cancer treatment.
Finally, in the series of down regulated biomarkers related to STAT1 and STAT3, we mention TWEAK. In CRC patients, increased tumor levels of TWEAK were found to be associated with statistically significantly higher overall survival 48 . In vitro invasion assays revealed TWEAK displaying an inverse relation with tumor invasive ability. Our study reveals significantly reduced TWEAK in tumor samples, which is negatively correlated in tumor tissues with MMP1 and MMP3 (Figs. 3, 4).
This study is a systematic investigation of a wide range of inflammatory factors involved in colon cancer related signaling pathways. We report the existence of a panel of underexpressed markers [namely IL-2, IL-10, IL-11, IL-12 (p40), IL-12 (p70), IL-19, IL-27 (p28), IL-28A/IFN-λ2, IL-29/IFN-λ1, IFN-α2] in CRC patients, which should be further studied to reveal their role in establishing a pro-tumoral inflammatory conditions. Our results bring into focus an overlooked association of markers which might be specific to CRC, namely diminished APRIL levels and high BAFF, as opposed to other tumor types which display both APRIL and BAFF downregulated levels. Together with IL-8 and MMP2, these two markers show strong correlation. Our results show that, APRIL, BAFF, IL8 and MMP2 expression levels are patient-wise inversely related to the immune infiltrate level. Moreover, in the case of these four biomarkers expressions in tumor, no statistically significant differences have been found across the staging groups. New treatments could be developed based on these four biomarkers related biosignature expression pattern.

Methods
Study design. This study was performed in accordance with the Declaration of Helsinki 1975, amended in 2013. All protocols and methods carried out were reviewed and approved by the Medical Ethics Committee of Elias University Emergency Hospital in Bucharest (no: 5939/2019). Prior to being included in the study, written informed consents have been signed by all participant patients.
Our patient cohort included 28 patients who underwent surgery in order to remove colorectal tumors at Elias University Emergency Hospital between January 2019 and May 2020. However, 6 of these patients have tested positive for SARS-CoV2, the day after surgery. Since at biological level, the inflammatory profile of these patients might have been influenced by this viral infection known to trigger the so-called "cytokine storm" 49 these 6 patients have been excluded from the study. The remaining n = 22 patients have been included in the study and submitted to the study protocol described below.
In Table 1, we present the demographic and clinical data of our study cohort. The cohort is formed 54.5% (12) male and 45.5 (10) female subjects. Body Mass Index (BMI) of the whole cohort has been categorized according to the World Health Organization classification of obesity 50 . Although both obesity and CRC incidence rates are increasing, it remains unclear what relationship, if any, exists between BMI, cancer recurrence and patient survival 51 . We have defined the following categories of tumor localization according to the colorectal anatomic segments: proximal region and distal region.
Within this definition 40.9% of patients (9) were subjected to surgery, concerning tumors localized in the distal region and 59.1% of patients (13) were subjected to surgery concerning tumors situated in the proximal region.
Prior to the surgery the patients have been submitted to standard arterial tension measurements. The mean value allowed the classification of these patients according to the European Society of Cardiology (ESC) and the European Society for Hypertension (ESH) criteria 52 in optimal, normal, high normal blood pressure (BP) and grade 1-3 hypertension. Table 1 shows that most of the patients 72.7% (16) had optimal or normal BP while 27.2% (6) of them were in the early stages of arterial hypertension (high normal BP and grade 1 stages).
All samples have been assessed by HE staining, in order to establish the staging and to evaluate the immune infiltrate. The cancer staging in Table 1 has been established using the ypTNM classification of tumors from histopathological samples by AJCC classification criteria grid 53 . We have classified our samples into three relatively broad qualitative categories: reduced, moderate and large immune infiltrate. For the purpose of treatment initiation, all patients have been submitted to MSI analysis through immunohistochemistry in private specialized laboratories outside of Elias University Emergency Hospital. ELISA assay for CEA and CA19-9. CEA and CA19-9 levels from tissue lysates were quantified by using CanAg CEA EIA kit (Fujirebio Diagnostics AB, Sweden) and CanAg CA19-9 EIA 120-10 (Fujirebio Diagnostics AB, Sweden), according to the manufacturer's protocol. The protein concentration of the tissue lysate for CEA assay was 1.5 mg/mL and for CA19-9 assay was 75 μg/mL. , mouse anti-human STAT1 monoclonal antibody (MCA3469Z, 1:350 dilution factor) and goat anti-human STAT3 polyclonal antibody (AHP1076, 1:500 dilution factor), were used. HRP conjugated secondary antibodies (STAR 120 and STAR 122 at 1:500 and 1:15.000 dilution factor respectively, Bio-Rad antibodies, Hercules, CA, USA) were used. For immunostaining of the membranes, they were incubated with the primary antibodies for 2 h and with the secondary antibodies for 1 h at room temperature, under constant homogenization. Blots were revealed using the Clarity Western ECL Substrate (Bio-Rad Laboratories, Hercules, CA, USA) and the chemiluminescence signal was detected using the ChemiDoc MP System. The target proteins expression was quantified using the Image Lab software version 5.2.1 (Bio-Rad Laboratories, Hercules, CA, USA), and normalized to the total proteins transferred onto the membrane (each protein band was normalized against the total proteins transferred in the corresponding lane) 56 . Cellular carbonylated protein detection was performed using the OxiSelect Protein Carbonyl Immunoblot Kit (Cell Biolabs, San Diego, CA, USA) with a post-transfer derivatization step of protein-bound carbonyl groups by treatment with a solution of dinitrophenylhydrazine (DNPH). The formed adducts were recognized by a primary rabbit anti-DNPH (diluted 1:1000), and HRP conjugated secondary antibody anti-rabbit IgG (diluted 1:1000) as previously described 54 .

Statistical analysis. Preliminary data validation. The biomarkers levels were analyzed with Bio-Plex
Manager software version 6.0 (Bio-Rad Laboratories, Hercules, CA, USA). The coefficient of variation (COV), percentage recovery and normality distribution of data was checked before proceeding to the statistical analysis. Data distribution for sample replicates has been assessed using COV with15% intra-assay precision and 20% inter-assay precision.
To ensure the data reliability quality controls have been run in parallel to the Bio-Plex sample assay. These assay controls are in fact samples containing analytes in known concentration provided by the assay manufacturer. Percentage recovery of 80-120% is considered acceptable to support accurate sample interpretation.
Statistical analysis steps. After the validation of COV and percentage step, statistical analysis was performed as follows: normality checks, log transformations for data normalization, comparisons control vs. tumor samples and correlations. Data analysis has been performed using the IBM SPSS Statistics 26 statistical analysis package.
Our first objective is to assess the statistically significant differences for all markers in control vs. tumor samples by applying the corresponding tests. The first step in choosing the comparison test is to check the normality distribution of the raw concentration and densitometry data.
The Shapiro-Wilks test (p > 0.05) and visual inspection of their histograms, normal Q-Q plots and box plots have been used at this stage to check the normality distribution for each biomarker concentrations data for both control and tumor sample. In the Shapiro-Wilks test the null hypothesis is the normal distribution of data. This hypothesis is rejected for p < 0.05 and it is accepted for p > 0.05.
Data sets which failed the initial normality test i.e. for which p < 0.05, have been normalized through log transformation i.e. y → log (y). The newly obtained data sets have been again submitted to the Shapiro-Wilks normality test.
Our samples were paired-related samples being collected from the same patients according to the above described protocol. The normally distributed data sets have been submitted to parametric means comparison testing (t Student tests, see Fig. 1) and the data sets which were not normalized post log transformation have been submitted to non-parametric comparison testing (Wilcoxon, see Fig. 2). The biomarkers concentrations were considered statistically significant distinct if p < 0.05. The biomarkers which failed the comparison tests have been excluded.
For the purpose of correlation, a new set of variables have been defined as the ratio between the biomarker level in tumor and that in control. The levels ratios normal distribution has been checked by Shapiro-Wilks testing. The ratios failing the first stage normality testing have been log-transformed and the Shapiro-Wilks test has been applied again.
Log-transformed variables and ratios that are normally distributed have been included in Pearson correlation matrices, while the remaining variables and the ratios that were not normally distributed have been included in Spearman correlation matrices. Since all the untransformed variables and ratios were not normally distributed, the Spearman correlation coefficients were also calculated for both tumor related data and ratios.
The correlation coefficients represent a measure of a monotonic relation between the paired data. The nearer this coefficient is to ± 1 the stronger the monotonic relationship is. If the correlation is significant at a 0.01 level (2-tailed) the coefficient value is marked by two stars (**) and if this correlation is significant at a 0.05 level the coefficient value is marked by a star (*).
All the correlation matrices between untransformed tumor related data, log-transformed normally distributed variables, log transformed normally distributed ratios and the untransformed ratios have been associated to a combined multi-type variables correlation matrix (see Fig. 3). The correlation coefficients have been colored specifically to identify their origin and a hierarchy has been established according to the following criteria: (1) Pearson correlations were considered stronger than Spearman correlations, (2) 0.01 (2 tailed www.nature.com/scientificreports/ stronger than the 0.05 (2 tailed) coefficients and (3) ratios related correlations were considered stronger than untransformed simple variable correlations. All normally distributed biomarkers and IL-8 which was not normally distributed but exhibited high upregulation in tumor with respect to control have been submitted to an Independent-Samples Kruskal-Wallis Test to check across the staging groups for any statistically significant differences in the biomarker expression in tumor tissue.