ENST00000535926 is an unfavorable prognosis-related and tumor-promoting transcript of the CHPF gene in luminal A and B breast cancer

Introduction

Glycosaminoglycans (GAGs) are long and unbranched polysaccharides, acting as extracellular matrix components (Afratis et al., 2012). GAGs are frequently dysregulated in the tumor microenvironment (TME) of breast cancer (Hampton et al., 2022; Kowitsch et al., 2017; Wei et al., 2020). Chondroitin sulfate is a major component of the GAG family and can covalently link to the core protein, forming chondroitin sulfate proteoglycans (CSPGs). CSPGs participate in diverse physiological events, such as cytokinesis, morphogenesis, and neuronal plasticity (Wei et al., 2020). Excessive chondroitin sulfate generation can exacerbate the aggressive behaviors of breast cancer cells by promoting growth, invasion, and metastasis (Cooney et al., 2011; Iida et al., 2015; Monzavi-Karbassi et al., 2007).

In humans, biosynthesis of chondroitin sulfate chains is initiated via the linkage of N-acetylgalactosamine (GalNAc) to a tetrasaccharide structure, glucuronic acid (GlcA)-GalNAc-GalNAc-xylose, on a core protein (Mikami and Kitagawa, 2013). Then, the chondroitin sulfate chain is elongated by chondroitin sulfate synthases (CHSY1, CHPF/CHSY2, CHPF2, and CHSY3), that have both β-1,3-GlcA and β-1,4-GalNAc transferase activity (Mikami and Kitagawa, 2013). CHPF upregulation is common in multiple types of cancer (Hou et al., 2019; Liao et al., 2021; Lin et al., 2021). CHPF promotes triple-negative breast cancer cell growth, metastasis, and accumulation of myeloid-derived suppressor cells (MDSCs) (Liao et al., 2021). Its prognostic significance in breast cancer was also confirmed by a recent bioinformatic analysis (Li et al., 2022).

However, breast cancer is a group of heterogeneous tumors with distinct clinical and genomic characteristics (Liu et al., 2022; Nielsen et al., 2014; Wang et al., 2022). A 50-gene quantitative polymerase chain reaction (qPCR assay; Prediction Analysis of Microarray 50, PAM50) of gene profiles can generally classify breast cancer into luminal A, luminal B, human epidermal growth factor receptor 2 (HER2)-enriched, basal-like (largely triple-negative), and normal-like breast cancer subtypes (Nielsen et al., 2014). In addition, CHPF is a gene with alternative splicing, with two isoforms (isoform 1: ENST00000243776 and isoform 2: ENST00000535926) expressing two proteins with different lengths (775 amino acids (aa) and 613 aa, respectively, Suppl. Fig. S1). The short isoform 2 (also called Klokin1) can facilitate the transport of parkin to the mitochondria, enhance cell viability, and protect neuronal cells against oxidative stress (Kuroda et al., 2012).

Considering the important tumor-promoting effects of CHPF in breast cancer, it is necessary to explore whether it has transcript-specific regulations. In this study, we explored the expression and prognostic significance of the two transcripts using data from The Cancer Genome Atlas (TCGA)-breast cancer (BRCA). In addition, the functional regulation of the two transcripts was also studied in MCF-7 and BT-474 cells.

Materials and Methods

Ethical statement

The bioinformatic part of this study is based on Genotype-Tissue Expression (GTEx), and TCGA, in which the data are anonymous and informed consent was appropriately handled. Ethical approval was not required, as no primary patient data were collected by any author in the current study.

Data retrieval and bioinformatic analysis

RNA-seq data of normal breast tissues in the GTEx project and tumor and tumor-adjacent normal in TCGA-BRCA were obtained as we previously described (Luo et al., 2020; Tian et al., 2020), using the UCSC Xena browser (https://xenabrowser.net/) (Goldman et al., 2020). RNA-seq data were obtained for 179 normal breast tissues in GTEx and 1210 cases in TCGA-BRCA (113 tumor-adjacent normal tissues and 1097 primary tumors). Among the 1097 primary tumor cases, 839 had PAM50 information, including 420 luminal A, 192 luminal B, 66 HER2+, 139 basal-like, and 22 normal-like. The transcript expression and isoform percentages were compared between GTEx and TCGA-BRCA. The survival data, including progression-free survival (PFS) and disease-specific survival (DSS), were extracted. Kaplan-Meier (K-M) survival curves were generated by median transcript or gene expression in PAM50 subgroups.

The top CHPF-correlated genes in TCGA-BRCA were identified using cBioPortal for Cancer Genomics (https://www.cbioportal.org/). Then, the correlation between these genes and the two CHPF transcripts was assessed separately by calculating Pearson’s correlation coefficients.

Cell culture and treatment

Since we observed that the dysregulation of certain CHPF transcripts might be associated with poor prognosis of luminal A and liminal B subtypes, we used luminal A representative tumor cell line MCF-7 and luminal B representative tumor cell line BT-474 for in vitro study. These two cell lines were cultured as described previously (Tian et al., 2020). Lentiviral ENST00000243776 or ENST00000535926 overexpression plasmids were generated based on the pLenti-CMV-puro backbone. Empty lentiviral plasmids (vector) served as a negative control. Lentivirus production and infection following the methods introduced previously (Luo et al., 2020; Tian et al., 2020).

Quantitative real-time polymerase chain reaction (qRT-PCR) assay

qRT-PCR assays were performed as we described previously (Tian et al., 2020). In brief, total RNA was extracted from MCF-7 and BT-474 cells with or without ENST00000535926 or ENST00000243776 overexpression using Trizol Reagent (Thermo Fisher Scientific, Waltham, MA, USA). Complementary DNA was reversely transcribed and used as the template for following PCR assays. Targeting mRNAs were amplified using SYBR Premix Ex Taq (TaKaRa, Dalian, China) with the ABI PRISM 7900HT Sequence Detection System (Applied Biosystems). We checked the data in Internal Control Genes (http://icg.big.ac.cn/index.php/Main_Page) (Sang et al., 2018). ACTB is among the appropriate reference genes for research and clinical analysis of breast tumor samples. The following primers were used: human CHPF, 5′-AACGCACGTACCAGGAGATCC-3′ (forward) and 5′-GGATGGTGCTGGAATACCCACG-3′ (reverse); human P3H1, 5′-GTCACTGGATGTGAGCAGACTG-3′ (forward) and 5′-ACTCGTGGTCAGAGATTACGCC-3′ (reverse); human RCN3, 5′-TCGAGAGGAGCTGACAGCCTTC-3′ (forward) and 5′-TCCTCCACCTGGACATAGCCAT-3′ (reverse); human EFEMP2, 5′-ACTCCTATGGGACCTTCCTGTG-3′ (forward) and 5′-CTCGTTGATGCAGCGGTACTGA-3′ (reverse); human ACTB, 5′-CACCATTGGCAATGAGCGGTTC-3′ (forward) and 5′-AGGTCTTTGCGGATGTCCACGT-3′ (reverse). The melting curve with a single peak was considered a pure and single amplicon and was used for the following data analysis. Representative melting curves were provided in Suppl. Fig. S2. Relative mRNA levels were normalized to ACTB and calculated using the 2−ΔΔCT method, based on three times of repeats in triplicate.

Western blot analysis

Western blotting was performed following the methods introduced previously (Tian et al., 2020). In brief, total proteins were extracted from MCF-7 and BT-474 cells, separated using 10% sodium dodecyl sulfate-polyacrylamide gel electrophoresis, and electro-transferred onto polyvinylidene difluoride membranes (Millipore, Billerica, MA, USA). Then, the membranes were blocked and incubated with primary antibodies against CHPF (1: 500, YT0924, Immunoway, Suzhou, China), P3H1 (1: 1000, 26691-1-AP, Proteintech, Wuhan, China), RCN3 (1: 1000, 27497-1-AP, Proteintech), EFEMP2 (1: 1000, 12004-1-AP, Proteintech), and β-actin (1: 2000, 20536-1-AP, Proteintech) overnight at 4°C. Then, the membranes were incubated with horse radish peroxidase-conjugated secondary antibodies for 1 h at room temperature. The protein band signals were developed using BeyoECL Star reagent (Beyotime, Shanghai, China) and an ImageQuant LAS-4000 imaging system (GE Healthcare, Piscataway, NJ, USA). The experiment was repeated three times.

Colony formation assays

Cells with or without ENST00000535926 and ENST00000243776 overexpression were seeded into 6-well plates (500 cells/well) and were cultured for about two weeks. Then, the cellular colonies were fixed with methanol and stained with 0.5% crystal violet. The plates were scanned, and the colonies were visually counted. The experiment was repeated three times in triplicate.

Estimation of tumor infiltration

The correlation between P3H1/LEPRE1 or RCN3 expressions and the infiltration of immune cells (CD8+ T cells, Tregs, MDSC, and cancer-associated fibroblasts) in luminal A and luminal B cases in TCGA were estimated using TIMER 2.0 (http://timer.cistrome.org/) (Li et al., 2020). This platform provides computational algorithms to infer immune infiltrates’ abundances from bulk tumor transcriptome profiles in TCGA. Estimation was performed in luminal A and luminal B cases separately.

Statistical analysis

Statistical analysis was performed using GraphPad Prism 8.1.2 (GraphPad Inc., La Jolla, CA, USA). One-way ANOVA and post-hoc Dunnett’s multiple comparisons test were performed for multiple group comparisons. Welch’s unequal variances t-test was conducted for a two-group comparison. The log-rank test was performed for survival analysis. A value with p < 0.05 was considered statistically significant.

Results

ENST00000535926 is the dominant CHPF transcript in breast cancer tissues

I expression profile of CHPF transcripts (ENST00000243776 and ENST00000535926) in normal and cancerous breast tissues was examined by comparing RNA-seq data in GTEx-breast (all normal tissues, n = 179) and TCGA-BRCA (PAM50 subgroups), in terms of transcript expression and isoform percentage. The wave diagram and plot chart showed that ENST00000535926 and ENST00000243776 have similar TPM values and percentages in normal breast tissues (Figs. 1A and 1B).

Figure 1: ENST00000535926 is the major CHPF transcript in breast cancer tissues. (A, B) Wave charts were generated to compare the expression (log2(TPM+0.001)) (A) and isoform percentage (B) of ENST00000535926 and ENST00000243776 in GTEx-normal breast tissues and primary breast tumors in TCGA. (C–F) Plot charts were generated to compare the expression (log2(TPM+0.001)) (C and E) and isoform percentage (D and F) of ENST00000535926 (C, D) and ENST00000243776 (E, F) in GTEx-normal breast tissues and the PAM50 subgroups in TCGA-BRCA. GTEx-B: GTEx-normal breast; adj. N: adjacent normal. One-way ANOVA and post-hoc Dunnett’s multiple comparisons test were performed. *p < 0.05; ***p < 0.001.

Compared to the GTEx-normal group, ENST00000535926 expression was significantly upregulated in tumor-adjacent normal tissues and all PAM50 tumor subgroups (Fig. 1C). In addition, its isoform percentage increased from approximately 40% in GTEx-normal breast tissues to over 70% in the cancer groups in TCGA (Fig. 1D). However, ENST00000243776 expression was significantly lower in tumor-adjacent normal tissues in TCGA-BRCA than in GTEx (Fig. 1E). Its expression was similar among GTEx-normal, luminal A, luminal B, HER2+, and basal tumor subgroups (Fig. 1E). In terms of isoform percentage, it dropped from around 50% in GTEx-normal breast to about 20% in the cancer groups in TCGA (Fig. 1F).

ENST00000535926 upregulation is associated with unfavorable survival in patients with luminal A and B tumors

Since the PAM50 subgroups have significant heterogeneity in the genomic background and clinical outcomes, we assessed the association between the expression of the two transcripts and survival in each PAM50 subgroup. Patients were categorized into two groups by median transcript expression. K-M survival analysis showed that high ENST00000535926 expression was associated with significantly shorter PFS and DSS in patients with luminal A and B tumors (Figs. 2A, 2B, 2E, 2F). However, no significant association was observed in HER2+ or basal subgroup (Figs. 2C, 2D, 2G, 2H). High ENST00000243776 expression was associated with significantly worse DSS in luminal B cases (Fig. 2N). However, no significant association was observed in PFS or DSS in other subgroups (Figs. 2I–2M, 2O, 2P).

Figure 2: ENST00000535926 upregulation is associated with unfavorable survival in patients with luminal A and B tumors (A–P). (K–M) Survival curves were generated to compare the difference in PFS (A–D, I–L) and DSS (E–H, M–P) in the PAM50 subtypes of breast cancer in TCGA. Patients were separated by median ENST00000535926 expression (exp.) (A–H) or ENST00000243776 expression (I–P). Log-rank test was conducted to assess the statistical difference between survival curves.

Chondroitin sulfate synthase isoform 2 protein stimulates the expression of P3H1 and RCN3 and promotes tumor cell colony formation

Since the two transcripts of CHPF have distinct prognostic significance, we tried to assess their functional differences in breast cancer. Using cBioPortal for Cancer Genomics, we characterized the top 19 genes correlated with total CHPF gene expression in luminal A and B cases in TCGA-BRCA. Then, the correlations between the expression of these genes and the two CHPF transcripts were assessed separately (Fig. 3A). By setting 0.3 as the cutoff for the difference in correlation coefficients, we found three genes (P3H1, RCN3, and EFEMP2) highly or moderately correlated with ENST00000535926 expression, but were only weakly correlated with or irrelevant to ENST00000243776 expression (Figs. 3A, 3C–3H). However, no significant correlation (|Pearson’s r| ≥ 0.2) was observed between the two isoforms and ESR1 expression in luminal A/B cases (Fig. 3B).

Figure 3: ENST00000535926 stimulates the expression of P3H1 and RCN3 and promotes tumor cell colony formation. (A, B) Summary tables showing the expressional correlation between the top CHPF correlated genes and the two CHPF transcripts (ENST00000535926 and ENST00000243776) (A) and between the expression of ESR1 and the two CHPF transcripts (B) in luminal A/B cases in TCGA-BRCA. (C–H) Plot charts were generated to show the expressional correlation between the gene and transcript pairs as indicated in luminal A and B cases in TCGA-BRCA. (I, J) Quantitative real-time polymerase chain reaction (qRT-PCR) assays were conducted to check the expression of CHPF mRNA 48 h after lentivirus-mediated ENST00000243776 (243776-OE) (I) or ENST00000535926 (535926-OE) (J) overexpression in MCF-7 and BT-474 cells. (K–M) qRT-PCR assays (K, L) and western blotting (M) were conducted to check the expression of P3H1, RCN3, and EFEMP2 at the mRNA (K, L) and protein (M) in MCF-7 and BT-474 cells 48 h after lentivirus-mediated ENST00000243776 (243776-OE) or ENST00000535926 (535926-OE) overexpression. (N–P) Representative images (N) and quantitation of (O, P) of colony formation of MCF-7 and BT-474 cells with lentivirus-mediated ENST00000243776 (243776-OE) or ENST00000535926 (535926-OE) overexpression. ##, comparison between 535926-OE and 243776-OE. *p < 0.05; ** and ## p < 0.01; ***p < 0.001.

Then, we overexpressed ENST00000535926 or ENST00000243776 in MCF-7 and BT-474 cells and detected the expression of these genes (Figs. 3I–3M). ENST00000243776 overexpression slightly increased the transcription of P3H1 (Fig. 3K). In comparison, ENST00000535926 overexpression increased the transcription of P3H1 and RCN3 (Fig. 3K). The results of western blotting revealed that ENST00000535926, but not ENST00000243776 overexpression, elevated the expression of P3H1 and RCN3 proteins levels (Fig. 3M). Colony formation assays confirmed that ENST00000535926 overexpression drastically enhanced colony formation of MCF-7 and BT-474 cells. These effects were much more potent than ENST00000243776 overexpression (Figs. 3N–3P).

P3H1 and RCN3 expression are negatively correlated with CD8+ T cell infiltration but positively correlated with cancer-associated fibroblast in luminal A/B tumors

Since the positive regulation of the CHPF isoform 2 on P3H1 and RCN3 expression was confirmed, we further analyzed the prognostic significance of these two genes. K-M survival analysis showed that although high P3H1 or RCN3 expression might be linked to poor PFS and DSS, the statistical analysis only confirmed a significant difference in DSS by median P3H1 separation (Figs. 4A–4H).

Figure 4: The survival significance of P3H1 and RCN3 expression in patients with luminal A and B tumors (A–H). (K–M) Survival curves were generated to compare the difference in PFS (A, C, E, and G) and DSS (B, D, F, and H) in the luminal A (A, B, E, and F) and luminal B (C, D, G, and H) subtypes of breast cancer in TCGA. Patients were separated by median P3H1 expression (A–D) or RCN3 expression (E–H). Log-rank test was conducted to assess the statistical difference between survival curves.

Considering that CHPF might regulate TME via its downstream effectors, we further assessed the correlation between P3H1 or RCN3 expression and tumor-infiltrating immune cells (CD8+ T cells, Tregs, myeloid-derived suppressor cells (MDSC), and cancer-associated fibroblasts) in luminal A and luminal B tumors, respectively. P3H1 showed a strong negative correlation with CD8+ T cell infiltration in luminal A tumors (Pearson’s r = −0.607) (Fig. 5A) and a moderate negative correlation (Pearson’s r = −0.549) in luminal B tumors (Fig. 5B). RCN3 showed a strong negative correlation with CD8+ T cell infiltration in both luminal A and B tumors (Pearson’s r = −0.653 and −0.67, respectively) (Figs. 5C and 5D).

Figure 5: P3H1 and RCN3 are negatively correlated with CD8+ T cell infiltration but were positively correlated with cancer-associated fibroblast in luminal A/B tumors (A–D). Plots charts were generated to show the correlation between P3H1/LEPRE1 (A, B) or RCN3 (C, D) expression and the infiltration of immune cells (CD8+ T cells, Tregs, MDSC, and cancer-associated fibroblasts) in luminal A (A and C) and luminal B (B and D) cases in TCGA. Data was obtained from TIMER2.0 (http://timer.cistrome.org/).

Discussion

RNA splicing is a process of RNA transformation in which the pre-mRNAs become mature mRNAs by excising introns and ligating exons (Yang et al., 2019). One direct consequence of RNA splicing is functional protein diversity (Birzele et al., 2008). Growing evidence demonstrates that alternative splicing is a molecular marker of human cancer and might provide new therapeutic targets (Zhang et al., 2021). Previous studies have characterized some important aberrant alternative splicing events in breast cancer (Hatami et al., 2013; Wang et al., 2016; Yang et al., 2019). Breast cancer type 1 (BRCA1) is a tumor suppressor gene in breast cancer. However, BRCA1-Δ11q (partial skipping of exon 11) can promote partial PARP inhibitor (PARPi) and cisplatin resistance (Wang et al., 2016). Compared to the full-length HER2, Δ16HER2 lacks exon 20 and encodes a small extracellular region. It can augment the critical oncogenic signals to promote HER2-driven tumorigenesis, cancer stemness, and drug resistance (Wang et al., 2016). Multiple protein isoforms of Kruppel-like factor 6 (KLF6) derived from alternative splicing have been observed in breast cancer, including KLF6-SV1, KLF6-SV2, and KLF6-SV3. The full-length KLF6 and KLF6-SV3 can be transferred to the nucleus using the nuclear localization signal (NLS). However, KLF6-SV1 and KLF6-SV2 localize to the cytoplasm due to the lack of NLS. KLF6-SV1 has been validated as a key driver of breast cancer metastasis by promoting epithelial-to-mesenchymal-like transition (Hatami et al., 2013).

Alternative splicing of CHPF and the functional differences of the isoforms have been reported in normal human neurons (Kuroda et al., 2012). The link between total CHPF expression and unfavorable survival of breast cancer patients was also observed (Li et al., 2022; Liao et al., 2021). However, it is not confirmed whether the link is transcript specific. In this study, using RNA-seq data from both TCGA-BRCA and GTEx, we compared the expression of the two CHPF transcripts and observed compositional alteration. ENST00000535926 expression was significantly upregulated in the tumor samples and became the dominant transcript. ENST00000535926, but not ENST00000243776 upregulation, was associated with significantly worse PFS and DSS in both luminal A and luminal B cases. However, no significant association was observed in PFS or DSS in other PAM50 subgroups. In addition, we validated that ENST00000535926 overexpression had a significantly stronger effect on promoting tumor cell colony formation of MCF-7 and BT-474 cells. These findings suggest that the prognostic significance of CHPF and its tumor-promoting effects might be transcript-specific in certain PAM50 subgroups.

Our correlation analysis revealed a significant difference in genes correlated with ENST00000535926 or ENST00000243776 expression. Molecular studies further confirmed that CHPF isoform 2 (encoded by ENST00000535926) can significantly elevate the expression of P3H1 and RCN3 at mRNA and protein levels in MCF-7 and BT-474 cells. P3H1 encodes Prolyl 3-Hydroxylase 1, a member of the collagen prolyl hydroxylase family. Its upregulation is associated with a poor prognosis of osteosarcoma (Huang et al., 2020). Its knockdown can reduce tumor cell proliferation, migration, and invasion (Huang et al., 2020). RCN3 encodes reticulocalbin 3, which acts as a molecular chaperone during protein biosynthesis and transport in the endoplasmic reticulum (Martinez-Martinez et al., 2017).

By performing immune-cell infiltration analysis, we confirmed that P3H1 and RCN3 expression are negatively correlated with CD8+ T cell infiltration but positively correlated with cancer-associated fibroblast infiltration in luminal A/B tumors. A recent single-cell RNA-seq data revealed that RCN3 is a fibroblast-specific biomarker of unfavorable prognosis in colorectal cancer (Zhou et al., 2020). Its upregulation can elevate the expression of immune-related factors, including TGF-β1, IL-10, and IL-6 in colorectal cancer cells (Ding et al., 2022). These cytokines and chemokines can help generate an immunosuppressive tumor microenvironment (Felcher et al., 2022; Huang et al., 2017). Typically, TGF-β1 can promote the transition of resident fibroblasts to cancer-associated fibroblasts and enhance cancer metastasis (Yoon et al., 2021). IL-10 expressed by tumor-infiltrating Tregs can drive the exhaustion of intratumoral CD8+ T cells (Damo and Joshi, 2019). Therefore, we infer that the CHPF-RCN3 axis might contribute to T cell exclusion or fibroblast recruitment to the breast cancer TME.

Patients with luminal A subtype have a very good prognosis, and systemic adjuvant therapy is most likely hormone therapy alone (Pellegrino et al., 2021). For luminal B cases, tumor-infiltrating lymphocytes are associated with prognosis, suggesting that immunotherapy might bring survival benefits (Hammerl et al., 2020). Therefore, targeting CHPF isoform 2 might be explored as a potential strategy to generate a TME favorable for immunotherapy for luminal B cases.

This study also has some limitations. First, although transcript-specific overexpression was used, we could not check how transcript-specific knockdown of CHPF influences the expression of the related genes. ENST00000535926 is a part of ENST00000243776, which means any shRNA targeting ENST00000535926 will lead to the knockdown of both transcripts. Second, only two ER+ cell lines were used for in vitro studies, which might influence the robustness of the findings. In the future, we will explore the specific regulatory mechanisms of ENST00000535926 on the tumor immune microenvironment and the malignant behaviors of other ER+ tumor cell lines.

Conclusion

This study revealed that ENST00000535926 is an unfavorable prognosis-related and tumor-promoting transcript of the CHPF gene in luminal A and B breast cancer. It stimulates the expression of P3H1 and RCN3. The expression of these genes is negatively correlated with CD8+ T cell infiltration but is positively correlated with cancer-associated fibroblast infiltration in both luminal A and B tumors.

Availability of Data and Materials: All data generated or analyzed during this study are included in this published article (and its supplementary information files).

Author Contributions: Study design: JL and CY; experiments carry out: JL and JPH; data collection: JL and JPH; analysis of results: JL, JPH, and YL; draft manuscript preparation: JL, JPH, and YL. All authors reviewed the final version of the manuscript.

Funding Statement: This work was supported by the Science and Technology Plan Project of Sichuan Province (Provincial Academy and Provincial University Cooperation Project) (2020YFSY0025).

Conflicts of Interest: The authors declare that they have no conflicts of interest to report regarding the present study.

Supplementary File

FIGURE S1: The protein sequences of Chondroitin sulfate synthase 2 (CHPF) isoforms.

FIGURE S2: The representative melting curves for quantitative real-time polymerase chain reaction analysis of CHPF expression.