Proteogenomic characterization of molecular and cellular targets for treatment-resistant subtypes in locally advanced cervical cancers

Abstract

We report proteogenomic analysis of locally advanced cervical cancer (LACC). Exome-seq data revealed predominant alterations of keratinization-TP53 regulation and O-glycosylation-TP53 regulation axes in squamous and adeno-LACC, respectively, compared to in early-stage cervical cancer. Integrated clustering of mRNA, protein, and phosphorylation data identified six subtypes (Sub1-6) of LACC among which Sub3, 5, and 6 showed the treatment-resistant nature with poor local recurrence-free survival. Elevated immune and extracellular matrix (ECM) activation mediated by activated stroma (PDGFD and CXCL1^high fibroblasts) characterized the immune-hot Sub3 enriched with MUC5AC^high epithelial cells (ECs). Increased epithelial-mesenchymal-transition (EMT) and ECM remodeling characterized the immune-cold squamous Sub5 enriched with PGK1 and CXCL10^high ECs. We further demonstrated that CIC mutations could trigger EMT activation by upregulating ETV4, and the elevation of the immune checkpoint PVR and neutrophil-like myeloid-derived suppressive cells (FCN1 and FCGR3B^high macrophages) could cause suppression of T-cell activation in Sub5. Increased O-linked glycosylation of mucin characterized adeno-LACC Sub6 enriched with MUC5AC^high ECs. These results provide a battery of somatic mutations, cellular pathways, and cellular players that can be used to predict treatment-resistant LACC subtypes and can serve as potential therapeutic targets for these LACC subtypes.

Keywords: And treatment resistance; Cancer subtyping; Locally advanced cervical cancer; Proteogenomics.

Fig. 1

Predominant alteration of keratinization-TP53 regulation axis in LACC. A Predominantly altered genes in LACC (PAG-LACC). Mutations per megabase in each patient are shown in the top panel. SMG mutations (2nd panel), alterations by HPV genome integration (3rd panel), and CNAs (amplification in red and deletion in blue; 4th panel) detected in each patient for the indicated genes are shown. Alteration frequencies are shown on the right, and the indicated clinical parameters for each patient are shown in the bottom. B-D Comparisons of frequencies of SMG mutations (B), alterations by HPV genome integration (C), and CNAs (D) between our LACC cohort and the two previous early-stage tumor-enriched cohorts. Red labels indicate the genes (PAG-LACC) having significantly (p < 0.05 by proportional test) higher alteration frequencies in our LACC cohort than in the early-stage tumor-enriched cohorts. *, p < 0.05; **, p < 0.01; ***, p < 0.001. n = 251, 289, and 102 patients for our, TCGA, and Huang et al. cohorts, respectively, for SMG mutations; n = 251, 178, and 45 patients for our, TCGA, and Huang et al. cohorts, respectively, for alterations by HPV genome integration and CNAs. E mRNAs, proteins, phosphorylated peptides up-regulated in tumors harboring the alterations (SMG, HPV integration, or CNA) in any of the indicated PAG-LACC. The colored bar represents the gradient of log₂-fold-changes of abundances for mRNAs, proteins, or phosphorylated peptides relative to their median levels. F Cellular pathways enriched by mRNAs (Gene) or proteins + phosphoproteins (Prot) correlated with alterations of the indicated PAG-LACC. The heat map shows the enrichment significance (p value from ConsensusPathDB) of each pathway by the mRNAs or proteins + phosphoproteins as –log₁₀(p). G Network model describing interactions among the mRNAs (blue node center), proteins (blue node border), or phosphoproteins (blue circled P) that have significant correlation with alterations of the PAG-LACC (red labeled) and are involved in keratinization and apoptosis/TP53 regulation. Gray nodes indicate molecules added to the pathways to increase connections among the nodes. Solid arrows, direct activations; dotted arrows, indirect activations; gray lines, protein–protein interactions; thick lines, plasma membrane

Fig. 2

Proteogenomic subtypes of LACC. A RNA1-4 clusters defined by mRNA signatures (rna1-4) in our LACC cohort. rna1-4 defining RNA1-4, respectively, are shown. Numbers of mRNAs in rna1-4 are indicated in parentheses. Clinical parameters for tumors are shown in the bottom, and the subtype bar plot shows subtypes predicted using molecular signatures defined by the TCGA [10]. B and C, Protein signatures (prot1-5 and phos1-4) defining Prot1-5 (B) and Phos1-4 (C) based on global proteome and phosphoproteome data. Numbers of proteins and phosphopeptides are indicated in parentheses. D Six subtypes (Sub1–6) identified via integrative clustering of mRNA, protein, and phosphorylation data. Each row in the heat map shows Z-values of the tumors with the corresponding molecular signature for each indicated cluster. Color bar, the gradient of Z-value representing distances from mean vectors with the signature. Colored bars at the top, Sub1-6; RNA1-4 clusters defined by rna1-4. E Disease-free survival of patients with tumors belonging to Sub1-6. n = 29, 15, 17, 23, 14, and 16 in Sub1-6, respectively. F Somatic mutations of frequently mutated genes enriched in Sub5 (red), Sub6 or Sub6 + ADC in Sub3 (green), or all Sub3, 5, and 6 (black). Fractions of tumors carrying mutations of each gene in Sub1-6 were shown in the right panel. G mRNAs, proteins, phosphorylated peptides up-regulated in tumors harboring ADC-enriched somatic mutations in the 11 genes (green labeled in F). The colored bar represents the gradient of log₂-fold-changes of abundances for mRNAs, proteins, or phosphorylated peptides relative to their median levels. H Cellular pathways enriched by mRNAs or proteins + phosphoproteins (Prot/Phos) correlated with the ADC-enriched somatic mutations. The heat map shows the enrichment significance (p value from ConsensusPathDB) of each pathway by the mRNAs or Prot/Phos as –log₁₀(p). I Network model describing interactions among the mRNAs (green node center), proteins (green node border), or phosphoproteins (green circled P) that have significant correlation with the ADC-enriched somatic mutations and are involved in O-linked glycosylation and apoptosis/TP53 regulation. PAG-LACC are labeled in red. Gray nodes indicate molecules added to the pathways to increase connections among the nodes. Solid arrows, direct activations; dotted arrows, indirect activations; gray lines, protein–protein interactions; thick lines, plasma membrane

Fig. 3

Characteristics of proteogenomic subtypes of LACC. A Cellular pathways represented by mRNAs (S1-G to S6-G) and proteins (proteins + phosphoproteins) (S1-P to S6-P) defining Sub1-6. The heat map shows the enrichment significance (p value from ConsensusPathDB) of pathways by the mRNAs or proteins defining Sub1-6 as –log₁₀(p). B, D, and F, Network models showing interactions between the mRNAs and proteins involved in immune and ECM pathways upregulated in activated stroma for Sub3 (B); EMT signaling and ECM remodeling pathways for Sub5 (D); and O-linked glycosylation (F). Node colors (center and boundary) indicate whether the corresponding mRNA and protein were selected as signatures for the corresponding subtypes (green for Sub3, yellow for Sub4, orange for Sub5, and dark green for Sub6). Presence of a circled P on a node indicates that the corresponding phosphoproteins had phosphorylated peptides that defined the corresponding subtypes. Kinases that could control the pathways in the network models in Sub3 and 5 were identified as previously described [–30] and included to the network models (blue-labeled nodes). Arrows, activation; inhibition symbols, inhibition; solid arrows, direct activation; dotted arrows, indirect activation; gray lines, protein–protein interactions. C and E Representative images for MMP2 (Sub3), ICOS (Sub4), and GPC4 (Sub6) (C) and BSG, VAV2, EGFR, and FOSL1 (Sub5) (E) and their quantification results in Sub1-6 from IHC analysis. Magnification 20 × . Scale bar = 100 μm. Blue color: Hematoxylin counterstaining; Brown color: DAB positive staining. n = 17, 11, 16, 20, 13, and 16 in Sub1-6, respectively, for MMP2 and ICOS; n = 17, 10, 15, 20, 13, and 16 in Sub1-6, respectively, for GPC4, BSG, VAV2, EGFR, and FOSL1. *, p < 0.05; **, p < 0.01 by comparing H-scores (grades or density) between a subtype of interest and the other subtypes using one-sided Fisher’s exact test (MMP2) or one-sided Student’s t-test (ICOS, GPC4, VAV2, EGFR, and FOSL1). G Representative images of CIC, ETS transcription factors (ETV4), and cytoskeletal proteins (CDH1/2, CTNNA1, and VIM) in #6507A (Sub5) and #6595 (Sub3) cells from Western blotting analysis. β-actin was used as a loading control. n = 2 independent experiments. H Quantification of indicated proteins relative to β-actin in the indicated cells. Band intensity from western blot images was quantified by ImageJ. Data are shown as mean ± s.e.m. n = 2 independent experiments. *, p < 0.05 by Student’s t-test. I, Representative immunofluorescence images of CDH1 (red; epithelial marker) and VIM (green; mesenchymal marker) in the indicated cells grown in monolayer culture system (n = 2 independent experiments). Magnification 40 × . Scale bar = 50 μm

Fig. 4

Cellular heterogeneity associated with LACC Sub1-6. A and G Uniform manifold approximation and projection (UMAP) plot showing subclusters of ECs (A) or myeloid cells (G). Associations of the individual subclusters with Sub1-6 are indicated in parenthesis. The stacked bar graph shows percentages of ECs belonging to individual subclusters. B Marker genes upregulated predominantly in each subcluster. Numbers in parenthesis represents the number of marker genes in the corresponding subcluster. The color bar represents the gradient of log₂-FC with respect to their mean values. C, H, and K Dot plot showing the enrichment p-value of marker genes for each subcluster of ECs (C), myeloid cells (H), or CAFs (K) in the indicated mRNA and protein signatures for Sub1-6 (see legend in bottom right). The associated subtype for each signature is denoted in parenthesis. Dot size represents the fraction of overlapping marker genes as shown in legend (bottom left). D, I, and L Signature scores of marker genes for each subcluster of ECs (D), myeloid cells (I), or CAFs (L) in Sub1-6. n = 33, 17, 27, 28, 18, and 23 in Sub1-6, respectively. E, J, and M Cellular pathways enriched by marker genes for each subcluster of ECs (E), myeloid cells (J), or CAFs (M). The color bar represents the gradient of –log₁₀(enrichment p-value). F UMAP plot showing ADC ECs with transferred labels of the subclusters of SCC ECs. The stacked bar graph shows percentages of ADC ECs having individual transferred labels

Fig. 5

Immune-suppressing pathways in treatment-resistant LACC Sub5. A PVR mRNA (left) and protein (right) expression patterns across Sub1-6. In violin plots, the line indicates the median value. B Representative images for PVR and its quantification results in Sub1-6 from IHC analysis. C Representative images for each subtype of 5-color mIHC staining scanned with PhenoImagerTM HT at 20 × magnification: PVR (red), TIGIT (yellow), NL-SC marker CD66B (cyan), cytokeratin (CK, white) and nucleus marker DAPI (blue). Scale bar = 200 μm. In the original and magnified images of spatial cell-to-cell distance analysis, TIGIT-expressing cells were represented as green dots while PVR-expressing tumor cells within the 30 μm radius from TIGIT-expressing cells were denoted as red dots. Box plots (right) showed the quantification results analyzed from 90 tumor tissues on TMA. D Quantification results in Sub1-6 from IHC analysis of CD66B. E Signature scores of mRNA/protein signatures of Sub5 involved in IL-17 signaling in the indicated myeloid subclusters. F IL-17R expressions on NL-SCs in bone marrow, spleen, blood, and tumor samples from #6507A (Sub5) tumor-bearing mice (n = 4). G Proportions of different IL-17R⁺ cell types in #6507A (Sub5) tumors. H Representative images showing migration of NL-SCs (black triangles) in control and IL-17-treated groups. Scale bar = 100 µm. I Number of migrated NL-SCs per field under indicated cytokine conditions. J Co-culture of NL-SCs and T-cells isolated from #6507A (Sub5) tumor models and naïve Balb/c mice, respectively. K Distribution of CFSE intensity for CD8⁺ T-cells measured at Day 4 after the co-culture. L MFI of CFSE intensity for CD8⁺ T-cells (n = 3/group). M Effector cytokine expression (GZMB: Granzyme B, Perforin, IFN-γ, TNF-α) in CD8⁺ T cells (n = 3/group). Data are shown as the mean ± SEM. *, p < 0.05; **, p < 0.01; ***, p < 0.001 from Student’s t-test (L and M) and one-way ANOVA with Tukey’s posthoc correction (F and I)