Meta-Analysis Reveals the Prognostic Relevance of Nuclear and Membrane-Associated Bile Acid Receptors in Gastric Cancer.

INTRODUCTION: Bile acids (BAs) arising from duodenogastric reflux are known to facilitate gastric cancer (GC) development. Although BAs traditionally contribute to carcinogenesis through direct cellular cytotoxicity, increasing evidence implicates nuclear and membrane BA receptors (BARs) as additional factors influencing cancer risk. Indeed, some BARs are already linked with GC, but conflicting evidence and lack of information regarding other endogenous BARs warrant further investigation. In this study, we meta-analyzed multiple data sets to identify clinically relevant relationships between BAR expression and prognosis, clinicopathology, and activity in GC.
METHODS: We collected transcriptomic data from the Gene Expression Omnibus and The Cancer Genome Atlas to analyze associations between BAR expression and GC prognosis, subtype, and clinicopathology. We also used Ingenuity Pathway Analysis to assess and predict functions, upstream regulators, and downstream mediators of membrane and nuclear BARs in GC.
RESULTS: BARs showed differential distribution in GC; membrane BARs (G protein-coupled BAR 1, sphingosine-1-phosphate receptor 2, and cholinergic receptor muscarinic 2) were enriched in diffuse-, genome-stable, and mesenchymal-type tumors, whereas nuclear BARs (pregnane-X-receptor, constitutive androstane receptor, and farnesoid-X-receptor) were enriched in chromosome instability and metabolic subtypes. High expression of all membrane but not nuclear BARs was associated with poor prognosis and unfavorable GC clinicopathologic features. Similarly, expression patterns of membrane but not nuclear BARs varied geographically, aligning with Helicobacter pylori infection and GC mortality rates. Finally, GC-related oncogenes, namely transforming growth factor β1, were associated with membrane BARs, whereas many metabolic-associated genes were associated with nuclear BARs. DISCUSSION: Through transcriptomic meta-analysis, we identified distinct expression profiles between nuclear and membrane BARs that demonstrate prognostic relevance and warrant further investigation.

INTRODUCTION

Gastric cancer (GC) is the fifth most common and third deadliest cancer worldwide, with a median overall survival of less than 1 year (1). In East Asia, for example, GC is the second leading cause of cancer-related deaths, the highest of any geographical region (2). Although overall GC incidence rates have been declining, incidence of non-cardia GC is rising in patients younger than 50 years (3). Helicobacter pylori infection represents the primary risk factor of developing GC in up to 78% of all cases (1), and eradication by a combination of antibiotics and proton pump inhibitors remains the mainstay of treatment (4). However, evidence has shown elimination of infection is only beneficial if no precancerous lesions are present (5) and, even then, only contributes to 33%–47% reduction in GC risk (6). Unfortunately, because up to two-thirds of GC patients initially present with late-stage GC (1), the long-term clinical utility of H. pylori eradication is uncertain. For example, studies show that atrophic gastritis with progression to adenocarcinoma might continue to occur even after H. pylori eradication (5) but the cause remains elusive and is suspected to be because of environmental and/or hereditary factors.

One likely determinant is duodenogastric reflux (DGR), which occurs when intestinal biliary contents such as bile acids (BAs) and bilirubin reflux through the pylorus and into the stomach (7). DGR is a known risk factor of GC (8–10) and is linked to active H. pylori infection (7,11) and gastric surgery but can also occur idiopathically. The BA component is believed to be the primary oncogenic source and contributes to neoplasia through direct gastric epithelial cell cytotoxicity (12–14). However, BAs can also signal directly through a subset of BA-sensitive nuclear (pregnane X receptor [PXR], constitutive androstane receptor [CAR], and farnesoid X receptor [FXR]) and membrane (G protein-coupled BAR [GPBAR] 1, sphingosine-1-phosphate receptor [S1PR] 2, and cholinergic receptor muscarinic [CHRM] 2) receptors (BARs), which are directly involved in the development of GC and other gastrointestinal (GI) cancers, as shown by previous studies (15–18). Nuclear BARs are primarily stimulated by hydrophobic BAs, such as CDCA, which are typically linked to chemical gastritis, followed by intestinal metaplasia (IM) and adenocarcinoma transformation. For example, activation of FXR promotes caudal type homeobox (CDX) 2-mediated IM (19,20), whereas PXR activity promotes GC cell multidrug resistance (21). On the other hand, membrane BARs have greater affinity for polar BAs and bile salts and are associated with a wider array of functions regarding GC tumorigenesis. The most well-known membrane BAR is GPBAR1, which has both tumor-suppressor and pro-oncogenic roles in GC (22–24). However, conflicting evidence among established BARs (FXR, PXR, and GPBAR1), combined with the fact that no information exists pertaining to the role of other BARs (CAR, S1PR2, and CHRM2) in GC despite normally being present in gastric tissue (25–27), and involved in other GI tumors (17,18) warrant both clarification and further investigation, respectively.

To address this gap in knowledge, we used a meta-analytic and bioinformatics-based approach to identify transcriptomic profiles of both nuclear (PXR, CAR, and FXR) and membrane (GPBAR1, S1PR2, and CHRM2) BARs in GC across a diverse set of publicly available data sets sourced from the Gene Expression Omnibus (GEO) and The Cancer Genome Atlas (TCGA). In addition, we assessed the relationships of each BAR gene with prognostic indicators such as survival, clinicopathology, and histologic and recently described molecular GC subtypes. Finally, using the Ingenuity Pathway Analysis (IPA) tool set, we probed potential upstream regulators, disease associations, and mechanistic networks of closely related genes common to nuclear and membrane BARs to elucidate potential pathways involved in mediating BAR-, and thus BA-, related changes in GC. Ultimately, building on results from previous studies with additional global BAR analyses will provide a greater understanding of the clinical role these receptors play in GC while providing potential avenues for therapeutic targeting for improving management and patient outcome.

METHODS

Data selection, acquisition, and processing

Data sets from RNA sequencing and microarray studies were sourced from GEO (https://www.ncbi.nlm.nih.gov/geo/) (28) and TCGA, respectively. Preprocessed RNA sequencing data were downloaded from the University of California Santa Cruz Xena browser (TCGA-STAD, https://xena.ucsc.edu/) (29,30) and cBioportal (TCGA-Molecular 2014, https://www.cbioportal.org/) (31,32). Data from both data sets were normalized as log-transformed transcripts per million (log2(TPM+1)). By contrast, GEO data sets were identified using the Search Tag and Analysis Resource for GEO by search terms and inclusion/exclusion criteria as outlined in Supplementary Figure S1 (see Supplementary Digital Content 1, http://links.lww.com/CTG/A486) (33). After additional manual curation, 9 studies (34–42) were found to be suitable for comparison between normal and tumor tissue, whereas 7 studies (34–36,43–45) were suitable for comparison between intestinal- and diffuse-type cancers. All GEO studies used in this investigation are detailed in Supplementary Table S1 (see Supplementary Digital Content 2, http://links.lww.com/CTG/A487). Raw data were subsequently downloaded from GEO, background corrected, normalized by probe, and log-transformed in R using robust multichip averaging (RMA) using the oligo package (Affymetrix-based arrays) or the lumi package (Illumina-based arrays) (46). Afterward, RMA-normalized data sets were collapsed according to max probe value, annotated by gene symbol in R, and prepared for downstream meta-analysis or data set merging with cross-platform batch removal.

Meta-analysis of GEO studies

Expression of nuclear and membrane BARs between normal and tumor tissues as well as histologic GC subtypes across multiple GEO studies were meta-analyzed. In the meta-analysis, the standardized mean difference (SMD) of gene expression in each study was first estimated using Hedge's g that is defined as the mean difference divided by the pooled SD (47). The summary SMD was then estimated by weighting each SMD by the inverse of the variance determined by random effects modeling (48). Results were represented as a Forrest plot. Summary statistics including Cochrane Q, heterogeneity (I2), and tau2 (T2), and meta-analysis Z score, were subsequently determined. Afterward, distribution patterns of each sample set were determined through plotting residuals on a normal quantile plot. Finally, publication bias (assuming a normal distribution pattern) was assessed using Egger Regression (49) and the test proposed by Begg and Mazumdar (50) and visually represented as a standard funnel plot in Supplementary Figures S2 and S3 (see Supplementary Digital Content 1, http://links.lww.com/CTG/A486).

Data set merging, batch correction, and cross-validation

Data set merging followed by interstudy and cross-platform batch correction represents a secondary means to assess and compare gene expression patterns across multiple studies (51). After RMA-normalized data sets were collapsed by probe and annotated, they were merged into a single data set according to matching gene symbols (see Supplementary Figure S1, Supplementary Digital Content 1, http://links.lww.com/CTG/A486). Interstudy batches were then identified by principle component analysis (PCA) and used to adjust the raw data through ComBat's Empirical Bayes approach provided in the sva R package (52). Specifically, the original parametric iteration of ComBat was used due to the significant batch heterogeneity to which the covariates of interest (normal, tumor, intestinal, and diffuse samples) were specified to preserve biological signatures postnormalization (53). Cross-platform normalization (CPN) was confirmed by PCA analysis and expression of nuclear and membrane BARs were determined, compared with the meta-analysis output, and validated with the TCGA-STAD data set. Finally, comparison of Spearman correlation coefficients of all BARs between the merged and TCGA-STAD data sets was used as a further validation metric and represented as a correlogram.

Analysis of BAR associations with survival, clinicopathology, and geographical region

Survival.

The impact of nuclear and membrane BAR expression on survival was analyzed by univariate cox-proportional hazard regression and plotted as either a hazard ratio (HR) Forrest plot or multicurve Kaplan-Meier plot for overall and histologic subtype-based differences, respectively. For assessing the impact of BAR expression on GC prognosis, disease-free survival (DFS) and overall survival (OS) probabilities were computed based on information obtained in GSE66229 and TCGA-STAD data sets. GSE66229 was chosen to compare against the validated TCGA-STAD study because it was the only study to provide information regarding survival status and time across the list of 9 studies. Because of the differences between both platform and baseline expression values, HRs and their corresponding 95% confidence intervals were computed based on P value optimization derived from comparisons between multiple expression cutoff values.

For analyzing the prognostic impact of BAR expression between histologic GC subtypes, OS probabilities were assessed over a 120-month period from GSE66229, GSE22377, and GSE15459 within the CPN-normalized merged data set. Expression cutoff values were determined from the median of the set due to the homogenous distribution of expression values between the studies. Thereafter, survival differences between low and high BAR expression cohorts were compared with each histological subtype cohort to identify patterns.

Clinicopathology.

Associations between BARs and clinicopathologic features of GC were assessed from the TCGA-STAD data set that contains the most robust list of clinical and pathologic descriptors. Expression of each BAR was determined for each factor, and results were computed as the log fold change (logFC) of expression between unfavorable and favorable characteristics. Afterward, the distribution of logFC values was visually represented as a heat map to identify potentially relevant patterns. Row and columns were clustered by pattern correlation. In addition, a heat map of Spearman correlation coefficients was similarly constructed for the comparison of BAR expression with continuous pathologic descriptors. Finally, a heat map of logFC values was constructed for comparison between BAR expression and molecular signatures across CpG island methylator phenotype subtypes described in the TCGA-Molecular 2014 data set (31). All P values for each independent comparison are described in Supplementary Tables S2 and S3 (see Supplementary Digital Content 2, http://links.lww.com/CTG/A487).

Geographical region.

For comparison of BAR distribution across different geographical regions, expression was compared and stratified according to region based on the descriptors provided by the TCGA-STAD data set (29). The impact of histologic subtype on the observed geographic distribution of BARs was determined using multifactorial analysis. In brief, the percent contribution of region to the observed variance was calculated from the 2 dimensions that captured the greatest variance in the data set and directly compared with histologic subtype. All P values for multiple comparisons between regions can be found in Supplementary Table S4 (see Supplementary Digital Content 2, http://links.lww.com/CTG/A487), and breakdown of ethnicity by region can be found in Supplementary Table S5 (see Supplementary Digital Content 2, http://links.lww.com/CTG/A487).

Ingenuity pathway analysis

To identify and predict shared functions and activities of BAR types in GC, a list of the top correlated genes or those whose expression patterns most similarly mirrored each BAR gene were analyzed using IPA (QIAGEN, https://www.qiagenbioinformatics.com/products/ingenuitypathway-analysis) (54). Expression similarity was first determined by Spearman correlation between BARs and all other genes within the TCGA-STAD data set. Significance was determined from BH-adjusted P values. Then, the top 7,500 and 2,500 significant results for nuclear and membrane BARs, respectively, were selected to ensure that the number of commonly overlapping genes fell within the recommended 1,000–2,000 range for downstream IPA analysis. Canonical pathways, disease and functions, upstream regulators, and mechanistic networks were subsequently reported as the predicted activity Z score.

Statistical Analysis

For meta-analysis, the SMD between conditions were compared and reported with 95% confidence interval values. Overall significance was computed from the meta-analysis Z score and reported as the 2-tailed P value. For comparisons of 2 groups, Welch unpaired t test assuming unequal variances was used. For comparisons of more than 2 groups, a 1-way ANOVA, followed by multiple t tests with P values adjusted by post hoc Tukey test. Univariate survival analyses were assessed using cox-proportional hazard regression with all reported P values being adjusted for multiple comparisons using the Benjamini-Hochberg (BH) method. A χ2 test was used for testing sampling bias and Fisher exact test was used for comparing sample distributions. Spearman correlation coefficients with corresponding significance P values were calculated for all correlation analyses. All statistical tests were performed in R, GraphPad Prism, or SPSS (IBM). Statistical significance was considered for *P < 0.05, **P < 0.01, and ***P < 0.001.

RESULTS

Comparative analysis reveals BAR expression is associated with GC clinicopathology and prognosis

To establish a framework for understanding the potential role BARs play in GC, we conducted a meta-analysis comparing expression patterns between normal (N = 335) and GC tissue (N = 590) across 9 independent studies sourced from GEO; for more details on each GEO study, see Supplementary Table S1, Supplementary Digital Content 1, http://links.lww.com/CTG/A487). Results of the meta-analysis reveal that the expression of every BAR, except PXR, was significantly reduced in tumor tissue (Figure 1a). Importantly, study results were found to be normally distributed and showed no significant publication bias (see Supplementary Figure S2, Supplementary Digital Content 1, http://links.lww.com/CTG/A486). To validate these results, we meta-analyzed the studies through data set merging, where data from each study are merged based on matching gene symbol followed by batch identification and CPN using ComBat. PCA of merged data sets pre-CPN and post-CPN revealed the identification of several distinct clusters (corresponding to individual studies) and their effective removal (Figure 1b). Analysis of the post-CPN merged data set showed expression patterns between nuclear and membrane BARs mirroring the original meta-analysis (Figure 1c, top) and, additionally, revealed intergene correlation patterns suggestive of a potentially common regulatory mechanism within but not between BAR types (Figure 1c, bottom). Furthermore, cross-comparison of the meta-analysis results with those of the TCGA-STAD RNAseq data set revealed similar findings, including reduced BAR expression, except PXR, in GC and common intergene correlation patterns, suggesting overall that expression of BARs are reduced during GC tumorigenesis (Figure 1d).

Figure 1.

Bile acid receptor (BAR) expression and prognostication varies in gastric cancer (GC). (a) Meta-analysis of BAR expression patterns between normal and gastric tumor tissue. Raw data from 9 studies sourced from the Gene Expression Omnibus (GEO) were background corrected, normalized, and log transformed using robust multichip averaging (RMA) and meta-analyzed using random effects modeling. Results are represented as a Forrest plot of the standardized mean difference (SMD), calculated as the difference of sample size–corrected average expression values between normal and tumor samples, and weighted by inverse variance. Statistical significance was based on the 2-tailed P value derived from the meta-analysis Z score. (b) Principle component analysis (PCA) plots showing the presence of distinct batches in the merged data set related to study (left) and their effective removal by cross-platform normalization (CPN) using ComBat (right). (c and d) Box plots showing expression patterns of nuclear and membrane BARs between normal and tumor samples in the post-CPN merged data set (c) and the TCGA-STAD validation data set (d, NNormal = 35, NTumor = 415) with corresponding gene-wide correlation patterns shown in the correlograms (below). Significance was calculated using Welch t test for the box plot and Pearson correlation coefficient for correlation; only significant associations were shaded. (e and f) Heat maps of expression (e) and Spearman correlation coefficients (f) across different clinicopathologic features of GC based on the TCGA-STAD data set. Row clustering was based on row-based correlation patterns. (g and h) Hazard ratio (HR) plots demonstrating the impacts of high nuclear and membrane BAR expression on overall survival (OS, g) and disease-free survival (DFS, h) compared between GSE66229 and TCGA-STAD. Expression cutoff values were calculated based on log-rank P value optimization through multiple comparisons and adjusted by the Benjamini-Hochberg method. *P < 0.05, **P < 0.01, ***P < 0.001.

Because findings from both the meta-analysis and TCGA contradicted our original belief that BARs were involved in promoting GC development, we explored the idea that the observed expression patterns could be explained by the heterogenous nature of GC. To address this issue, we first analyzed the impacts of high BAR expression on both GC OS and DFS from data provided by studies containing information on such parameters including TCGA-STAD and GSE66229. When comparing survival outcomes between cancer patients with high and low BAR expression, a disconnect between HRs and expression patterns were noted because many BARs, especially membrane types, were associated with reduced OS (Figure 1e) and DFS (Figure 1f) in both data sets. This revelation then prompted us to probe the relationship between BAR expression and clinicopathology as a factor potentially explaining this divergence. Results from each comparison were represented as a heat map, and when both rows and columns were clustered according to correlation patterns, it was determined that nuclear and membrane BARs associated more closely with each other and demonstrated reduced expression (logFC < 0) for nuclear BARs and enhanced expression (logFC > 0) for membrane BARs with unfavorable clinicopathology such as advanced tumor stage, grade, pathologic stage, treatment failure, and metastasis. In addition, we found that membrane BARs were correlated directly with tumor size and leukocyte infiltration and inversely with cancer age of onset, unlike nuclear BARs, which showed little association with any of these clinical parameters (Figure 1h; see Supplementary Table S2, Supplementary Digital Content 1, http://links.lww.com/CTG/A487). Taken together, our results show that, despite apparent downregulation in tumor tissue, high expression of BARs, and membrane BARs, in particular, is mainly associated with poor prognosis and unfavorable clinicopathologic features of GC.

Comparative analysis of BAR expression and prognosis in histologic GC subtypes

To explain the divergence between BAR expression patterns with prognosis and clinicopathology, we asked whether the observed differences could be potentially explained by the heterogenous nature of GC. For example, GC tumors differ substantially in phenotypic and molecular composition, the most well-known stratification system classifying tumors based on histologic characteristics and include intestinal- and diffuse-type cancers. Disease course and prognosis differs between these histologic subtypes, and results from our analysis of clinicopathology such as age of onset, treatment failure, grade, and pathologic stage indeed support this notion. Thus, we next examined expression patterns of both nuclear and membrane BARs between intestinal- and diffuse-type GC through similar meta-analyses. Data from 7 GEO studies containing information about histologic subtype were collected and meta-analyzed as before (Figure 2a; see Supplementary Table S1, Supplementary Digital Content 1, http://links.lww.com/CTG/A487). We found that expression of nuclear BARs favored intestinal-type cancers, whereas expression of all membrane BARs was significantly associated with diffuse-type cancers (Figure 2a). These expression patterns were further validated after data set merging and CPN (Figure 2b), which showed similar results, especially for membrane BARs (Figure 2c), with the initial meta-analysis and the TCGA-STAD data set (Figure 2d), despite a slight difference in the proportion of sampled histologic subtypes between the merged and TCGA-STAD data set (Figure 2e).

Figure 2.

Bile acid receptor (BAR) expression and prognostication varies across gastric cancer (GC) histologic subtypes. (a) Meta-analysis of BAR expression patterns between histologic subtypes (intestinal-type and diffuse-type cancers). Raw data were collected, processed, and meta-analyzed in the same way as mentioned previously from 7 studies that included descriptions about cancer subtypes. Significance was based on the 2-tailed P value calculated from the summary meta-analysis Z scores. (b) PCA analysis performed in the same way as previously mentioned, showing the presence of distinct interstudy batch effects that were effectively identified and removed by CPN using the ComBat function. (c and d) Box plots showing expression patterns of nuclear and membrane BARs between intestinal-type and diffuse-type cancers in both the post-CPN merged data set (c) and the TCGA-STAD validation data set (d, N = 235). Statistical significance was calculated using Welch t test. (e) Stacked bar plot showing the comparable distribution of histologic subtypes between the merged and TCGA-STAD data sets. Distribution similarity was assed using Fisher exact t test. (f and g) Kaplan-Meier plots comparing overall survival between histologic subtypes for nuclear BARs (f) and membrane BARs (g). Survival data were sourced from studies within the post-CPN, which contained information regarding survival time and status; this included GSE66229, GSE22377, and GSE15459. Expression cutoff values were determined individually by population median. Significance for each comparison was calculated using Cox regression and represented within the color chart inset. *P < 0.05, **P < 0.01, ***P < 0.001.

Regarding the prognosis between histologic GC subtypes cancers, we next assessed the prognostic value of BARs between these tumor types using data combined from 3 of the 7 GEO studies (GSE66229, GSE22377, and GSE15459). Kaplan-Meier curves of PXR and FXR showed little association between expression and OS, whereas high CAR expression was associated with reduced survival in both histologic subtypes (Figure 2f). By contrast, high expression of membrane BAR genes, especially GPBAR1, was consistently associated with reduced survival in diffuse-type cancers when compared with that of intestinal-type cancers (Figure 2g).

Overall, these results suggest that the expression of membrane but not those of nuclear BARs are differentially distributed among histologic subtype. Specifically, this data also point to membrane BARs being implicated in diffuse-type GC tumorigenesis, which is further supported by positive correlations of each receptor with genetic signatures associated with diffuse-type cancers, as reported by Zhang et al. 2016 (55) (see Supplementary Figure S4, Supplementary Digital Content 1, http://links.lww.com/CTG/A486). Importantly, these results might provide a potential explanation underlying the difference between expression and prognosis because diffuse-type cancers are known to be more aggressive and associated with prognostically unfavorable tumor characteristics (56).

Relating BAR expression with molecular alterations and genetic signatures in GC

Genetic and epigenetic modifications have been implicated as primary drivers of GC tumorigenesis (57,58). Because certain molecular alterations are strong prognostic factors and cluster with different GC subtypes (31), we assessed the distribution of nuclear and membrane BAR expression patterns across molecular subtypes using data from the TCGA-Nature 2014 and GSE15459 data sets.

Distinct clusters are apparent within microsatellite-instable (MSI) and genome-stable (GS) signatures (Figure 3a). For example, all membrane BARs showed inverse associations with MSI subtype molecular signatures, including ARID1A mutation, MLH1 silencing, and hypermutation. By contrast, expression patterns across all membrane BARs were positively association with GS subtype molecular signatures, including microsatellite stability (MSS), CLDN18–ARHGAP fusion, and RHOA mutations. These patterns were either less apparent or not present for nuclear BARs. Expression of membrane BARs in chromosomal-instable (CIN) and Epstein-Barr virus subtypes revealed general reduction, whereas nuclear BAR expression in CIN subtypes were elevated; however, these expression patterns were considered intermediate. Taken together, these findings are in agreement with the fact that GS tumors comprised mostly of the diffuse subtype, whereas MSI tumors comprised mostly of the intestinal subtype (Figure 3b), in line with our previous results (31) (for individual P values, see Supplementary Table S3, Supplementary Digital Content 2, http://links.lww.com/CTG/A487).

Figure 3.

Bile acid receptor (BAR) expression is associated with distinct gastric cancer (GC) molecular signatures and subtypes. (a) Heat map of expression across the 4 TCGA-Molecular GC subtypes: microsatellite instability (MSI), Epstein-Barr virus (EBV), chromosomal instability (CIN), and genomic stable (GS) along with their associated signatures. (b) Distribution of histologic subtypes across each TCGA-Molecular cancer subtypes. (c) Box plots showing BAR expression stratified by Lei-Molecular subtypes: proliferative (N = 70), metabolic (N = 40), and mesenchymal (N = 51), from GSE15459 for both nuclear BARs (left) and membrane BARs (right). Significance was computed by 1-way ANOVA with multiple t tests corrected for multiple hypothesis testing via post hoc Tukey's test. (d) Distribution of histologic subtypes across each Lei-Molecular cancer subtypes. *P < 0.05, **P < 0.01, ***P < 0.001.

We next investigated BAR expression across GC tumors classified by the recently developed Lei-Molecular system from GSE15459. Unlike the TCGA-Molecular system, which categorizes tumors based solely on genetic and epigenetic alterations (31), the Lei-Molecular system has greater therapeutic utility and classifies tumors into 3 major subtypes based on molecular signatures as a function of pharmacologic sensitivity (59). We found that expression of all nuclear BARs was greatest in metabolic tumors and lowest in mesenchymal tumors, especially for PXR (Figure 3c). In contrast to nuclear BARs, expression of all membrane BAs was lowest in the proliferative subtype, whereas GPBAR1 and S1PR2 expression was highest in mesenchymal tumors and S1PR2 in metabolic tumors (Figure 3c). Importantly, these expression patterns also closely mirror the distribution of histologic subtypes; the proliferative subtype is enriched in intestinal-type cancers, whereas the mesenchymal subtype is enriched with diffuse-type cancers (Figure 3d).

Geographical distribution of BAR expression in GC

Because regional GC prevalence, incidence, and mortality rates are closely tied to H. pylori infection (1), combined with our results indicating BAR expression relates to H. pylori infection (Figure 1g), we next assessed whether BAR expression was geographically distributed. When categorized by region, expression of nuclear BARs showed no trends (Figure 4a), whereas expression of membrane BARs progressively increased from western to eastern hemisphere regions (Figure 4b). Specifically, expression of membrane BARs corresponded with GC mortality rate, being highest in the Eastern European and Asian cohorts (highest mortality rates) and lowest in the North American and European cohorts (lowest GC mortality rates) (1) (for comparison between regions P values, see Supplementary Table S4, Supplementary Digital Content 2, http://links.lww.com/CTG/A487). Because ethnicity is also an important risk factor of GC development and survival outcomes, we assessed whether regional BAR expression patterns could be recapitulated by ethnicity (see Supplementary Table S5, Supplementary Digital Content 2, http://links.lww.com/CTG/A487). We found that there were no significant differences in expression of nuclear BARs between ethnic groups but that expression of membrane BARs was lower in the Black race (see Supplementary Figure S5, Supplementary Digital Content 1, http://links.lww.com/CTG/A486). However, little variation was noted among other groups, suggesting that the observed trends in expression are more a consequence of region, not ethnicity.

Figure 4.

Bile acid receptor (BAR) expression varies by gastric cancer sample source region. (a and b) Box plots showing the expression of (a) nuclear and (b) membrane BARs in South America (N = 62), North America (N = 55), Europe (N = 51), Eastern Europe (N = 135), and Asia (N = 21). (c) Multifactorial analysis (MFA) of the top 2 dimensions identifying distribution patterns between region (left) and histologic subtypes (right), including samples labeled as not otherwise specified (NOS). (d) Bar charts showing the relative contributions to observed variance in the expression of BARs across regions. (e) Stacked bar chart showing the distribution of subtypes according to region. (f) Modified pie charts showing the distribution of different countries within the included regions. All data were sourced from the TCGA-STAD data set because it contains the most robust list of samples by region.

To confirm that expression differences were not biased, we assessed the proportions of histologic subtypes for each region. Through multifactorial analysis, we found that the top 2 dimensions were able to explain 23.4% and 18.5% (41.9% total), respectively, of the observed variance with distinct clustering among both region and histologic subtypes, revealing unique expression patterns and magnitudes among these groups (Figure 4c). Moreover, we found that region contributed more to the observed variance in both dimensions (>40%) compared with that by histologic subtype (30%–40%) and the membrane BARs themselves (15%–25%) (Figure 4d). When further stratifying histologic subtypes based on region, bias was detected (χ2 P = 0.03) in the South America cohort when only comparing intestinal- and diffuse-type cancers but not samples labeled “not otherwise specified”. An oversampling of intestinal-type cancers was found in this cohort, likely contributing to the low expression of membrane BARs (Figure 4e). However, no sampling bias (χ2 test, P = 0.87) was found between the remaining regions, indicating an even distribution of histologic subtypes among these regions.

It is worthwhile to note that some regions were overrepresented by certain countries. For example, all samples in the South America group were from Brazil, whereas Vietnam contributed approximately 95% of the samples in the Asia group (Figure 4f). Although an analysis of a larger cohort of samples originating from a more diverse pool is required to fully assess the validity of these results, these initial findings suggest that expression of membrane but not nuclear BARs might be associated with region-based differences in GC.

Analysis of BAR-associated upstream regulators and functional pathways in GC by IPA

Previous findings have found conflicting results regarding the role nuclear and membrane BARs play in GC development. For example, a few studies have indicated a pro-oncogenic role for GPBAR1 (15) while others show a tumor suppressor–like role in GC cells and tumors (60). In addition, little is known about what roles other BARs, such as CAR, S1PR2, and CHRM2 play in GC, despite endogenous expression. Because our data suggests common roles for nuclear and membrane BARs based on conserved expression and distribution patterns, we ascertain which functions, regulators, and mechanisms are shared between each BAR type that could provide further insights into how BAs promote GC development at the cellular level. We accomplished this by using the TCGA-STAD data set to elucidate a set of genes that significantly correlated with each BAR, followed by determining the degree of overlap and submitting the list of common differentially correlated genes (DCGs) to IPA. For a complete list of genes and IPA results, see Supplementary Digital Content 3 (http://links.lww.com/CTG/A488) for nuclear BAR DCGs and Supplementary Digital Content 4 (http://links.lww.com/CTG/A489) for membrane BAR DCGs.

Analysis of DCGs revealed that there was a lower degree of overlap between DCGs for nuclear BARs (Figure 5a) than that for membrane BARs (Figure 5b). In fact, to reach the recommended 1,000–2,000 gene threshold for IPA analysis, 7,500 DCGs were needed for nuclear BARs, whereas only 2,500 DCGs were needed for the membrane BARs, supporting previous findings that gene similarity based on expression patterns was greater between the latter than the former. IPA analysis of each common DCG list for predicted upstream regulators revealed that genes involved in lipid metabolism such as NR5A1 and PPARA were predicted to be active, whereas hallmark tumor-suppressor genes including BRCA1 and TP53 were predicted to be downregulated (Figure 5c). By contrast, upstream regulators predicted to be active for membrane BARs were primarily oncogenes including transforming growth factor β 1 (TGFβ1) (61), ETS-related gene (ERG) (62), transglutaminase 2 (TGM2) (63), and TEA domain family member 1-4 (TEAD1-4) (64), whereas regulators predicted to be inhibited were other oncogenes such as SAM pointed domain containing ETS transcription factor (SPDEF) (65), zinc finger protein 217 (ZNF217) (66), and non-POU domain containing octamer binding (NONO) (62) (Figure 5d). Of note, multiple genes promoting epithelial-mesenchymal transition (EMT) were predicted active with membrane BARs such as hypoxia-inducible factor 1-alpha (HIF1A), twist family BHLH transcription factor 1 (TWIST1) and vascular endothelial growth factor A (VEGFA), whereas those opposing EMT, such as CDH1 (67), were predicted to be inhibited (see Supplementary Digital Content 4, http://links.lww.com/CTG/A489).

Figure 5.

Ingenuity pathway analysis (IPA) of common genes associated with nuclear and membrane bile acid receptors (BARs). (a) Venn diagram showing the overlap of genes between each nuclear BAR from a list of the top 7,500 differentially correlated genes (DCGs). Positively and negatively correlated with each nuclear BAR were identified based on Spearman correlation coefficient, and the top 7,500 significantly correlated genes (false discovery rate [FDR] < 0.05) were chosen to obtain a robust list of genes common between pregnane X receptor (PXR), constitutive androstane receptor (CAR), and farnesoid X receptor (FXR). (b and c) IPA analysis of the top 20 upstream regulators (b) and disease and functions (c) predicted from the common gene list. (d) Pathway analysis of the top mechanistic network based on data from both common genes, predicted upstream regulators, and predicted disease and functions. (e) Venn diagram showing the overlap of genes between each membrane BAR from a list of the top 2,500 DCGs identified as described earlier. (f–i) IPA analysis of the common membrane BAR gene list regarding predicted upstream regulators (f), disease and functions (g), and summary top mechanistic network (i).

Regarding disease and functions, common DCGs of nuclear BARs were predicted to be positively associated with lipid and carbohydrate metabolism while being negatively associated with more traditional cancer pathways such as tissue tumorigenesis and epithelial neoplasm (Figure 5e). On the other hand, no disease and function pathways were predicted to be inhibited for membrane BAR DCGs but, overwhelmingly, supported the notion that membrane BARs are associated with oncogenic pathways, especially those related to tumor cell and leukocyte migration/invasion in addition to angiogenesis (Figure 5f). Remarkably, the predicted association between membrane BAR DCGs and immune cell trafficking is in direct agreement with our initial findings of leukocyte infiltration of tumor samples (Figure 1e), suggesting potential involvement of BARs in coordinating inflammatory responses in GC tissue.

Finally, we assessed potential mechanisms by which nuclear and membrane BARs might promote the listed disease and function by analyzing top functional networks of common DCGs. For nuclear BARs, the top functional network predicted the activity of terpenoid metabolism and transport of lipids and other molecules based primarily on the stimulation of NR5A1 by DCGs including various CYP and APO proteins (Figure 5g). For membrane BARs, a more extensive functional network was predicted, which pointed to the primary activation of TGFβ1, TEAD1–4, and myocardin Related Transcription Factor B (MRTFB) and inhibition of SPDEF, overall resulting in enhanced tumor cell and leukocyte migration and inflammation and nucleotide synthesis (Figure 5h). When taken together, these results point to a protumorigenic role for membrane BARs and a weak tumor-suppressor/metabolic role for nuclear BARs, the opposition of which being exemplified in the comparison of predicted canonical pathways as seen in Supplementary Figure S6 (see Supplementary Digital Content 1, http://links.lww.com/CTG/A486).

DISCUSSION

DGR is a known risk factor of developing GC (8–10). Refluxed BAs might induce oncogenesis through direct injury to gastric epithelial cells or stimulation of BARs, the latter being increasingly recognized as novel drivers of GI carcinogenesis. Although a few nuclear and membrane BARs have been investigated in GC, their overarching importance remains unclear. Through meta-analysis of multiple publicly available data sets, we report distinct transcriptional profiles of nuclear and membrane BARs in GC (see Visual Abstract, http://links.lww.com/CTG/A490).

Initially, we found global downregulation of BARs in tumor tissue, but further analysis identified conserved expression patterns among histologic and molecular subtypes, which closely corresponded to clinicopathology and prognosis. Traditionally, GC is classified as either intestinal or diffuse type based on histologic features. Although both subtypes share similar risk factors, tumor behavior and pathophysiology are distinct. For example, intestinal-type cancers are more common and follow a standard progression model starting from atrophic gastritis, resulting in IM and eventual adenocarcinoma transformation. Indeed, IM typically results from inflammation downstream of chronic H. pylori infection but is also known to be caused by BA-associated chemical gastritis. Compared with diffuse-type cancers, intestinal-type cancers are usually less aggressive, are of lower grade, and present at earlier stages (68). Indeed, BA-associated chemical gastritis is known to promote IM through cytotoxicity-induced inflammation (13) and stimulation of CDX2 signaling through FXR. In this analysis, we found that expression of nuclear BARs was associated with lower grade and tumor stage. In addition, we found that nuclear BAR expression patterns were more closely associated with intestinal-type cancer than with membrane BARs. Importantly, we found that many shared nuclear BAR-related DCGs included hallmark genes of IM and adenocarcinoma such as CDX1 and CDX2, FABP1 and FABP2, CDH17, cyclin D1, MyoA1, ERBB2, and Villin-1 (69). IPA analysis of these common DCGs identified upstream regulators that suggested a role in intestinal-type carcinogenesis. For example, the predicted activation of hepatocyte nuclear factor-4 alpha (HNF4α) (70) and inhibition of TP53-mediated signaling (71) is indicative of a role in IM and adenocarcinoma transformation, respectively. In addition, the predominant metabolic signature associated with nuclear BARs supports the high expression of these receptors found in metabolic-type GC, which is linked with a type of gastric IM termed spasmolytic polypeptide-expressing metaplasia (59). Furthermore, analysis of both disease and functions and top networks show that nuclear BAR activity centers around physiologic roles more typical of enterocytes (i.e., transport and metabolism of lipids and carbohydrates) than gastric epithelial cells.

On the other hand, diffuse-type cancers are less common, have no known precursor lesions, and do not have a defined mechanism (72). These cancers are more aggressive, occur at younger age, carry a poor prognosis, and typically present at later stages (68). Similar to intestinal-type cancers, BAs have also been implicated in the pathogenesis of these tumors. For example, gastric stump cancers at anastomotic sites, which occur because of BA reflux arising from Billroth 2 construction postgastrectomy (73), are more histologically and phenotypically related to diffuse-type GC (74,75). Indeed, meta-analysis revealed elevated expression of all membrane BARs in diffuse-type cancers that was consistently associated with poor prognostic and clinicopathologic indicators, most notably being younger age of onset. In addition, we found that all membrane BARs were positively correlated with both hallmark genetic signatures (DTNA, IDS, GPR161, RHOQ, and TSHZ2) and phenotypic signatures (tumor cell migration, invasion, and angiogenesis) of diffuse-type GC (55). In lieu, IPA analysis of this list identified multiple EMT-related genes including TCF4, ZEB1, MMP2, and FGFR1 as top contributors to the functional network (67). Most importantly, the predicted activation of TGFβ1 and TGFβ2, which are key mediators driving tumor cell proliferation and metastasis (76), and inhibition of CDH1 (77), which is intimately linked with loss of cell–cell adhesion typical of this tumor type, strongly suggests a role for membrane BARs in diffuse-type GC development.

The dichotomy between nuclear and membrane BARs with intestinal- and diffuse-type GC was further reinforced when expression patterns were additionally stratified by recently described molecular subtypes. Among CpG island methylator phenotype-classified cancers, we found that membrane BAR expression was elevated in GS-type tumors and positively associated with GS-related molecular signatures such as MSS and RHOA mutations and CLDN18–ARGHAP fusion event. Compared with the other molecular subtypes, GS-type tumors are the most aggressive, carry the poorest prognosis (56), and are closely associated with diffuse-type tumors (78,79). We additionally found that membrane BAR expression was lowest in MSI-type tumors and negatively associated with MSI-related signatures including MLH1 hypermethylation, genomic hypermutation, and ARID1A mutation (80,81). In contrast to GS tumors, MSI-type tumors are associated with more favorable prognoses and are enriched in intestinal-type cancers (31,68,82). Moreover, analysis of nuclear and membrane BAR expression between Lei molecular subtypes revealed that nuclear BARs were enriched in metabolic-type tumors, whereas GPBAR1 and CHRM2 were enriched in mesenchymal-type tumors. These expression patterns are important because they indicate relative chemotherapy sensitivity and provide insights into tumor phenotype. For example, metabolic-type tumors primarily resemble premalignant GC lesions, are intestinal subtype enriched, and indicate potential sensitivity to 5-FU treatment (59). Similarly, mesenchymal-type tumors are diffuse-type enriched, associated with enhanced stem cell properties including migration, invasion, and adhesion, and are potentially sensitive to PI3K-mTOR inhibitors (59).

Because atrophic gastritis resulting from chronic H. pylori infection is the most common cause of both intestinal- and diffuse-type GC (1), combined with the fact that DGR might be caused by gastric dysmotility arising from gastritis-associated somatostatin inhibition (11), we further investigated this relationship. First, we noted that GPBAR1, CHRM2, and FXR expression is elevated in tumors during concurrent infection. Second, we also found that membrane but not nuclear BARs were enriched in GC samples originating from regions with high GC mortality rates such as Eastern Europe and Asia, compared with regions with low GC mortality rates including North America and Europe, independent of ethnicity. These results are notable because H. pylori infection rates correspond with regional GC prevalence, incidence, and mortality rates (2). For example, elevated BAR expression in the Asia cohort is of note because GC is the second leading cause of cancer-related mortality due to greater H. pylori seroprevalence (83) and virulence (2). Although no mechanistic link between H. pylori and BARs has yet been proposed, H. pylori-associated IM and adenocarcinoma progression involves many of the same pathways targeted by FXR (20), GPBAR1 (84–87), and S1PR2 (88). It is possible that membrane BARs either facilitate or sustain gastritis-associated inflammation. For example, previous reports show a prominent role for membrane BARs in mediating inflammatory processes (89,90) and coordinating leukocyte trafficking (91). This notion is not only congruent with our IPA results, indicating a strong immunologic and chemotactic component, but also might explain the conserved positive association between membrane BAR expression and GC tissue leukocyte infiltration.

Despite the inclusion of a robust list of data sets following a strict analytic pipeline, a major limitation of this study comes from the reliance on publicly available data with no follow-up in vitro or in vivo data. Hence, conclusions from this study are not intended to delineate mechanisms or impute clinical or prognostic significance but, instead, build on what is known and provide rationale to spur further investigations into the role BARs play in GC tumorigenesis. Collectively, through meta-analysis and bioinformatics of various public data sets, we identified distinct transcriptomic profiles of BARs in GC. To our knowledge, this is the first study to take a global approach to predicting the differential role these BARs play in GC and their relation to prognosis and clinicopathology. Of note, the potential implications for membrane BARs in diffuse-type tumors are of particular interest because molecular mechanisms underlying diffuse-type tumor progression are lacking. Based on the results of this study, we propose that, in response to refluxed BAs, stimulation of nuclear BARs by more hydrophobic BAs might promote IM and adenocarcinoma transformation, whereas stimulation of membrane BARs by more hydrophilic BAs might promote diffuse-type GC tumorigenesis. Ultimately, this analysis supports the growing consensus that BARs are prognostically relevant biomarkers and important mediators of GC tumorigenesis.

CONFLICTS OF INTEREST

Guarantor of the article: Deborah A. Altomare, PhD.

Specific author contributions: M.R.: designed and conceptualized the study. M.R., J.A., S.L., and S.P.N.: facilitated acquisition and analysis of data. M.R., J.A., S.L., and S.P.N.: interpreted data. M.R.: performed statistical evaluation. M.R. and T.R.R.: prepared figures. M.R., T.R.R., and D.A.A.: wrote and revised manuscript text. T.R.R., D.H., and D.A.A.: critically evaluated the manuscript. All authors contributed to draft revisions and have approved the final version of the manuscript.

Financial support: Research reported in this publication was supported by the National Cancer Institute of the National Institutes of Health under Award Number F30CA257492. The content is solely the responsibility of the authors and does not necessarily represent the official views of the National Institutes of Health.

Potential competing interests: None to report.

WHAT IS KNOWN

✓ Gastric cancer (GC) is the fifth most common and third deadliest cancer worldwide.

✓ Risk factors of GC stem are heterogenous and include duodenogastric reflux (DGR).

✓ DGR-mediated bile acid (BA) reflux is known to simulate chemical gastritis, increasing risk of GC.

✓ BAs facilitate tumorigenesis through direct cytotoxicity and interaction with nuclear and surface receptors.

✓ Bile acid receptors (BARs) are known oncogenes in gastrointestinal cancers and might underlie DGR-associated GC.

WHAT IS NEW HERE

✓ Membrane, more so than nuclear, BARs are important prognostic markers of GC.

✓ High expression of membrane BARs associates with advanced clinicopathology and poor prognosis.

✓ Membrane BAR expression is enriched in diffuse, genome-stable, and mesenchymal GC subtypes.

✓ Nuclear BAR expression is enriched in intestinal and metabolic GC subtypes.

✓ Geographical expression patterns of membrane BARs follow H. pylori infection and mortality rate.

✓ Nuclear BAR activity is tied with enterocyte-like functions such as lipid and carbohydrate transport.

✓ Membrane BAR activity is predicted to signal primarily through TGFβ1 and 2 and other oncogenes.

TRANSLATIONAL IMPACT

✓ Assessing global BAR expression patterns in GC will enhance our understanding of the role BAs plays in GC tumorigenesis and their influence on tumor behavior and prognosis.