GLUTag cells were a kind gift from Prof Daniel Drucker (Toronto). Standardized cell culture protocols for GLUTag cells used in this study are described in detail in Cyranka et al (Cyranka et al., 2019).

Active GLP-1 secretion from GLUTag cells was measured by fluorescence resonance energy transfer (FRET)-based ELISA (62GLPPEG; Cisbio) as described in Cyranka et al. (2019). To facilitate a high open probability for NMDA receptors, all secretion experiments were done in Mg²⁺-free Krebs buffer (basal buffer) and was prepared as follows: 120 mM NaCl, 5 mM KCl, 2 mM CaCl₂, 24 mM NaHCO₃, 15 mM HEPES, 10 µM Glycine, 0.1% bovine serum albumin (BSA, fatty acid free, A6003, Sigma), pH 7.4 at 37°C, similar to the buffer previously used by Marquard et al. (2016). The following additives were used when indicated in the experiments: 10 µM dizocilpine (MK-801), 100 µM glutamic acid, 300 µM glycine, 5.6 mM glucose, and 10 mM glucose. All experiments were done on ice to minimize GLP-1 degradation. Cell viability in the presence of different buffer conditions and additives were confirmed by Neutral Red cell viability assay (Abcam, ab234039), performed according to the manufacturer’s instructions.

An overnight culture of GLUTag cells (0.6 × 10⁵ cells/well) in 24-well plates coated with 0.4% Matrigel™(Corning) was transfected with siRNA targeting Grin1 (probe ID: 62062; AM16708, Ambion). The siRNA complex was prepared according to manufacturer’s protocol with Lipofectamine RNAiMAX (Life Technologies). The amount of siRNA used per well was 600 nM (15 pmol), and 1.5uL of Lipofectamine. The negative control was non-targeting siRNA sequence Silencer Negative Control No. 1 (AM4611, Life Technologies). After transfection, cells were kept at 37°C, 5% CO₂ for up to 72 h s.

Total protein in whole cell lysates of GLUTag cells was quantified by BCA (Bicinchoninic Acid) Protein Assay (Pierce™ BCA Protein Assay Kit - Thermo Fisher Scientific) according to manufacturer’s protocols, and final protein concentrations were calculated using a bovine serum albumin (BSA) standard curve supplied with the kit. Lysate samples were stored at −80°C until Western blot analysis.

SDS-PAGE gel electrophoresis and Western blot analysis were performed according to standard laboratory protocols. The anti-NMDAR1/Grin1 primary antibody was AB9864, rabbit mAB from Merck-Millipore, visualized by A16035 donkey anti-rabbit IgG:HRP from Life Technologies. Western blots were quantified using Fiji (NIH) and normalised to the loading control for each lane.

For cell culture studies, data are presented as mean ± standard deviation. Data was analysed by repeated measures one-way ANOVA or multiple t-tests (Prism) using a Holm-Šídák’s post-hoc analysis multiple comparison correction. A p value ≤ 0.05 for a 95% confidence interval was regarded as statistically significant. Statistics were performed using GraphPad Prism6 (Graphpad Software, CA, United States). GLP-1 secretion assays are shown as three experimental data points or more, with each experimental data point representing the mean of three technical repeats, unless otherwise indicated in the legends to the figures. Experiments were performed on cells from different frozen stocks and different cell passages to control for interbatch variation. A total of 6–24 data points from the technical repeats were pooled and analysed for normal distribution using Prism, and normality of the data was confirmed.

GLUTag cells were cultured in low glucose (1g/L) DMEM (Cat.no. 11880028, Life Technologies) supplemented with GlutaMAX (Cat.no. 35050061, Life Technologies) and 10% in-house heat inactivated FBS (Cat.no. F7524, Sigma). Cells were seeded at 1 × 10⁶cells/well in 6-well dishes coated with 0.4% Matrigel (BD, Cat.no. 354234) and grown for 48 h s until the cell harvest. mRNA was isolated using a Dynabeads mRNA DIRECT Kit (Cat.no. 61012, Invitrogen™). Extracted mRNA was used to prepare Illumina compatible sequencing libraries with QuantSeq 3′mRNA-Seq Library Prep Kit-FWD (Cat.no. 15, Lexogen) based on the manufacturer’s instructions. The final sequencing libraries were pooled and sequenced in Illumina NextSeq 500 system.

Publicly available RNA-sequencing datasets of human and murine enteroendocrine cells (Roberts et al., 2019) as well as primary ileal organoid cultures (Goldspink et al., 2020) were retrieved from NCBI’s Gene Expression Omnibus (GEO) (Edgar et al., 2002) under the respective accession numbers GSE114853, GSE114913 and GSE148224. In the dataset GSE148224, cells cultured in both IF and IF* media were included in the analysis due to the similarity of their transcriptome (Goldspink et al., 2020). The FACS-sorted cells from each study were renamed and arranged into groups for the purpose of this paper, as detailed in Table 1.

Quality control and adaptor trimming of the sequencing data were performed with FastQC (v. 0.11.9) and Trim Galore! (v. 0.6.4_dev). Salmon (v. 1.5.1) (Patro et al., 2017) was used to align reads to decoy-aware indexes built using transcript sequences and primary assembly genomes from Gencode (human data: GRCh38, release 38; mouse data: GRCm38, release M23). For the QuantSeq 3′ mRNA sequencing data, the --noLengthCorrection argument was used. Alignment statistics were collected with MultiQC (v. 1.11) (Ewels et al., 2016). The average mapping rate of the newly sequenced GLUTag dataset was 66% (standard deviation (s.d.) 4%), leading to an average number of mapped reads of 27.3 million (M) (s.d. 5.8 M). For the publicly available datasets GSE114853, GSE114913 and GSE148224, these statistics were, respectively, 23% (s.d. 8%), 2.5 M (s.d. 1.4 M); 62% (s.d. 3%), 6.2 M (s.d. 4.3 M); 29% (s.d. 7%), 8.5 M (s.d. 2.8 M). Transcript quantification import, transcript-per-million (TPM) transformation (whenever relevant) and gene-level summarizations were performed by the R package Tximeta (v. 1.10) (Love et al., 2020), using the countsFromAbundance = “no” argument.

In order to compare the human and mouse transcriptomes, datasets were merged using one-to-one orthologues pairs using Ensembl data (release 102) (Howe et al., 2021) from BioMart (v. 2.46.3) (Smedley et al., 2009). To minimize the batch effect of sequencing technologies and independent studies, the ComBat_seq() function of the R package sva (v. 3.38.0) (Leek et al., 2012; Zhang et al., 2020) was used, defining the datasets as batches. Subsequently, the corrected count data was normalized by library size and variance-stabilizing transformed using the vst() function of the R package DESeq2 (v. 1.30.1) (Love et al., 2014). Normalized counts were averaged within each group to be compared and the Pearson correlations’ R-squared values calculated between each group (All R values were positive before squaring). The results were then plotted as a heatmap using the default parameters (complete linkage hierarchical clustering using the Euclidean distances) of the pheatmap function of the R package ComplexHeatmap (v. 2.6.2) (Gu et al., 2016). BioMart annotations and the IUPHAR/BPS Guide to PHARMACOLOGY complete “target and family” list (v. 2021.3) were used to restrict the analyses to specific gene types and families (Harding et al., 2022). For within-dataset comparisons (gene expression plots), non-batch corrected transcript-per-million (TPM) values were used (corresponding to count-per-million (CPM) for the QuantSeq 3′ mRNA sequencing data).

Data processing was performed in R (v. 4.0.5) using the tidyverse packages (v. 1.3.1). Plots were generated using combinations of the packages ComplexHeatmap (v. 1.42.3), ggplot2 (v. 3.3.5), ggpubr (v. 0.4.0), lemon (v. 0.4.5), RColorBrewer (v. 1.1–2), cowplot (v. 1.1.1) and ggplotify (v. 0.1.0).

The GLUTag-related sequencing data generated and discussed in this publication has been deposited and is available under the GEO accession number GSE193866. The R code used to generate the results is available in the github repository https://github.com/novonordisk-research/GLUTag_analysis.

The addition of the non-competitive NMDA receptor antagonist MK-801 (10 μM) for 2 h to GLUTag cells cultured in buffer containing glycine with glutamate showed a 37% increase in GLP-1 release relative to control samples (Figure 1A, open diamonds, 1,045 ± 142 pg GLP-1/ml vs. filled diamonds, 1,655 ± 382 pg GLP-1/ml), while buffer containing glycine alone showed a 34% increase in GLP-1 release (Figure Fig1A, open triangles, 812 ± 12 pg GLP-1/ml vs. filled triangles, 1,223 ± 80 pg GLP-1/ml) when compared to controls. The stimulatory effect of MK-801 was not significant in the presence of basal buffer in the absence of Mg²⁺ ions (Figure 1A, filled circles), glutamate (Figure 1A, hexagons) or in the presence of 5.6 mM glucose (Figure 1A, squares). The data are consistent with the effect of NMDA receptor antagonism being most prominent in Mg²⁺-free buffer in the presence of both the NMDA receptor agonists, glutamate and glycine, as these buffer conditions would ensure the highest channel open probability for NMDA receptor. GLUTag cell viability were not affected by the buffer conditions or additives used in this study (Supplementary Figure S1) and GLP-1 secretion under basal buffer conditions was not affected by the presence (460 ± 102 pg GLP-1, n = 6) or absence (505 ± 17 pg GLP-1/ml, n = 6) of Mg²⁺. Unfortunately, these studies could not be extrapolated to ex vivo primary gut crypt cultures as these cultures were not compatible with Mg²⁺-free buffer and established crypts detached rapidly in the absence of Mg²⁺ ions.

In order to control for any possible off-target effects of MK-801—which is known to also target acetylcholine receptors (nAChR), and to unequivocally confirm a direct role of NMDA receptor activity in the regulation of GLP-1 release from GLUTag cells, the essential NMDA receptor subunit, GluN1 (also called NMDAR1 or zeta-1), was knocked down by siRNA targeting of the gene Grin1, (Figure 2; Supplementary Figure S2). Western blot analysis showed a significant decrease in GluN1 protein expression levels in siRNA-treated GLUTag cells (Figure 2A, WB inserts and Figure 2B, 32% ± 10.5%). This decrease in protein expression levels coincided with up to 31% increased GLP-1 release from siRNA-treated GLUTag cells (Figure 2, squares, 746 ± 79 pg GLP-1/ml, 1,385 ± 444 pg GLP-1/ml, and 1812 ± 166 pg GLP-1/ml for Mg²⁺-free buffer, 1 mM glucose, and 10 mM glucose, respectively) when compared to control siRNA probe (Figure 2, triangles, 476 ± 13 pg GLP-1/ml, 1,049 ± 257 pg GLP-1/ml, and 1,002 ± 300 pg GLP-1/ml for Mg²⁺-free buffer, 1 mM glucose, and 10 mM glucose, respectively). This effect was most prominent in the presence of 10 mM glucose (1812 ± 166 pg GLP-1/ml Grin1 siRNA vs. 1,002 ± 300 pg GLP-1/ml siRNA control).

As ex vivo gut crypt cultures were not compatible with buffer conditions that favor a high open probability for NMDA receptors, NMDA receptor expression in GLUTag cells and in both murine and human L-cells were investigated to ascertain if GLUTag cells would be a good model system to investigate NMDA receptor activity. To this effect, RNA-sequencing of GLUTag cells was performed and their transcriptome compared to murine and human cells. Publicly available datasets of FACS-sorted in-vitro primary ileal organoid cells (Goldspink et al., 2020) as well as FACS-sorted in vivo gastrointestinal cells from both mouse and human gut tissue (Roberts et al., 2019) were downloaded and processed. The different sample types from each study are detailed in Table 1, alongside the related abbreviations used in this paper.

Using Pearson correlations and clustering of the R ² values, the datasets were used to determine which enteric cell transcriptomes were the most strongly correlated with the transcriptome of GLUTag cells. Merging of the datasets, using variance-stabilized normalized data, led to a clustering by species and datasets (data not shown). This is partially due to the different sequencing technologies and protocols used for each independent dataset. To correct for the batch effects due to the different methods utilized (using GEO datasets as batch), the software Combat-Seq (Zhang et al., 2020) was used to observe transcriptome correlations between different tissues/cell types rather than between studies/technologies (as in (Gilad and Mizrahi-Man, 2015)). A principal component analysis showed that the main variance in the resulting data was driven by differences in cell types (Supplementary Figure S3). Moreover, genes contributing the most to the first principal component are relevant to the cells’ functions, such as their secretory phenotype (e.g., Gcg, Gip, Pyy). After this correction, all sample groups were clustered using genes from three subsets of the transcriptomes: protein coding genes (n = 15779), genes coding for ligands and target proteins (using the IUPHAR/BPS Guide to PHARMACOLOGY, “an expert-curated resource of pharmacological targets and the substances that act on them” (Harding et al., 2022)) (n = 2,621) or genes coding for G-protein-coupled receptors (GPCRs) and ion channels (n = 584). In all three cases, the transcriptome of GLUTag cells was observed to be most similar to entero-endocrine cells (regardless of species/tissues) compared to samples containing non-EEC cells (Figure 3A), although the GLUTag cells seem to preferentially cluster with human ileal cells rather than the mouse duodenal cells when restricted to GPCRs and ion channel genes (Figure 3A; Supplementary Figure S4).

Next, the expression data of the NMDA receptor subunit genes were specifically analysed. In GLUTag cells, the two main NMDA receptor subunits expressed are Grin3a and Grin1 (Figure 3B). After rescaling the expression of these genes within L-cells and L-cell models (GLUTag and GLUVenus POS organoid cells), GLUTag cells clustered with human L-cells due to the similarly high expression of Grin1 relative to other NMDA receptor subunit genes.