CILP acted as a late event in BRCA with poor prognosis
The Function module of LinkedOmics was used to analyze mRNA sequencing data from 1093 BRCA patients in TCGA. Generally, CILP mRNA levels were reduced with higher grade of clinical stage (Figure 1a). Thus, it was lower in the M1 stage than in the M0 stage (Figure 1b). As shown in Figure 1c, the mRNA expression of CILP different across the normal and Luminal, Her2 positive, and TN intrinsic subclasses of BRCA. The expression of CILP was significantly lower in HER2 positive and TNBC than luminal subclass. To further examine the prognostic potential of CILP in BRCAs, Kaplan-Meier plotter database was used to evaluate the CILP prognostic value. Interestingly, the poor prognosis in BRCA was shown to correlate with lower CILP expression (Figure 1d). However, CILP expression had no effects on survival in early BRCA (stage I to II), and showed a worse overall survival (OS) in advanced BRCA (stage III to IV) of RNA-seq data (Figure 1e).
From the above findings, CILP may serve as a late event occurring in BRCA with poor prognosis.
Identification of markedly different expressed genes between breast cancer tissue and brain metastases tissue from breast cancer
To identify DEGs between BRCA tissue and BM tissue from BRCA, data from 51 samples from two independent mRNA expression arrays of GSE100534 and GSE125989 datasets (25, 26) were downloaded from the GEO database and normalized and merged using the “oligo” and “sva” packages in RStudio (Figure 2a). According to the thresholds set (adjusted P-value<0.05 and log2FC ≥ 1), 90 DEGs (62 downregulated and 28 upregulated) were identified (Figure 2b). Subsequently, we visualized the 90 DEGs in the form of heatmap including all samples (Figure 2c).
CILP was the hub gene significantly associated with BCBM identified by construction of weighted co-expression network (WGCNA)
We used the merged data to construct the gene co-expression networks, which included a total of 51 samples complete with clinical features data for the WGCNA. Raw data were preprocessed identically using RStudio for background correction and normalization. Finally, we obtained a total 11,786 genes, and all were chosen for WGCNA. WGCNA is a systematic biological approach used to analyze the expression patterns of multiple genes in different samples, which forms clusters or modules containing genes with the same expression pattern (27). If certain genes are located in the same module, they are likely to possess the same biological functions and the association between modules and sample characteristics such as clinical traits can be studied (28, 29). The selection of the soft-thresholding power is an important step when constructing a WGCNA. We performed the network topology analysis for thresholding powers from 1 to 20 and identified the relatively balanced scale independence and mean connectivity of the WGCNA. In this study, we selected the power of β=8 (scale free R^2=0.83) as the soft-thresholding to achieve a scale-free network (Figures 3a, 3b). As a result, we set the cut height as 0.25 to merge similar modules (Figure 3c). Twenty gene-co-expression modules were identified using merged dynamic tree cut (Figure 3d). The yellow module showed the highest correlation with BM phenotype of invasive BRCA (Figure 3e). There were 399 genes in the yellow module, data for yellow module were selected for further analysis. Subsequently, we identified the module containing CILP, that is., the yellow module, as a key module accounting for the highest correlation with BRCA brain metastasis phenotype (R2=−0.6, P=3e-06) (Figure 3e). A total of 24 genes, including the CILP gene, were identified as candidate hub genes for the yellow module (Figure 3f). It is well-known that CILP is an inhibitor of transforming growth factor (TGF)-beta signaling. Interestingly, TGF-beta can enable cells to acquire a migratory and infiltrating phenotype, thus enabling the spread of cancer cells.
We further performed enrichment analysis to explore the biological processes and pathways in which the yellow module was involved. GO biological processes and the KEGG functional enrichment analyses were conducted using the R package “clusterProfiler” (21). GO terms such as extracellular matrix organization, collagen fibril organization, cell-substrate adhesion, epithelial cell proliferation, transmembrane receptor protein serine/threonine kinase signaling pathway, response to TGF-beta, collagen metabolic process, were mainly enriched in the “biological process” category (Figure 3g). The results of the KEGG enrichment analysis revealed that focal adhesion, complement and coagulation cascades, proteoglycans in cancer, ECM-receptor interaction, relaxin signaling pathway, TNF signaling pathway, apoptosis, TGF-beta signaling pathway, and PI3K-AKT signaling pathway, and many other related KEGG biological pathways were both important for the development of carcinoma (Figure 3h). Moreover, we constructed a network of protein-protein interaction (PPI) for all genes in the yellow module by Cytoscape, which consisted of 56 nodes and 366 edges according to the weight of the edge (≥0.07) (Figure 3i).
In summary, we may conclude that the CILP gene is significantly associated with the BM phenotype of BRCA.
CILP Expression was correlated with immune infiltration levels in breast cancer favoring BCBM
The composition of significant tumor-infiltrating immune cells was assessed using the CIBERSORT algorithm for the BRCA and BCBM tissues. For the merged GEO dataset, the heatmap (Figure 4a) and the histogram (Figure 4b) indicated that plasma cells, macrophages M2, T cells CD8, memory B cells, and follicular helper T cells accounted for a large proportion of the BCBM immune cell infiltration. While the heatmap (Figure 4c) showed M0 macrophages, macrophages M1, macrophages M2, CD8+ T cells, and resting mast cell accounted for a large proportion of BRCA immune cell infiltration. As shown in Figure 4d, the violin plot suggested that resting mast cells (P=0.010) showed obvious differences in the immune cell fractions between brain metastasis and primary BRCA via the Wilcoxon rank-sum test. Moreover, the Wilcoxon rank-sum test suggested that memory CD4+ T cells (P<0.001), follicular helper T cells (P<0.001), regulatory T cells (Tregs) (P<0.001), resting NK cells (P=0.001), monocytes (P<0.001), M0 macrophages (P<0.001), M1 macrophages (P<0.001), M2 macrophages (P<0.001), resting mast cells (P=0.003), and activated mast cells (P=0.001) were significantly different among the immune cell fractions between BRCA tissues and adjacent normal tissues (Figure 4e). The level of infiltrating resting mast cells was low in BCBM, and the relationship between resting mast cells and BCBM is worthy of further study. BCBM occurs during stage IV of neoplasms. We further analyzed the association between infiltrating immune cells and pathological staging. We concluded that M0 macrophages (Staging: P=0.007), M2 macrophages (Staging: P=0.018), activated memory CD4+ T cells (Staging: P=0.019), resting mast cells (Staging: P=0.001), monocytes (Staging: P<0.001), and plasma cells (Staging: P=0.007) all showed significant correlations with clinical features (Figures 4f). It was obvious that resting mast cells showed low expression in stage IV of BRCA. Subsequently, we performed Kaplan-Meier survival analysis, and found that only samples enriched in M2 macrophages had a poor prognosis (P=0.005) in BRCA (Figure 4g). Subsequently, we analyzed the relationship between CILP and resting mast cells, and both were involved in BCBM (Figures 3e-f and 4d). Further studies (Figures 4h) showed that the expression of CILP positively correlated with resting mast cells (R=0.33; P<0.001), resting memory CD4+ T cells (R=0.3; P<0.001), resting dendritic cells (R=0.16; P<0.001), gamma delta T cells (R=0.14; P<0.001), and naive B cells (R=0.13; P<0.001), instead negatively correlated with activated dendritic cells (R=-0.21; P<0.001), memory B cells (R=-0.11; P=0.0014), M0 macrophages (R=-0.14; P<0.001), and follicular helper T cells (R=-0.18; P<0.001), suggesting that CILP was associated with the immune microenvironment of BRCA. Of these, resting mast cells involved in brain metastasis had a significantly weak positive correlation with the expression of CILP. In the Timer database, we also found a weak positive association between CILP and mast cells (Figure 4i). Essentially, resting mast cells showed a lower expression of stage IV of BRCA, and the expression of CILP was also found lower in stage IV of BRCA, which is consistent with the positive association between CILP and resting mast cells.
Thus, it is reasonable to speculate that CILP may influence the progression of BRCA favoring BMs through the infiltration of resting mast cells.
GSEA analysis identified the potential pathways involved in BCBM
To identify the potential function of CILP genes in BRCA samples from TCGA database, GSEA was used to identify the potential pathways involved in brain metastasis processes. Nine enriched pathways were detected, including the B cell receptor signaling pathway; Hedgehog signaling pathway; lysosome, ubiquitin mediated proteolysis; DNA double-strand break response; antigen activation of the B cell receptor (BCR) leading to generation of second messengers; CD22-mediated BCR regulation; creation of C4 and C2 activators; and FCERI-mediated MAPK activation; all of which were closely related to tumor occurrence, development, and metastasis (P-value<0.05, FDR<0.25) (Figure 5).