We downloaded the gene expression profiling dataset created by high-throughput sequencing (GSE37815) and the microarray-based gene methylation profiling dataset (GSE37817) from the Gene Expression Omnibus database (GEO, https://www.ncbi.nlm.nih.gov/geo/). In total, five normal bladder mucosae and 18 primary BCa samples were included in GSE37815 (platform: GPL6102 Illumina human-6 v2.0 expression beadchip). As for the DNA methylation datasets, GSE37817 included six normal bladder mucosae and 18 primary BCa samples based on the GPL8490 platform (Illumina HumanMethylation27 BeadChip).

To identify the potential prognostic hub genes of the MeDEGs, we performed GEO2R (http://www.ncbi.nlm.nih.gov/geo/geo2r/) to compare two or more groups of samples in a GEO Series to screen genes that are differentially expressed across specific experimental conditions. In the present study, GEO2R was used to identify the differentially expressed genes (DEGs) and differentially methylated genes (DMGs). |t| >2 and P <0.05 were considered statistically significant. Furthermore, hypomethylation-high expression genes were obtained after the overlap of upregulated and hypomethylated genes, and hypermethylation-low expression genes were obtained after the overlap of downregulated and hypermethylated genes. The hypomethylation-high expression genes and hypermethylation-low expression genes were then identified as methylation-regulated differentially expressed genes (MeDEGs).

To obtain the functional annotations of hub gens, we utilized the Database for Annotation, Visualization, and Integrated Discovery (DAVID, https://david.ncifcrf.gov/) is a straightforward web tool that can provide integrative and systematic annotation for users to unravel the biological interactions of multiple genes. It was utilized to perform functional and pathway enrichment analyses. Gene ontology (GO) analysis, including the biological process (BP), cellular component (CC), molecular function (MF), and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis, were conducted for the selected MeDEGs by DAVID (28, 29). P <0.05 was considered statistically significant.

In this study, STRING (http://string-db.org; version 11.0) was adopted to describe protein co-regulation of hypomethylation-high expression genes and hypermethylation-low expression genes, respectively, and measure functional interactions among nodes (30). The interaction specificity score above 0.4 (the default threshold in the STRING database) was considered statistically significant. Cytoscape (version 3.6.0) was used to visualize interaction networks obtained from STRING (31). MCODE (version 1.4.2) of Cytoscape is a plug-in to cluster a given network to identify densely connected regions based on topology (32). It was utilized to find the most related module network with selection threshold as follows: MCODE scores >5, degree cutoff = 2, node score cut-off = 0.2, Max depth = 100 and k-score = 2.

Gene Expression Profiling Interactive Analysis (GEPIA, http://gepia.cancer-pku.cn/) is an online tool that can provide customizable functionalities based on data from The Cancer Genome Atlas (TCGA; https://tcga-data.nci.nih.gov/tcga/) and the Genotype-Tissue Expression project (GTEx; https://www.gtexportal.org/home/index.html) (33). GEPIA performs survival analysis based on gene expression levels, using a log-rank test for the hypothesis evaluation. The horizontal axis (x-axis) represented the time in days, and the vertical axis (y-axis) showed the probability of surviving or the proportion of people surviving. The cut-off value was defined via median value or using “survminer” R package. The lines presented the survival curves of the two groups.

Based on data from the TCGA database, the GSEA tool (version 2.10.1 package) was used to predict associated up- and down-regulated genes and their significantly involved hallmarks pathways (34). The student’s t-test statistical score was implemented in consistent pathways, and the mean of the DEGs was calculated for each analysis. A permutation test 1000 times was utilized to recognize the significantly involved pathways. The adjusted P using Benjamini and Hochberg (BH) and false discovery rate (FDR) method by default were used to correct for the occurrence of false positive results. Significantly related genes were defined with an adjusted P <0.01 and FDR <0.25.

BSAS methylation next-generation sequencing of BCa cell lines was conducted by GeneChem Biotechnology Co., Ltd., Shanghai (GSGC0257632). Microsoft Office Excel software and Methylation Plotter software were used to examine the results. The Kruskal-Wallis test is a nonparametric test of three or more groups of data and was used when the ANOVA test could not be utilized. Bisulfite PCR amplification was performed using the High Pure PCR Template Preparation Kit (Roche) with forward and reverse primers provided. DNA methylation libraries were developed using the VAHTS Turbo DNA Library Prep Kit for Illumina^® (ND102-0102).

GEO2R was adopted to identify the DEGs and DMGs, respectively. For DEGs of gene expression microarray, 2425 overlapping up-regulated genes and 2563 overlapping down-regulated genes were screened ( Figures 2A, B ). A total of 1114 overlapping hypermethylation genes and 10,024 overlapping hypomethylation genes were discovered for DMGs of gene methylation microarray. The study identified 1162 hypomethylated, upregulated genes and 203 hypermethylated, downregulated genes after integrating the DEGs and DMGs ( Figures 2C, D ).

The PPI network of hypomethylation-upregulated genes and hypermethylation-downregulated genes was visualized using Cytoscape (version 3.6.0) [28]. A Cytoscape plug-in called MCODE (version 1.4.2) clusters a given network to select densely connected regions based on topology [29]. The results are presented in Figures 2E, F . As a result, in the hypomethylation-upregulated genes module, CXCL10, CXCL16, CX3CR1, SSTR1, C3AR1, CNR1, ADCY2, BDKRB2, POMC, GPR55, GNG12, GNG11, GNGT2, CCL4, C5AR1, GNB4, CXCL13, GNG7, GPR18, APP, FPR1, GPR17, and GPR18 were confirmed as hub genes. While in the hypermethylation-downregulated genes module, PTGDR, ADCY9, OXER1, GALR1, ADCYAP1, ADCY8, ADCY7, VIPR2, GNG10, GNB4, CCR10, GNG4, and CXCR3 were confirmed as hub genes.

For hypomethylation-upregulated genes, changes in biological processes were mostly enriched in angiogenesis, signal transduction, aging, and immune response. The hypermethylation-downregulated genes were primarily enriched in extracellular matrix organization, signal transduction, cell adhesion, cAMP-mediated signaling, and cellular response to glucagon stimulus. Moreover, the study found that the hypomethylated, upregulated genes were associated with extracellular exosome, plasma membrane, and extracellular region. Whereas the hypermethylated, downregulated genes were associated with the proteinaceous extracellular matrix, plasma membrane, and extracellular matrix in the cellular component group. For hypomethylated, upregulated genes, changes in molecular function were primarily enriched in protein binding, heparin-binding, actin filament binding, and extracellular matrix structural constituent. On the other hand, for hypermethylated, downregulated genes, changes were significantly enriched in collagen binding, phosphorus-oxygen lyase activity, and extracellular matrix binding. Furthermore, pathway enrichment was performed using KEGG. The study revealed that hypomethylated genes predominantly participated in morphine addiction, retrograde endocannabinoid signaling, and cholinergic synapse. For hypermethylated genes, the most significantly enriched pathways involved focal adhesion, pathways in cancer, and the PI3K-Akt signaling pathway.

Significant survival outcomes of hub genes in the PPI network are displayed in Figure 3 . According to the expression of each gene, overall survival for SKCM patients was acquired. The study found that high mRNA expression of ADCY2 (P = 0.047) was significantly associated with worse prognosis for SKCM as well as APP (P = 0.022), BDKRB2 (P = 0.029), FPR1 (P = 0.025), GNB4 (P = 0.013), GNG11 (P = 0.011), ADCY9 (P = 0.036), and ADCYAP1 (P = 0.014). Significant genes and pathways were obtained using GSEA.

After screening and analysis, this study found that ADCY2 was significantly associated with prognosis. A total of 34 cancers were then selected to observe the difference in ADCY2 expression between tumor and normal tissues. In BCa, the expression level in tumor tissues was significantly lower than that in normal tissues (n = 435, Figure 4A ). Hence, a separate survival analysis of BCa patients and male and female subgroups was performed. The findings revealed that the high expression of ADCY2 predicted a worse prognosis, and it was significant in female patients, which might be due to a higher incidence of BCa in males than in females ( Figures 4B, C ).

Simultaneously, the significantly down-regulated genes were listed, including MFAP4, LMOD1, CNN1, COMP, SFRP4, and so on ( Figure 5A ). The correlation trend of these genes in populations with high and low ADCY2 expression was observed ( Figure 5B ).

In the KEGG pathways, ADCY2 was found to be primarily enriched with protein absorption, PI3K-Akt signaling pathway, and focal adhesion pathway. An enrichment of ADCY2 was also observed in the BCa pathway, which showed that ADCY2 is involved in the development of BCa. Meanwhile, in GO term, ADCY2 was mainly enriched in extracellular structure organization and extracellular matrix organization pathways ( Figure 5C ).

Several immune-related cytokines were screened, mainly from five families, for correlation analysis with ADCY2 in pan-cancer. ADCY2 was found to be associated with most cytokines and presented a tumor immune microenvironment dominated by MHC and chemokine in most cancer types ( Figure 6A ). Notably, an immune infiltration analysis was performed in BCa using EPIC and CIBERSORT algorithms. T cell CD4+ memory resting, B cell naive, B cell memory, and macrophage M2 were significantly clustered in the ADCY2^high group, while uncharacterized cells were significantly increased in the ADCY2^low group ( Figures 6B, C ).

Next, the study screened for immune checkpoint molecules. The association between ADCY2 and some immune checkpoints in pan-cancer was first explored, including those that promote immunotherapy and inhibit the efficacy of immunotherapy ( Figure 7A ). The immune molecules, such as HAVCR2, PDCD1LG2, and TIGIT, which were significantly clustered in the high and low ADCY2 groups of BCa, were then studied independently ( Figure 7B ). However, the difference between the two groups was not significant. Finally, TIDE analysis was performed in the two BCa sample groups. The high score of the ADCY2^high group suggested a high expression of ADCY2 associated with tumor heterogeneity, indicating a worse immunotherapy effect and prognosis for BCa patients ( Figure 7C ).

In human BCa 5637 and MGH1 cells, the study developed and verified the effectiveness of ADCY2 primers using BSAS technology. Figure 8A displays the three pairs of primers for ADCY2 methylation, where the CG site of the amplified fragment is indicated in bold and red. The schematic diagram of the average methylation information of all sites in the 5637 and MGH1 cells was developed using the Methylation plotter ( Figure 8B ). Besides, the methylation levels of all sites in the grouped samples are shown in boxplots and dendrograms ( Figures 8C–E ). The findings revealed that the expression of ADCY2 is highly regulated by the methylation of the promoter regions.