Extensive genomic alterations of GSDM genes in cancer
Genomic alterations were defined as gene amplifications, deep deletions, fusion, truncating mutations, in-frame mutations, missense mutations or splice mutation using the cBioPortal definitions [15, 16]. Copy number alterations (CNAs) were classified as gene amplifications, gain, shallow deletion and deep deletions. We performed mutation and CNAs analysis of GSDM genes to identify genomic alterations across the Pan-Cancer cohort.
Oncoprint representation from cBioPortal revealed the distribution of GSDM genomic alterations in the Pan-Cancer (Fig. 1A). As a result, GSDM family genes were altered in 1588 (14%) of 10,953 patients in 32 TCGA studies, and the overall alteration frequency of each gene in the GSDM family was 1.4–8%. GSDMC (8%) and GSDMD (6%) had the overall higher gene alteration frequency than GSDMA (3%), GSDMB (3%), GSDME (2%) and PJVK (1.4%). CNAs played a predominant role in gene alteration frequency in GSDMD, GSDMA, GSDMB and GSDMC, while it played an approximate equal role as mutation in GSDME. Noticeably, PJVK had more mutations other than CNAs, setting it apart.
The overall average somatic mutation frequency of each gene in the GSDM family was 0.5–1.3% (Fig. 1B). The majority of mutations were almost distributed evenly in the coding sequence of GSDM genes. The most common mutations of GSDM genes were missense and truncating mutations. GSDMC (1.3%) had the highest mutation rate among the six members, while GSDMD (0.5%) had the lowest. In GSDMC, there were 135 Missense, 21 Truncating and 7 splice mutations and the hotspot mutations (R40Q/*, R34H/C, R7C/H) were all located within the Gerdermin domain. In GSDMD, there were 54 Missense, 8 Truncating and 2 splice mutations and the hotspot mutations (E15K/Q) were also located within the Gerdermin domain.
Of the 57 cancer types, 53(93%) had genomic alterations in GSDM family genes (Fig. 1C). Genomic alteration frequencies varied widely across cancer types, from no alterations in embryonal carcinoma, synovial sarcoma, paraganglioma and pleural mesothelioma (epithelioid), to 44.83% alteration frequency in esophageal adenocarcinoma. Genomic alteration type also varied greatly across tumor types, from no alterations or single type in acute myeloid leukemia, pheochromocytoma, chromophobe renal cell carcinoma, mixed germ cell tumor, seminoma, myxofibrosarcoma, pleural mesothelioma (biphasic) and endocervical adenocarcinoma, to multiple type alteration in most other cancers. Among the 57 types, esophageal adenocarcinoma (44.83%), serous ovarian cancer (40.41%) and uterine serous carcinoma (36.7%) had the highest frequency of GSDM genes alteration, while chromophobe renal cell carcinoma (1.54%), pheochromocytoma (1.36%) and acute myeloid leukemia (0.5%) had the lowest. Serous ovarian cancer had the highest amplification frequency (38.01%). Cutaneous melanoma had the highest GSDM genes mutation frequency (15.99%). GSDMA mutation was more likely associated with intrahepatic cholangiocarcinoma, uterine mixed endometrial carcinoma and cutaneous melanoma. GSDMB mutation was more likely associated with cutaneous melanoma. GSDMC and GSDMD mutations were more likely associated with mucinous stomach adenocarcinoma. GSDME and PJVK mutation were related to uterine endometrioid carcinoma and uterine mixed endometrial carcinoma, respectively.
In most cancer types, amplification was the major contributor to the GSDM gene alteration (Fig. 1C). Esophageal adenocarcinoma had the highest frequency of GSDMA (20.69%) and GSDMB (21.84%) amplification. Serous ovarian cancer had the highest frequency of GSDMC (31.51%) and GSDMD (26.71%) amplification. Mucinous carcinoma and serous ovarian cancer had the highest frequency of GSDME (5.88%) and PJVK (2.91%) amplification, respectively. Intrahepatic cholangiocarcinoma, cutaneous melanoma, uterine endometrioid carcinoma and uterine mixed endometrial carcinoma had the highest GSDMA (6.67%), GSDMB (2.7%), GSDME (6.52%) and PJVK (4.76%) mutation frequency, respectively. Mucinous stomach adenocarcinoma had the highest GSDMC (9.09%) and GSDMD (4.55%) mutation frequency.
Expression of GSDM genes extensively changed in cancer
We explored the expression levels of the six GSDM family members in all 33 cancer types available in TCGA pan-cancer data. In pan-cancer, GSDMD had the relatively highest expression level while GSDMC had the lowest (Fig. 2A). Overall, GSDM family genes tended to be upregulated compared with their adjacent normal tissue (Fig. 2B, D–I). However, a noticeable intra- and inter-cancer heterogeneity regarding the expression levels of the corresponding genes was unveiled for all six GSDM members. There was no intrinsic unified pattern in every single GSDM gene expression. In addition, expression levels of GSDM family members were significantly (P < = 0.001) correlated with each other in pan-cancer based on Pearson's correlation tests (Fig. 2C). However, the absolute value of the correlation coefficient ranged from 0.03 to 0.38, testifying a weak or negligible correlation. These findings revealed the intrinsic differences in the expression of GSDM genes between different cancer types or between different GSDM family members. Whether a specific GSDM gene was an oncogene or an anti-oncogene cannot be decided without choosing a cancer type. The complexity of expression spectrum necessitated a further study of each individual GSDM gene member as an entity.
Genetically, the copy number aberrations (CNAs) might be a considerable molecular mechanism leading to the disturbance of gene expression of GSDM family (Additional file 3: Fig. S3). Overall, the expression level of GSDMs was upregulated in gain and amplification but downregulated in shallow deletion and deep deletions in pan-cancer. Regarding the widely amplifications and relatively few deep deletions in GSDM genes, the overall upregulated expression of corresponding gene seemed to be reasoned.
Epigenetically, we identified the effects of promoter methylation on GSDM family gene expression in 33 types of cancer by integrating methylation levels and expression profile data (Fig. 2J). 14 (42.4%) of 33 types of cancer were significantly (FDR < = 0.05) differential methylated in GSDM gene promoters. Overall, GSDM promoters were hypomethylated in variety of cancer and the expression of GSDM genes was negatively correlated with promoter methylation level (Fig. 2K). This finding corresponded to the overall up-regulated GSDM genes expression.
GSDM genes involving widely in cancer-related pathways and drugs sensitivity
We investigated the correlation between the expression of GSDM genes and ten famous cancer related pathways, including TSC/mTOR, RTK, RAS/MAPK, PI3K/AKT, Hormone ER, Hormone AR, EMT, DNA Damage Response, Cell Cycle and Apoptosis pathways. As a result, the expression of GSDM family genes was correlated with the activation or inhibition of a variety of cancer related pathways (Fig. 3A). However, the activating or inhibiting role on the cancer related pathways of GSDM genes varied depending on types of cancer, also, types of pathways. Noticeably, the high expression of GSDME was significantly (FDR < = 0.05) correlated with EMT activation in 41% of cancer types.
We further investigated the correlation between drug sensitivity and GSDM family genes expression level in CellMiner [18]. We found that the expression level of GSDM family genes was significantly correlated (P < 0.05) with drug sensitivity in 309 (41.3%) of 748 (under FDA approved or clinical trial) chemical compounds and natural products (Fig. 3C, Additional file 8: Data S1). GSDMA, GSDMB, GSDMC, GSDMD, GSDME and PJVK were correlated with sensitivity of 34, 31, 49, 39, 88 and 37 drugs, respectively. Obviously, the expression of GSDME was more likely to be associated with drug sensitivity. The extent or the direction of the correlation was not unified regarding every GSDM gene and every drug. For example, the high expression of GSDME increased drug sensitivity of Telatinib, P–529 and IDH–C227, but increased drug resistance of Vinorelbine and Des-fluoro-TAK-960. Interestingly, GSDMA and GSDMB generally went to the same direction of correlation when associated with drug sensitivity of a specific drug (Fig. 3B, Additional file 9: Data S2), but the opposite direction of GSDME. For example, the drug sensitivity of PI-103 was negatively correlated with the expression of GSDMA (r = − 0.35) and GSDMB (− 0.43), but it was positively correlated with the expression of GSDME (r = 0.31). The drug sensitivity of AM-5992 negatively correlated with GSDME (r = − 0.41), but it positively correlated with the expression of GSDMA (r = 0.39) and GSDMB (r = 0.36). These results implied the possibility that GSDMA and GSDMB were involved in some pathways together with GSDME, but they played a different role from GSDME. In addition, GSDMC exhibited a great difference in spectrum of correlation with drug sensitivity from the rest of the GSDM family members, indicating a different role GSDMC might play in regulation of drug sensitivity.
We next examined the protein–protein-interaction (PPI) network of GSDM family genes. Through the visualized networks, we found no interaction within GSDM genes, but GSDM genes interacted with several proteins known to be involved in regulation of cell proliferation, transcription and drug resistance (Fig. 3D). These findings implied that the GSDM family genes were widely involved in the regulation of biologic behaviors of cancer and reactions to drugs. However, there was a nonnegligible heterogeneity in the specific mechanism across cancer types or GSDM gene types.
Patient survival was associated with the expression of GSDM genes
Since the extensive alteration and differential expression of GSDM genes among cancers, we further investigated the association of GSDM alteration and expression with patient overall survival (OS) and progression free survival (PFS). We found that altered GSDM genes was significant associated with worse OS and PFS in pan-cancer (Fig. 4A). Univariate Cox proportional hazards model demonstrated that the expression level of GSDM genes widely associated with the survival risk (Fig. 4B, Additional file 10: Data S3). However, whether a specific GSDM gene was a high or low survival risk gene varied depending on types of cancer. For instance, high expression of GSDME increased survival risk in KIRC, LIHC and UCEC, but decreased survival risk in ACC. According to KM survival curves (Fig. 4C), OS benefited with high expression of GSDMB in BLCA and SKCM, but with low GSDMB expression in KIRC. High GSDMC expression associated with survival advantage in COAD and LGG, but with survival disadvantage in KIRC, KICH and PAAD. High expression of GSDMD predicted a better OS in BLCA and SKCM, but a worse OS in ACC, LGG and UVM. High expression of PJVK also associated with a better OS in LAML, MESO, PAAD and SARC, but with a worse OS in KIRC and COAD. Interestingly, the expression of GSDMA did not affect OS significantly. In brief, the expression of GSDM genes was extensively involved in patient survival, but it varied by cancer types.
GSDM genes were associated with immune subtypes and tumor microenvironment (TME)
Since the GSDM genes were a family of pore-forming proteins implicated in the immune response, we investigated the correlation between GSDM genes and immune infiltrates in cancer to understand how each of the GSDM family members was associated with immune components. A global transcriptomic immune classification of solid tumors had identified six immune subtypes (ISs) (C1-C6) [19]. In pan-cancer, patients classified into type C3 (inflammatory) and C5 (immunologically quiet) associated with significant survival advantage while patients characterized into type C4 (lymphocyte depleted) and C6 (TGFβ dominant) had a significant survival disadvantage (Additional file 4: Fig. S4). Further analysis demonstrated the expression of GSDM family genes was significantly (P < 0.001) correlated with immune subtypes (Fig. 5A). High expressions of GSDMA, GSDMB and GSDMC were associated with subtype C1 (Wound healing), C2 (INF-r dominant) and C6, indicating a cancer promoter role of these three members. High expressions of GSDME and PJVK were associated with subtype C5, indicating these two genes might mainly play a cancer suppressor role. However, GSDMD might play a dual role, since high expression of GSDMD correlated significantly with patients characterized into C3 as well as C4 and C6.
We further investigated the correlation of GSDM genes with ImmuneScore, StromalScore and EstimateScore which estimated the level of immune cell infiltrate, stromal cells infiltrate and tumor purity respectively using the ESTIMATE program [20] (Fig. 5B, Additional file 11: Data S4, Additional file 12: Data S5). The higher ImmuneScore or StromalScore were, the more extensive the immune or stromal components infiltrated in TME. The EstimateScore, negative correlation with tumor purity, was the sum of StromalScore and ImmuneScore. As a result, the expressions of all GSDM genes but PJVK were mainly positively correlated with ImmuneScore. Specifically, GSDMD was strong correlated with LGG (r = 0.68) and PCPG (r = 0.66), and GSDMA was strong correlated with ImmuneScore in ACC (r = 0.69) and SARC (r = 0.60) (Additional file 11: Data S4). PJVK was negatively correlated with ImmuneScore significantly in LAML (r = − 0.5). Regarding StromalScore, the expressions of all GSDM genes but GSDMB and PJVK were mainly positively correlated with StromalScore. To be specific, GSDMD was strong correlated with StromalScore in PCPG (r = 0.76) and LGG (r = 0.64). The most striking was that the expression of GSDME was very strong correlated with StromalScore in COAD (r = 0.84) and PRAD (r = 0.83) and it had a strong correlation in READ (r = 0.72), KICH (r = 0.66) and BRCA (r = 0.60) with StromalScore. In spite of the relatively small coefficient, the expression of GSDMB (|r|< 0.37) and PJVK(|r|< 0.47) were mainly negatively correlated with StromalScore. Not surprisingly, the high expressions of GSDMA, GSDMC, GSDMD and GSDME were observed correlated with low tumor purity (Additional file 5: Fig. S5). Meanwhile, the high expression of PJVK indicated a high tumor purity in most types of cancer. Stromal and immune cells were proposed to play an important role in cancer growth, metastasis and drug resistance [20], indicating GSDM family genes might have a role in regulating tumor behavior by interacting with TME.
GSDM genes were associated with tumor cell stemness
To evaluate the relationship between GSDM genes and tumor stemness, we investigated the correlation between the expression of GSDM genes and stemness score (RNAss and DNAss). As a result, the association between GSDM genes expression and tumor stemness varied depending on the types of cancer and types of GSDM genes (Fig. 5C, Additional file 13: Data S6). To be specific, RNAss was strongly correlated with GSDMD in LGG (r = -0.66) and PCPG (r = − 0.67), strongly correlated with GSDME in PRAD (r = − 0.75). DNAss was found strongly correlated with GSDMB in OV (r = − 0.71). Besides these negative correlations, GSDMD was found positively associated with RNAss in THYM (r = 0.61), with DNAss in OV (r = 0.64). In spite of the different results estimated by RNAss and DNAss based on their algorithm differences, GSDM genes were found to be correlated with tumor stemness in varying degrees.
GSDM genes were associated with colorectal cancer cell invasion and metastasis
Next, we investigated the genes associated (P < 0.05) with GSDM genes in cell lines in CCLE using Venn diagrams [21]. Interestingly, in colorectal cancer cell lines, we found the genes positively correlated with GSDMA had no intersection (none gene) with the genes negatively correlated with GSDMB, but it had intersection (1022 genes) with the genes positively correlated with GSDMB (Fig. 6A). Also, the genes negatively correlated with GSDMA had almost no intersection (only one gene) with the genes positively correlated with GSDMB, but it had intersection (672 genes) with the genes negatively correlated with GSDMB. When it came to GSDME, we found a similar relationship. The genes positively correlated with GSDME had almost no intersection with the genes positively correlated with GSDMA (none gene)/GSDMB (two genes), but it had intersection with the genes negatively correlated with GSDMA (609 genes)/GSDMB (863 genes). In turn, the genes negatively correlated with GSDME had almost no intersection with the genes negatively correlated with GSDMA (four genes)/GSDMB (none gene), but it had intersection with the genes positively correlated with GSDMA (601 genes)/GSDMB (1002 genes). And that, the genes positively correlated with GSDMD had interaction (291 genes) with the genes negatively correlated with GSDME, but had almost no interaction (18 genes) with the genes positively correlated with GSDME. In brief, we suspected that GSDMA, GSDMB, GSDMD and GSDME shared some biofunctions and the direction of regulation of these biofunctions by GSDMA and GSDMB was same to GSDMD but the opposite to GSDME.
We further performed GO enrichment analysis on the GSDMA− GSDMB− GSDMD− GSDME+ gene list (negatively correlated with GSDMA, GSDMB and GSDMD but positively correlated with GSDME) and GSDMA+ GSDMB+ GSDMD+ GSDME− gene list (positively correlated with GSDMA, GSDMB and GSDMD but negatively correlated with GSDME). The GO analysis results revealed that the former list significantly (p-value < 0.05 and q-value < 0.05) enriched in biological processes related to axonogenesis and morphogenesis (shape and size of cancer cells), while the latter gained none enrichment in biological process. Also, the majority of these genes in the two lists were located at the leading edge or projection of cancer cells. Since cell morphological changes were the prerequisite for cell migration and abnormal cell migration underlined invasion and metastasis of cancer cells [22], we suspected that GSDMA, GSDMB, GSDMD and GSDME might be involved in colorectal cancer cell invasion and metastasis, and the former three might played as negative regulators while the last one played as a positive regulator.
To examine the role of GSDME in cell migration and invasion, we stably expressed GSDME in HCT116 and SW480 cells (Fig. 7A). The transwell assays revealed that HCT116-GSDME and SW480-GSDME group exhibited more invasive (SW480: t = 4.218, p = 0.0135; HCT116: t = 9.627, p = 0.0007) and migrating (SW480: t = 5.003, p = 0.0075; HCT116: t = 6.354, p = 0.0031) than the control group (Fig. 7B, C). The wound-healing assay indicated that over expression of GSDME in HCT116 (24 h: t = 4.312, p = 0.0125; 48 h: t = 6.394, p = 0.0031, 72 h: t = 15.24, p = 0.0001) and SW480 cells (24 h: t = 3.726, p = 0.02; 48 h: t = 4.141, p = 0.014, 72 h: t = 12.55, p = 0.0002) increased the migratory rate (Fig. 7D).
GSDME was involved in angiogenesis in colorectal cancer (CRC)
During the drug sensitivity analysis, we found that the some GSDME related drugs were angiogenesis inhibitors, such as Telatinib (r = 0.61) and Lenvatinib (r = 0.45) (Additional file 8: Data S1). In addition, during the tumor microenvironment estimation, we found that expression of GSDME was positively correlated with StromalScore in CRC(r = 0.84) (Additional file 11: Data S4). Based on this, we suspected that GSDME might related to angiogenesis in CRC tissue. Therefore, further investigation was taken into the correlation between expression of GSDME and endothelial cell (EC) markers using GEPIA [23]. As a result, the mRNA level of all the EC markers showed positively correlation with the mRNA level of GSDME in CRC tissue (Fig. 6B, Additional file 5: Fig. S5A). Gene Set Enrichment Analysis (GSEA) [24, 25] highlighted that the GO Endothelium Development Gene Set upregulation was related to increased expression of GSDME (Fig. 6C).
We doubted whether GSDME promoted the expression of vascular endothelial growth factor (VEGF) in CRC cells. We next analyzed the association between expression of GSDME and the VEGF in Cancer Cell Line Encyclopedia (CCLE). As expected, the expression of GSDME was positively correlated with VEGFB (r = 0.46) and VEGFC (r = 0.40) (Additional file 6: Fig. S6). Hence, we suspected that GSDME was positively correlated with the production of VEGF by cancer cell. Also, the expression level of GSDME associated (Fvalue = 2.913, Pr = 0.034) with stages of CRC (Fig. 6D). And as a contrast, we also analyzed GSDMC and found a positive correlation between the mRNA level of GSDMC and mRNA level of EC marker, but no significant correlation has been found between GSDMC and VEGF in CCLE (Additional file 7: Fig. S7). Since VEGF promoted cancer angiogenesis and metastasis [26], we suspected that GSDME might be involved in CRC growth and metastasis via VEGF related pathways in CRC.
To determine the expression of GSDME and angiogenesis, IHC assay was performed using 16 CRC tissue (Fig. 7E). The number of micro vessels per unit area (micro vessel density, MVD) was visualized by CD34 and revealed that angiogenesis was promoted in higher GSDME expression CRC tissue (t = 2.278, p = 0.039).