An Analysis Regarding the Association Between Connexins and Colorectal Cancer (CRC) Tumor Microenvironment

Introduction

Colorectal cancer (CRC) is a tumor of the gastrointestinal system, ranking third in incidence among malignancies and second in lethality among cancers.¹ Its pathogenesis appears to involve multiple factors and steps, including genetics, diet, and lifestyle habits.^2,3 Statistics show that the tumor frequently recurs and metastasizes, spreading especially to the liver, lung, bone, and ovaries, and the five-year survival rate for such patients is approximately 10%.^4,5 Treatment methods for CRC include chemotherapy, surgery, radiotherapy, and molecular-targeted therapy.^6,7 Molecular-targeted therapy for CRC has become a hot topic.⁸ The number of molecularly stratified treatment options is increasing, as is the use of biomarkers to guide prediction and treatment decisions.^9,10 Therefore, the identification of novel markers and targets for treating CRC is necessary for improving the disease prognosis.

Gap junction protein (Connexin) is a four-pathway transmembrane protein involved in the assembly of gap junctions among adjacent cells, thereby facilitating gap junction cell-to-cell communication (GJIC), which allows small molecules and ions to diffuse freely directly from the cytoplasm of neighboring cells without encountering the extracellular environment, resulting in a rapid and coordinated cellular response to internal and external stimuli.^11,12 It is well established that connexins are often considered tumor suppressive due to reduced expression or altered cytoplasmic localization of connexin in primary tumor cells, resulting in GJIC deficiency.¹³ However, enhanced expression levels of connexins have also been reported in metastatic lesions, including connexin 43 (Cx43, GJA1) and connexin 26 (Cx26, GJB2), which may function to mediate the promotion of cancer cell migration and adhesion.^14,15 An interesting study in 2016 found that Cx43 localized on fibroblasts could help tumor cells expel excess acid to maintain a favorable ecological environment, which in part reinforces the complex and variable role of connexins.¹⁶ However, to date, studies on the role of connexins in CRC are still very scarce. The present study seeks to fill this lacuna.

In this study, human colorectal cancer specimens were analyzed by immunohistochemistry and immunofluorescence methods. Meanwhile, the potential functions and mechanisms of action of connexins, especially GJA4, were initially explored by using bioinformatics software and methods.

Materials and Methods

Reagents

The reagents and antibodies used are listed in the Supplementary Materials (Supplementary Table 1). The concentrations used were based on those used in the previous research or on the manufacturer’s recommendations. Experimental procedures are also provided in the Supplementary Material.

Data Sources and Process

The workflow of our research is shown in Supplementary Figure 1. A total of 620 CRC samples, including gene expression profiles and matching clinical data as well as mutation data, were downloaded through the University of California, Santa Cruz (UCSC) browser. The gene expression profiles (FPKM values) of The Cancer Genome Atlas (TCGA) – Colorectal adenocarcinoma (COAD&READ) dataset was processed through R software to transform into Transcripts Per Kilobase million (TPM), which is closer to the microarray data.¹⁷ We then selected the Gene Expression Omnibus (GEO) dataset GSE44861 and GSE10950 for expression and immune cell infiltration analysis.^18–20 Gene expression data were extracted using R software (version 4.1.1), and data matrix was constructed for further analysis.

Establishment of Consistent Clustering Based on Connexins

Consensus clustering algorithm was performed to estimate the number of unsupervised classes in TCGA-COADREAD. The TCGA-COADREAD cases were classified into 2 clusters based on the expression profile of the connexins by “consensusclusterplus” R package. The procedure was repeated 500 times to ensure the stability and reproducibility of the classification. Principal component analysis (PCA) was used to determine the validity of clustering.²¹

Least Absolute Shrinkage and Selection Operator (LASSO) and COX Analysis

The LASSO regression analysis was conducted through the R package “glmnet” to determine the research objective.²² We conducted the univariate and multivariate Cox regression analysis to access the prognostic values of connexins coding genes. The forest plots representing the P-value, hazard ratio (HRs), and 95% confidence interval (CIs) of all variables were visualized with the “forestplot” R package.

Expression and Survival Analysis

The relationship between RNA sequences and clinical characteristics in TCGA was visualized through the “survival” R package and MEXPRESS database, mainly including Overall Survival (OS), Disease Specific Survival (DSS), Progression Free Interval (PFI), Progression Free Survival (PFS), Disease Free Survival (DFS), and “T (Tumor)”, “N (Node)”, “M (Metastasis)” stages.²³ And then, TCGA expression profile and Gene Expression Omnibus (GEO) data were used to assess the differential expression. In addition, the Human Protein Atlas (HPA) and our own immunohistochemical data were used to validate the expression of GJA4 at the protein level.²⁴

Ethics Statement and Specimen Collection

The study’s protocol was approved by the ethics committee of the Jiangsu Province Hospital of Chinese Medicine, and informed consent was obtained from clinicians and patients (2020NL-107-01). This study complied with the Declaration of Helsinki. The CRC tissue and the surrounding non-tumorous tissue (margin, 5 cm) were collected during surgery from 30 previously treatment-naïve patients with CRC at the Jiangsu Provincial Hospital of Traditional Chinese Medicine. Inclusion criteria: (1) All patients were first pathologically diagnosed with CRC. (2) The diagnosis was referenced to the latest (8th edition) American Joint Committee on Cancer (AJCC) CRC staging system.²⁵ (3) All patients had complete clinical information. After extraction, tissue specimens were washed with cold phosphate-buffered saline and immediately placed in liquid nitrogen. They were then transferred and stored at −80°C until further examination. Immunohistochemistry (IHC) and Tissue Section Immunofluorescence (IF) staining were performed as previously described.^26,27 The precise details are given in Supplementary Materials.

Enrichment Analysis

The “limma” R package was used to obtain the Differentially Expressed Genes (DEGs) between C1 and C2 clusters and genes co-expressed with GJA4. The “clusterProfiler” R package and GSCALite web tool were used to perform the functional enrichment analysis. The protein–protein interaction (PPI) was constructed by STRING online database and visualized by Cytoscape.^28,29 The Cytoscape cytoHubba plugin was used to search for hub genes by degree order, with default parameters being used as parameter settings.³⁰

The gene set enrichment analysis (GSEA) was performed through R 4.1.1.³¹ The gene set “subset of Gene Ontology (GO)” and “subset of hallmark” were downloaded from the Molecular Signatures Databases and were used as reference gene sets.³² The False Discovery Rate (FDR) <0.1 was considered to be with statistical significance.

Single Cell Analysis

First, we conducted a single cell analysis based on GSE146771 by Tumor Immune Single Cell Hub (TISCH) database to explore in which cell types GJA4 may express.^33,34 Then, we downloaded an independent dataset from scTIME Portal, which consists of 1 CRC patient’ s cells. All samples were imported into Seurat V3 and visualized in UMAP after a standardized quality control process. Next, fibroblasts were distinguished from all cells. The “FindAllMarkers“ function (min.pct = 0.25, logfc.threshold= 2, tes.use = ”wilcox”) was set to find markers for all cell types. ‘Dotplot’ function was used to plot dot plots, and ‘Vlnplot’ function was used to plot violin plots.

Genetic Analysis

We conducted the “maftools” package to assess the frequency of mutations in 21 connexins and the Gene Set Cancer Analysis (GSCA) database to calculate the relationship between connexin genetic alterations and expression levels.³⁵

Immune Analysis

The CIBERSORT and ssGSEA algorithms were used for the overall assessment of immune infiltration score, and the “estimate” package was used to calculate immune and stromal scores as well as tumor purity.^36,37 The Tumor Immune Estimation Resource (TIMER) database was used for reliable assessment of fibroblast infiltration, and pairwise comparisons were quantified by Spearman’s rank correlation, with P < 0.05 representing significance. All results were processed using the ‘ggplot2’ and ‘pheatmap’ packages in R.³⁸

Immunohistochemical and Immunofluorescence Staining

Immunohistochemical and Immunofluorescence Staining and H-sore procedures were performed as described previously.^39,40 A detailed description is shown in Supplementary Materials.

Results

The Genetic Characteristics and Transcriptional Variations of 21 Connexins

We collected the currently available research data on 21 connexins including GJA1, GJA3, GJA4, GJA5, GJA8, GJA9, GJA10, GJB1, GJB2, GJB3, GJB4, GJB5, GJB6, GJB7, GJC1, GJC2, GJC3, GJD2, GJD3, GJD4, and GJE1. Figure 1A demonstrates that gap junctions among tumor cells are potent regulators of EMT, while little has been reported about the role of gap junctions among non-tumor cells. We first showed various relationships, including interactions, colocalization, and common pathways among connexins (Figure 1B). In CRC tissue, most connexins activated the EMT signaling pathways but inhibited the cell cycle and DNA damage response (Figure 1C). Based on the mRNA data of paired tumor and normal samples present in TCGA-COADREAD, several genes were found to be highly expressed in CRC tissues: GJA3, GJA4, GJB1, GJB3, GJB4, GJB5, GJB6, GJC3, and GJD2 (Figure 1D). We then confirmed that in CRC, genetic variation and DNA methylation were important factors affecting the expression levels of connexins, and that Copy Number Variations (CNVs) and mRNA expression levels were positively correlated, while gene methylation levels and mRNA expression levels were negatively correlated (Figure 1E). In addition, for some of the connexins, these two factors were also able to influence the prognosis of CRC patients (Figure 1F). In the TCGA database, there was a significant co-expression relationship between GJA1 and GJA4, GJA5, GJC1, and GJD3 (P < 0.05, Figure 1G). Out of 544 samples, 96 samples had mutated connexins, with a frequency of 17.65%. Among these 96 samples, GJA8 had the highest frequency of mutations, mainly missense mutations (Figure 1H). These results strongly suggested that imbalance of connexins can lead to CRC.

Figure 1 Expression variation of connexins. (A) pattern map shows the molecular regulatory mechanism of connexins and their potentially important regulatory role in EMT of tumor/non-tumor cells. (B) Connexins with neighboring genes showing physical interactions, co-expression, co-localization, predicted common pathways, genetic interactions, and common protein domains. (C) The correlation between the 11 connexins in CRC and important cancer associated signaling pathways. The solid line represents activation and the dashed line represents inhibition. (D) Differential expression of connexins between CRC and paired normal tissues. (E) Relationship between CNVs, DNA methylation and connexins expression. Red indicates positive correlation; blue indicates negative correlation. The deeper color indicates a larger correlation index. The bubble size indicates the FDR. (F) The correlation between CNVs, DNA methylation of the 21 connexins and CRC patients’ survival including DFI, DSS, OS, PFS. Red shows a positive correlation and blue shows a negative correlation. The darker color indicates a larger correlation index. Bubble size indicates the FDR. (G) Correlation between expression levels of connexins, red represents positive correlation, blue represents negative correlation, and shades of yellow indicates the P values. (H) Mutation frequency of 21 connexins in 620 patients with CRC in the TCGA-COADREAD cohort. Each column represents an individual CRC patient. (*P < 0.05, **P < 0.01, ***P < 0.001). Abbreviation: NS, not significant.

Identification and Analysis of Connexin Molecular Clusters of CRC

To identify potential connexin molecular clusters of CRC, we classified 620 CRC cases in the TCGA cohort using consensus clustering analysis. Based on the “ConsensusClusterPlus” R package, all tumor samples were classified into k (k = 2–6) different clusters based on the expression levels of 21 connexins coding genes in GC. According to the results of cluster analysis, we chose k = 2 being the number of clusters, indicating that CRC patients were accurately classified into two clusters (C1 and C2 clusters) (Figure 2A–C). The principal component analysis (PCA) showed that the classification was effective in distinguishing between two clusters of CRC patients (Figure 2D), and Kaplan–Meier survival analyses showed that C1 patients had worse OS, DFS, and PFS relative to C2 patients (Figure 2E–G). In addition, Figure 2H shows the connexins expression profiles based on both clusters and the correlation between the distribution of clinicopathological parameters including T/N/M and pathological stage with the two clusters.

Figure 2 The landscape of connexins in CRC and the characteristics of two clusters. (A) The cumulative distribution function (CDF) curves in consensus cluster analysis. CDF curves of consensus scores by different subtype numbers (k = 2, 3, 4, 5, and 6) were displayed. Relative change in area under the CDF curve for k = 2–6 was also shown. (B) The consensus score matrix of all samples when k = 2–6. (C) The PCA distribution of TCGA-COADREAD samples by expression profile of connexins. Each point represents a single sample; different colors represent the C1 and C2 clusters respectively. (D–F) Survival analysis of (D) OS, (E) PFS, and (F) DFS based on two clusters. (G) Expression distribution of 21 connexins between two clusters, connexins with 0 expression in most samples have been removed. (H) Distribution of clinical characteristics of two subtypes, including ‘tumor-node-metastasis” (TNM) stages and pathological stage. (I) The thermogram shows variations in mRNA expression of chemokines, interleukins, alterons, and other cytokines between two clusters. (J) Representative pictures of pathological HE staining of two clusters. (K–L) The differences of (K) TME infiltrating cells and (L) common immune checkpoints between the two connexins phenotypes (Kruskal–Wallis test). (*P < 0.05, **P < 0.01, ***P < 0.001).

We then analyzed the effect of connexins on the tumor immune microenvironment. We observed that chemokines, interleukins, and alterons had higher expression in cluster C1 than in cluster C2, including CXCL10, CCL7, VEGFB, VEGFC, and TGFBR1 (Figure 2I). In addition, we also found that m6A and ferroptosis-related genes were differentially distributed between cluster C1 and C2 (Supplementary Figure 2A and B). TCGA pathology slide confirmed that there was more infiltration of immune cells in the CRC nests of Cluster C1 patients but less infiltration of immune cells in the CRC tissue of Cluster C2 patients (Figure 2J). Then CIBERSORT algorithm based on TCGA-COADREAD demonstrated that patients in cluster C1 have more M2 macrophage infiltration than those in cluster C2 (Figure 2K) (P < 0.001). We finally compared the differences of immune checkpoints (ICPs) between the 2 clusters and most ICPs including CD274, CTLA4, HAVCR2, LAG3, PDCD1, PDCD1LG2, and TIGIT were significantly higher in C1 than in C2 (Figure 2L) (P < 0.001).

Identification of Functional Differences Between Two Clusters

We first identified the 20 genes with the highest mutation frequencies in CRC (Figure 3A) and then compared the differences in mutation frequencies of these 20 genes between the two clusters (Figure 3B). The differences in gene expression between subtypes of connexin and the characteristics of the associated signaling pathways were then explored. We observed that in CRC, especially GJC1, GJA5, GJA4 showed a 100% level of activation in the EMT signaling pathway, but consistent inhibition of the cell cycle and DNA damage response (Figure 3C). A total of 305 DEGs were screened based on the analysis by CRC gene expression profiles; 276 were up-regulated and 29 were down-regulated in cluster 1 relative to cluster 2 (Figure 3D). Then, all DEGs were used for the assembly of PPI, and based on the degree value from cytoHubba, 10 hub genes including FN1, DCN, POSTN, COL1A1, THBS2, MMP2, COL1A2, COL3A1, COL5A1 and BGN were identified (Figure 3E), all of which have strong activation effect on EMT in CRC (Figure 3F). GSEA showed that EMT, angiogenesis and macrophage activation were significantly activated in the cluster C1 (Figure 3G).

Figure 3 Analysis of connexins-related functions. (A) Mutation annotation of the top 20 frequently mutated genes in the TCGA-COADREAD cohort.(B) Comparison of the difference in mutation between C1 and C2 clusters. (C) The heatmap shows the correlation between the expression level of the 21 connexins in important cancer signaling pathways. The global percentage of cancers in which a gene has an effect on the pathway in COAD and READ, is shown as the percentage: (number of activated/inhibited cancer types/2 *100%). “Pathway activate” (red) represents the percentage of cancers in which a pathway may be activated by given genes, inhibition in a similar way shown as “pathway inhibit” (blue). (D) Volcano plot of Differentially expressed genes (DEGs) induced by alterations in connexins levels. red dots, with statistical significance; blue dots, not statistically significant; abscissa, expression differences (log2 fold change); ordinate, significance of differences (−log10 P value). (E) Network of DEGs. Darker colors and larger sizes indicate higher degrees. (F) The correlation between 10 hub genes from the network and important cancer signaling pathways. Red stands for COAD and blue stands for READ; the solid line represents activation and the dashed line represents inhibition. (G) GSEA of the status of special biological pathways in two connexins phenotypes. (H) LASSO model fitting. Each curve represents a gene. The profiles of coefficients were plotted versus log(λ). Vertical lines indicate the positions of seven genes with coefficients greater than 0 determined by 10-fold cross-validation. λ was determined from 10-fold cross-validation. The x-axis represents log(λ); the y-axis represents binomial deviance. Optimal values calculated from minimum criteria and one standard error of the criteria are indicated by the dotted vertical lines. (I) Univariate Forest plot and multivariate forest plot showing association between 4 candidate genes expression and OS, PFI, DSS in CRC. (J) Differential expression of GJA3, GJA4, and GJB6 between normal tissue and CRC tissue by using public datasets from Gene Expression Omnibus (GSE10950 and GSE44861). (K) Differential expression of GJA3, GJA4, and GJB6 between C1 and C2 clusters based on TCGA-COADREAD data. (L) Relationship between GJA3, GJA4, GJB6 levels and immune cell infiltration based on ssGSEA. The size of the bubble represents the strength of the correlation and the color of the bubble represents the P-value. (*P < 0.05, ****P < 0.0001). Abbreviation: NS, not significant.

Subsequently, 21 connexins were used for LASSO regression to identify robust markers for follow-up studies. Cross-validation was performed in ten rounds to prevent overfitting (Figure 3H). Sequential Cox regression (both univariate and multivariate) indicated that GJA3, GJA4, GJB6 were significantly correlated with CRC prognosis (Figure 3I). The calculation results based on two independent datasets (GSE44861 and GSE10950) showed that GJA4 was overexpressed in CRC tissues (Figure 3J), consistent with the results from TCGA (Figure 1D). Finally, these three genes acting as potential cancer-promoting factors were shown to be significantly highly expressed in cluster C1 (Figure 3K), and all of them, especially GJA4, were significantly associated with the degree of immune cell infiltration (Figure 3L).

Clinical Analysis of GJA4

The results of the HPA database suggested that GJA4 was significantly highly expressed at the protein level in CRC tissues relative to normal tissues (GJB6 was not included) (Figure 4A and B). Then, the role of GJA4 in CRC was generally assessed. In cases from MEXPRESS categorized based on different clinical factors, including colon polyps present, lymphatic invasion and sample type, GJA4 displayed the differences in expression patterns (Figure 4C and D). The results of independent calculations based on TCGA-COADREAD showed that the GJA4 expression level was correlated with pathologic stage and T stage (Figure 4E). Kaplan–Meier survival curves were then constructed to validate the prognostic value of GJA4 and the results showed that patients with higher GJA4 expression level had shorter PFI and disease-specific survival (DSS) (Figure 4F) (P < 0.05). Finally, an IHC H-score based on our own samples showed that GJA4 was expressed at significantly higher levels in CRC tissues than in paracancerous tissues (Figure 4G and H) (P < 0.05).

Figure 4 Clinical value of GJA4. (A) Representative immunohistochemistry images of GJA4 and GJA3 in CRC and noncancerous tissues derived from the HPA database. (B) Staining strengths were annotated as not detected, low, medium, and high. Bar plots indicate the number of samples with different staining strengths. (C and D) Correlations between GJA4 level and clinicopathological characteristics in (C) COAD and (D) READ from MEXPRESS. (E) Association of GJA4 mRNA expression with T/N/M stages and pathological stage in CRC patients. (*P < 0.05). (F) OS, DSS, and PFI about GJA4 based on TCGA-COADREAD. (G and H) Representative images of different immunohistochemical staining intensities for GJA4 based on our own samples and statistical comparison of GJA4 expression levels (H-SCORE) in (G) paracancerous and (F) CRC tissue (n = 30). (*P < 0.05, ***P < 0.001). Abbreviation: NS, not significant.

Single-Cell Level Analysis of GJA4

Herein, we attempted to localize GJA4 at the single cell level to investigate its potential function. We observed that in normal intestinal epithelium, GJA4 was predominantly expressed on enteroendocrine cells (Figure 5A). It is interesting to note that in CRC tissues, we observed that GJA4 was predominantly expressed on fibroblasts but not on CRC cells, like GJA1, an oncogene acting as an acid transmitter previously shown to be expressed on fibroblasts (Figure 5B). The results of the TIMER web tool also showed that the expression levels of GJA1 and GJA4 were highly positively correlated with the abundance of fibroblasts (Figure 5C and D). We then included an independent sample for validation and found that GJA4 was consistent with the expression of fibroblast markers such as FAP, PDPN, COL1A2, DCN, COL3A1, COL6A1, ACTA2 (Figure 5E–G). The co-expression pattern of GJA4 and ACTA2 which was not observed on the normal intestinal epithelium (Figure 5H) was confirmed by the IF staining on the CRC specimen (Figure 5I). These results suggested that in CRC tissues, GJA4 was not expressed in tumor cells, and fibroblasts overexpressed GJA4.

Figure 5 Single cell analysis and immunofluorescence of GJA4. (A) Distribution of GJA4 in normal intestinal epithelium at single level. (B) Single-cell analysis of GJA4 and GJA1 based on single cell RNA‐seq dataset GSE146771. Uniform Manifold Approximation and Projection (UMAP) plots showing different CRC cell types after quality control, reduction of dimensionality, and clustering; Enrichment scores of genes from the Hallmark EMT signaling gene set for each cell, from gene set variation analysis; UMAP plots showing expression of GJA4 and GJA1 clusters; The correlation between GJA1 and GJA4 was also calculated based on TCGA-COADREAD and GEPIA, respectively. (C and D) Dot plot showing the correlation between (C) GJA4, (D) GJA1 and fibroblasts based on the TIMER web tool. (E) Fibroblasts are distinguished from all other cells based on UMAP dimensionality reduction analysis. The t-distributed Stochastic Neighbour Embedding (t-SNE) is also used for dimensionality reduction of the data. (F) Heatmap demonstrating each cluster and its gene marker. (G) Dot plots demonstrate the expression distribution of GJA4 and seven fibroblast marker genes. (H and I) Double immunofluorescence staining images of GJA4 and ACTA2 in the normal and CRC tissue. The representative view of the co-staining of GJA4 and ACTA2 is shown in the enlarged images view below. Scale bars, 100 and 20 mm (enlarged images). Nuclei (DAPI) in blue.

Immune Microenvironment Analysis of GJA4

Considering fibroblasts as one of the powerful regulators of the immune microenvironment, we further analyzed the relationship between GJA4 and immune cells in this section. Through ESTIMATE analysis, we found that GJA4 was negatively correlated with tumor purity (P < 0.001, R = −0.609; P = 0.002, R = −0.457) and positively correlated with stromal score (P < 0.001, R = 0.536; P = 0.002, R = 0.600; Figure 6A) in the GSE44861 and GSE10950 datasets, respectively. Cellular communication from the scTIME Portal revealed complex interactions among fibroblasts, myeloid cells and malignant cells (Figure 6B). Schematic sketch of intercellular communication in the development of stromal cells is shown in Figure 6C. Fibroblasts have been identified to secrete cytokines to induce macrophages to polarize toward the M2 phenotype to exert inflammatory suppressive effects; therefore, we calculated the correlation between macrophage abundance and the GJA4 level. We found a positive correlation between the degree of infiltration of M2 macrophages and GJA4 levels using GSE10950 and GSE44861 (Figure 6D and E), which was then confirmed by IHC staining (Figure 6F). Finally, we searched for genes positively associated with GJA4 based on its median expression level (Figure 6G) and performed enrichment analysis (Figure 6H), which showed that GJA4 was closely associated with EMT-related signaling pathways such as “Focal adhesion”, “extracellular structure organization”, “extracellular matrix organization”, and “ECM–receptor interaction”.

Figure 6 The role of GJA4 in the immune microenvironment. (A) The correlation between GJA4 mRNA expression with stromal/immune score and tumor purity in the GSE44861, GSE10950 and TCGA datasets. (B) Cellular communication from scTIME Portal. (C) Schematic diagram demonstrating the driving effect of fibroblasts on macrophage polarization. (D) Proportion of 22 types of Tumor-initiating cells (TICs) in the CRC tumor samples based on GSE10950 and GSE44861. (E) Correlation of the GJA4 expression levels with macrophage abundance based on GSE10950 and GSE44861. (F) Correlation between GJA4 and the M2 macrophage marker CD163 based on immunohistochemical H-score calculation. (G) Volcano map showing the genes co-expressed with GJA4. red dots, positively correlated genes; grey dots, not statistically significant; abscissa, expression differences (log2 fold change); ordinate, significance of differences (−log10 P value). (H) Functional enrichment analysis of the GJA4-related genes including GO and KEGG.

Discussion

Currently, metastatic colorectal cancer remains the second leading cause of cancer death and places a heavy burden on the health-care system.⁴¹ As a critical step in the process of tumor metastasis, activation of the EMT process is necessary for metastatic growth.⁴²

Gap junction is one of the most widely distributed cell junctions and cell-to-cell communication channels in animal tissues, providing a way for the diffusion of low molecular weight substances between cells. In the past few decades, it has been found that it has an inseparable relationship with EMT.⁴³ There is no doubt that the loss of direct intercellular communication is one of the essential features of cancer cells, a fact that somehow implies the cancer-suppressive properties of gap junctions. Cx43 reduces the growth of rat gliomas in vitro and in vivo,⁴⁴ and Cx32 has also been observed to show a reduction in the growth of human liver cancer cells when cells are injected into mice.⁴⁵ As research has become more advanced and expanded, there is now a wealth of evidence that connexins are broadly involved in cell proliferation, apoptosis, drug resistance, migration and invasion, but not all studies point to a tumor suppressive effect.^46,47 Interestingly, many of the connexins for which substantial evidence exists as tumor suppressors have been progressively found to be overexpressed and may lead to more aggressive tumor characteristics in some cases. As early as 2000, it was discovered that Cx26 plays a role in tumor cell intra/extravasation by forming heterologous gap junctions with endothelial cells.⁴⁸

Considering such a complex and variable role of connexins, we systematically evaluated the differential expression of 21 connexins-encoding genes at the transcriptome level and found that most of them were overexpressed in tumor tissues compared to normal tissues, in addition, members including GJA4, GJA5, and GJA8 were found to be indicators of poor prognosis. We then classified TCGA-COADREAD cases into two clusters based on the all connexins expression levels by Consensus Clustering, and survival analysis showed that patients in cluster C1 had shorter OS, PFS and DFS than those in cluster C2. Subsequently, we found significant differences in the levels of immune-related factors between the two clusters, and not only that, we also observed that M2 macrophages had higher infiltration levels in cluster C1 relative to cluster C2, whereas macrophage M1, NK cell resting and B cell plasma showed the opposite trend, suggesting that cluster C1 showed a distinct immunosuppressive profile relative to cluster C2. Notably, the sensitivity of tumors to treatment targeting ICPs appears to be governed by the quality and quantity of tumor immune infiltration, and one mechanism by which M2 macrophages promote an immunosuppressive microenvironment appears to be through enhanced expression of ICPs.^49–51 We then calculated that ICPs levels were significantly higher in the cluster C1 than in the cluster C2. In addition, we obtained genes with significantly different expression levels between the two clusters and performed enrichment analysis, which showed that EMT as well as angiogenic phenotypes were significantly activated in the cluster C1, and that most of the hub genes obtained based on the MCC algorithm were also key players in EMT, results that are consistent with the now widely confirmed close association between connexins and EMT. The consistent association between EMT-related gene expression and immune cell infiltration suggests an inextricable link between EMT and influencing the response to ICP blockade.

The above results suggested that changes in the expression profile of connexins are closely associated with EMT activation/inhibition. To further identify reliable diagnostic and prognostic indicators, we combined LASSO and COX regression to identify GJA3, GJA4, and GJB6 as independent prognostic factors, and all three were significantly highly expressed in the cluster C1. To determine the sole study target, we further used transcriptomic data from GEO database and immunohistochemical data from HPA database as well as our own patients and found substantial evidence supporting the high expression of GJA4 in tumor tissue. Notably, GJA4 levels correlate with some clinical features of CRC patients, especially high expression in patients with advanced disease. With increasingly refined clinical observations in large samples, there is now a tendency to believe that individual connexins may exhibit dual effects, such as acting as tumor suppressors in primary tumor initiation, but having the opposite effect of promoting cancer progression in advanced disease.

It has been shown that the proliferation of insulinoma cells is significantly slowed upon induced expression of Cx37, a result that is associated with a specific conformation of Cx37.⁵² In addition, it has also been reported that the Cx37 C1019T polymorphism may contribute to tumor cell growth.⁵³ In conclusion, like all other connexin members, Cx37 may have a dual role in tumor development. Unlike Cx43, which has been extensively studied, Cx37 has been reported very sparsely in tumors and has hardly even been reported in CRC. We first observed Cx37 expression in normal intestinal epithelial tissue, mainly on endocrine cells, which are a type of intestinal epithelial cells that secrete granule-releasing peptides and amine-like hormones. However, for CRC tissue, Cx37 was found to be abundantly expressed on fibroblasts, a fraction of cells also enriched for EMT-related genes. Alzbeta Hulikova et al found that the intercellular transmission of acid through gap junctions composed of Cx43 involvement was high between fibroblasts but not between CRC cells.¹⁶ Notably, loss of Cx43 expression was associated with a significant reduction in relapse-free and overall survival of CRC patients and was able to negatively regulate the Wnt signaling pathway thereby enhancing CRC cell apoptosis.^54–56 We speculated that this may reveal an interesting fact that the significance of connexins expression in the parenchymal and stromal components of tumors may be completely opposite. Since gap junctions are mainly composed of 2 connexons, each connexon is composed of 6 subunits (connexins) to form functional unit. Our following calculations revealed a strong positive correlation between Cx37 and Cx43, suggesting a co-activation relationship between the two that may be jointly involved in constituting gap junctions. Results based on four algorithms of the TIMER web tool also indicated that both Cx37 and Cx43 expression levels were significantly positively correlated with fibroblasts abundance. The results of immunofluorescence directly showed that Cx37 was barely expressed on fibroblasts in normal intestinal epithelial tissues, whereas it appeared dramatically on fibroblasts in CRC tissues. While gap junction is a lethal factor for CRC cells, the fact that fibroblasts can help tumor cells metabolize potentially toxic amounts of acid through gap junctions to maintain a pH that is favorable for cancer cell proliferation and invasion results in CRC cells not only excreting excess acid but also minimizing apoptosis and energy burden. We suggest that this may be a microcosm of tumor-stromal interactions.

Since the expression of Cx37 was consistent with the abundance of fibroblasts, the most prominent component of the tumor stromal microenvironment, we continued to investigate the function of Cx37 in the tumor microenvironment. We calculated the relationship between GJA4 levels and stromal scores, immune scores based on three independent datasets and observed a significant positive correlation. It is well established that high immune/stromal scores tend to imply low tumor purity, which not only represents unfavorable survival and high metastatic recurrence rates but also correlates with the immune evasion phenotype. Tumor-associated macrophages usually acquire an immunosuppressive state that suppresses antitumor immunity, which can be induced by CAFs secreting cytokines, the best known of which is stromal-derived growth factor-1.^57–59 Further calculations showed that GJA4 was significantly positively correlated with CD163 both at the protein level and at the transcriptome level, suggesting a role of GJA4 in monocyte trans-differentiation toward the M2 phenotype. Complex malignant crosstalk among cancer cells, cafs and M2 macrophages can synergistically promote EMT of tumor cells to enhance stemness and ultimately contribute to tumor metastasis. Enrichment analysis based on GJA4 expression levels also implied that GJA4 was closely associated with EMT.

Conclusions

Overall, this study provides a new perspective to understand gap junctions. There is a heterogeneity of gap junctions in CRC, and different gap junctions subtypes showed different prognosis and immune infiltration patterns. How connexins function in tumor progression may be closely related to their localization in the complex tumor microenvironment. We preliminarily recognize Cx37 expressed on fibroblasts as a potential tumor promoter and is associated with a state of immunosuppression. The underlying complex molecular mechanisms and the prognostic and diagnostic value of Cx37 still need to be explored and validated in more experiments and clinical trials with large samples.