Mechanistic Clues Provided by Concurrent Changes in the Expression of Genes Encoding the M1 Muscarinic Receptor, β-Catenin Signaling Proteins, and Downstream Targets in Adenocarcinomas of the Colon

Muscarinic receptors (MRs) in the G protein-coupled receptor superfamily are recipients and mediators of parasympathetic neural transmission within the central and enteric nervous systems. MR subtypes, M1R–M5R, encoded by CHRM1-CHRM5, expressed widely throughout the gastrointestinal (GI) tract, modulate a range of critical, highly regulated activities in healthy tissue, including secretion, motility, and cellular renewal. CHRM3/M3R overexpression in colon cancer is associated with increased cell proliferation, metastasis, and a worse outcome, but little is known about the role of the other four muscarinic receptor subtypes. To address this gap in knowledge, we queried the NCI Genomic Data Commons for publicly available TCGA-COAD samples collected from colon tissue. RNA-seq data were collected and processed for all available primary adenocarcinomas paired with adjacent normal colon. In this unbiased analysis, 78 paired samples were assessed using correlation coefficients and univariate linear regressions; gene ontologies were performed on a subset of correlated genes. We detected a consistent pattern of CHRM1 downregulation across colorectal adenocarcinomas. CHRM1 expression levels were positively associated with those for APC and SMAD4, and negatively associated with CTNNB1, the gene for β-catenin, and with coordinate changes in the expression of β-catenin target genes. These findings implicating CHRM1/M1R as an important deterrent of colon cancer development and progression warrant further exploration.


INTRODUCTION
Despite meaningful advances in colon cancer prevention and treatment in the United States, metastatic disease remains a highly prevalent cause of morbidity and mortality (Siegel et al., 2021). To address this concern, more must be learned about the molecular mechanisms that modulate progression from primary colon cancer to locoregional and distant metastases. An emerging area, and the focus of our current research, is the role muscarinic receptors (MRs) and alterations in post-MR signaling play in cancer progression (Ali et al., 2021;Schledwitz et al., 2021a;Tolaymat et al., 2021). The five MR subtypes, M 1 R-M 5 R, respectively, encoded by CHRM1-CHRM5, are G protein-coupled receptors comprised of seven transmembrane helix domains with an extracellular acetylcholine (ACh)-binding region, and extracellular allosteric binding regions for a variety of non-ACh ligands (Maeda et al., 2019;Tolaymat et al., 2021). Of these subtypes, M 1 R and M 3 R, expressed in healthy and neoplastic colon epithelial cells, appear to play prominent roles in disease initiation and progression (Harrington et al., 2010;Cheng et al., 2017;Schledwitz et al., 2021b;Tolaymat et al., 2021). M 3 R activation results in pro-proliferative downstream signaling, and its activity is most clearly associated with colon cancer progression. In mice, functional M 3 R deficiency, whether achieved by ablating Chrm3 or more recently by selective M 3 R antagonism, reduces tumor number and attenuates cancer cell invasion (Raufman et al., , 2011Hering et al., 2021). In contrast, M 1 R activation may protect against colon cancer, although its actions are understudied relative to M 3 R . In mouse models of colon cancer, Chrm1/M 1 R deficiency does not reduce tumor size or number. Notably, in mice with combined Chrm1/M 1 R and Chrm3/M 3 R deficiency, M 1 R deficiency outweighs the effects of M 3 R deficiency; that is mice with combined M 1 R and M 3 R deficiency have as many colon tumors as control mice (Cheng et al., 2014). The molecular biology underlying the differential actions of M 1 R and M 3 R in colon cancer remains obscure.
β-catenin signaling is a critical player in colon cancer pathogenesis; over 90% of sporadic human colon cancers exhibit mutations in at least one protein in the β-catenin pathway -∼10% exhibit activating mutations in β-catenin and ∼80% have inactivating mutations in APC, a facilitator of cytoplasmic β-catenin degradation (Muzny et al., 2012;Cheng et al., 2019). Our previous work supported the presence of functional interactions between MR and β-catenin signaling (Raufman et al., 2011). The Apc Min/+ mouse model generally recapitulates the effects of human β-catenin pathway mutations, although tumors predominate in the small intestine rather than colon (Boivin et al., 2003). Compared to Apc Min/+ littermate control mice, we observed that Chrm3/M 3 R-deficient Apc Min/+ mice had ∼70% fewer small intestinal adenomas (Raufman et al., 2011). Moreover, nuclear staining for β-catenin was significantly reduced in adenomas from Chrm3/M 3 Rdeficient Apc Min/+ mice (Raufman et al., 2011). The mechanistic underpinnings of this interaction between MR and β-catenin signaling were unclear.
To shed light on the role of MR gene expression in colon cancer, we explored a subset of publicly available TCGA RNA-seq data for colon adenocarcinoma and matched normal tissue samples. We found CHRM1 expression was downregulated in adenocarcinomas and CHRM1 expression levels correlated with changes in the expression of β-catenin signaling pathway genes and their downstream targets. Gene Ontology analysis applied across all genes associated with CHRM1 suggested the effects of CHRM1 downregulation in colon cancer are not subtle, supporting a conceptual framework in which CHRM1 plays a protective role against colon cancer development.
All experiments were approved by the Institutional Animal Care and Use Committee of the University of Maryland, Baltimore, and by the Research and Development Committee at the VA Maryland Health Care System.

RNA-Seq Data Accession and Processing
Data publicly available on the National Cancer Institute (NCI) Genomic Data Commons were downloaded from the TCGA-COAD analysis of tissue samples of colorectal cancer (Grossman et al., 2016). A subset of these data comprised of paired adenocarcinoma and adjacent normal colon samples were selected for further analysis. Fragments per kilobase of transcript per million mapped reads (FPKM) normalized RNA-seq count data and metadata were collected for 78 primary adenocarcinoma samples with an adjacent normal control (Supplementary Figure 1A). Count tables and metadata were pooled into an expression set and the following processing steps performed using the "hpgltools" package in R (Belew and Hughitt, 2018). Expression data were log2 transformed, quantile normalized, and underwent filtering for low count genes (using the default cbcb filter) prior to SVAseq batch correction (Supplementary Figure 1B). Principal component analysis (PCA) was used to assess appropriateness of the pipeline, and clustering based on condition was present following the pipeline described above (Supplementary Figure 1C).

Differential Expression and Model Analysis
Differential expression profiles were generated for the processed adenocarcinoma and matched normal tissue counts using basic log2 fold change calculations, edgeR, and limma-voom analyses for each set of tissues via the "hpgltools" package (McCarthy et al., 2012;Ritchie et al., 2015;Belew and Hughitt, 2018). An adjusted p-value using Benjamini Hochberg FDR was applied (Benjamini and Yosef, 1995). Univariate mixed linear models comparing the final processed counts of various genes and CHRM1, controlling for participant ID, were generated using the nlme package in R (Pinheiro et al., 2021). R 2 -values and residual standard error (RSE) were assessed for all models to ensure generation of a reliable model; residuals were assessed to ensure they were homoscedastic and independent and identically distributed (i.i.d.). The "MuMIn" package was used to generate R 2 -values for the models (Bartoń, 2020). For the initial models described in Figure 1, the conditional R 2 s ranged from 0.00926 to 0.2330, where a conditional R 2 indicates the variation captured by the model (i.e., CTNNB1 and CHRM1 levels, accounting for variation among samples), while the marginal R 2 describes only the variance explained by the fixed effects (i.e., CTNNB1 and CHRM1 levels only). Three of the six models produced conditional R 2 s suggesting small relationships despite the statistical significance: conditional R 2 s of 0.00926, 0.0317, and 0.115 (Axin1, GSK3b, and FZD1, respectively), resulting in no further analysis, and thus in three final models for APC, SMAD4, and CTNNB1 (Equations 1-3): √ Residuals, RSE, and R 2 of the models described in Table 3 were probed to assess model fit. Residuals were found to be independent, normally distributed, and homoscedastic. Individual Pearson correlation coefficients were all ≥ 0.30, and all models had R 2 within a reasonable range given the complex interplay and network of genes that is interacting at any given time. Although RSEs were not low, the associative (vs. predictive) intent of our models and the relatively small denominators altogether with the other results suggests the models are reasonable.
The RSEs of these ranged from 35.74 to 36.62%, where an RSE describes the typical variance from a model, altogether suggesting Equations 1-3 are predicative but imperfect models ( Table 1). "Heatmap.2" was used to generate the heatmap in Figure 2A (Warnes et al., 2020) while the "ggplot2" package was used for all other portions of Figures 1-4 with the exception of Figures 4B,C (Wickham, 2016). The "stats" package in R was used for t-tests and to generate Pearson correlation coefficients (R Core Team, 2021). Gene Ontologies were performed using the "gprofiler2" package in R (Kolberg et al., 2020).

Statistical Analysis for Animal Data
Animal data are presented as means ± SE and were analyzed by two-tailed unpaired Student's t-test. P < 0.05 were considered statistically significant.
FIGURE 1 | Changes in CHRM1 mRNA expression levels are associated with changes in the expression levels of APC, CTNNB1, and SMAD4. We used univariate mixed models with participant ID as a random effect to model processed (low count filtered, quantile normalized, log2 transformed, and batch corrected) CHRM1 expression data as a function of key genes in colorectal cancer. CHRM1 expression levels was significantly associated with APC (A), CTNNB1 (encoding β-catenin) (B), and SMAD4 (C) expression levels, which showed differential clustering of values based on cancer status.
Frontiers in Physiology | www.frontiersin.org FIGURE 3 | CHRM1 expression is associated with multiple β-catenin signaling target genes. We modeled CHRM1 processed expression data as a function of the relevant β-catenin signaling and target genes using mixed (univariate) models, with participant ID as a random effect. β-catenin target genes previously associated with colon cancer and having a significant positive (blue symbols) or negative (red symbols) association with CHRM1 are shown.
FIGURE 4 | Gene pathways differentially expressed in colon adenocarcinomas associated with changes in CHRM1 expression levels. (A) We used a series of univariate mixed models (participant ID as random effect) to assess the relationship between expression levels of CHRM1 and other genes in the dataset, with 4,724 genes demonstrating statistically significant (p-value ≤ 0.01) negative associations and 4,463 demonstrating significant positive associations with CHRM1. (B) We detected 3,921 genes that had a model with a surrogate effect magnitude ≥ 0.25 with KEGG data available. Seventy-seven GO terms were upregulated in adenocarcinoma samples-the top 10 are displayed. (C) Twenty-five GO terms were downregulated in adenocarcinoma samples.

Chrm3 Gene Ablation-Induced Attenuation of β-Catenin Signaling and Intestinal Neoplasia in Mice Identifies Crosstalk Between Muscarinic Receptor and β-Catenin Signaling
Previous evidence supporting a functional interaction between MRs and β-catenin signaling consisted of finding fewer and smaller adenomas in the small intestines of Apc Min/+ mice deficient in Chrm3/M 3 R (Raufman et al., 2011). Moreover, we observed diminished nuclear staining for β-catenin in these small intestinal adenomas (Raufman et al., 2011). Although Apc Min/+ FIGURE 5 | M 3 R deficiency attenuates colon tumor formation in Apc min/+ mice with aberrant β-catenin signaling. In Apc min/+ mice with aberrant β-catenin signaling resulting in exuberant intestinal polyposis, M 3 R deficiency attenuates colon tumor formation. Apc min/+ mice with varying M 3 R expression were created as described in Methods. We euthanized mice at 14 weeks of age, harvested colons, and counted and analyzed tumors. Each symbol represents one Apc min/+ mouse; bars show means. mice are generally considered an animal model of genetic colon cancer, β-catenin signaling is almost equally important in sporadic colon cancer, where approximately 90% of cancers harbor mutations in genes involved in β-catenin signaling. Collectively, these previous findings suggested an important role for MR signaling in the translocation of β-catenin, a key transcription co-factor in colon cancer, to its site of activity in the cell nucleus. Activated β-catenin signaling is known to upregulate the transcription of numerous genes that promote colon neoplasia (Cheng et al., 2019;Chen et al., 2021).
To determine if M 3 R deficiency also attenuated neoplasia in Apc Min/+ mouse colons, we analyzed colon tumor data recorded during our previous study (Raufman et al., 2011). This previously unpublished analysis revealed that besides impacting the production of small intestinal adenomas, in Apc Min/+ mice M 3 R deficiency also robustly attenuated colon neoplasia; we detected colon tumors in all six Apc Min/+ mice expressing M 3 R but in only one of eight mice with homozygous M 3 R deficiency (p = 0.005; Figure 5). These findings provide strong support for the concept that functional crosstalk between MR and β-catenin signaling robustly impacts the progression of colon neoplasia.

CHRM1 Expression Is Attenuated in Human Colorectal Adenocarcinomas
Comparing muscarinic receptor subtype expression levels in adenocarcinoma vs. normal colon samples, we consistently observed reduced CHRM1, CHRM2, and CHRM4 and increased CHRM3 RNA levels (log2fc normal relative to cancer tissue; Figure 6). Following low count filtration, CHRM5 expression data were removed from our analysis. Only CHRM2 and CHRM4 levels were found to have an adjusted p ≤ 0.01 for log2fc across all three analytical methods, while CHRM1 had an adjusted p ≤ 0.01 for all methods except edgeR ( Table 2). Since normal tissue samples were derived from full thickness colon tissue samples and M 2 and M 4 muscarinic receptor subtypes are predominantly localized in intestinal smooth muscle, we were not surprised to observe reduced CHRM2 and CHRM4 expression levels in samples of adenocarcinomas which derive from intestinal epithelial cells; in fact, these findings provided FIGURE 6 | Downregulated expression of CHRM1, CHRM2, and CHRM4 and upregulated expression of CHRM3 mRNA levels in colon cancer. Log2fc was assessed for CHRM1-4 levels in adenocarcinoma compared to normal colon tissues using basic fold change, limma-voom, and edgeR normalization. Changes in CHRM2 and CHRM4 expression were statistically significant (adjusted p ≤ 0.05) across all three assessment methods, while changes in CHRM1 were significant when assessed by basic fold change and limma models; changes in CHRM3 were only significant using basic log2fc. ***p ≤ 0.001, *p ≤ 0.05. an internal validation of our methodology. Since M 1 R and M 3 R are expressed by both normal and neoplastic intestinal epithelial cells, we considered CHRM1 and CHRM3 more disease-relevant candidates for further exploration. Based on these considerations, we explored the relationship between reduced expression of CHRM1, increased expression of CHRM3, and the expression levels of genes commonly mutated in colorectal cancer.
CHRM1 Expression Levels Are Associated With APC, SMAD4, and Inversely Associated With CTNNB1 Expression Levels We first explored the relationship between CHRM1 levels and APC, KRAS, TP53, and SMAD4 RNA levels by generating a series of univariate linear mixed models with participant ID as a random factor. The relationship describing CHRM1 expression values as a function of KRAS changes was not statistically significant. Whereas the model of CHRM1 expression as a function of TP53 technically showed a significant association, the scatterplot revealed a bimodal distribution of cancer samples along the axis for TP53, suggesting the relationship between CHRM1 and TP53 expression is inconsistent in colorectal cancer. Thus, even if a relationship exists, it appears to be independent of cancer status.
In contrast, the model using APC expression levels yielded a statistically significant relationship, as did the model using SMAD4 (p ≤ 0.01) ( Table 3, Equations 1 and 2, Figures 1A,C), prompting us to investigate the relationship between CHRM1 levels and other genes relevant to β-catenin signaling, particularly since SMAD4 activates this pathway (Hussein et al., 2003;Romero et al., 2008;Freeman et al., 2012;Du et al., 2020). We detected a statistically significant linear relationship between CHRM1 and Axin1, CTNNB1 (gene encoding β-catenin; Figure 1B), FZD1 (gene encoding the Wnt receptor), and GSK3b RNA levels ( Table 3) using the stringent p-value of effect size (beta coefficient) but did not detect statistically significant relationships between CHRM1 and WNT, GSK3a, CK1α, DVL1, LRP5, or LRP6. A combination of R 2 -values, RSE, and assessment of residuals as well as simple Pearson correlation used to assess the models, identified three genes of interest-APC, SMAD4, and CTNNB1 (Equations 1-3).
Performing the same analysis using CHRM3, we detected statistically significant linear relationships between increasing Log2fc was assessed for CHRM1-4 levels in adenocarcinoma vs. solid normal colorectal tissue samples using basic fold changes, limma-voom, and edgeR normalization, respectively, employing the processing scheme illustrated in Supplementary Figure 1B. LRP5 and Axin1 RNA levels, decreasing GSK3b RNA levels, and upregulated CHRM3 expression. However, none of the models were i.i.d. with conditional R 2 s ≥ 0.15 and correlation coefficients of ≥ 0.30. As a result of this analysis, moving forward we focused on gene associations with CHRM1 expression levels.

Downstream Genes in the β-Catenin Signaling Pathway Predict Cancer Status
To clarify the relationship between CHRM1 and components of the β-catenin pathway and gauge the potential for functional MR-induced effects on β-catenin signaling, a subset of 72 genes previously identified to be downstream of β-catenin signaling and curated in our dataset were analyzed in normal colon vs. primary colon adenocarcinoma tissues (Herbst et al., 2014). Statistically significant fold changes (adjusted p ≤ 0.01) were observed for 58 genes; 23 were downregulated and 35 were upregulated in cancer (Supplementary Table 1). We generated a heatmap of expression levels for those β-catenin-related genes which demonstrated clustering by cancer status (normal vs. adenocarcinoma tissue) using unsupervised clustering, confirming this gene set predicted cancer status (Figure 2A). To examine the potential variation across each individual gene, we generated a plot of the normalized counts for each sample and each gene ( Figure 2B). Whereas the variation in expression of some genes was greater in normal compared to cancer samples, most showed a wider range Univariate linear mixed models with participant id incorporated a random factor generated to model the relationship between CHRM1 and APC, CTNNB1, and SMAD4 levels were statistically significant. FPKM normalized and processed (log2 transformed, low count filtered, quantile normalized, and SVAseq batch corrected) expression levels were used for model generation.
of variation in expression in the cancer samples, suggesting some genes were upregulated and playing a role in a subset of the cancers, but there is variability across samples. We used the paired Student's t-test to analyze these findings, comparing the ranges of normalized values; this revealed a mean increase in variation of 0.7348 in adenocarcinoma (95% CI = 0.4672-1.0023).

Association of CHRM1 Expression Levels With Those of Genes in the β-Catenin Signaling Pathway
Using a series of univariate regressions, we explored the relationship between expression levels for the subset of the aforementioned 72 genes and CHRM1. Statistically significant (adjusted p-value ≤ 0.01) increased RNA expression for 10 (AC2, AC9, ITF2, AC6, CASL, ALEX1, TIAM-1, AC5, FAAR, and FLVI2/BMI1) and decreased RNA expression for 16 genes (BMP4, FOSL1, BCL1, B2T, CMD1, AUTS9, ADCADN, CD87, API4, CCL28, VEGF, LEF-1, FAN1, NMA, CFND, and SFRS3) were associated with increasing CHRM1 levels, adjusting for participant ID with a regression coefficient of ≥ 0.1. Regression coefficients were used as surrogate effect sizes given the univariate nature of the mixed model which only controlled for participant ID. Expression levels of an additional seven genes (MRTL, ASH2, DKK-4, CLD1, MMP-7, ATF, and FEX) were more modest decreased in association with increasing CHRM1 levels (Figure 3 and Supplementary Table 2). After expanding the univariate regression analysis to all 36,987 genes, we found 9,187 genes were significantly related to CHRM1expression levels-of these, 5,776 had a surrogate effect size with a magnitude ≥ 0.25 (Figure 4A and Supplementary  Table 3). Several are pseudogenes, or long noncoding (lnc)RNAs. To identify pathways associated with downregulated CHRM1 expression in colon cancer we performed a Gene Ontology (GO) analysis on the 5,776 gene subset which also had an adjusted p-value of ≤ 0.01 and had available KEGG data (3,921 genes). Of 25 GO terms downregulated in cancer, the following were prominent: phospholipid binding, S/T kinase activity, nucleoside triphosphatase regulator activity, GTPase regulator and activator activity, guanyl nucleotide exchange factor activity, Ras guanyl nucleotide exchange factor activity, and PI and PIP binding ( Figure 4C). Of 77 GO terms upregulated in cancer, the following were prominent: catalytic activity on DNA/RNA, histone binding, nuclease and methyltransferase activity, DNA/helicase activity, and single stranded DNA and RNA as well as damaged DNA binding ( Figure 4B).

DISCUSSION
Previous studies implicated MRs in the progression and spread of many cancer types, including colon cancer (Shah et al., 2009). To date, no studies attempted to discern the relationship between MR subtype expression and key players in colon cancer. We sought to fill this gap in knowledge using an unbiased approach to explore expression of MR subtypes in publicly available RNA-seq data. We found CHRM1 levels were downregulated in colon adenocarcinoma, while, as expected, CHRM3 was upregulated. Further, we identified an imperfect but linear relationship between CHRM1 and APC and CTNNB1 mRNA levels, supporting a relationship between CHRM1/M 1 R and β-catenin signaling in colon cancer. APC plays a key role in the β-catenin destruction complex to ensure that, in the absence of activation by Wnt ligands, cytosolic β-catenin is ubiquitinated and degraded (Parker and Neufeld, 2020). Notably, Ras degradation is also regulated by GSK3β, a member of the β-catenin destruction complex (Lee et al., 2018). As M 1 R activation results in downstream Ras/MAPK activation (Qian et al., 1995), this illuminates a potential node of intersection between M 1 R and β-catenin signaling.
Although M 3 R activation results in similar downstream cascades to those modulated by M 1 R activation, M 3 R expression is upregulated in colorectal adenocarcinoma and promotes cancer progression. The divergent effects of M 1 R vs. M 3 R activation on colon cancer were described previously  and remain poorly understood. Notably, other differences between the expression and actions of these MR subtypes are reported. For example, M 1 R and M 3 R differ in their propensity to engage in cAMP signaling via non-canonical pathways. Following treatment with a muscarinic agonist, Chinese hamster ovary (CHO) cells expressing only M 1 R maximally accumulate four times as much cAMP as CHO cells expressing only M 3 R, despite similar receptor density, pharmacokinetics, and efficiency of signaling via phospholipase C, inositol triphosphate, and calcium. Investigation of the underlying mechanism suggested post-M 1 R stimulation of cAMP accumulation involves activation of G s (Burford et al., 1995;Burford and Nahorski, 1996). As β-catenin is a target of phosphorylation by PKA, which is activated by cAMP (Katoh and Katoh, 2017), reduced expression of CHRM1/M 1 R in colon cancer would attenuate β-catenin signaling mediated by this non-canonical M 1 R/cAMP/PKA/β-catenin axis. This differential impact of M 1 R and M 3 R activation on cAMP signaling adds additional complexity and potential fine tuning to the impact of MR expression and activation on β-catenin signaling.
It is possible that differential localization of MRs within cells comprising the tumor microenvironment also plays a role in discrepant effects; M 3 R is more consistently expressed throughout multiple tissue layers, including circular and longitudinal muscle, myenteric nerve cell bodies, and mucosal epithelial tissue, while M 1 R expression is largely restricted to myenteric and submucosal nerve cells and epithelial cells (Harrington et al., 2010). In support of this, the wide variation in identified expression levels in CHRM3 in colon cancer, specifically across a range of experimental methods and cell lines, suggests tissue localization may play a role in the observed changes (Felton et al., 2018;Ali et al., 2021). Further, unlike the other muscarinic receptor subtypes, both M 1 R and M 3 R undergo pre-coupling with G i/o G-proteins, that are their non-preferential G proteins. Differential pre-coupling in different tissues could result in variable downstream effects following MR activation (Jakubík et al., 2011).
Our exploration of the relationship between CHRM1 and genes known to be differentially expressed downstream of β-catenin (Figure 3) was similarly thought-provoking. Despite great variation across expression levels of those genes, particularly in cancer tissue, they remained predictive of cancer status. By exploring the relationship between those genes and CHRM1 levels, we identified a subset of genes with a statistically significant association. Univariate analyses identifying concurrent changes in gene expression suggested functional interactions. Of particular interest were positive correlations with ITF2 and AC2, which had the highest surrogate effect sizes. ITF2 is implicated as both a tumor suppressor and proto-oncogene in colon cancer (Davidsen et al., 2018); relevant to the current analysis, in human colon cancer cells and tissue, ITF2 is reported to prevent activation of the β-catenin-TCF4 complex and transcription of β-catenin gene targets (Shin et al., 2014). Decreased levels of AC2 and the general class of adenylyl cyclases was previously demonstrated in human colon cancer cell lines (Nelson and Holian, 1988;Yi et al., 2018;Fan et al., 2019). These changes in AC2 expression may also be relevant to the non-canonical M 1 R/cAMP/PKA/β-catenin axis described above.
Additionally, we identified a negative correlation with SFRS3, implicated previously in colon cancer (Kuranaga et al., 2018;Zhou et al., 2020). We detected more modest inverse correlations with BMP4 and FOSL1 (but with higher conditional R 2 -values and high p-values); these genes are also implicated in colon cancer cell proliferation and metastasis (Diesch et al., 2014;Yokoyama et al., 2017;Liu et al., 2021). Notably, MMP7, has a well-documented role in colon cancer progression and metastasis (Kitamura et al., 2009;Koskensalo et al., 2011;Polistena et al., 2014), as does VEGF via vasculature remodeling in tumoradjacent tissues (Ellis et al., 2000;Liu et al., 2011;Ahluwalia et al., 2014;Yang et al., 2015). Notably, our group previously observed that M 3 R activation induces MMP7 expression in human colon cancer cells (Xie et al., 2009). Hence, the inverse relationship between CHRM1 and MMP7 expression levels is of particular interest. While we anticipated detecting genes previously implicated in colon cancer, the relative relationship to changes in CHRM1 expression suggests a common relationship between the concurrent changes observed in CHRM1 and APC/CTNNB1, and these downstream β-catenin target genes.
Lastly, we used the relationship between expression levels of genes with a significant relationship to CHRM1 expression to identify gene ontology (GO) terms and pathways of importance to cancer progression. Several GPCR regulationand activation-related GO terms were downregulated, as were several potentially CHRM1 cascade-related GO terms such as PI and PIP binding and guanyl nucleotide exchange factor activity. Although these genes were specifically stratified based on altered levels correlating with changes in CHRM1 expression, it remains interesting that changes in CHRM1 and related genes had sufficient impact to associate with downregulated pathways. While the sheer number of genes with changing expression levels paralleling CHRM1 was high, the fact that the most concurrently downregulated pathways were CHRM1-related while the most concurrently upregulated GO terms were growth-and proliferation-related, implies the strong correlations were unlikely due to chance and reflect the impact of CHRM1 downregulation in colon cancer. This further validates that CHRM1/M 1 R downregulation observed in colon cancer is not only non-random, but very likely involved in important aspects of development and/or progression of disease.
Our findings suggest an important relationship exists between CHRM1/M 1 R and β-catenin signaling in colorectal cancer. Nonetheless, we acknowledge important limitations. Our analytical approach cannot identify a specific mechanism (s) that underlies directionality of gene changes nor a specific cause and effect relationship. Moreover, our results do not exclude the possibility that additional players are involved; additional studies are required to elucidate the processes driving these changes. For example, RNA-seq experiments in Chrm1/M 1 R knockout mice may elucidate the temporal dynamics and confirm the directionality of such changes in gene and protein expression, thus better delineating the mechanisms whereby these pathways intersect.
Collectively, these patterns of changes in gene expression levels and cancer-related pathways support a conceptual framework wherein CHRM1/M 1 R expression contributes to protection against the development and progression of colorectal adenocarcinoma. Further elucidation of these novel mechanistic insights regarding intersection of CHRM1/M 1 R and β-catenin signaling, beyond the scope of the current project, warrant additional experimentation.

DATA AVAILABILITY STATEMENT
The datasets analyzed for this study can be found in the Genomic Data Commons, at https://portal.gdc.cancer.gov/, under the TCGA-COAD project.

AUTHOR CONTRIBUTIONS
MA and J-PR conceptualized the project. KC performed the experiments and analyzed the results. MA, AS, and J-PR wrote the initial draft, proofread, edited, contributed additional material, and completed the final draft. MA, AS, KC, and J-PR approved the final manuscript. All authors contributed to the article and approved the submitted version.