Hallmarks of glycogene expression and glycosylation pathways in squamous and adenocarcinoma cervical cancer

Patient samples

The Ethics Committee approved the present study of the Instituto Mexicano del Seguro Social (R-2012-785-061). All patients signed an informed consent form according to the guidelines of the Human Ethics Committee. Biopsies were obtained at the Hospital de Especialidades, General Manuel Ávila Camacho, Instituto Mexicano del Seguro Social from women between the ages of 20 and 73. For the CC group, inclusion criteria included adult women with SC diagnosis with no previous treatment (mean age 42). In contrast, for the normal cervix group, samples were obtained from adult women with uterine myomatosis diagnosis (mean age 35). A total of 23 CC samples and 15 biopsies of normal cervical tissue were included in the analysis.

Microarray assay

Normal cervix tissues and CC samples were collected and maintained in RNAlater solution (Qiagen, Hilden, Germany) at −80 °C. Total RNA was extracted using the RNeasy Plus system (Qiagen, Hilden, Germany) following the manufacturer’s protocol. RNA samples were pooled into two groups (normal and CC). Microarray analysis was performed at the Microarray Unit of the Instituto de Fisiología Celular, Universidad Nacional Autónoma de México (UNAM). A total of ten micrograms of total RNA were used to synthesize labeled cDNA employing the First-Strand cDNA labeling kit (Invitrogen, Waltham, MA, USA), incorporating dUTP-Alexa555 or dUTP-Alexa647. Equal quantities of labeled cDNA were hybridized using IniHyb hybridization solution (TeleChem International INC., Los Altos, CA, USA). The microarray included 10,000 genes from the NCB and GenBank databases. Acquisition and quantification of array images were performed using a Scan Array 4,000 instrument and ScanArray 4,000 software (Packard BioChips,Waltham, MA, USA). The Alexa555 and Alexa647 density and background mean values were calculated with ArrayPro Analyzer software from Media Cybernetics. GenArise1.38.0 software developed in the Computing Unit of the Cellular Physiology Institute of UNAM was used for the data analysis and statistics. Analyzed data were submitted to the NCBI-Gene Expression Omnibus database (accession number GSE159976).

Glycogene expression analysis

We selected glycogenes with a z-score ≥ 2 to identify the glycogenes with expression changes between the CC and normal cervix. A total of 401 glycogenes reported to date were used for glycogene expression analysis; these were obtained from the published reports of Venkitachalam et al. (2016), and Aco-Tlachi et al. (2018). Glycogenes included glycosyltransferases and glycosidases, sugar-nucleotide synthases, sugar-nucleotide transporters, chaperones, and some genes related with energetic metabolism.

Protein-protein enrichment analysis

Glycogenes identified as altered were submitted to the STRING database version 10.5 (http://string-db.org/) to evaluate predicted protein-protein interactions (Szklarczyk et al., 2019). We considered the following settings for the analysis: text mining, experiments, databases, coexpression, neighborhood, gene fusion, and cooccurrence as interaction sources, with no more than five interactors and a minimum interaction score of 0.4 as the confidence level. Biological process (GO) and KEGG pathways were chosen as functional enrichments of the network. The results exclusively included a protein-protein enrichment p-value of at least ≤ 0.05 and an FDR of at least ≤ 0.05.

RNA-seq data analysis

RNA-seq data of 401 glycogenes from TCGA and Firehose Legacy were obtained from the cBioPortal for Cancer Genomics (Gao et al., 2013). We analyzed glycogene expression considering the mRNA expression of genes with a z-score threshold ≥ 2.0 relative to all samples (log RNA Seq V2 RSEM). The population included women from 20 to 88 with a diagnosis of SCC and AC, four samples of cervical adenosquamous carcinoma, and all the group stages from I to IVB (AJCC staging system).

Hierarchical clustering

Clinical data of patients obtained from the cBioPortal for Cancer Genomics and their respective RNA-seq data of glycogene expression were organized and submitted in the online free software Morpheus (https://software.broadinstitute.org/morpheus/) to perform the heat map analysis and hierarchical clustering. One minus Pearson correlation was used as the metric, and the average was used as the linkage method. Venn diagram of glycogenes clusters were obtained using the tool http://bioinformatics.psb.ugent.be/webtools/Venn/.

Pathway enrichment analysis

Altered glycogenes were analyzed using DAVID (Database for Annotation, Visualization, and Integrated Discovery) software (http://david.abcc.ncifcrf.gov/). Based on gene ontology (GO), genes were classified according to their function in the KEGG pathway with at least a p-value < 0.05 and FDR < 0.05.

Comparison of glycogene expression between cervical cancer and healthy cervical tissue

We performed a microarray expression assay to identify changes in glycogene expression in CC compared to a normal cervix. From the 178 glycogenes contained in the microarray, results indicated that 15 glycogenes were altered in CC (z-score threshold ≥ 2.0). The glycogenes PIGC, PIGN, SMPD3, CHI3L1, HEXB, and ST6GAL1 were upregulated (Table S1-column A), and GYS1, CHST12, CTBS, HAS3, HPSE, HYAL2, GALNT11, UGT2B4, and UGT2B28 were downregulated (Table S2-column A). Bibliographic research indicates that most altered glycogenes are reported to have expression changes in other types of cancer (Tables 1 and 2).

Then, to further map the 15 altered genes onto the global glycogene network, we performed a protein-protein interaction (PPI) network functional analysis in the STRING database (Szklarczyk et al., 2019) using the 401 glycogenes (Fig. S1A). Results showed a complex interaction among the 401 glycogenes, since some of them were involved in several glycosylation pathways (Figs. S1B, S1C, S1F). Thus, we refine the search for identification of specific glycosylation pathways by performing a PPI analysis, including only the upregulated and downregulated glycogenes. The glycogenes SMPD3, CHI3L1, HEXB, and ST6GAL1, were involved in the immune system response (p-value ≤ 0.005), while PIGC and PIGN participate in the synthesis of GPI anchors (Fig. 1A). In contrast, results revealed that HYAL2 and HAS3 were involved in the hyaluronan (HA) metabolic process. To further identify putative targets and additional glycosylation pathways not displayed in the previous analysis, we included the upregulated or downregulated glycogenes and five more interactors in the subsequent analysis. The results again showed that PIGC and PIGN were involved in synthesizing GPI anchors (Fig. 1C). In contrast, CHI3L1 and HEXB, together with FUCA2 and HEXA, were engaged in neutrophil degranulation (p-value ≤ 0.001) (Fig. 1C). Concerning the downregulated glycogenes, UGT2B4 and UGT2B28, along with UGP2, were implicated in pentose and glucuronate interconversions (p-value ≤ 0.001) (Fig. 1D). Other results showed that CHST12 could interact with CHST3 and CHST7 in chondroitin sulfate (CS)/dermatan sulfate synthesis (p-value ≤ 0.001) (Fig. 1D).

Comparison between adenocarcinoma and squamous cervical

Next, to identify glycogene signatures in CC, we analyzed the expression of 401 glycogenes in 297 samples using the RNA-seq public data obtained from the cBioPortal for Cancer Genomics using the samples that corresponded to the histological types of SC and AC (Table 3). Unsupervised hierarchal clustering analysis revealed that 25 genes exhibited low expression in almost all samples compared with the rest of the glycogenes (Table S2-column B), suggesting a hallmark of glycogene expression among all CC types.

Then, to test whether AC and SCC could display distinctive glycogene expression, we compared the gene expression between them. The results showed differences in glycogene expression. Samples of SCC displayed a heterogeneous pattern of gene expression (Fig. 2A). In contrast, AC displayed characteristic glycogene expression that included 107 glycogenes grouped into three different clusters: AC-A, AC-B, and AC-C (Fig. 2A). Cluster AC-A included 51 genes with a high expression (Fig. 2A; Table S1-column B), while AC-B and AC-C displayed low gene expression of 23 and 33 glycogenes (Table S2-column C and column D), respectively, compared with the rest of the CC population (Fig. 2A). The results of the functional analysis revealed that glycogenes in the AC-A cluster participate in N-glycan biosynthesis, mucin-type O-glycan biosynthesis, glycosphingolipid biosynthesis (lacto and neolacto series), GAG biosynthesis keratan sulfate (KS), glycosphingolipid biosynthesis (globo and ganglio series), other glycan degradation, and other pathways, such as amino sugar and nucleotide sugar metabolism, starch and sucrose metabolism, and carbohydrate digestion and absorption (Fig. 2B). In comparison, genes with a low expression in AC-B are involved in GAG biosynthesis (CS, dermatan sulfate (DS)) and heparan sulfate/heparin (HS), mucin-type O-glycan biosynthesis, and metabolic pathways (Fig. 2C). Similarly, glycogenes in cluster AC-C were implicated in metabolic pathways and biosynthesis of CS, DS, HS but also GPI-anchor biosynthesis, glycosphingolipid biosynthesis (lacto and neolacto series), galactose metabolism, GAG degradation, and other pathways, such as amino sugar and nucleotide sugar metabolism, starch and sucrose metabolism, and insulin signaling pathways (Fig. 2D).

We further compared the glycogene expression of adenosquamous carcinoma with AC and SCC since the histological type includes both adenoid and squamous cells from the origin.

Identification of hallmark glycogenes in squamous cervical cancer

To further identify a glycogene expression hallmark in SCC, we removed the 20 glycogenes with low expression in all the CC samples and performed an unsupervised hierarchal clustering analysis exclusively in SCC samples (n = 253). The results indicate that almost all SCC samples exhibit a low expression of CHIA, GLT6D1, SPACA3, FUT9, B3GALT1, GALNTL5, NEU2, UGT2B15, UGT2A3, UGT2B7, and LYZL4 (Table S2-column E). Moreover, the results showed that two clusters of samples, SCC-1 and SCC-2, display specific glycogene expression (Fig. 3A). SCC-1 comprised 29 samples that showed a cluster of 30 glycogenes with high expression (SCC-1A; Table S1-column C) and a cluster of 19 genes with low expression (SCC-1B) (Fig. 3 and Table S2-column F). In contrast, SCC-2 formed by 36 samples exhibited a cluster of 31 genes with low expression (Fig. 3A and Table S2-column G).

SCC-1 included samples with a diagnosis from stage I to stage III, but most of the samples were stage I, and patients displayed (Table 4) an 89% survival rate. Samples included moderately differentiated (G2) and poorly differentiated (G3) tumors (Table 4). Enrichment analysis showed that genes were impl1icated in N-glycan biosynthesis, other types of O-glycan biosynthesis, GPI-anchor biosynthesis, lysosomes, GAG degradation, biosynthesis of CS/DS and HS, and other processes, such as processes in the endoplasmic reticulum and amino sugar and nucleotide sugar metabolism (Fig. 3C). Glycogenes in cluster SCC-1B were involved in N-glycan biosynthesis, CHS/DS biosynthesis, sphingolipid metabolism, and protein processing in the endoplasmic reticulum (Fig. 3C).

Regarding cluster SCC-2 and the clinicopathological characteristics, most of the samples corresponded to stage I; however, two samples with a diagnosis of stage IV and comprised of differentiation grades from G1 (well-differentiated) to G3 were also included. Overall survival analysis showed that patients displayed an 88% survival rate (Table 4). The glycogenes were involved in CS/DS biosynthesis, GAG degradation, other glycan degradation, and glycosphingolipid biosynthesis (ganglio, globo, lacto and neolacto series), among others (Fig. 3C). These results indicate that at least two types of SCC samples display a hallmark of glycogene expression. However, we cannot conclude a specific association with tumor characteristics or patient survival because not all the samples contained complete information.

Then, we evaluated individually whether cancer stage, molecular classification, tumor differentiation grade, or overall survival status were associated with a glycogene signature; however, the results did not show any glycogene expression pattern. These data indicate that glycogene expression is not individually associated with any of the abovementioned clinical features. To refine the analysis, we exclusively analyzed SCC samples that contained complete information: cancer stage, molecular classification, tumor differentiation grade, and overall survival status (n = 158). The results showed that from 158 samples, 67 were organized into four clusters that displayed a specific glycogene expression pattern (Fig. 4A). The first cluster included 25 patient samples (SCC-3); the second cluster included 23 samples (SCC-4); the third cluster included ten samples (SCC-5); and the fourth cluster included nine samples (SCC-6) (Fig. 4A). Analysis of the clinical characteristics of clustered patient samples (Table 4) showed that SCC-3 tumors were stage I (56%), with a nonkeratinizing molecular classification (68%), a G3 differentiation grade (52%), and a 68% survival rate. Cluster SCC-5 was formed by samples from stage I (90%), with nonkeratinizing tumors in 70%, a G3 differentiation grade in 80% of the cases, and patients displayed a 50% survival rate. Concerning SCC-6, the cluster comprised only nonkeratinizing tumors with stages I to IV and an 88% survival rate (Table 4).

SCC-3 displayed a high expression of 6 glycogenes related to HS/heparin biosynthesis (Fig. 3B; Table S1-column D; Fig. S2C). SCC-5 showed a high expression of 9 glycogenes implicated in CS/DS biosynthesis and glycosphingolipid biosynthesis (lacto and neolacto series) (Fig. 4C; Table S1-column F; Fig. S2C). In contrast, SCC-6 exhibited a high expression of 42 glycogenes, some of which were involved in glycosphingolipid biosynthesis (ganglio and globo series) (Fig. 4E; Table S1-column G; Fig. S2C). Regarding SCC-4, the results showed two clusters of glycogenes, SCC-4A and SCC-4B, that corresponded to the same glycogenes in clusters SCC-1A and SCC-1B (Tables S1-column E and S2-column F; Fig. S2A). A characteristic of SCC-4 was the low expression of genes involved in glycosphingolipid biosynthesis (lacto and neolacto series) and the high expression of genes involved in GAG degradation (Fig. 4D). In summary, all the data suggest that there are at least five types of SCC with distinguishable glycogene expression, where SCC-4 probably corresponds to SCC-1 since they share some samples and glycogenes (Table 5; Fig. S2B).

Cervical cancer displays upregulation of glycogenes implicated in GPI synthesis and downregulation of genes related to hyaluronan metabolism compared to healthy tissue

The human genome contains approximately 700 gene-encoding enzymes, including transporters and chaperones required for the cellular glycosylation machinery, glycan modifications, and degradation (Schjoldager et al., 2020). In this study, we analyzed the expression of 401 glycogenes, previously reported in Venkitachalam et al. (2016) and Aco-Tlachi et al. (2018), to obtain the glycome scenario in CC. We first compared the gene expression between CC and healthy tissue, and then we reached the glycogene expression in AC and SCC. The first analysis allowed us to understand the general glycosylation changes during cell transformation. The second analysis allowed us to identify the glycogene signatures in AC and SCC.

A comparison between CC and normal tissue indicated that CC displays upregulation of glycogenes related to GPI synthesis and downregulation of genes associated with the metabolism of HA. Regarding GPI-anchor biosynthesis, glycolipids act as anchors to specific cell surface proteins called GPI-anchored proteins. In breast, bladder, and gastric cancer, the high expression of other GPI-anchor biosynthesis glycogenes, such as PIGU, PIGT, and PIGX, is associated with oncogenesis, poor prognosis, and tumorigenesis (Guo et al., 2004; Zhao et al., 2012; Gamage & Hendrickson, 2013; Nakakido et al., 2016). Here, for the first time, we describe that in comparison with healthy tissue, CC also displays upregulation of two GPI-anchor biosynthesis glycogenes; interestingly, a high expression of some GPI-anchored proteins has been reported in CC (Jing et al., 2012; Liu et al., 2017; Jöhrens et al., 2019). To know whether the increase of GPI-anchored proteins is due to the rise of GPI biosynthesis should be tested in the future. Concerning HA, CC displayed downregulation of HAS3 and HYAL2 in comparison with healthy tissue. In HA metabolism, hyaluronan synthases 1 and 2 synthesize high molecular weight HA, while the hyaluronan synthase 3 (encoded by HAS3) synthesizes low molecular weight HA (Itano et al., 1999; Itano & Kimata, 2002); in contrast, HA degradation is induced by a family of hyaluronidases including Hyal2 (Csoka, Frost & Stern, 2001; Kaneiwa et al., 2012). Some evidence in SCC indicates that patients with a good prognosis display abundant levels of HA compared with those with the worst prognosis (Sano & Ueki, 1987). The accumulation of HA is typical in several types of cancers, and its levels directly correlate with increased malignancy and a poor prognosis (Passi et al., 2019). For example, in endometrial cancer, HYAL1 and HYAL2 are downregulated compared to healthy tissue, and this phenotype correlates with the accumulation of hyaluronan (Nykopp et al., 2010). Studies in animals indicate that HYAL2 or HYAL1 inhibits tumor growth and may control intercellular interactions (Wang et al., 2008). Examination of whether downregulation of HYAL2 could lead to HA accumulation can be interesting in CC. Concerning HAS3, it is unknown whether the low expression of HAS3 leads to a specific decrease in low molecular HA synthesis. Still, the low expression is associated with a poor prognosis in urothelial carcinoma (Chang et al., 2015). Characterization of the content levels of high-or low-molecular-weight HA in CC can be interesting, especially because their effects are different in cancer cell behavior and chemotherapy resistance (Price, Lokman & Ricciardelli, 2018; Tavianatou et al., 2019a; Tavianatou et al., 2019b). In summary, the results suggest that CC may display an increase in the synthesis of GPI anchors and dysregulation in the metabolism of HA.

The results obtained by the microarray assay allowed us to understand the global changes from the normal cervix to a cancer stage, and some of them show some coherence with the RNAseq analysis. Thus, the low expression of UGT2B4 and UGT2B28 in CC, analyzed in the microarray assay, was also found in all CC samples in the RNAseq analysis (Table S2-column B). Moreover, the low expression of HAS3, HPSE, and GYS1 was found in the AC-B and AC-C clusters (Table S2-column C and column D, respectively). In the same way, the high expression of ST6GAL1 was also found in the cluster AC-A (Table S1-column B). Despite these results, one limitation of our study is to validate the glycogene expression by RT-PCR. However, the validation can be difficult, since implies to have a previously characterized glycogene expression of the samples, considering the different glycogene signatures found in the SCC. Hence, the microarray assay of pooled samples allowed us to understand global changes of CC but not the molecular diversity of CC.

Adenocarcinoma displays unique glycogene expression, but squamous cancer shows at least five types of glycogene signatures

Genomics studies suggest that CC is not a unique entity from molecular and genetic viewpoints (The Cancer Genome Atlas Research Networkk (TCGA), 2017; Srinivasan, 2019). Molecular characterization of CC includes three subgroups: keratin-low squamous, keratin-high squamous, and adenocarcinoma-rich (The Cancer Genome Atlas Research Networkk (TCGA), 2017). Regarding the glycome, our result suggests that AC displays a unique glycogene expression signature distinguishable from SCC, which in turn displays at least five types of glycogene signatures. Results shown in this study could explain some of the results where certain glycogenes in CC display high expression while other are exhibit low expression (Varchalama et al., 2009), and even where some of them are associated with specific cancer characteristics associated with patient prognosis (Varchalama et al., 2009; Wang et al., 2003).

First, we showed that AC displays a glycogene signature characterized by 42 glycogenes (Table 5; Table S3), where some of them, such as ALG11, ALG9, and B3GALT5, can be relevant for cancer patients (Yu et al., 2020; Liu et al., 2021; Liao et al., 2021). In comparison, SCC was characterized by a heterogenicity of glycogene expression that can be distinguished by a selective glycogene expression (Fig. S2). The results suggest that the composition of characteristics such as tumor stage, differentiation grade, molecular characterization, and global patient survival comprise a global value related to the glycogene expression pattern. SCC-2 was characterized only by the low expression of glycogenes (Table S2-column G). Also, we identified four kinds of SCC. Type SCC-1/SCC-4 displayed, specifically, the expression of 38 glycogenes, while six glycogenes recognized SCC-3. In comparison, SCC-5 showed the expression of eight glycogenes, and SCC-6 could be distinguished by 36 glycogenes (Table 5; Table S3). The glycogene signatures of SCC-1/4, SCC-2, and SCC-6 were associated with 88–89% of overall survival. In contrast, the glycogene signatures of SCC-3 and SCC-5 displayed 68% and 50% of overall survival, respectively. These results suggest that glycogenes CHST2, COLGALT1, DSE, EXT1, EXT2, and GALNT18 in SCC-3 can be helpful as a biomarker of prognostic. In addition, detecting high levels of CHST11, CHST15, and GFPT2 can also be useful as prognostic biomarkers since their increased expression is associated with 50% of overall survival (in SCC-5). In comparison, their low expression is associated with 88−89% of overall survival (SCC-2). Further characterization of protein expression levels characterization and their association with survival will be necessary to consider some key glycogenes as valuable biomarkers.

Regarding the different expression of glycogenes in SCC, the whole molecular context can affect the gene expression in CC: miRNA, methylation, signaling pathway signatures, structural aberration (Yang et al., 2009; Farkas et al., 2013; Vojta et al., 2016; The Cancer Genome Atlas Research Networkk (TCGA), 2017) and the integrated HR-HPV genome (The Cancer Genome Atlas Research Networkk (TCGA), 2017; Aco-Tlachi et al., 2018; Cisneros-Ramírez et al., 2020; Gagliardi et al., 2020), which we did not explore.

Glycogene expression diversity in CC and putative consequences in their glycosylation pathways

Glycosylation changes in cancer can be due to changes in the availability and abundance of the sugar nucleotide donors and cofactors, altered enzyme activity, and to dysregulation at the transcriptional level of the glycogenes (Pinho & Reis, 2015; Schjoldager et al., 2020). The evidence shown here suggests that CC may display a glycome diversity due to variability in the glycogene expression, whether AC or certain SC types. Regarding AC, results indicate that it is characterized by the increase of genes implicated in KS biosynthesis, N-glycan biosynthesis, and glycosphingolipids genes related to ganglio and globo series. In contrast, AC displayed low expression of glycogenes involved in the CS/DS and HS biosynthesis, GAG degradation, and GPI-anchor biosynthesis (Fig. 5). Regarding SC, KEGG analysis revealed that the SCC-2 displayed low expression of glycogenes implicated in the CHS biosynthesis, GAG degradation, and glycosphingolipids biosynthesis (Fig. 5).

In comparison, SCC-1/SCC-4 displayed an increase in the glycogenes implicated in HS/heparin biosynthesis and GAG degradation (Fig. 5). Type SCC-3 exhibited high expression of glycogenes implicated in HS/heparin biosynthesis (Fig. 5). SCC-6 was distinguished by the high expression glycogenes involved in several glycosylation pathways including in HS and KS biosynthesis and glycosphingolipid biosynthesis of the ganglio and globo series (Fig. 5). Finally, SCC-5 displayed the expression of eight glycogenes involved in CS/DS biosynthesis, glycosphingolipid biosynthesis of lacto and neolacto series, and mucin-type O-glycan biosynthesis (Fig. 5). All the data suggest a putative glycome diversity in CC, and future studies can be directed to characterize glycoconjugates accompanied by a glycogene signature.