RNA-Seq Analysis of Colorectal Tumor-Infiltrating Myeloid-Derived Suppressor Cell Subsets Revealed Gene Signatures of Poor Prognosis

Elevated levels of myeloid-derived suppressor cells (MDSCs), including polymorphonuclear MDSCs (PMN-MDSCs) and immature MDSCs (I-MDSCs), are usually associated with disease progression in cancer patients, including colorectal cancer (CRC). However, biological mechanisms and molecular pathways regulated by MDSC subpopulations in the CRC tumor microenvironment (TME) have not been fully investigated. In this study, we performed transcriptomic analysis of tumor-infiltrating I-MDSCs and PMN-MDSCs isolated from tumor tissues of six CRC patients, compared to antigen-presenting cells (APCs). We also compared the transcriptomic profiles of tumor-infiltrating PMN-MDSCs to I-MDSCs. Our results showed different molecular pathways regulated by each MDSC subset, potentially reflecting their phenotypical/molecular/functional characteristics in the CRC TME. Moreover, we identified gene signatures in PMN-MDSC and I-MDSC of poor overall survival (OS) and disease-free survival (DFS) using the Cancer Genome Atlas (TCGA) dataset from patients with colon adenocarcinoma (COAD). However, functional studies are required to validate these findings.

The contribution of MDSCs to cancer pathogenesis and progression is well-established (3,6). Increased level of MDSCs has been associated with poor prognosis and short survival periods in CRC patients (7)(8)(9). MDSCs within the tumor microenvironment (TME) exert their suppressive activity on T cells to inhibit their anti-tumor activities (6,10). MDSCs mediate immunosuppression by expressing co-inhibitory ligands, such as PD-L1, which induces T cell dysfunction upon the interaction with its receptor PD-1, and by expressing suppressive molecules, such as arginase-1 (ARG1), inducible nitric oxide synthase (iNOS), interleukin-10 (IL-10), and transforming growth factorb (TGF-b) (11,12). Furthermore, MDSCs promote tumorigenesis via other means, such as the induction of angiogenesis and tumor growth/metastasis, and activation of cancer-associated fibroblasts (CAFs) (12,13). Up to date, the biological mechanisms and signaling pathways regulated by MDSC subpopulations have not been fully explored. Thus, further insights into these mechanisms and pathways could result in the identification of potential therapeutic targets for cancer.
Previously, we reported increased number of PMN-MDSCs and I-MDSCs in CRC tumor tissues, compared to paired-normal tissues (14), implicating the importance of MDSC role in CRC tumorigenesis and immunosuppression (7). Additionally, we reported the transcriptomic profiles of CRC tumor-infiltrating I-MDSCs and PMN-MDSCs from tumor tissues of only two patients, compared to antigen-presenting cells (APCs) (14). In this study, we extended our investigation and included more CRC patients.
We found that immune response-mediated pathways associated with dendritic cell (DC) maturation (15), triggering receptor expressed on myeloid cells 1 (TREM1) signaling (16), nuclear factor of activated T cells (NFAT)-mediated regulation of immune response (17), and Fcg receptor-mediated phagocytosis (18) were commonly downregulated in PMN-MDSCs and I-MDSCs, compared to APCs. We also compared the transcriptomic profiles of PMN-MDSCs vs. I-MDSCs and found that pathways supporting tumor growth and survival, related to metabolism, lipid biosynthesis, stress response, increased production of glucose and interaction between DCs and natural killer (NK) cells were different in tumor-infiltrating PMN-MDSCs and I-MDSCs. Therefore, these results may reflect the molecular profile and functional characteristics of each MDSC subset in the CRC TME. We also validated RNA-Seq data by confirming the mRNA expression of selected genes in the different myeloid subpopulations using qRT-PCR. We also utilized the Cancer Genome Atlas (TCGA) from patients with colon adenocarcinoma (COAD) to analyze our RNA-Seq data and identify gene signatures for PMN-MDSC and I-MDSC to predict overall survival (OS) and disease-free survival (DFS).

Sample Collection and Storage
Tumor tissues (TT) were obtained from six CRC patients (#05, 07, 08, 09, 44, and 53) who underwent surgery at Hamad Medical Corporation, Doha, Qatar. Demographical details and clinicopathological features of study population are shown in Table 1. These patients were treatment-naïve prior to surgery and provided written informed consent prior to sample collection. This study was performed under ethical approvals from Hamad Medical Corporation, Doha, Qatar (protocol no. MRC-02-18-012) and Qatar Biomedical Research Institute, Doha, Qatar (protocol no. 2018-018). Tissue specimens were frozen in 1 ml of freezing medium (10% dimethyl sulphoxide (DMSO; Sigma-Aldrich, Missouri, USA), 50% fetal calf serum (FCS; HyClone, GE Healthcare Life Sciences, Utah, USA), and 40% RPMI-1640 medium (Life Technologies, New York, USA)), then stored in liquid nitrogen to be used in batches for subsequent analyses. Tissue specimens were thawed and processed, as previously described, followed by cell staining/ sorting (14,(19)(20)(21).
The flow chart for the experimental design and tissue processing is shown in Figure 1. All experiments were performed in accordance with relevant guidelines and regulations.

Library Preparation
cDNA libraries were generated using Exome TruSeq Stranded mRNA Library Prep Kit (illumina, San Diego, USA) following the manufacturer's protocol, as previously described (22).
Libraries that passed quality control were subjected to clustering using TruSeq PE Cluster Kit v3-cBot-HS (illumina). Sequencing of clustered samples was performed on an illumina HiSeq 4000 instrument using HiSeq 3000/ 4000 SBS kit (illumina).

RNA Sequencing Data and Functional Annotation Analyses
As previously described, 150 bp depth paired-ends reads were trimmed and aligned to the hg19 human reference genome in CLC Genomics Workbench-12 (Qiagen) using default settings (14,22). The abundance of gene expression is determined by the score of TPM (Transcripts Per Million), mapped reads in CLC Genomics Workbench 12. Hierarchical clustering, principal component analysis (PCA) and differential gene expression analyses were performed, as previously described (23), using 2.0-fold change and P value <0.05 cutoffs. Ingenuity Pathways Analysis (IPA) software (Ingenuity Systems; www.ingenuity.com) was utilized to perform functional and pathway analyses on differentially expressed genes, as described previously (14,22). Raw data comparing myeloid subsets are shown in Supplementary Table  1. The flow chart for the bioinformatic tools used for RNA-Seq analysis is shown in Figure 1. Protein-protein interaction (PPI) networks among the significantly up/downregulated genes were determined by web-based online tool, STRING V11.0 (http:// string-db.org). Quantitative Real-Time Reverse Transcription PCR QuantStudio 6/7 Flex Real-time PCR system (Applied Biosystems, California, USA) was used to perform qRT-PCR for genes including CSF2, IL1B, PRF1, GZMA, IFNG, IL2RA, CD40, and b-actin with PowerUp SYBR Green Master Mix (Applied Biosystems). The relative mRNA expression was determined by the normalization to b-actin, and represented as the mean (log 10 ) ± standard error of the mean (SEM).
Supplementary Table 2 lists the sequences for the primers used.
The Cancer Genome Atlas Analysis for RNA Sequencing Data APC were selected for overall survival (OS) and disease-free survival (DFS) analyses using the GEPIA2 database, on a cohort of 269 patients with colon adenocarcinoma (COAD) from the Cancer Genome Atlas (TCGA) dataset. The 200 genes were subjected to univariate survival analysis and genes exhibiting poor survival and significant Log-rank test P value (≤0.05) were identified. The refined gene list was then subjected to forward combined gene survival analysis and genes that improved the performance (lower Log-rank test P value) were retained, while those that did not were dropped from the signature. Analysis were conducted employing the GEPIA2 algorithm, as detailed previously (24). Patients were divided into high and low groups based on median gene expression; top 50% was designated as high and bottom 50% was designated as low, and the Log-rank test was used for curve comparison. The survival signature score is calculated by mean value of log 2 (TPM + 1) of each gene.

Statistical Analyses
Statistical analyses were performed using GraphPad Prism 8 software (GraphPad Software, California, USA). For samples that passed the Shapiro-Wilk normality test, paired t-tests were performed, while those which did not pass the normality test were subjected to Wilcoxon signed-rank tests. Statistically significant P values are represented as follows; ***P < 0.001, **P < 0.01, *P < 0.05.

Hierarchical Clustering and Comparisons of Differentially Expressed Genes in Colorectal Tumor-Infiltrating PMN-MDSC, I-MDSC-MDSC, and APCs
We have previously reported that levels of PMN-MDSCs, I-MDSCs, and APCs in TT are higher than those in NT from the same CRC patients (7,14). In this study, we performed RNA-Seq to characterize the differential gene expression in these myeloid subsets. The gating strategy for sorting myeloid subsets was previously described (14) and as shown in Figure 1. CD33 + HLA-DR -/low CD14 -CD15 + were identified as PMN-MDSCs, CD33 + HLA-DR -/low CD14 -CD15were identified as I-MDSCs, and CD33 + HLA-DR + CD14 + were identified as APCs.
The hierarchal clustering shows a distinct cluster of genes which are differentially regulated in the tumor-infiltrating PMN-MDSC and I-MDSC, compared with APCs, from six CRC patients ( Figure 2A). We found that a total of 3,133 upregulated genes in PMN-MDSC and 824 downregulated genes, compared to APCs. A total of 1,206 genes were upregulated in I-MDSC and 940 genes were downregulated in I-MDSC, compared to APCs. Additionally, we found that 788 genes were upregulated and 253 genes were downregulated in PMN-MDSC, compared to I-MDSC. Principal component analysis (PCA) showed that datasets for the myeloid populations (PMN/I-MDSC) were clustered distinctly from the control APCs suggesting the similarity of gene expression patterns in MDSC subsets, compared to APCs ( Figure 2B). MDSCs clustered distinctly from APCs within the first two principal components, accounting for approximately 48.6% of the observed variation ( Figure 2B). Volcano plot shows the genes that were upregulated (shown in red), downregulated (shown in green), or remained unchanged (shown in grey) when comparing PMN-MDSC vs. APC ( Figure  2C), I-MDSC vs. APC ( Figure 2D) and PMN-MDSC vs. I-MDSC ( Figure 2E). Only genes that were significantly up/downregulated, with a fold change > 2 and P value < 0.05 cutoffs were selected for subsequent analyses.

Functional Annotation Analyses of Colorectal Tumor-Infiltrating PMN-MDSCs and I-MDSC-MDSCs
APCs are responsible for antigen processing and presentation to activate adaptive immune responses, while MDSCs are known to suppress immune responses (25). In agreement with this, our RNA-Seq data confirmed these functional characteristics of myeloid subpopulations. We analyzed the differentially expressed genes in colorectal tumor-infiltrating PMN-MDSCs and I-MDSCs, compared to APCs ( Figures 3A, B). Functional annotation analyses for top significantly affected transcripts, with a fold change > 2 and P value < 0.05 cutoffs, showed that genes related to immune system processes, positive regulation of immune response, and defense responses were downregulated in PMN-MDSCs, compared to APCs ( Figure 3A). Additionally, genes involved in DC maturation, Fcg receptor-mediated phagocytosis and NFAT-mediated regulation of immune response were significantly downregulated in PMN-MDSCs ( Figure 2A). On the other hand, genes related to ketogenesis and potentially associated with resistance to cancer immunotherapy (anti-PD-1/PD-L1) were upregulated in PMN-MDSCs ( Figure  3A). We also found that genes involved in the positive regulation of immune response, immune system processes and defense responses were downregulated in I-MDSCs ( Figure 2B). Furthermore, genes related to IL-8, IL-6, CXCR4 and IL-3 signaling pathways, Toll-like receptor signaling pathways and DC maturation were significantly downregulated in I-MDSCs ( Figure  3B). These data showed the functional characteristics of PMN-MDSCs and I-MDSCs as compared to APCs, and the different pathways associated with immune response regulation, which were downregulated in tumor-infiltrating PMN-MDSCs and I-MDSCs from CRC patients.
We found that genes related to oxidative phosphorylation, histamine degradation and NF-E2-related factor 2 (NRF2)mediated oxidative stress response were upregulated in PMN-MDSCs, compared to I-MDSCs ( Figure 3C). However, genes associated with the crosstalk between DCs and natural killer (NK) cells were significantly downregulated in PMN-MDSCs, compared to I-MDSCs ( Figure 3C). These findings suggest that pathways related to metabolism, stress response and interaction between DCs and NK cells vary in colorectal tumor-infiltrating I-MDSCs vs. PMN-MDSCs.
Next, we performed protein-protein interaction (PPI) and enrichment network analysis, using STRING web-based tool, to show the interactions of proteins within the significantly affected pathways within myeloid subsets. For this analysis, we selected deregulated genes from PMN-MDSC vs. APC (Supplementary  Figures 1-3). These data further confirm the significance of identified networks form IPA pathway analyses.

Up/Downregulated Canonical Pathways in Colorectal Tumor-Infiltrating PMN-MDSCs
Next, differentially expressed genes in colorectal tumorinfiltrating myeloid subpopulations were subjected to Ingenuity Pathway Analysis (IPA). Based on our analysis, we identified canonical and signaling pathways that were differentially upregulated or downregulated in tumor-infiltrating PMN-MDSCs and I-MDSCs (Figure 4). We found that TREM1 signaling, NFAT regulation of the immune response, Fcg receptor-mediated phagocytosis in macrophages and monocytes and DC maturation were downregulated in PMN-MDSCs, compared to APCs (-5.0 < Z score > -3.0, Figure 4A). On the other hand, pathways related to arginine biosynthesis, cancer immunotherapy and cell cycle were upregulated in PMN-MDSCs (2.0 < Z score > 3.0, Figure 4A).

Up/Downregulated Canonical Pathways in Colorectal Tumor-Infiltrating PMN-MDSCs
Next, we compared pathways that were differentially upregulated or downregulated in PMN-MDSCs, compared to I-MDSCs. We found that pathways related to the crosstalk between DCs and NK cells were downregulated, while pathways related to oxidative phosphorylation, ketogenesis, glycolysis, glucogenesis, cholesterol biosynthesis and NRF2-mediated oxidative stress response were all upregulated in PMN-MDSC, compared with I-MDSCs (-3.0 < Z score > 6.0, Figure 4C). These findings suggest that pathways related to metabolism, stress response, histamine degradation and interaction between DCs and NK cells are different in tumor-infiltrating I-MDSCs and PMN-MDSCs. Thus, they provide novel insights into the functional pathways regulated by each MDSC subset in the colorectal TME.

qRT-PCR Validation of Selected Gene Expression in Colorectal Tumor-Infiltrating PMN-MDSC and I-MDSC
We validated the expression of selected genes from RNA-Seq data in tumor-infiltrating PMN-MDSCs, I-MDSCs and APCs (the latter used as a control) by qRT-PCR. The selected genes were amongst the top significantly affected transcripts with a fold of change >2 and P value cutoff <0.05. We validated that IL1B, IL2RA, and CD40 genes were downregulated in PMN-MDSCs ( Figure 5A); these genes are important for APC and T cell activation (26)(27)(28). We validated that CSF2 and GZMB gene expressions were upregulated in PMN-MDSC ( Figure 5A); these genes encode GM-CSF and granzyme B, respectively. Upregulated expression of GM-CSF and granzyme B in PMN-MDSC could be associated with their capability of recruiting neutrophils in the TME, which potentially increases tumorassociated neutrophils, and releasing cytolytic molecules to kill immune cells that positively regulate anti-tumor immunity (29). We also validated the downregulation of IL1B and the  upregulation of IFNG, GZMB, CSF2, and PRF1 genes in I-MDSCs, compared to APCs ( Figure 5B). Although the upregulated genes, IFNG, GZMB and PRF1 (PRF1 encodes perforin 1), have been implicated in CD8 + T cell activation and function (30), they could be also associated with the promotion of I-MDSC suppressive function to enhance PD-L1 expression and kill reactive CD8 + T cells and favor tumor growth (31)(32)(33). Finally, we validated the upregulation of IL1B and the downregulation of IFNG and IL2RA genes in PMN-MDSCs, compared to I-MDSCs ( Figure 5C). This latter finding suggests that PMN-MDSCs and I-MDSCs can mediate different pathways and induce different effects on the immune response and tumor microenvironment.  (Figures 6C, D). It is worth noting that top upregulated genes in PMN-MDSC or I-MDSC vs. APC had little impact on prognosis and survival rates, hence, they were not included in these analyses.

DISCUSSION
Increased MDSC numbers in CRC tumor tissues, compared to paired-normal tissues have been previously reported, implicating the importance of MDSC function in tumorigenesis and immunosuppression (14). We have previously reported the transcriptomic profiles of CRC tumor-infiltrating I-MDSCs and PMN-MDSCs, compared with APCs (14). We found that pathways related to JNK and Wnt signaling, and SNARE complex activation were upregulated in I-MDSCs. Meanwhile, pathways related to CRC progression, cell migration, and MDSC recruitment were upregulated in PMN-MDSCs (14). These findings were obtained from RNA-Seq data analyses in two CRC patients (14). In this study, we included more patient samples in order to identify additional variations and emergence of different signaling pathways. Moreover, it is noteworthy that sorting few cells of different myeloid subsets; therefore, it was challenging to perform RNA extraction from few cells and additional amplification steps were required. Few studies have reported that cryopreservation/thawing procedure can reduce the proportion of MDSC subsets in the circulation (peripheral blood mononuclear cells; PBMC) and the expression of arginase-1 (34,35). However, the use of cryopreserved PBMC samples is acceptable as no significant changes were seen in the number of the most suppressive MDSC subset, PMN-MDSC, in fresh versus cryopreserved PBMC (36). We have also shown that expression of arginase-1 mRNA in CD33 + HLA-DRmyeloid cells, presumably suppressive subsets of myeloid cells, in cryopreserved PBMC is higher than that of CD33 + HLA-DR + APCs (37). Additionally, we have sorted different MDSC subsets from tumor tissue, which could be less susceptible to be affected by freeze-thaw procedure, unlike PBMC. Hierarchical cluster analysis and PCA showed that variations within myeloid subsets (I-MDSCs and PMN-MDSCs) are much less than those compared between myeloid subsets and APCs. These data confirm the differences in the functional characteristics of MDSCs versus APCs. In support of this, we found that genes associated with DC maturation, positive regulation of immune response, immune system processes and defense response were downregulated in PMN-MDSCs, compared with APCs. DCs, which also function as APCs, play an important role in the activation of anti-tumor immunity (15). DC maturation is very important for the activation of adaptive immunity and T cell responses (15). One of the mechanisms by which cancer cells evade anti-tumor immunity is by compromising APC function and DC maturation (38). Genetic, epigenetic and cell-mediated mechanisms can contribute to APC dysfunction and impaired DC maturation (39). Some of these cell-mediated mechanisms are driven by the action of MDSCs. Tan et al. reported that PMN-MDSCs, induced by modified vaccinia TianTan in mesothelioma mouse model, suppress DC function by releasing IL-10, resulting in the impaired induction of anti-tumor cytotoxic T cell response (40). In another study, PMN-MDSCs were shown to suppress DC maturation and T cell proliferation in autoimmune arthritis mouse model (41). In addition to DC maturation, we found that other immune response-related pathways, such as TREM1 signaling (16), NFAT-mediated regulation of immune response (17) and Fcg receptor-mediated phagocytosis (18), were downregulated in PMN-MDSCs, compared to APCs. On the other hand, PMN-MDSCs showed an upregulation of pathways related to arginine biosynthesis, which potentially lead to increased production of arginine and its consumption by tumor cells to maintain their survival and growth (42,43), and pathways potentially associated with resistance to cancer immunotherapy primarily anti-PD-1/anti-PD-L1 therapy (44).
In contrast to tumor-infiltrating APCs, pathways related to the positive regulation of immune and defense response, including IL-8/IL-6/IL-3/CXCR4 signaling, TREM1 signaling and leukocyte extravasation pathways were downregulated in I-MDSCs. Several cytokines signaling pathways, such as IL-8, IL-6, and IL-3 have been implicated in the positive regulation of immune response and the activation of T cell responses. IL-8 (45) and IL-6 (46) signaling pathways have a double-sword function in anti-tumor immunity. IL-8 signaling is known to recruit leukocytes to the TME via a chemoattractant gradient; these leukocytes play a role in innate immunity and could promote anti-tumor immune responses (45). It has been demonstrated that IL-6 production by TLR-activated APCs drives the activation of CD4 + T cells and B cell antibody response, suggesting a role for IL-6 signaling in the activation of immune responses (47). NFAT-mediated regulation of immune response and DC maturation pathways were found to be downregulated in both PMN-and I-MDSCs. NFAT is a key regulator for several key pathways, including IL-3 production, which is responsible for both lymphoid and myeloid differentiation (48). We found that IL-3RA gene was downregulated in both PMN-and I-MDSCs, compared with APCs (Figures 2A, B). Together, these data suggest the potential relationship between NFAT and IL-3 pathways, and may indicate that downregulation of NFAT-IL-3 pathways interfere with the maturation of myeloid cells in the CRC TME.
We also compared the transcriptional profiles of CRC tumorinfiltrating PMN-MDSCs with I-MDSCs. RNA-Seq data, canonical pathway and IPA analyses showed that genes and pathways related to the crosstalk between DCs and NK cells were downregulated in PMN-MDSCs, while metabolic pathways related to oxidative phosphorylation, ketogenesis, glycolysis, glucogenesis, cholesterol biosynthesis and NRF2-mediated oxidative stress response were all upregulated in PMN-MDSC. These pathways could lead to the increased consumption of glucose, cholesterol metabolism, reduced mitochondrial respiration, and increased generation of reactive oxygen species, which all promote tumor growth and survival (49). PMN-MDSCs could be responsible for enhancing the immunosuppressive environment within the tumor by mediating glycolysis and oxidative phosphorylation, leading to lipid accumulation (50). In turn, lipid uptake by tumor-infiltrating MDSCs enhances their immunosuppressive activities, which they exert on T cells, thereby promoting tumorigenesis (50). These results may reflect the phenotypical, molecular and functional characteristics of each MDSC subset in the CRC TME.
Based on TCGA analysis for the downregulated genes from our RNA-Seq data, high gene signature comprising of TREM2, CATSPER1, NPL, and PRAM1 in tumor-infiltrating PMN-MDSC predicted poor OS rates in COAD patients. Triggering receptor expressed on myeloid cells 2 (TREM2) gene is important for the activation of macrophages and DCs, and secretion of inflammatory cytokines, while PML-RARA Regulated Adaptor Molecule 1 (PRAM1) gene could be associated with neutrophil function, myeloid differentiation and T cell activation (51)(52)(53)(54). Therefore, their downregulation in PMN-MDSC could negatively influence anti-tumor immune responses and favor disease progression. However, the function of CATSPER1 and NPL in myeloid cell biology and function remains unclear. Additionally, high gene signature, with only thromboxane A synthase 1 (TBXAS1), in PMN-MDSC was found to be the best predictor for short DFS in patients. The role of TBXAS1 has been implicated in DC maturation and Th response, inflammatory response in myeloid cells and lipid metabolism (50,55). Hence, low levels of TBXAS1 could favor tumor escape from anti-tumor immunity and may also support PMN-MDSC suppressive functions. In tumor-infiltrating I-MDSC, high gene signature comprising of TNNT1, NPL and CATSPER1 predicted poor OS rates in COAD patients. Moreover, high signature of TNNT1, IL3RA and INHBA in I-MDSC predicted short DFS in COAD patients. Downregulation of IL-3 signaling pathway is associated with impairment of DC maturation and activation of adaptive immune response (56), and hence could be associated with diminished ability of tumor eradication, and poor survival rates. Like NPL and CATSPER1, the roles of TNNT1, and INHBA (gene encodes a member of TGF-b family and act as a growth/differentiation factor) in myeloid cell biology and function are not clear. Together, our findings highlight different molecular and functional pathways regulated by PMN-MDSCs and I-MDSCs in CRC TME, compared to APCs, and provide novel insights into gene signatures for each subset, which could predict poor OS and DFS in CRC patients. Some of these genes/pathways could be targeted in MDSCs and could have a clinical benefit in CRC patients. Notably, this study included six CRC patients and there were no functional assays used to validate the findings presented in the manuscript.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: BioProject ID PRJNA664006.

ETHICS STATEMENT
The studies involving human participants were reviewed and approved by Hamad Medical Corporation, Doha, Qatar (protocol no. MRC-02-18-012) and Qatar Biomedical Research Institute, Doha, Qatar (protocol no. 2018-018). The patients/ participants provided their written informed consent to participate in this study.