Transcriptomic Immune Profiles Can Represent the Tumor Immune Microenvironment Related to the Tumor Budding Histology in Uterine Cervical Cancer

Tumor budding (TB) histology has become a critical biomarker for several solid cancers. Despite the accumulating evidence for the association of TB histology with poor prognosis, the biological characteristics of TB are little known about in the context related to the tumor immune microenvironment (TIME) in uterine cervical cancer (CC). Therefore, this study aimed to identify the transcriptomic immune profiles related to TB status and further provide robust medical evidence for clinical application. In our study, total RNA was extracted and sequenced from 21 CC tissue specimens. As such, 1494 differentially expressed genes (DEGs) between the high- and low-TB groups were identified by DESeq2. After intersecting the list of DEGs and public immune genes, we selected 106 immune-related DEGs. Then, hub genes were obtained using Least Absolute Shrinkage and Selection Operator regression. Finally, the correlation between the hub genes and immune cell types was analyzed and four candidate genes were identified (one upregulated (FCGR3B) and three downregulated (ROBO2, OPRL1, and NR4A2) genes). These gene expression levels were highly accurate in predicting TB status (area under the curve >80%). Interestingly, FCGR3B is a hub gene of several innate immune pathways; its expression significantly differed in the overall survival analysis (p = 0.0016). In conclusion, FCGR3B, ROBO2, OPRL1, and NR4A2 expression can strongly interfere with TB growth and replace TB to stratify CC patients.


Introduction
Uterine cervical cancer (CC) is one of the most common cancers in women [1]. Although progress has been made in CC prevention and management, it still has a poor outcome. According to Globocan 2020, 3218 new CC cases and 1014 CC deaths occur annually in Korea [1]. Therefore, novel markers for CC diagnosis, prognosis, and treatment have long been of great interest.
Recently, tumor budding (TB) histology has become a critical biomarker for several solid cancers, including colorectal cancer, head and neck cancer, pancreatic cancer, etc. [2][3][4][5][6]. In addition, we previously showed that TB status is a potential independent prognostic factor and can strongly interfere with CC treatment strategies [7]. However, the main disadvantage of this histopathological marker is that it can only be established after surgery. In this scenario, insights into the biological characteristics of TB can provide a helpful direction to investigate potential tools for CC management. The tumor microenvironment (TME), referred to as a "house" of cancer cells, comprises non-resident components (tumor-derived cells and infiltrating leukocytes) and resident components (blood vessels, nerve fibers, mesenchyma, and structural components) [8][9][10]. Some researchers suggested six TME subgroups, including hypoxic niche, immune microenvironment, metabolism microenvironment, acidic niche, innervated niche, and mechanical microenvironment [8]. The tumor immune microenvironment (TIME) is an important aspect because the theory of TIME can provide a powerful account of how to develop effective anti-tumor therapies. It is revealed that the immune response in TME can sometimes block cancer development [11]. Unfortunately, in many cases, tumor cells often activate the immunosuppressive mechanisms and lead to tumor escape from the host immune response [8][9][10].
Until now, TB has been little known in relation to the TIME. Nearchou (2019) and Dawson (2020) found that CD8+/CD3+ T-cell location and density are correlated with TB status in colorectal cancer [12,13]. However, there have not been any published studies about the relationship between tumor buds and TIME in CC. Therefore, this study aimed to identify the transcriptomic immune profiles related to TB status and further provide robust medical evidence for clinical application, such as stratification, prognosis, and immune-targeted therapeutic strategy.

Sample Collection
This study was conducted between 2011 and 2018. Twenty-one tissues from early stage and locally advanced CC patients were obtained after radical hysterectomy. The exclusion criteria included a history of preoperative chemotherapy, radiotherapy, and synchronous malignancies. The 2009 International Federation of Gynecologic Obstetrics (FIGO) staging for carcinoma of the cervix was used to stage the patients [14]. The flowchart in Figure 1 represents the overall process of this study.
prognostic factor and can strongly interfere with CC treatment strategies [7]. However, the main disadvantage of this histopathological marker is that it can only be established after surgery. In this scenario, insights into the biological characteristics of TB can provide a helpful direction to investigate potential tools for CC management.
The tumor microenvironment (TME), referred to as a "house" of cancer cells, comprises non-resident components (tumor-derived cells and infiltrating leukocytes) and resident components (blood vessels, nerve fibers, mesenchyma, and structural components) [8][9][10]. Some researchers suggested six TME subgroups, including hypoxic niche, immune microenvironment, metabolism microenvironment, acidic niche, innervated niche, and mechanical microenvironment [8]. The tumor immune microenvironment (TIME) is an important aspect because the theory of TIME can provide a powerful account of how to develop effective anti-tumor therapies. It is revealed that the immune response in TME can sometimes block cancer development [11]. Unfortunately, in many cases, tumor cells often activate the immunosuppressive mechanisms and lead to tumor escape from the host immune response [8][9][10].
Until now, TB has been little known in relation to the TIME. Nearchou (2019) and Dawson (2020) found that CD8+/CD3+ T-cell location and density are correlated with TB status in colorectal cancer [12,13]. However, there have not been any published studies about the relationship between tumor buds and TIME in CC. Therefore, this study aimed to identify the transcriptomic immune profiles related to TB status and further provide robust medical evidence for clinical application, such as stratification, prognosis, and immune-targeted therapeutic strategy.

Sample Collection
This study was conducted between 2011 and 2018. Twenty-one tissues from early stage and locally advanced CC patients were obtained after radical hysterectomy. The exclusion criteria included a history of preoperative chemotherapy, radiotherapy, and synchronous malignancies. The 2009 International Federation of Gynecologic Obstetrics (FIGO) staging for carcinoma of the cervix was used to stage the patients [14]. The flowchart in Figure 1 represents the overall process of this study.

Pathological Process
Specimens were selected from the tumor area and stained with hematoxylin and eosin. All tissue slices (ranging from 8 to 25 per patient) were carefully assessed for histopathological features. TB is described as a single neoplastic cell or cell cluster of up to four neoplastic cells at the invasive front of the tumor [15]. High TB was defined as ≥5 buds/high-power field (Supplementary Figure S1). Detailed examination methods and high-TB criteria have been described previously [7,16].

Clinical Parameters and Follow-Up
Clinicopathological parameters included age, tumor stage (early stage < IIb and late stage ≥ IIb), histological subtype, TB status (high and low), overall survival (OS), and relapse-free survival (RFS). All patients were followed-up every three months for the first two years, every six months for the next five years, and then annually [7,16].

RNA Extraction and Sequencing
For each patient, total RNA was extracted from 2 to 4 sections (5 µm each) of a block using the ReliaPrep TM FFPE RNA Miniprep System (Promega, Madison, WI, USA), according to the manufacturer's protocol. Next, we conducted library preparation using the TruSeq RNA Exome Kit and RNA sequencing with NovaSeq 6000 (Illumina, San Diego, CA, USA). The quantity and quality were checked using a Qubit 4 Fluorometer (Thermo Fisher Scientific, Waltham, MA, USA), and Bioanalyzer (Aligent Technologies, Santa Clara, CA, USA).

Bioinformatic Analysis
Raw sequences in the FASTQ format from the sequencer were assessed for read quality using FASTQC (https://www.bioinformatics.babraham.ac.uk/projects/fastqc/, accessed on 8 September 2021). All low-quality reads and sequencing adapters were removed using Trimmomatic [17,18]. Kallisto was used to quantify RNA-seq data [19]. DESeq2 package was used to compare gene expression between high-and low-TB groups [20]. The p value and |log2-fold-change| thresholds were 0.05 and 1.0, respectively. A volcano plot of differentially expressed genes (DEGs) was created using the "ggplot2" package [21].

Functional Analysis
The IR-DEGs were uploaded to the DAVID website (https://david.ncifcrf.gov/tools. jsp, accessed on 19 April 2022) for analyzing the Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways [23]. Benjamini-and-Hochberg-adjusted p < 0.05 was considered statistically significant.

Immune-Cell-Type Analysis
Gene-length-normalized expression was imported into the TIMER2.0 website (http: //timer.cistrome.org/, accessed on 21 April 2022) to estimate the proportions of different immune cell types using the TIMER (TIMER immune cells) and CIBERSORT (CIBERSORT immune cells) algorithms [24][25][26]. The fraction of immune cells between the high-and low-TB groups was compared using the Wilcoxon rank-sum test and a boxplot was drawn with the "ggpubr" package [27]. Benjamini-and-Hochberg-adjusted p < 0.05 was considered statistically significant.

Candidate Gene Analysis
We performed Least Absolute Shrinkage and Selection Operator (LASSO) regression, a machine learning method in the "glmnet" package, to select hub genes from IR-DEGs [28]. Then, the correlation analysis (Spearman's method in the "corrplot" package) between the hub genes and immune cell types was used to select the significantly correlated genes [29]. The candidate genes were the overlapping genes between CIBERSORT and TIMER. The p value threshold was 0.05.

Receiver Operating Characteristic (ROC) Curve and Survival Analysis
The "pROC" package was used to determine the diagnostic accuracy of the candidate genes between the high-and low-TB groups in CC [32]. Next, the optimal cut-off point for survival analysis was determined using maxstat (maximally selected rank statistics) in the "survminer" package [33]. The different expression levels of candidate genes in OS and RFS were explored using the log-rank test in the "survival" package [34,35].

Identification of IR-DEGs
The volcano plot indicates 1464 DEGs (929 upregulated and 535 downregulated genes) between the high-and low-TB groups (Figure 2 and Supplementary Table S1). Subsequently, 106 IR-DEGs that overlapped between the DEGs and immune-related gene database ( Figure 3A) were selected. A heatmap was then used to visualize the hierarchical clustering of the identified IR-DEGs ( Figure 3B and Supplementary Table S2).
The volcano plot indicates 1464 DEGs (929 upregulated and 535 downregulated genes) between the high-and low-TB groups (Figure 2 and Supplementary Table S1). Subsequently, 106 IR-DEGs that overlapped between the DEGs and immune-related gene database ( Figure 3A) were selected. A heatmap was then used to visualize the hierarchical clustering of the identified IR-DEGs (Figure 3B and Supplementary Table S2).

Functional Enrichment Analysis
The functions of 106 IR-DEGs were investigated by GO enrichment and KEGG pathway analyses. As shown in Figure 4A and Supplementary Tables S3 and S4, the 50 most genes) between the high-and low-TB groups (Figure 2 and Supplementary Table S1). Subsequently, 106 IR-DEGs that overlapped between the DEGs and immune-related gene database ( Figure 3A) were selected. A heatmap was then used to visualize the hierarchical clustering of the identified IR-DEGs ( Figure 3B and Supplementary Table S2).

Correlation Analysis
We first investigated the proportion of immune cells using the CIBERSORT and TIMER algorithms (Supplementary Table S5). Figure 6A (CIBERSORT) and 6B (TIMER) indicate that each immune cell type was not significantly different between the low-and high-TB groups.
Spearman's method was used to identify the correlation between hub genes and immune cell types generated by CIBERSORT and TIMER (Figure 7). Eleven ( Figure 7A) and five genes ( Figure 7B) were significantly correlated with CIBERSORT and TIMER immune cells, respectively. Finally, the Venn diagram ( Figure 7C) identified four candidate genes (overlapping genes), including one upregulated (FCGR3B) and three downregulated (ROBO2, OPRL1, and NR4A2) genes.

Correlation Analysis
We first investigated the proportion of immune cells using the CIBERSORT and TIMER algorithms (Supplementary Table S5). Figure 6A (CIBERSORT) and 6B (TIMER) indicate that each immune cell type was not significantly different between the low-and high-TB groups.
Spearman's method was used to identify the correlation between hub genes and immune cell types generated by CIBERSORT and TIMER (Figure 7). Eleven ( Figure 7A) and five genes ( Figure 7B) were significantly correlated with CIBERSORT and TIMER immune cells, respectively. Finally, the Venn diagram ( Figure 7C) identified four candidate genes (overlapping genes), including one upregulated (FCGR3B) and three downregulated (ROBO2, OPRL1, and NR4A2) genes.

PPI Network Analysis
A network of 106 IR-DEGs was constructed using STRING and displayed in Cytoscape (Figure 8). The confidence score was 0.4 and the PPI enrichment p value was 1 × 10 −16 . FCGR3B, ROBO2, OPRL1, and NR4A2 were significantly enriched in seven KEGG pathways and 51 GO terms (Tables 2, 3, S6 and S7). Interestingly, FCGR3B was the hub gene for the natural-killer-cell-mediated cytotoxicity and phagosome pathways.

PPI Network Analysis
A network of 106 IR-DEGs was constructed using STRING and displayed in Cytoscape ( Figure 8). The confidence score was 0.4 and the PPI enrichment p value was 1 × 10 −16 . FCGR3B, ROBO2, OPRL1, and NR4A2 were significantly enriched in seven KEGG pathways and 51 GO terms (Tables 2, 3, S6 and S7). Interestingly, FCGR3B was the hub gene for the natural-killer-cell-mediated cytotoxicity and phagosome pathways. "Red" and "blue" represent the upregulated and downregulated genes, respectively. The gray lines (edges) indicate interactions between connected nodes. Star highlights the candidate gene. The green, purple, and orange circles represent the natural-killer-cell-mediated cytotoxicity, phagosome, and axon guidance pathway, respectively.  "Red" and "blue" represent the upregulated and downregulated genes, respectively. The gray lines (edges) indicate interactions between connected nodes. Star highlights the candidate gene. The green, purple, and orange circles represent the natural-killer-cell-mediated cytotoxicity, phagosome, and axon guidance pathway, respectively. Regulation of multicellular organismal process 42 2.67 × 10 −08

ROC Curve and Survival Analysis
ROC analysis was conducted to investigate the accuracy of the four candidate genes as diagnostic biomarkers for low-and high-TB levels. The area under the curve (AUC) of FCGR3B, ROBO2, OPRL1, and NR4A2 was 81.7%, 92.2%, 82.2%, and 94.4%, respectively ( Figure 9). However, only the FCGR3B expression level significantly differed in OS (p = 0.0016; Figure 10).

Discussion
Recent studies have noted the high correlation between TB status and cancer progression [5][6][7]16]. Because tumor progression results from the interaction between tumor cells and the TME, insights into biological TB can provide potential tools for cancer management [5]. However, the relationship between tumor bud and TME, particularly TIME, has been unclear. Therefore, this is one of the first studies identifying the transcriptomic immune profiles of TB in CC.
In this study, 1464 DEGs were obtained by comparing the gene expression between high-and low-TB groups. By intersecting these DEGs and 2483 public immune-related genes, we selected 106 IR-DEGs. Functional analysis showed that these IR-DEGs mainly involve innate immune pathways, such as the natural killer (NK)-cell-mediated cytotoxicity, cytokine-cytokine receptor interaction, chemokine signaling, JAK-STAT, and Tolllike receptor signaling pathways. Most current reports consider CD8+T cells, FOXP3+T cells, and CD68+ macrophages to be the primary innate immune response in the tumor microenvironment [5,36]. Nevertheless, NK also plays a crucial role in the TIME. For example, Garcia-Iglesias et al. reported that NK-activating receptors and the cytotoxic

Discussion
Recent studies have noted the high correlation between TB status and cancer progression [5][6][7]16]. Because tumor progression results from the interaction between tumor cells and the TME, insights into biological TB can provide potential tools for cancer management [5]. However, the relationship between tumor bud and TME, particularly TIME, has been unclear. Therefore, this is one of the first studies identifying the transcriptomic immune profiles of TB in CC.
In this study, 1464 DEGs were obtained by comparing the gene expression between high-and low-TB groups. By intersecting these DEGs and 2483 public immune-related genes, we selected 106 IR-DEGs. Functional analysis showed that these IR-DEGs mainly involve innate immune pathways, such as the natural killer (NK)-cell-mediated cytotoxicity, cytokine-cytokine receptor interaction, chemokine signaling, JAK-STAT, and Toll-like receptor signaling pathways. Most current reports consider CD8+T cells, FOXP3+T cells, and CD68+ macrophages to be the primary innate immune response in the tumor microenvironment [5,36]. Nevertheless, NK also plays a crucial role in the TIME. For example, Garcia-Iglesias et al. reported that NK-activating receptors and the cytotoxic activity of NK cells significantly decrease in high-grade squamous intraepithelial lesions and CC [37].
To investigate the candidate genes from the 106 IR-DEGs, multiple analyses were conducted. Briefly, the LASSO regression algorithm was implemented and 16 hub genes were identified. Then, the correlation between hub genes and immune cell types was analyzed to select candidate genes. Finally, one upregulated (FCGR3B) and three downregulated (ROBO2, OPRL1, and NR4A2) genes were identified. PPI network analysis revealed the functional insights of the four candidate genes in TB status.
First, high FCGR3B expression induces TB progression by suppressing NK-cell-mediated cytotoxicity and phagocytosis. FCGR3B encodes a low-affinity immunoglobulin γ Fc region receptor III-B protein (FcγRIIIb), a member of the IgG Fc receptor family (FcγRs) [38]. FcγRs play crucial roles in cancer immunotherapy via antibody-dependent cellular cytotoxicity and phagocytosis [38].
Second, low ROBO2 expression may lead to a high-TB status through the axon guidance signaling pathway. The ROBO2 gene encodes roundabout guidance receptor two protein that functions as an axon guidance receptor by binding to secreted SLIT ligands [45][46][47]. The ROBO family is frequently downregulated in several cancers and considered an antioncogene [48]. The SLIT/ROBO signaling can inhibit tumor progression through some mechanisms, such as preventing cell migration and angiogenesis, enhancing cell-cell adhesion, and blocking endothelial cell proliferation [49]. Interestingly, this signaling can regulate macrophage immune responses by inducing cytoskeletal changes in macrophages, preventing macrophage spreading and inhibiting macropinocytosis [50]. Although the ROBO2 expression level in this study was not significantly related to patient survival prognosis because of the small sample size, ROBO2 can still be a promising marker for predicting TB status (AUC = 92.2%).
Lastly, NR4A2 (nuclear receptor subfamily 4 group A member 2) and OPRL1 (opioidrelated nociceptin receptor 1) were significantly enriched in several important biological processes, such as signal transduction, signaling receptor binding, signaling receptor activity, cell-cell signaling, and regulation of cell death. These findings suggest that OPRL1 and NR4A2 downregulation may enhance TB progression in CC. In other words, these genes function as tumor suppressors. With substantial diagnostic accuracy from ROC analysis, OPRL1 and NR4A2 are promising candidates for predicting TB status. Although there are several conflicting data regarding the role of NR4A2 and OPRL1 in cancer, our results are in accordance with current observations: NR4A2 can trans-activate Foxp3, involved in the differentiation, maintenance, and function of regulatory T cells, and plays a significant role in cancer cell development and survival [51][52][53]. Furthermore, Inamoto et al. reported that NR4A2 is a tumor suppressor in human bladder cancer tissues [51]. Regarding OPRL1, it belongs to the Aγ family of G protein-coupled receptors, which are involved in various diseases, including cancer [54,55]. Bedini et al. revealed that OPRL1 acts as a tumor inhibitor in U87 glioblastoma cells by blocking lipopolysaccharide [56]. In addition, OPRL1 can activate markers on the surface of T cells, enhance CD4+ T and CD8+ T-cell proliferation, and alter cytokine secretion, which is closely involved in tumor progression [57].
This study lacks experimental validation. Therefore, further studies should be conducted to clarify the mechanisms of these candidate genes in TB progression.

Conclusions
Four immune-related genes (one upregulated (FCGR3B) and three downregulated (ROBO2, OPRL1, and NR4A2) genes) that can impact TB formation and development were detected through a multiple-step analysis. These genes may be potential biomarkers for replacing TB status to stratify CC patients. Furthermore, they can provide novel ideas for immune-targeted therapeutic strategies in the future.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/genes13081405/s1. Table S1: Significant differentially expressed genes between the low-and high-tumor-budding group. Table S2: Immune-related differentially expressed genes. Table S3: Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways. Table  S4: Gene Ontology (GO) terms. Table S5: Proportion of immune cell types using the TIMER and CIBERSORT algorithms. Table S6: GO-term of 106 IR-DEGs with STRING enrichment. Table S7: GO-term of 106 IR-DEGs with STRING enrichment. Figure S1: Representative microscopic features of high tumor budding (TB) histology.