Single-cell RNA Sequencing Reveals How the Aryl Hydrocarbon Receptor Shapes Cellular Differentiation Potency in the Mouse Colon

The large intestine (colon) is made up of diverse cell types with distinct cellular differentiation programming and differentiation trajectories (1, 2). Typically, stem cells replenish the intestinal epithelium every 3 to 5 days, and a constant pool of Lgr5⁺ stem cells is required for intestinal homeostasis (3). This is noteworthy because Lgr5⁺ crypt stem cells are the cell-of-origin of colon cancer, and a stem cell/progenitor cell hierarchy is maintained in early neoplastic lesions (4). Recently, it has been demonstrated that dietary and microbial cues regulate intestinal tumorigenesis in mouse models by targeting the aryl hydrocarbon receptor (Ahr; refs. 5–8). This has been linked to the antagonism of Wnt signaling (6, 9) and the Ahr–FoxM1 axis (6), which mediate colonic stem/progenitor cell behavior. Collectively, these findings suggest that Ahr signaling regulates the intestinal stem cell niche both intrinsically and extrinsically. However, precisely how Ahr and its dietary/microbial ligands interact in terms of stem cell homeostasis in the colonic crypt is still under investigation.

Single-cell analysis is rapidly becoming a valuable tool to dissect cellular heterogeneity and define cell identity in complex systems (10, 11). For example, single-cell analyses have revealed conserved populations and signaling mechanisms associated with colonic epithelial diversity in health and the regenerating intestine (12–15). Thus, we performed single-cell RNA sequencing (scRNAseq) on colonic crypts from wild-type (WT) and Ahr knockout (KO) mice to further elucidate the effects of Ahr on the signaling pathways that are integral to the maintenance and differentiation of epithelial adult stem cells. As part of this effort, single-cell entropy (16, 17) and RNA velocity (18, 19) analyses were used to assess crypt cell overall differentiation potential (potency) and entropy-based measures. In addition, quantitative inference and analysis of intercellular communication networks was performed. Herein, we report that deletion of Ahr elevates differentiation potency, cellular differentiation trajectories (velocity length), and perturbs intercellular signaling cross-talk in most colonic crypt cell types. These results support our premise that Ahr is a potential therapeutic target to recalibrate remodeling of the intestinal stem cell niche.

Animals were housed under conventional conditions, adhering to the guidelines approved by the Institutional Animal Care and Use Committee at Texas A&M University (College Station, TX). Stem cell–targeted Lgr5-EGFP-IRES-Cre^ERT2, Ahr^f/f and tdTomato^f/f mouse strains have all been previously described (5). The mouse genotypes used in this study were Lgr5-EGFP-Cre^ERT2 X Tomato^f/f (WT, control), and Ahr^f/f X Lgr5-EGFP-Cre^ERT2 X Tomato^f/f (Ahr KO). Male mice were fed ad libitum an AIN-76A semipurified diet (Research Diets, D12450B) and housed on a 12-hour light–dark cycle. Littermate controls were cohoused with the KO mice. Mice (n = 5 per genotype, 8–10 weeks of age) were injected intraperitoneally with 2.5 mg of tamoxifen (Sigma, T5648) dissolved in corn oil (25 mg/mL) once a day for 4 consecutive days.

Two weeks after tamoxifen injection, the large intestine was removed, washed with cold PBS without calcium and magnesium (PBS^−/−), everted on a disposable mouse gavage needle (Instech Laboratories), and incubated in 15 mmol/L EDTA in PBS^−/− at 37°C for 35 minutes as previously described (5). Following transfer to chilled PBS^−/−, crypts were mechanically separated from the lamina propria by vigorous vortexing. After dissociation with trypsin, epithelial cells were subsequently filtered through a 40-μm mesh and Tomato-expressing cells (includes GFP⁺/Tomato⁺ as well as GFP⁻/Tomato⁺) were collected using a MoFlo Astrios Cell Sorter (Beckman Coulter), using DAPI to exclude dead cells. Since Tomato⁺ cells represent colonic stem cells and their progeny, we were able to examine the effects of Ahr knock-out on stem cells and all other cell types originating from the Ahr knocked out stem cells. Samples were processed using the 10x Genomics scRNAseq pipeline described below. A total of 62,741 cells from 10 mice were sequenced. These included 34,889 sorted colonocytes from the WT and 27,852 from the KO mice. The average number of genes detected per sample was 20,141. From all sequenced cells, 40,690 (21,263 from WT and 19,427 from KO samples) were removed using criteria developed by the scRNAseq quality control procedure (20). Typically, excluded cells had either a high proportion of mitochondrial reads (greater than 10%) or exhibited an extremely large or small library size.

Single-cell sample preparation was conducted according to Sample Preparation Protocol provided by 10x Genomics as follows: a cell suspension (1 mL) from each mouse genotype was pelleted by centrifugation (400 g, 5 minutes). The supernatant was discarded and the cell pellets resuspended in 1× PBS with 0.04% BSA, followed by two washing procedures by centrifugation (150 × g, 3 minutes). Cells were resuspended in 500 μL 1× PBS with 0.04% BSA followed by gently pipetting 10 to 15 times and enumerated using an Invitrogen Countess automated cell counter (Thermo Fisher Scientific) and the viability of cells was assessed by trypan blue staining (0.4%). Subsequently, single-cell Gel beads in Emulsion (GEM) and sequencing libraries were prepared using the 10x Genomics Chromium Controller in conjunction with the single-cell 3′ kit (v3). Cell suspensions were diluted in nuclease-free water to achieve a targeted cell count of 5,000 for each sample. cDNA synthesis, barcoding, and library preparation were carried out according to the manufacturer's instructions. Libraries were sequenced in the North Texas Genome Center facilities using a NovaSeq6000 sequencer (Illumina, San Diego). For the mapping of reads to transcripts and cells, sample demultiplexing, barcode processing, and unique molecular identifier (UMI) counts were performed using the 10x Genomics pipeline CellRanger v.2.1.0 with default parameters. Specifically, for each library, raw reads were demultiplexed using the pipeline command “cellranger mkfastq” in conjunction with “bcl2fastq” (v2.17.1.14, Illumina) to produce two fastq files: the read-1 file containing 26-bp reads, consisting of a cell barcode and a UMI, and the read-2 file containing 96-bp reads including cDNA sequences. Sequences were aligned to the mouse reference genome (mm10), filtered and counted using “cellranger count” to generate the gene-barcode matrix.

Dimension reduction of expression matrices and cell clustering was performed using t-SNE and k-means clustering algorithms, respectively. Cell type assignment was performed manually using the “SC_SCATTER” function of scGEAToolbox (20). Cell cycle phase assignment was made using the “CellCycleScoring” function in the Seurat R package (21), which utilizes phase-specific marker genes generated by the “cc.genes” dataset (22). Cell differentiation potency was computed using Correlation of Connectome and Transcriptome (CCAT; 16, 17). In addition, differential gene expression was performed using MAST (23) from the Seurat R package (21). Briefly, cells for all the samples from each experimental group were concatenated, normalized using the library size of 10,000 as a scaling factor, and log-transformed as by default in Seurat (21). Labeled cell types were compared across experimental groups to quantify the differences in the level of expression. For each cell type, all the genes expressed in a minimum of 5% of the cells were tested. Following the default in Seurat (21), genes with |log2(FC)| > 0.25 and P value < 0.05 after multiple test correction were considered as differentially expressed. Expression profiles of differentially expressed genes in 10 different cell-type groups were computed. Subsequently, the concatenated list of genes identified as significant was used to generate a heatmap. Genes were clustered using hierarchical clustering. The dendrogram was then edited to generate two major groups (up- and down-regulated) with respect to their change in the KO samples. Identified genes were enriched using Enrichr (24). We subsequently performed an unbiased assessment of the heterogeneity of the colonic epithelium by clustering cells into groups using known marker genes as previously described (25, 26).

Single-cell potency was measured for each cell using the CCAT—an ultra-fast scalable estimation of single-cell differentiation potency from scRNAseq data. CCAT is related to the Single-Cell ENTropy (SCENT) algorithm (27), which is based on an explicit biophysical model that integrates the scRNAseq profiles with an interaction network to approximate potency as the entropy of a diffusion process on the network.

To estimate the RNA velocities of single cells, two count matrices representing the processed and unprocessed RNA were generated for each sample using “alevin” and “tximeta” (28). The python package scVelo (19) was then used to recover the directed dynamic information by leveraging the splicing information. Specifically, data were first normalized using the “normalize_per_cell” function. The first- and second-order moments were computed for velocity estimation using the “moments” function. The velocity vectors were obtained using the velocity function with the “dynamic” mode. RNA velocities were subsequently projected into a lower-dimensional embedding using the “velocity_ graph” function. Finally, the velocities were visualized in the precomputed t-distributed stochastic neighbor embedded (t-SNE) embedding using the “velocity_embedding_stream” function. All scVelo functions were used with default parameters. To compare RNA velocity between WT and KO samples, we first down-sampled WT cells from 12,227 to 6,782 to match the number of cells in the KO sample. The dynamic model of WT and KO was recovered using the aforementioned procedures, respectively. To compare RNA velocity between WT and KO samples, we calculated the length of velocity, that is, the magnitude of the RNA velocity vector, for each cell. We projected the velocity length values with the number of genes using the prebuilt t-SNE plot. Each cell was colored with a saturation chosen to be proportional to the level of velocity length. We applied the Kolmogorov–Smirnov test on each cell type, statistically verifying differences in the velocity length.

Cellular communication analysis was performed using the R package CellChat (29) with default parameters. WT and KO single cell data sets were initially analyzed separately, and two CellChat objects were generated. Subsequently, for comparison purposes, the two CellChat objects were merged using the function “mergeCellChat”. The total number of interactions and interaction strengths were calculated using the “compareInteractions” function. Significant signaling pathways were identified using the “rankNet” function based on the difference in the overall information flow within the inferred networks between WT and KO cells. The enriched pathways were visualized using the “netVisual_aggregate” function.

The data generated in this paper are publicly available in Gene Expression Omnibus (GEO) at GSE167595. The source code for data analyses is available at https://github.com/chapkinlab.

Colons were removed 2 weeks following the final tamoxifen injection. At this timepoint, loss of Ahr potentiates FoxM1 signaling to enhance colonic stem cell proliferation, resulting in an increase in the number of proliferating cells per crypt, compared with WT control (5). In order to define the effects of Ahr deletion on colonic crypt cell heterogeneity, scRNAseq was performed on 19,013 cells, including 12,227 from WT (Lgr5-EGFP-Cre^ERT2 X tdTomato^f/f) and 6,786 from KO (Lgr5-EGFP-IRES-Cre^ERT2 x Ahr^f/f x tdTomato^f/f) mice. Single cells from colonic crypts were sorted using FACS of Cre recombinase recombined (tdTomato⁺) cells (Fig. 1A). Tomato gene expression was detected in approximately 1.8% of cells (Supplementary Fig. S1). As a measure of scRNAseq data quality control, we used a customized mitochondrial DNA threshold (%mtDNA) to filter out low-quality cells by selecting an optimized Mt-ratio cutoff (ref. 30; Supplementary Fig. S2). Numbers of cells obtained from samples before and after quality control filtering of scRNAseq data are shown in Supplementary Fig. S3.

The transcriptomic diversity of data was projected onto two dimensions by t-SNE. Unsupervised clustering identified 10 clusters of cells. Based on known cell-type markers (Supplementary Table S1), these cell clusters were assigned to distinct cell types, namely noncycling stem cell (NSC), cycling stem cell (CSC), transit-amplifying (TA) cell, enterocyte (EC), enteroendocrine cell (EEC), goblet cell (GL, type 1 and 2), deep crypt secretory cell (DCS, type 1 and 2), and tuft cell (Fig. 1B). We observed two distinct subclusters for GL and DCS. Relative proportions of cells varied across clusters and differed between WT and KO samples (Fig. 1C). Notably, the relative abundance of CSC in the KO samples (15.2%) was only approximately half that in the WT samples (28.7%). This apparent discrepancy with previous findings (5) may be attributed to the known GFP mosaicism associated with the Lgr5-EGFP-IRES-CRE^ERT2 model (5) and the initial isolation of tdTomato⁺ cells used in this study. The annotated cell types were also independently defined using cluster-specific genes, i.e., genes expressed specifically in each cluster. Figure 1D demonstrates the 2-D t-SNE plots of WT and KO samples. Figure 1E shows examples of these cluster-specific genes. Some of these cluster-specific genes served as marker genes, which were used for cell-type annotation. For example, Lgr5 was found to be highly expressed in CSCs and NSCs (Fig. 1F).

Samples for each experimental group (WT, n = 5; and KO, n = 5) were pooled to compare the expression level of genes from the same cell type across experimental groups. We used MAST (23) and the Seurat R package (21) to identify genes with |log2(FC)| > 0.25, where FC is fold change, and adjusted P value < 0.05 following multiple test correction. A total of 115 genes exhibited a significant expression change in at least one cell type. As most of these genes showed the same directional change in different cell types, their profiles were concatenated and analyzed jointly. For each of the 115 genes, the log2(FC) values between KO and WT expression across different cell types were assessed. Using the FC profile (i.e., according to whether genes were expressed higher or lower in the KO samples relative to the WT samples), genes were clustered and divided into two major groups: KO upregulated (n = 40; Fig. 2A) and KO downregulated (n = 75; Fig. 2B). No genes were significantly KO upregulated in one cell type, and significantly KO downregulated in another cell type, or vice versa. Enrichment analysis based on Enrichr (24) revealed that the Ahr knockout in colonic crypts induced the overexpression of ribosomal genes or genes related to translation (Rps28, Rps27, Rps29, Sec61g, Rpl37a, Rpl38, Pabpc1, Rpl39, and Rps21; FDR = 4.13e–9), as well as the MAPK/TRK pathway (Egr1 and Fos; FDR = 4.00e–2). Consistent with previous studies (31, 32), most of the known Ahr target genes were modulated in the KO samples (Supplementary Fig. S4). KO upregulated genes included Fos and Hspa1a (Fig. 2C), both targets of Foxm1, suggesting an effect of Ahr deletion on Foxm1-regulated genes. This is consistent with the ability of the Ahr-FoxM1 axis to mediate oncogenic activation (5, 33, 34). The list of KO downregulated genes was enriched with several functions, including cholesterol homeostasis (Lgals3, Fdps, Sqle, Hmgcs1 and Ethe1; FDR = 1.21e–4), oxidative phosphorylation (Ndufb8, Ndufb7, Ndufs7, Cox4i1, Mgst3, Cox5b and Cox5a, FDR = 1.21e–4), and the p53 pathway (FDR = 0.75e–2). The downregulation effect on the p53 pathway is consistent with the ability Ahr to attenuate oncogenic activation (5, 33, 34). In contrast, cytochrome P450 genes, e.g., Cyp1a1 and Cyp1b1, were not affected.

In general, pluripotent stem cells are endowed with the capacity to differentiate into all major cell lineages and thus have a higher entropy/differentiation potency (16). To identify novel stem-or-progenitor cell phenotypes in our scRNAseq data, we utilized the CCAT computational method (16, 17). This approach measures global signaling entropy and can estimate a cell's differentiation potential. Thus, CCAT was applied to measure the stemness of all cell types in an unbiased manner (Fig. 3A). By comparing the potency level across different cell types, we found that NSC, CSC, and TA cells had a significantly higher potency than the other cell types (all P values < 1.05e–10, the Kolmogorov–Smirnov tests between the three high-value cell types versus the other cell types). We subsequently compared the potency levels between different cell types in the WT and Ahr KO samples. The comparisons were performed independently for each of the cell types. Across all cell types, cells in the KO samples tended to have a higher potency compared with WT (Fig. 3B). The differences between the WT and KO samples were highly significant for all cell types except tuft cells (P values of Kolmogorov–Smirnov tests described in Table 1). The same pattern was observed when only G1-phase cells were included in the analysis (Fig. 3C and D; Table 1), ruling out cell cycle as a confounding factor. These findings suggest that the deletion of Ahr elevates differentiation potency in most colonic crypt cell types in the KO samples.

In order to further assess the effects of Ahr KO on cellular differentiation trajectories, we performed RNA velocity analysis. RNA velocity is a time derivative of an individual cell's expression state, which can be used to predict the future state of single cells (18, 19). In the RNA velocity analysis, the ratio of unspliced to spliced mRNA abundance is used to determine the velocity of each cell. For each cell, a velocity vector is computed by combining velocities across genes. The direction of the vector points to the future state of the cell; the length of the vector (or velocity length) indicates the rate of change in global mRNA abundance during the dynamic process of cell differentiation. Using the RNA velocity analysis tool, scVelo (19), we constructed a velocity field map to highlight the cell trajectories that give rise to different cell types (Fig. 4A). The field map, shown as a streamline plot, depicts the dynamics of cell transition from NSC to enterocytes.

More specifically, projections of velocity vectors form a strong directional flow originating from NSC, passing through CSC, and ending at enterocytes. This dynamic process agrees with the lineage relationship that quiescent intestinal stem cells (NSCs) give rise to active CSCs (31), and that CSCs further differentiate into mature differentiated cells such as enterocytes, enteroendocrine cells, and tuft cells (32).

Next, to detect the difference in the cell differentiation rate following Ahr deletion, we compared the average velocity length between single cells from WT and KO samples (Fig. 4B). In order to control for the uneven numbers of cells in the two groups, we subsampled the same number of WT (n = 6,782) and KO cells. An equal number of cells in the two groups ensured that dynamic parameters could be estimated in a comparable fashion. The estimated results showed that, across different cell types, cells from KO samples tended to have significantly greater velocity length than their WT counterparts (all P values < 0.05, Kolmogorov–Smirnov test; Supplementary Table S2). Among them, the differences were found to be most pronounced in NSCs and CSCs (Fig. 5A). The phase portraits derived from the learned dynamics for two representative genes, Notch2 and Ezr, are shown (Fig. 5B and C). The phase portrait of a gene is a scatter plot of inferred unspliced and spliced numbers of the gene across all cell types. Each point in the scatter plot is a cell. The purple dashed line corresponds to the estimated “steady-state” and the purple curve represents the learned dynamics. Cells located in positions that deviate from the steady-state line have either higher or lower mRNA abundance changing rate than expected when cells are in a steady state. By comparing cell position in the two gene phase portraits, we observed that cells in the KO samples tended to be more deviated from and above the steady-state line compared with those in WT samples. This is consistent with an increase in the expression of Notch2 and Ezr in the KO samples, which were upregulated by 19.8% and 13.2%, respectively. Notch2 encodes a protein member of the Notch family that plays a role in a variety of developmental processes. It functions as a receptor for membrane-bound ligands to regulate cell-fate determination, affecting the implementation of differentiation, proliferation, and apoptotic programs such as intestinal stem remodeling, homeostasis (35), and asymmetric division (36). Ezr encodes ezrin, a cytoplasmic peripheral membrane protein that functions as a protein-tyrosine kinase substrate in microvilli. This protein promotes stem cell properties (37), plays a key role in cell surface structure adhesion, migration, and organization, and has been implicated in colon cancer initiation (38).

Among its related pathways are proteoglycans in cancer and cell adhesion_ECM remodeling. Collectively, these data provide additional mechanistic insight into how Ahr modulates colonic crypt homeostasis.

We also analyzed putative ligand-receptor interactions in the colon and modeled cell–cell communication between different cell types, which is critical for maintaining cellular and tissue homeostasis (39). For this purpose, cellular communication analysis was assessed using CellChat (29). Upon comparison of the WT and Ahr KO groups, 129 and 122 significant ligand-receptor pairings were identified, respectively (Supplementary Tables S3 and S4). Global interactions between ligand-receptor pairs were distributed between different cell types and comparisons between WT and KO groups examined with respect to the differential interaction strength between the two groups (Fig. 6A and B).

Interactions showing increased strength in KO cells are given in Fig. 6A. In comparison, interactions showing decreased strength in KO cells are given in Fig. 6B. CellChat allows grouping of multiple individual ligand-receptor pairs into pathways, so we configured detected interactions into signaling pathways and compared the sum of interaction intensity among all ligand-receptor pairings of cell groups (also known as the information flow) between the WT and KO groups (Fig. 6C). While signaling pathways such as stem cell factor receptor c-Kit (KIT), granulin (GRN), epidermal growth factor (EGF), and galectin-9/TIM-3 (GALECTIN) were detected in both WT and KO cells, some pathways were exclusively detected in either WT [e.g., guanylate cyclase activity (GUCA)] or KO [e.g., visfatin signaling pathway (VISFATIN), class 3 semaphorin proteins signaling pathway (SEMA3), and macrophage migration inhibitory factor (MIF)]. For example, three representative pathways, GUCA (WT only), MIF (KO only), and EGF (strength increased in KO), were uniquely impacted by Ahr status. The GUCA pathway was detected only in WT cells (Supplementary Fig. S5A), which suggests in KO cells, the deletion of Ahr downregulated the activity of the GUCA pathway. This is noteworthy, because GUCY2C is a key receptor that mediates homeostatic intestinal barrier function, and suppresses colitis and tumorigenesis (40). GUCA2A and GUCA2B are two predominant ligands in this pathway. In intestinal cancer, the loss of GUCA2A and GUCA2B suppresses GUCY2C signaling early in transformation (41). Interestingly, the MIF pathway was exclusively detected in KO cells (Supplementary Fig. S5B). This is consistent with previous findings indicating that Ahr suppresses pathogenic inflammatory activity (42). During intestinal inflammation, the CD74 signaling receptor for cytokine macrophage migration inhibitory factor is strongly activated (43). Finally, with respect to EGF, Ahr is known to modulate the EGF pathway directly (44). Our results indicate that following Ahr deletion, increased EGFR interactions involving enterocytes were detected (Supplementary Fig. S5C), suggesting a compensatory response. This is noteworthy, because hyperactivation of the EGFR signaling axis is sufficient to drive tumorigenesis (45).

Ahr, a ligand-activated transcription factor, controls the maintenance and differentiation of intestinal stem cells and integrates dietary and microbial cues to modulate crypt homeostasis and colon cancer risk (5, 6, 9). Mounting evidence suggests that enhancement of extrinsic dietary and intrinsic microbial-derived ligands can favorably modulate Ahr signaling and, thus, should be part of the colon cancer prevention armamentarium. Modulation of Ahr signaling is also associated with many chronic diseases, including inflammatory bowel diseases where Ahr expression/activation is protective (46–48). In this study, we provide additional mechanistic evidence demonstrating how the loss of Ahr augments colonic Lgr5⁺ stem cells and nonstem cell differentiation potency and cell fate transitions. The phenotypic plasticity of single cells, defined as the ability to adopt an alternate cell fate in response to perturbation, was estimated in silico from their RNA sequencing (RNA-seq) profile using signaling entropy. As expected, NSC, CSC, and TA cells had a significantly higher potency than the other well differentiated cell types because these cells are largely uncommitted, or undifferentiated (16). Interestingly, intestinal Ahr deletion elevated single-cell entropy (a measure of differentiation potency or cell stemness) in both Lgr5⁺ stem cells (noncycling, cycling) and differentiated cells, e.g., goblet cells and enterocytes. This suggests that Ahr is directly capable of regulating the capacity of committed cells to dedifferentiate into stem cells and potentially promote the regeneration of epithelial cells (49). These findings have broad implications for cancer biology because the accumulation of undifferentiated stem cells is preferentially primed for transformation and often serve as the cells of origin for cancer (50). We also provide evidence of an Ahr-dependent underlying physiologic form of cell plasticity that could be coopted by dedifferentiation and acquisition of stem cell-like properties to induce intestinal tumorigenesis (51). This is consistent with recent studies indicating that Ahr signaling plays a protective role in carcinogen-induced colon cancer, colitis-associated colon tumorigenesis and Apc-dependent mouse models (5–9, 52).

Comparison of RNA velocity in colonic crypt single cells was also used to trace Ahr-driven remodeling of the stem cell niche. RNA velocity has facilitated the study of cellular differentiation in single-cell RNA-seq data. It describes the rate of gene expression change for an individual gene at a given time point based on the ratio of its spliced and unspliced mRNA (18, 19). Interestingly, nearly all cell types, e.g., Lgr5⁺ stem cells, EC, goblet cells, EEC, and tuft cells, had a significantly greater velocity length relative to their WT counterparts. We observed both higher expression levels and a higher rate of change in transcriptional kinetics. For example, Notch2 and Ezr both exhibited a higher expression level and increased transcriptional rate in the KO samples. These findings are consistent with previous studies demonstrating that loss of Ahr signaling augments features of stemness, i.e., colonic stem cell and nonstem cell progenitor cell self-renewal, clonogenic capacity, and organoid growth (5, 6, 9). Similarly, Ahr KO also inhibits the differentiation of colonic stem cells toward goblet cells and enterocytes (5, 9). It is worth noting that the RNA velocity comparison analysis we adapted helped reveal the changes in transcriptional rate in many key genes, which were undetectable when only a steady expression comparison analysis was carried out.

Here, we further probed the role of Ahr in regulating stem cell proliferation. Ahr KO upregulated Fos and Hspa1a expression, both targets of Foxm1, suggesting an effect of Ahr deletion on Foxm1-regulated genes. This is consistent with our previous findings indicating that Ahr acts as a transcriptional repressor of FoxM1, a master driver of cell-cycle progression (5, 53). Collectively, these findings indicate that modulation of the Ahr-FoxM1 axis, in part, controls colonic stem cell/progenitor cell proliferation. This is noteworthy because the lifetime risk of cancer is highly correlated with the total number of stem cell divisions (54, 55). Additional work is needed to determine whether Ahr-Foxm1 can serve as a potential target for cancer chemoprevention. Interestingly, in complementary systematic analyses assessing cell–cell communication patterns, we also documented for the first time, the ability of Ahr to mediate cross-talk via soluble and membrane-bound factors within the context of the colonic crypt.

With respect to the translational relevance of our findings, previous studies demonstrate the importance of the Ahr and its ligands in colonic stem cell growth and colon carcinogenesis. For example, loss of the Ahr in mouse models enhances development of colon cancer in genetic APC^min mouse models (5–9). In addition, loss of the Ahr in Lgr5⁺ colonic epithelial cells increases colon stem cell growth (5, 9). Ligands such as plant-derived indole-3-carbinol decrease colon cancer development and growth in genetic and carcinogen-induced mouse models (7, 8) and Ahr ligands also decrease Lgr5⁺ colonic stem cell growth (5, 9). Our recent study provides evidence that roasted coffee extracts are Ahr-active and decrease Lgr5⁺ colonic stem cell growth in cells expressing the Ahr but not Ahr KO cells (56). Thus, dietary and possibly microbial derived Ahr ligands play important chemoprotective roles in colon carcinogenesis and the contributions of Ahr regulated Wnt, Foxm1, and other genes/signaling pathways identified in this study are currently being investigated.

In conclusion, our novel single cell RNA-seq data further highlight the molecular function of Ahr in modulating putative stem cell driver genes, cell potency, lineage decisions, and intercellular communication networks in vivo. These findings support the feasibility of using dietary and gut microbial-derived Ahr ligands to modulate the stem cell niche in order to reduce oncogenic signaling and cancer risk.

L.A. Davidson reports grants from Texas A&M University during the conduct of the study. No disclosures were reported by the other authors.

Y. Yang: Data curation, software, formal analysis, investigation, visualization, methodology, writing–original draft, writing–review and editing. D. Osorio: Conceptualization, resources, data curation, software, formal analysis, validation, investigation, visualization, methodology, writing–original draft. L.A. Davidson: Conceptualization, resources, formal analysis, supervision, funding acquisition, investigation, writing–original draft, project administration, writing–review and editing. H. Han: Resources, investigation, methodology. D.A. Mullens: Resources, formal analysis, investigation. A. Jayaraman: Resources, supervision, project administration. S. Safe: Resources, supervision, project administration. I. Ivanov: Resources, formal analysis, supervision, methodology. J.J. Cai: Conceptualization, resources, software, supervision, funding acquisition, investigation, visualization, methodology, writing–original draft, project administration, writing–review and editing. R.S. Chapkin: Conceptualization, resources, data curation, software, formal analysis, supervision, funding acquisition, validation, investigation, visualization, methodology, writing–original draft, project administration, writing–review and editing.

We would like to thank Kerstin Landrock for maintenance of the mouse colonies, Andrew Hillhouse for single-cell transcriptome library generation, Gus A. Wright for cell sorting, and Rachel Wright for graphical art. We also thank Hans Clevers for providing the Lgr5^CreERT2 reporter mice. D.A. Mullens is a recipient of a Hagler Institute for Advanced Study Fellowship.

Funding was provided by Texas AgriLife Research, Hagler Institute for Advanced Study, Allen Endowed Chair in Nutrition & Chronic Disease Prevention (to R.S. Chapkin), Sid Kyle Chair Endowment (to S. Safe), Cancer Prevention Research Institute of Texas (RP160589; to R.S. Chapkin, A. Jayaraman, S. Safe), and the NIH R01-ES025713 (to R.S. Chapkin, A. Jayaraman, S. Safe), R01-CA202697 (R.S. Chapkin, A. Jayaraman, S. Safe), R01-AT010282 (R.S. Chapkin, A. Jayaraman, S. Safe), R35-CA197707 (R.S. Chapkin) and P30-ES029067 (R.S. Chapkin, A. Jayaraman, S. Safe).

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.