Hypermethylation and Downregulation of UTP6 Are Associated With Stemness Properties, Chemoradiotherapy Resistance, and Prognosis in Rectal Cancer: A Co-expression Network Analysis

Background To identify the hub genes associated with chemoradiotherapy resistance in rectal cancer and explore the potential mechanism. Methods Weighted gene co-expression network analysis (WGCNA) was performed to identify the gene modules correlated with the chemoradiotherapy resistance of rectal cancer. Results The mRNA expression of 31 rectal cancer patients receiving preoperative chemoradiotherapy was described in our previous study. Through WGCNA, we demonstrated that the chemoradiotherapy resistance modules were enriched for translation, DNA replication, and the androgen receptor signaling pathway. Additionally, we identified and validated UTP6 as a new effective predictor for chemoradiotherapy sensitivity and a prognostic factor for the survival of colorectal cancer patients using our data and the GSE35452 dataset. Low UTP6 expression was correlated with significantly worse disease-free survival (DFS), overall survival (OS), and event- and relapse-free survival both in our data and the R2 Platform. Moreover, we verified the UTP6 expression in 125 locally advanced rectal cancer (LARC) patients samples by immunohistochemical analysis. The results demonstrated that low UTP6 expression was associated with worse DFS and OS by Kaplan-Meier and COX regression model analyses. Gene set enrichment and co-expression analyses showed that the mechanism of the UTP6-mediated chemoradiotherapy resistance may involve the regulation of FOXK2 expression by transcription factor pathways. Conclusion Low expression of the UTP6 was found to be associated with chemoradiotherapy resistance and the prognosis of colorectal cancer possibly via regulating FOXK2 expression by transcription factor pathways.

Methods: Weighted gene co-expression network analysis (WGCNA) was performed to identify the gene modules correlated with the chemoradiotherapy resistance of rectal cancer.
Results: The mRNA expression of 31 rectal cancer patients receiving preoperative chemoradiotherapy was described in our previous study. Through WGCNA, we demonstrated that the chemoradiotherapy resistance modules were enriched for translation, DNA replication, and the androgen receptor signaling pathway. Additionally, we identified and validated UTP6 as a new effective predictor for chemoradiotherapy sensitivity and a prognostic factor for the survival of colorectal cancer patients using our data and the GSE35452 dataset. Low UTP6 expression was correlated with significantly worse disease-free survival (DFS), overall survival (OS), and event-and relapse-free survival both in our data and the R2 Platform. Moreover, we verified the UTP6 expression in 125 locally advanced rectal cancer (LARC) patients samples by immunohistochemical analysis. The results demonstrated that low UTP6 expression was associated with worse DFS and OS by Kaplan-Meier and COX regression model analyses. Gene set enrichment and co-expression analyses showed that the mechanism of the UTP6-mediated chemoradiotherapy resistance may involve the regulation of FOXK2 expression by transcription factor pathways.

INTRODUCTION
Preoperative chemoradiotherapy (CRT) and radical surgery have become the standard of care for stage II/III rectal cancer patients (Heald et al., 1982). The benefits of this multimodality therapy have been well documented, including tumor downsizing and downstaging, increased radical resection, and reduced local recurrence (Sauer et al., 2004(Sauer et al., , 2012van Gijn et al., 2011). However, the treatment may increase the perioperative mortality rates by 48% (Rahbari et al., 2013). Moreover, 15% to 45% of the rectal cancers would are estimated to develop resistance to CRT, and might be exposed to CRT-related toxicitiestoxicity without oncological benefit (Ha et al., 2017). Therefore, the identification of valid biomarkers for colorectal cancer (CRC) with resistance to CRT has become imperative.
Currently, high-throughput sequencing is commonly used to screen and identify differentially expressed genes (DEGs). Conventional molecular biology methodology evaluates the variations in and functions of genes independently. Weighted gene co-expression network analysis (WGCNA) is an unbiased systematic biological approach, which clarifies the transcriptome function at a system-level, determines gene-gene correlations, and identifies gene modules with a high correlation across microarray data. It also bridges the gap between individual genes and tumorigenesis and progression (Zhang and Horvath, 2005;Langfelder and Horvath, 2008;Tian et al., 2017). WGCNA can facilitate the network-based gene screening approaches that screen key biomarkers associated with clinical traits in various cancers. However, this efficient bioinformatics approach has not yet been adopted to explore network-centric genes associated with CRT resistance in rectal cancer patients.
In this context, this study aimed to screen the relevant hub genes for CRT in our expression profile using WGCNA. Then, the hub genes were verified using testing data sets and patient tissue samples.

Subjects and Collection
In total, 31 locally advanced rectal cancer (LARC) patients receiving preoperative CRT between March 2016 to December 2016 in Fujian Medical University Union Hospital, China were enrolled as the internal dataset in this study, which was used for the gene-chip analysis. The CRT, treatment, and follow-up protocol were described in a previous study (Zhang et al., 2017). All patients received preoperative CRT and underwent R0 resection of primary colorectal tumors after neoadjuvant CRT. And the samples were collected at diagnosis and without any treatment before biopsy from colonoscopy. The chemoradiotherapy protocol was as follows: preoperative radiotherapy consisted of 45 Gy to the pelvis Abbreviations: WGCNA, weighted gene co-expression network analysis; CRT, chemoradiotherapy; CRC, colorectal cancer; LARC, local advanced rectal cancer; GEO, Gene Expression Omnibus; TCGA, the cancer genome atlas; PPI, proteinprotein interaction; GSEA, Gene set enrichment analysis; ROC, receiver operating characteristic curve.
for 5 weeks (180 cGy/25 fractions) and a tumor boost of 5.4 Gy. Concomitant chemotherapy was administered with oral capecitabine (825 mg/m 2 twice daily from day 1 to day 14 per cycle).
Moreover, to further validate the UTP6 expression in the LARC patients' colonoscopic samples before CRT. Consecutive LARC patients who underwent CRT and radical resection between 2011 and 2014 were identified. A total of 125 patients were enrolled in the present study. This study was approved by the Institutional Review Board of Fujian Medical University Union Hospital (2013051).

RNA Extraction, Quality Control, Labeling, Array Hybridization, and Data Analysis
Total RNA was extracted from the tissue mentioned previously using Trizol reagent (Invitrogen), according to the manufacturer's protocol. RNA quantity and quality were measured by a NanoDrop ND-1000 and RNA integrity was assessed by standard denatured agarose gel electrophoresis. Sample labeling and array hybridization were performed according to the human GeneChip Analysis protocol (Affymetrix, Santa Clara, CA). The scanned data obtained from each microarray were normalized to correct for small differences in the amounts of each of the cRNA probes applied to the microarray and were processed for signal values using Affymetrix software (LIMS 5.0). The differentially expressed mRNAs were identified through fold-change filtering. Hierarchical clustering was performed using Agilent Gene Spring GX software (version 11.5.1). The upregulated and downregulated mRNA groups were defined as the mRNA expression in the CRT-responsive group compared with the CRT-non-responsive group. The Gene Ontology (GO) functional analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis were performed using standard enrichment computation method.

Data Collection From the Gene Expression Omnibus and the Cancer Genome Atlas Database
Datasets were downloaded from the Gene Expression Omnibus (GEO) database 1 on December 19, 2019. The GSE35452 dataset included mRNA expression gene microarray of 46 rectal cancer patients receiving preoperative CRT and biopsy specimens were collected before preoperative CRT, which was used as the external validation dataset. Cancer stem cells are regarded as having self-renewal and differentiation properties and are one of the integral factors mediating the response to CRT (Baumann et al., 2008;Wahab et al., 2017). GSE14773 and GSE24747 were used to determine the hub gene expression in CRC stem cells. The GSE14773 dataset included one cell line with markers of CRC stem cells, high CD44, and CD166 expression, and the other one was a parental control CRC cell. The GSE24747 dataset included two groups, CD133 + and CD133-, CACO-2 cells. DNA methylation is considered an important component in gene expression regulation, and high levels of DNA promoter region methylation result in transcriptional silencing. Additionally, promoter regions methylation is reported to be associated with the CRT-resistance and considered as a predictor for CRT (Wilting and Dannenberg, 2012;Ivanova et al., 2013). The GSE104271 dataset was utilized to explore whether the hub genes were regulated by DNA methylation. As described by Huang et al. (2018) previously, FIBP knockdown cell lines were defined as cancer cells compared to cancer stem cells. The methylation data of the CRT-resistant genes CRC and normal colon tissues were obtained from The Cancer Genome Atlas (TCGA) in the UALCAN database 2 (Chandrashekar et al., 2017). A flow diagram of the present study is shown in Figure 1.

Definitions
Tumor response to CRT was graded according to the American Joint Committee on Cancer pathological tumor regression grade 2 http://ualcan.path.uab.edu/ (AJCC TRG) (Amin et al., 2017;Benson et al., 2018); that is, TRG 0, no residual tumor cells; TRG 1, single cells or small groups of cells; TRG 2, residual cancer with the desmoplastic response; and TRG 3, minimal evidence of tumor response. Pathological complete response (pCR) was defined as the absence of viable tumor cells in the resected specimen, either at the primary site or in the lymph nodes. Venous blood samples were obtained within 1 week before CRT.

Co-expression Network Construction
The WGCNA algorithm was previously described in detail (Zhang and Horvath, 2005). Briefly, first, we identified the qualification of the profile data. The co-expression network was constructed using the "WGCNA" package in R software (Horvath and Dong, 2008;Mason et al., 2009). Next, the correlation matrix was established and the soft threshold power was determined by analysis of the network topology. Then, the topological overlap matrix (TOM) was established (Yip and Horvath, 2007;Ivanova et al., 2013;Botía et al., 2017). Based on the phenotypic data of the groups, we calculated each module p-value by the t-test gene significance.

CRT-Resistance Modules and Hub Gene Identification
To explore the relevant modules, we examined the association between module eigengenes (MEs) and CRT-resistance using Pearson's correlation analysis. To identify the hub genes, the CRT-resistance model with the highest correlation coefficient (P < 0.05) in the data set was chosen, and the module also had the highest specific weight of all of the modules. The hub genes in the module were defined by module connectivity as measured by thean absolute value of the >0.3 in Pearson's correlation >0.3 analysis. The protein-protein interaction (PPI) network was constructed by all genes in the CRT-resistance module, nd PPI network analysis was performed to screen the hub genes by Cytoscape. And the top 10 values were included for further analysis.
Statistical analysis was performed using SPSS software (version. 23 SPSS Inc., Chicago, IL, United States) and R software (version. 3.4.1). The continuous variables were reported as means and standard deviation from the analysis of variance test. The survival outcomes were assessed using Kaplan-Meier and logrank analyses. Receiver operating characteristic (ROC) curve analysis was performed. Finally, the prognostic significance of the hub genes in CRC patients was analyzed using the R2: Genomics Analysis and Visualization Platform, 3 on January, 03, 2020. The optimal cutoff points for the expression of UTP6 and FTSJ3 were calculated and determined using Cut-off Finder 4 Visit and download data, 2019.12.26). Budczies et al. (2012), a new bio-informatics tool for biomarker assessment and outcomebased cut-point optimization, which identified the cut-off with the minimum p-values from log-rank χ2 statistics in terms of disease-free survival (DFS). A P-value of <0.05 was considered statistically significant.

Gene Set Enrichment Analysis and Co-expression Gene Analysis
To explore the potential function of UTP6 in LARC patients, gene set enrichment analysis (GSEA) was performed in the patients in our previous datasets. A P-value of <0.05 and an enrichment score (ES) of >0.3 were set as the cutoff criteria. To further explore the correlation between UTP6 and relevant genes in the GSEA results, we analyzed the co-expressed genes in the TCGA data in the UALCAN dataset (Visit and download data, 2020.02.01).

Immunohistochemical Analysis of UTP6 in the LARC Patients' Sample
The protein expression of UTP6 in specimens obtained before and after CRT in 125 LARC patients was assessed using the immunohistochemical streptavidin-biotin complex method (Zhang et al., 2018). Phosphate-buffered saline (PBS) was used as the negative control and the image of the positive control from GE Healthcare Life Sciences. Immunoreactivity was scored by semi-quantitative analysis, and the fields were randomly selected in five directions (up, center, down, left, and right) under high magnification ( * 400). The color was determined based on the intensity score as follows: 0 (no staining), 1 (light yellow), 2 (brown), and 3 (deep brown). The percentage of positive cells was scored as 0 (<5%), 1 (5-25%), 2 (25-50%), 3 (50-75%), and 4 (>75%). The mean value was calculated for each case with the aforementioned scoring methods and the final score was obtained by multiplying these two scores. The score between 0 and 4 was defined as the low expression and >4 was defined as high expression. All analyses were performed in a doubleblind manner.

Cluster Analysis
A gene chip array was used to examine the gene expression profiles in primary tumor cells. A supervised hierarchical cluster analysis of the gene expression profiling data showed that the two groups had a clustering trend (Figures 2A,B). A total of 18419 genes were detected in the gene chip array. The SAM for the DEGs revealed that the two groups significantly differed in genes related to tumor cell biology, including 798 upregulated genes and 450 downregulated genes.

GO Enrichment and KEGG Analysis
GO enrichment analysis was performed to investigate the molecular mechanism of the differently expressed genes involved in the resistance to CRT in LARC patients. We detected the top significant GO-enriched terms and KEGG terms in both the significantly upregulated and downregulated genes in rectal cancer patients (Figures 2C-F). The results showed that the top three significant GO terms for the upregulated genes were related to pancreatic A cell differentiation, muscle contraction, and type B pancreatic cell differentiation involved in the immune response. Among the downregulated genes, the top three significant GO terms were related to the ER-associated ubiquitindependent protein catabolic process, mRNA processing, and RNA splicing. Moreover, the top three significant KEGG terms in the upregulated genes were related to neuroactive ligand-receptor interaction, lysine degradation, and olfactory transduction involved in the immune response. In the downregulated genes, the top three significant GO terms were related to protein processing in the endoplasmic reticulum, amino sugar and nucleotide sugar metabolism, and the proteasome.

Development of Weighted Co-expression Network and Identification of Key Modules
To identify the hub genes, a weighted co-expression network was utilized to analyze our data ( Figure 3A). A total of 31 modules were identified, as shown in Figure 3C. We further analyzed the relationship between CRT-resistance and the modules ( Figure 3C). Among these modules, the module shown in dark turquoise ( Figure 3C) had the strongest negative association with CRT-resistance (r = −0.45, P = 0.01), while the light cyan module showed the highest positive association with CRT-resistance (r = 0.25, P = 0.18). Then, we selected the dark turquoise module as the hub module. GO and KEGG enrichment analyses were conducted to explore the functions of the dark turquoise module. The results demonstrated that the dark turquoise module was mainly enriched for translation, DNA replication, and the androgen receptor signaling pathway ( Figure 3B).

Hub Gene Identification
Through co-expression analysis, 92 genes in the dark turquoise module were considered genes with high module connectivity. Then, the genes were analyzed by the PPI network and 59 genes in the dark turquoise module were identified as hub genes in the co-expression network. Ten hub genes were analyzed by the degree of PPI and the correlation between clinical traits and module connectivity was analyzed? (Figure 3D and Table 1). Finally, we chose the most associated genes, UTP6 and FTSJ3, as the actual hub genes.

Hub Gene Validation
To validate the hub genes, we examined the expression of UTP6 and FTSJ3 in rectal cancer tissues of CRT-resistant and CRT-sensitive cases in our datasets. In the internal testing data sets, the relative expression of UTP6 (6.46 ± 0.40 vs. 6.07 ± 0.40, P = 0.02; Figure 4A) and FTSJ3 (6.21 ± 0.32 vs. 5.75 ± 0.47, P < 0.01, Supplementary Figure 1A) was significantly increased in the CRT-sensitive tissues. The ROC curve demonstrated that both UTP6 (P = 0.02, AUC = 0.76, Figure 4D) and FTSJ3 (P < 0.01, AUC = 0.79, Supplementary Figure 1D) could efficiently discriminate CRT-resistant from CRT-sensitive rectal cancer cases. Moreover, Pearson's analysis was conducted to determine whether the hub gene expression was associated with the TRG grade. The results demonstrated that UTP6 (r = −0.35, P = 0.02, Figure 4C) and FTSJ3 (r = −0.37, P = 0.04, Supplementary Figure 1C) had significant associations with the TRG grade. We further explored the prognostic impact of UTP6 and FTSJ3 on the survival of LARC patients. Cut-off Finder was used to identify the cutoff values for UTP6 and FTSJ3. Cut-off Finder identified 6.18 and 5.97 as the cutoff values for UTP6 and FTSJ3 expression, respectively (Supplementary Figure 2). After a median follow-up period of 42 months (range, 4 -47 months), low expressions of UTP6 and FTSJ3 were associated with significantly worse disease-free survival (DFS) compared to high expression (UTP6 low expression VS. high expression: 50.0% vs. 100.0%, P < 0.001; FTSJ3 low expression VS. high expression: 54.5% vs. 100.0%, P < 0.001), as shown in Figure 4F and Supplementary Figure 1F. Similarly, a high expression of UTP6 and FTSJ3 were correlated with better overall survival (OS) compared to low expression (UTP6 high expression VS. low expression: 100.0%, vs. 50.0%, P < 0.001; FTSJ3 high expression VS. low expression: 100.0% vs. 57.7%, P < 0.001), as demonstrated in Figure 4G and Supplementary Figure 1G. (C) Pearson's correlation analysis between the AJCC TRG grade and UTP6 expression. ROC curves and AUC statistics to evaluate the predictive efficiency of the UTP6 in our data and external data to distinguish CRT-resistance from CRT-sensitive CRC cases from (D) our data and (E) GSE35452. The disease-free survival (F) and overall survival (G) between low and high expression of UTP6. (H) UTP6 expression in the pCR and non-pCR group in the LARC patients' sample. (I) ROC curves and AUC statistics to evaluate the predictive efficiency of the UTP6 in LARC patients' samples to distinguish pCR from non-pCR LARC cases. The overall survival (J) and disease-free survival (K) between low and high expression of UTP6 in LARC patients' sample.

Hub Gene Validation in External Data
To validate the hub genes in our data set, we examined the expression level of UTP6 and FTSJ3 by comparing the rectal cancer tissues of CRT-resistant and-sensitive cases using an external dataset. In the external GSE35452 database, UTP6 expression was significantly higher in the CRT-sensitive tissues compared to the CRT-resistant tissues (0.45 ± 1.51 vs. −0.47 ± 1.45, P = 0.042; Figure 4B), but the expression of FTSJ3 was not statistically different between the two groups (0.11 ± 0.40 vs. −0.05 ± 0.53, P = 0.246; Supplementary Figure 1B). Thus, we chose UTP6 as the "real" hub gene for further analysis. ROC analysis revealed that UTP6 could efficiently discriminate CRTresistant from CRT-sensitive cases (P = 0.045, AUC = 0.672; Figure 3E). However, FTSJ3 could not efficiently distinguish CRT-resistant from CRT-sensitive cases (P = 0.12, AUC = 0.63; Supplementary Figure 1D).
FIGURE 5 | High UTP6 expression was associated with a better event-, disease-, and relapse-free survival. (A-H) High UTP6 expression was associated with a significantly better event-, disease-and relapse-free survival (both P < 0.05).

UTP6 Validation in LARC Patients' Sample
To further verify UTP6 expression in the tumor samples, colonoscopy samples were collected from LARC patients before CRT. The immunohistochemical analysis demonstrated that higher UTP6 scores were associated with better prognosis in LARC patients following CRT. The 3-year OS rate of the low-UTP6 group was significantly lower than that in the high-UTP6 group (60.5% vs. 89.5%; P < 0.01, Figure 4J). Lower UTP6 scores were correlated with better DFS (Figure 4K). The 3-year DFS rate for the low-UTP6 group was significantly higher than that in the high-UTP6 group (42.4% vs. 87.1%; P < 0.01). Moreover, we compared the UTP6 expression in the pCR and non-pCR group. The results revealed that the relative expression of UTP6 was significantly increased in the pCR group (P < 0.01; Figure 4H). The ROC curve demonstrated that UTP6 could efficiently discriminate pCR from non-pCR LARC cases (P < 0.001, AUC = 0.77, Figure 4I). In addition, we detected the UTP6 expression in the rectal cancer tissues and adjuvant cancer tissues, the result demonstrated that UTP6 expression value was similar in rectal cancer and adjacent cancer tissues (P = 0.571, Supplementary Figure 3).

Association of UTP6 Expression With Patient Characteristics and Perioperative Clinicopathological Parameters in the LARC Patients Following CRT
No significant differences were observed between UTP6 low expression group and UTP6 high expression group in terms of gender, age, American Society of Anesthesiologists (ASA) grade, interval time between CRT and surgery, distance from the anal verge, clinical T stage, clinical N stage, pre-CRT CEA level, and pre-CRT CA199 level, as shown in Table 4.
No significant differences were observed between UTP6 low expression group and UTP6 high expression group in terms of pathological type, postoperative complication, circumferential resection margin (CRM) involvement, tumor differentiation, perineural invasion, vascular invasion, and organ preservation procedure (Table 5). Compared to the UTP6 low expression group, the UTP6 high expression group was associated with an increased metastasis to the lymph nodes (0.45 ± 1.4 vs. 3.1 ± 7.1, P < 0.001), pathological T stage, pathological N stage (all P < 0.01), and poorer TRG grade (P < 0.01). In the training dataset, in the high-risk score group a larger tumor size was seen (2.7 ± 1.2 vs. 3.5 ± 1.6, P = 0.002), more lymph nodes were retrieved (11.2 ± 8.7 vs. 15.7 ± 13.8, P = 0.028), and poorer histopathology was observed (P = 0.026) compared with UTP6 low expression group.

Analysis DNA Promoter Region Methylation of UTP6 in CRC Tissues and Expression, Promoter Methylation of UTP6 in CRC Stem Cells
To explore whether UTP6 was modulated by methylation, we evaluated data on the methylation of the UTP6 promoter region in CRC and normal colon tissues. Based on the ULACAN database, methylation of the UTP6 promoter was decreased in normal colon tissue compared to CRC tissue (0.0409 ± 0.0013 vs. 0.0441 ± 0.0005, P = 0.017, Figure 6A). Further, the expression and promoter methylation levels of UTP6 were analyzed in CRC stem and non-stem cells. The relative UTP6 expression was significantly reduced in CRC stem cells compared to CRC nonstem cells (GSE14773: 10.33 ± 0.13 vs. 10.80 ± 0.14; P = 0.04; GSE24747: 3.44 ± 0.09 and 3.81 ± 0.00, P = 0.03; Figures 6C,D). We analyzed promoter methylation in CRC stem cells in the GSE104271 dataset. The results demonstrated that the UTP6 promoter was hypermethylated at cg10893370 (stem cell vs. non-stem cell: 0.06 ± 0.00 vs. 0.09 ± 0.00, P = 0.008, Figure 6B) and cg13453082 (0.05 ± 0.01 vs. 0.08 ± 0.00, P = 0.03, Figure 6B) sites in CRC stem cells.

GSEA and Co-expression Gene Analysis
GSEA was conducted to investigate the potential mechanism of the UTP6-mediated CRT-resistance in CRC. Our data demonstrated that the negatively correlated KEGG pathways were enriched for the ABC transporter signaling pathway (Figures 7A-C). The positively correlated pathways included transcription factor pathways and cell cycle signaling pathways.
To further explore the mechanism of the UTP6-mediated CRT response, the genes co-expressed with UTP6 in the TCGA colon tumor data from ULACAN were analyzed. The results demonstrated that FOXK2 was significantly correlated with UTP6 (r = 0.55, P < 0.01; Figure 7D).

DISCUSSION
To date, reliable molecular markers for CRT-resistance in LARC patients are still unavailable. In this study, a gene chip array was performed to detect the gene expression in LARC patients who received CRT. Then WGCNA, an advanced methodology for multigene analysis, was conducted to identify the gene coexpression modules associated with CRT resistance. UTP6 was identified and validated as a hub gene correlated with CRT resistance. High UTP6 expression was correlated with better survival in CRC patients. UTP6 was hypermethylated in CRC tissues, especially in the CRC stem cell subpopulation. The underlying mechanism of UTP6 in maintaining CRC stemness might involve transcription factor pathways, especially FOXK2. Currently, microarray expression profiling has been utilized to screen biomarkers in patients with CRT-resistant rectal cancer. Conventional molecular biology methodology identifies DEGs, but it is difficult to correlate biological information provided by gene names with biological functions independently. WGCNA has emerged as an effective method to discover the relationship between networks, genes, phenotypes, and biological sample information, thus avoiding the shortcomings of the traditional methods (Bakhtiarizadeh et al., 2018;Gao et al., 2016;Magani et al., 2018).
Herein, we performed WGCNA to identify the actual hub gene in our patient cohort. The results demonstrated that the dark turquoise module was chosen as the hub module and enriched for translation, DNA replication, and androgen receptor signaling pathways. Those pathways are involved in multi-drug resistant CRT tumors, indicating that the dark turquoise module was significantly associated with CRT-resistance. As described previously (Liu et al., 2017;Wang et al., 2018), PPI network analysis and clinical trait correlation were performed to screen for the actual hub gene. The results suggested that the UTP6 gene had the highest correlation coefficient. We further analyzed the expression of UTP6 in our previous data and testing data set (GSE35452), which showed that UTP6 was the relevant gene for CRT-sensitivity. ROC analysis showed that UTP6 could effectively discriminate between CRT-resistant and -sensitive cases. To our knowledge, this was the first study to identify and verify UTP6 as an effective new marker for the prediction of CRT response.
Having shown that UTP6 was associated with CRT-resistance and survival in CRC patients, the underlying mechanism remained unclear. Cancer stem cells are widely known as tumorinitiating cells due to stem cell properties. Increasing evidence has shown that cancer stem cells are involved in the resistance to conventional cytotoxic therapies, radio-and chemotherapy (Baumann et al., 2008;Wahab et al., 2017). Herein, we evaluated the relative UTP6 expression in two independent data sets including CRC stem cells and parental cells (GSE14773 and GSE24747). The results demonstrated that the expression of UTP6 was significantly decreased in CRC stem cells compared to CRC non-stem cells.
DNA methylation, a main epigenetic modification in the mammalian genome (Baylin et al., 2001;Herman and Baylin, 2003), may inactivate genes and affect the CRT response (Wilting and Dannenberg, 2012;Ivanova et al., 2013). Several studies have reported that DNA methylation could promote drug or radiation resistance (Su et al., 2011), and thus, could be used to predict treatment response. In the present study, we found hypermethylation in colon cancer tissue and the cg10893370 and cg13453082 sites located in the UTP6 promoter in CRC stem cells. These results indicated that the hypermethylation of UTP6 may result in a decreased expression of UTP6 in CRC stem cells and that the downregulation UTP6 expression in cancer stem cells may be due to hypermethylation of the promoter sites. Several studies have reported that DNA methylation could affect gene expression in the transcriptome, which is involved in CRT-resistance in various cancer stem cells (Kurscheid et al., 2015;Kostopoulou et al., 2018). Importantly, DNA methylation can be reversed by anti-tumor drugs, such as decitabine (Ariazi et al., 2017;Flotho et al., 2018). The findings suggest that the combination of demethylation therapy targeting UTP6 and cancer stem cells during CRT might be beneficial to rectal cancer patients. UTP6 (also known as HCA66) is an essential component of the UTPB, a large complex composed of six proteins Utp1/Pwp2, Utp6, Utp12/Dip2, Utp13, Utp18, and Utp21 (Grandi et al., 2002;Krogan et al., 2004). Ferraro et al. (2011) reported that Apaf1 could bind with UTP6 and participate in the regulation of centrosomal microtubule nucleation, spindle assembly, cell migration, and mitochondrial network organization. Meanwhile, UTP6 has been proposed as a promising prognostic indicator in bladder cancer (Shivakumar et al., 2017). A previous study reported that UTPB was the consist of the UTP6, which is an important complex in regulating ribosome synthesis (Hunziker et al., 2016;Zhang et al., 2016). Herein, we found that UTP6 was not only related to the ABC transporter signal pathway, but also transcription factor pathways, which are associated with CRT resistance (Shinto et al., 2014;McCoy et al., 2015). Pearson's correlation analysis showed that the expression of UTP6 was correlated with that of FOXK2, which is a transcription factor known as an anti-oncogene (Shan et al., 2016;de Moraes et al., 2019). The results indicated that the anti-chemoradiotherapy resistance effect of UTP6 may be related to FOXK2 expression.
Unfortunately, the precise mechanism of UTP6 involved in CRT sensitivity in rectal cancer patients is still unclear. There were several limitations to our study. Although the underlying mechanism of UTP6 was explored by bioinformatics analysis, the results warrant further confirmation by in vitro or in vivo experiments. To our knowledge, this was the first study to explore CRT-resistance-related genes through WGCNA using a patient cohort with complete clinicopathological and follow-up data. The present study provided insight into the role of UTP6 in CRT-sensitive in rectal cancer patients.

CONCLUSION
Through WGCNA, UTP6 was identified and validated as a hub gene and a valid novel predictor for CRT-sensitivity. High expression of UTP6 associated with better survival in LARC patients. UTP6 was hypermethylated in CRC tissues, especially in CRC stem cells subpopulation. GSEA demonstrated that the mechanism underlying UTP6 maintains CRC stemness might be involved in the transcription factors pathway especially FOXK2. Further research is needed to confirm the above findings.

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: https://www.ncbi.nlm. nih.gov/geo/, GSE145037.

ETHICS STATEMENT
Studies relative to humans in this article were approved by the Ethics Committee of The Fujian Medical University Union Hospital. The patients/participants provided their written informed consent to participate in this study. Written informed consent was obtained from the individual(s) for the publication of any potentially identifiable images or data included in this article.

AUTHOR CONTRIBUTIONS
YZ, QG, YW, GG, and XL designed and performed the experiments, analyzed the data, and wrote the manuscript. YP, JZ, and YY performed the experiments. All authors read and approved the final manuscript.