Distinct microbial and immune niches of the human colon.

Gastrointestinal microbiota and immune cells interact closely and display regional specificity; however, little is known about how these communities differ with location. Here, we simultaneously assess microbiota and single immune cells across the healthy, adult human colon, with paired characterization of immune cells in the mesenteric lymph nodes, to delineate colonic immune niches at steady state. We describe distinct helper T cell activation and migration profiles along the colon and characterize the transcriptional adaptation trajectory of regulatory T cells between lymphoid tissue and colon. Finally, we show increasing B cell accumulation, clonal expansion and mutational frequency from the cecum to the sigmoid colon and link this to the increasing number of reactive bacterial species.

The colon, as a barrier tissue, represents a unique immune environment where immune cells display tolerance towards a diverse community of microbes - collectively known as the microbiome. The microbiome is critical for many aspects of health and an imbalance of commensals and pathogenic microbes is linked with many disease states [1]. Thus, understanding what constitutes a healthy, homeostatic relationship between host immune cells and the microbiome of the human colon is of critical importance.

The composition of the microbiota at any location in the intestines is determined by the availability of nutrients and oxygen, the transit rate of luminal content and compartmentalized host immune activity, and as such, is spatially distinct [2]. Regional differences are most evident when comparing the small intestine and distal colon in humans or other mammals [2]. Within the colon, increasing bacterial diversity from proximal (including the cecum, ascending and transverse colon) to distal regions (comprising the descending and sigmoid colon connecting to the rectum) has been reported [3].

The intestinal immune system has a symbiotic relationship with the microbiome and is central to the maintenance of epithelial barrier integrity. The lamina propria and associated lymphoid tissues contain one of the largest and most diverse communities of immune cells - including both lymphocytes and myeloid cells [4]. There is marked regional variation in immune cells along the gastrointestinal tract, with T helper (TH) 17 cells decreasing in number from duodenum to colon, and T regulatory (Treg) cell numbers being highest in the colon [5]. Immune cells can respond to environmental cues including the microbiota. Mouse studies have demonstrated that specific bacterial species can fine-tune intestinal immune responses, including TH17 [6,7], Treg [8], or TH1 [9,10] and B cell activation [11]. However, the extent to which there is regional variation in the mucosal microbiome within an individual, and how this might influence local immune cell niches along the colon, has not been investigated to date.

Here, we catalogued the mucosal microbiome in different regions of the human colon, a gastrointestinal organ with the most diverse and dense microbiome content and region-restricted disease states [2]. In parallel, we applied single-cell RNA-seq (scRNA-seq) to make a census of steady-state immune cell populations in the adjacent tissue and in draining mesenteric lymph nodes (mLN), results of which are available at www.gutcellatlas.org. We demonstrate previously unappreciated changes in the proportions and activation status of T and B cells in distinct regions of the healthy human colon from proximal to distal, and relate these differences to the changing microbiota.

Results

Microbiome composition differs along distinct colon regions

To create a map of bacterial composition at the mucosal surface of the colon, we performed 16S ribosomal RNA (16S rRNA) sequencing of swabs from the mucosa surface of the cecum, transverse colon and sigmoid colon of twelve disease-free Caucasian deceased transplant donors (Methods, Figure 1a and Supplementary Table 1). The major gut phyla - Bacteroidetes, Firmicutes, Proteobacteria and Actinobacteria, were present throughout the gut of each donor (Figure 1b and Supplementary Table 2). While diversity of operational taxonomic units (OTUs) was consistent across the colon, and considerable variability existed between donor as previously reported [12] (Supplementary Figure 1a), we did observe changes in the composition of the microbiome. Most notably at the level of phyla, Bacteroidetes was more prevalent in sigmoid colon (Figure 1c and Supplementary Figure 1b). This was mostly attributable to an increase in Bacteroides, which dominates the colonic microbiome of individuals on high protein and fat diets typical in Western countries [13] (Figure 1c and Supplementary Figure 1c). Additionally, Enterococcus was more prevalent in the proximal colon, and Coprobacillus and Escherichia/Shigella were more abundant in the distal colon, although these proportions varied considerably between donors (Figure 1b,c and Supplementary Figure 1c).

Figure 1

Variation in the microbiome from proximal to distal colon.

a) Workflow for 16S ribosomal sequencing of matching mucosal microbiomes and scRNA-seq profiling of immune cells from mesenteric lymph node (mLN), and lamina propria of cecum, transverse colon and sigmoid colon. b) Phylogenetic tree representing diversity and mean abundance of bacterial species in the cecum, transverse colon and sigmoid colon. Mean abundance was calculated as the percentage of operational taxonomic units (OTUs) for each species from total as determined by 16S rRNA sequencing and averaged for twelve donors (black scale). Unassigned OTUs are shown as black branches. Bacteria groups of interest are highlighted. c) Relative abundances of OTUs at genus level of bacteria species in colon regions as in (b).

Past studies characterizing the colonic microbiome typically rely on stool samples, which do not accurately recapitulate the composition of bacteria at the mucosal surface [14]. Our catalog of mucosal bacteria throughout the colon demonstrates heterogeneity in the microbiome at the mucosal surface from proximal to distal colon, and reveals specific genera with preference for colonizing certain colon regions.

Immune cell heterogeneity in steady-state colon

Next, we sought to determine whether the heterogeneity we observed in the colonic microbiome was accompanied by differences in the adjacent host immune cells. To this end, we generated high-quality transcriptional data from over 41,000 single immune cells from the mesenteric lymph nodes (mLN) and lamina propria of cecum, transverse colon and sigmoid colon (Figure 1a). We acquired tissue biopsies from five deceased transplant donors (Methods and Supplementary Table 1). Samples were dissociated to release cells from the mLN and lamina propria of colonic tissue. Immune cells were enriched either by flow-sorting of CD4+ T cells (live CD45+CD3+CD4+) and other immune cells (live CD45+CD3−CD4−), or by CD45+ magnetic bead selection or Ficoll gradient (Methods and Supplementary Figure 2a). Each fraction was then subjected to scRNA-seq (Figure 1a). Despite the enrichment for immune cells, we captured epithelial cells and fibroblasts and these were computationally removed from further analysis (Methods).

Pooled analysis and visualization with Uniform Manifold Approximation and Projection (UMAP) of all four tissues from all donors revealed distinct clusters in the lymphoid and colonic tissues (Figure 2a and Supplementary Figure 2b). We identified 25 cell types and states in the intestinal lamina propria and mLN (Figure 2a, Supplementary Table 3 and Supplementary Table 4), consistent with the immune populations described in recent reports [15]. Among these were follicular and memory B cells, IgA+ and IgG+ plasma cells, effector and memory CD4+ T cells, T regulatory (Treg) cells, CD8+ T cells, γδ T cells, innate lymphoid cells (ILCs), natural killer (NK) cells, mast cells and myeloid cells (Figure 2a). Sub-clustering of myeloid cells showed two distinct populations of conventional dendritic cells: cDC1 expressing XCR1, CADM1, CLEC9A, BATF3 and IDO1, and cDC2 expressing CLEC10A and CD86 (Figure 2a,b). In addition, we identified monocytes expressing CD14 and CD68, macrophages expressing FCGR3A (gene encoding CD16), LYVE1+ macrophages [16] and plasmacytoid DCs (pDCs) expressing IRF4 and SELL. B cells, DCs and γδ T cells had MKI67+ cycling populations suggesting higher rates of proliferation compared to other colonic immune cell populations (Figure 2a,b).

Figure 2

Profiling immune cells along the steady-state colon.

a) UMAP illustration of pooled scRNA-seq data of immune cells of mLN, cecum, transverse colon and sigmoid colon from five donors colored by cell type annotation (left) and tissue of origin (right). b) Heatmap of mean expression of marker genes used to annotate cell types in (a). Point size shows the fraction of cells with nonzero expression. c) Relative percentages of CD4+ T subtypes within all CD4+ T cells for each tissue as in (a). d) Relative percentages of cell types within all non-CD4+ T immune cells for each tissue as in (a), with B lineage cells shown in the left panel and all other cell types in the right panel.

To determine how immune cells differ along the colon from proximal to distal regions, we investigated the relative proportions of cell types within mLN and colon regions (Figure 2c,d). Since CD4+ T cells and all other immune cells were sorted separately during initial tissue processing, our analyses here were also kept separate. As noted from our visual inspection of the UMAP plot (Figure 2a), there were major differences in the cell types present between mLN and the colon. In particular, there were marked differences in the activation and memory status of T and B cells between the colon and lymph nodes, suggesting that these cell types are molded by their environment (Figure 2c,d). In the mLN, CD4+ T cells were typically CXCR5+, ICOS follicular helper cells, and SELL+ (encodes CD62L), CCR7+ central memory cells (Figure 2b,c). In contrast, colonic CD4+ T cells had a more effector phenotype, expressing high levels of the tissue residency marker CD69 [17], falling into the TH17 (CCR6+, IL22+, CCL20) or TH1 (CXCR3, IFNG+) subtypes. There was an inverted gradient in the relative proportion of TH17 and TH1 cells, with the cecum dominated by TH17 cells that reduced in frequency in the transverse colon, and still further in the sigmoid colon, and TH1 cells following the opposite trend, being more abundant in the sigmoid colon (Figure 2c). This distinct distribution of colonic TH1 and TH17 cells is concordant with spatial variation of the microbiome, hinting at a relationship between the two populations.

B cells in the mLN were predominantly CD19, MS4A1 (encodes CD20), CD40+, TNFRSF13B+ (encodes TACI), CD38 memory or CXCR5+, TCL1A+, FCER2+ (encodes CD23) follicular B cells (Figure 2b,d). In contrast, in the three regions of the colon, the main population of B cells were SDC1+ (encodes CD138), CD38+, plasma cells. Plasma cells were enriched in the sigmoid colon relative to both the cecum and transverse colon, whereas the proportion of memory B cells was lower in the sigmoid colon (Figure 2d). This suggests that conditions in the sigmoid colon may favor the generation of plasma cells rather than memory B cells from germinal center responses in this region, or that the tissue contains more plasma cell niches.

T helper cells disseminate through the colon and adopt region-specific transcriptional profiles

We next investigated CD4+ T effector cells across the colon. These were annotated above based on expression of functional markers (IL17A, IL22 and CCL20 versus IFNG and CXCR3) (Supplementary Figure 3a) rather than transcriptional regulators (RORC and TBX21) that are lowly expressed. Correlation analysis of TH1 and TH17 cells between different colon tissues revealed high transcriptional similarity between these effector cell subtypes (Figure 3a), with mLN versus peripheral tissue signature accounting for the greatest amount of variability (Spearman’s corr = 0.88). Within the effector T cells of the colon, transverse colon and cecum cells cluster by T helper subtype. TH1 and TH17 cells of the sigmoid colon did not cluster with their respective effector subtypes from other regions. Differential gene expression analysis between the effector cells in the sigmoid colon versus those in cecum and transverse colon revealed higher expression of activation-related molecules including TANK (TRAF family member associated NF-kB activator) (adjusted P < 10−10), CD83 (adjusted P < 10−10) and PIM3 (adjusted P < 10−8) (Figure 3b). Expression of CCL20, encoding the ligand for CCR6 that is expressed by epithelial and myeloid cells more highly in small intestines than colon [18], was also slightly increased by T helper cells of the proximal colon (Figure 3b). Although this is likely due to higher abundance of CCL20+ TH17 cells at this site (Figure 2c and Supplementary Figure 3a). Conversely, sigmoid colon effector T cells showed higher expression of KLF2 (adjusted P < 10−14) (Figure 3b) that encodes a transcriptional factor that transactivates the promoter for Sphingosine-1-phosphate receptor 1 (S1PR1) and is critical for T cell recirculation through peripheral lymphoid tissue [19], LMNA (adjusted P < 10−47) that encodes a molecule that is reported to promote TH1 differentiation [20] and EEF1G (adjusted P < 10−13) that encodes a driver of protein synthesis.

Figure 3

Dissemination of T helper cells in colon and region-determined transcriptional profiles.

a) Correlation matrix of mean transcriptional profiles of TH1 and TH17 cells from cecum, transverse colon, sigmoid colon and mLN (n=5 donors). b) Mean expression level of differentially expressed genes of pooled TH1 and TH17 between cecum, transverse colon and sigmoid colon as in (a). Point size shows the fraction of cells with nonzero expression. c) UMAP projection of Smartseq2 profiled flow-sorted T cells annotated as TH1 and TH17 cells of the cecum, transverse colon and sigmoid colon (n=1 donor). Colored lines connect cells sharing the same CDR3 sequence. d) Heatmap of numbers of members within clonal families in TH1/TH17 subsets (left) and colon region (right) as in (c).

Next, we looked into the clonal relationships between T helper cells. We performed plate-based Smartseq2 on TCRα/β+ flow-sorted cells from colon regions and mLN of a sixth donor to capture paired gene expression and TCR sequences from individual T cells. Clonal groups were shared between TH1 and TH17 subtypes, supporting the notion that effector fate of CD4+ T cells is determined after their initial activation (Figure 3c,d) [21]. Additionally, clonal expansion was observed by TH1 cells of the sigmoid colon (Figure 3d), in line with greater abundance in this tissue. Likewise, clonal expansion of TH17 cells was greatest in the cecum matching accumulation seen with the droplet-based scRNA-seq analysis (Figure 3d). Several TH1 and TH17 clonal sisters were shared between clonal sites (Figure 3d), evidence that T helper clones disseminate to distant regions of the colon.

Together these data demonstrate region-specific transcriptional differences relating to activation and tissue migration in TH1 and TH17 cells of the proximal and sigmoid colon. Identification of clonal sharing between these colon regions supports the idea that these observed transcriptional differences are due to cell-extrinsic rather than intrinsic factors.

Activation trajectory of colonic CD4+ T regulatory cells

Treg cells are known to play a role in balancing the immune activity of other CD4+ T cell subsets. Documented Treg cell activation by Clostridium spp.[8] and previous descriptions of tissue-specific transcriptional profiles [22] inspired us to interrogate these cells in greater detail. Firstly, we noted that the relative proportion of Treg cells did not change significantly from proximal to distal colon (Figure 2c). We then investigated whether the transcriptional profiles of Treg cells from different compartments indicated distinct activation states or functionalities.

As previously observed in the mouse [22], sub-clustering of Treg cells from the mLN revealed major populations of central Treg cells and effector Treg cells (Figure 4a). Central Treg cells were defined by highest expression of SELL, while effector Treg cells were characterized by genes associated with the TNFRSF-NF-κB pathway (TNFRSF9 and TNF) (Supplementary Figure 4a). In addition, we observed a population, previously termed non-lymphoid tissue-like Treg cells (NLT-like Treg). These have characteristics of non-lymphoid tissue Treg cells, including high expression of FOXP3, PRDM1, RORA, IL2RG, IL2RA and CTLA4 (Figure 4a & Supplementary Figure 4a). A fourth population, termed Treg-like cells, lacked the conventional markers of Treg cells - FOXP3 and IL2RA - but clustered more closely with Treg cells than conventional T cells (Figure 4a and Supplementary Figure 4b). These expressed the highest levels of PDCD1 (gene encoding PD-1) among Treg cell populations and, uniquely for this cell type, the transmembrane protein-encoding gene MS4A6A (Supplementary Figure 4a). This population of Treg-like cells could represent a Treg population that transiently loses FOXP3 expression [23].

Figure 4

Treg activation pathway from lymphoid to peripheral tissue.

a) UMAP visualization of Treg subtypes in mLN (left) and pooled from cecum, transverse colon and sigmoid colon (right) (n=5 donors). b) Relative proportions of Treg subsets within all Treg cells from mLN and colon tissue regions as in (a). Bars show mean proportion across all donors (circles). c) Density of Treg subclusters as in a across ‘pseudospace’ (top) and expression kinetics of genes contributing to pseudospace smoothed into 100 bins (filtered by qval < 0.001 and expression in >15 cells). Top bar shows the most represented tissue within each bin. Various dynamically expressed immune-related molecules are annotated, with key genes colored red.

In the colon, Treg cells clustered into three populations as previously described in mouse (Figure 4a) [22]. KLRB1+ (also known as CD161) Treg cells were characterized by expression of LAG3, IL2RA, CTLA4, KLRB1, ICOS and FOXP3 (Supplementary Figure 4a) suggesting a robust regulatory function, analogous to CD161+ Treg cells described in human colon tissue [24]. Non-lymphoid tissue Treg (NLT Treg) cells express IKZF2, GATA3 and DUSP4 (Supplementary Figure 4a) consistent with the profile of thymic-derived Treg cells [25]. The third population exhibited a profile reminiscent of lymphoid tissue Treg cells with expression of SELL, CCR7, TCR7, CXCR5 and RGS2, and were termed Lymphoid Tissue-like Treg (LT-like Treg) cells (Supplementary Figure 4a). LT-like Treg cells are likely newly arrived in the colon tissue from mLN. A small number of Ki67 Treg cells were also identified in the colon and mLN (Figure 4a and Supplementary Table 3). The proportions of these subsets varied by donor, but were mostly consistent between colon regions (Figure 4b).

The presence of NLT-like Treg cells in the mLN and LT-like Treg cells in the colon, both with profiles suggesting migration, led us to recreate this migration pathway in silico. We ordered all Treg cells along “pseudospace” using Monocle2. This gave rise to a smooth pseudospace trajectory from resting central Treg cells in the mLN to highly regulatory Treg cells in the colon (Figure 4c). As seen in the mouse [22], NLT-like Treg cells and LT-like Treg cells blended in the middle of the trajectory, in accordance with these cells representing transitioning and migratory populations between lymphoid and peripheral tissues. In order to understand which gene signatures drive the migration and tissue adaptation of Treg cells in human tissues, we determined the genes changing along the previously calculated pseudospace (Figure 4c). Genes expressed at the beginning of pseudospace included SELL, CCR7 and CXCR4, permitting entry into lymph nodes (Figure 4c). At the end of pseudospace, the most highly expressed genes were FOXP3, IL2RA, CTLA4, IL10 and LAG3. These genes encoding suppressive molecules were co-expressed with TNF receptor genes (TNFRSF4, TNFRSF18, TNFRSF1B) indicating a reliance on the TNFRSF-NF-κB axis. Chemokine receptor-encoding genes, CXCR3, CXCR6, CCR6 and CCR4, were also expressed by Treg cells in the periphery, matching previous reports of TH1- and TH17-like Treg cells in the colon (Figure 4c)[22].

Together, these results highlight heterogeneity in Treg cell states in mLN and colon, and reveal a possible FOXP3-transiently absent population. We also infer a continuous activation trajectory of these Treg cell states between draining lymph nodes and colon, and highlight genes regulating Treg cell migration between tissues and their adoption of Th-like profiles.

B cells display a proximal-to-distal colon activation gradient

Following from our observations in CD4+ T cells across the colon, we next focused on humoral responses by performing a more in-depth analysis of B cells in different colon regions. We compared transcriptional profiles of plasma cells between different colonic regions. This analysis revealed CCL3 and CCL4 as highly enriched in cecal plasma cells (log fold change of 0.61 and 0.70 and adjusted P <10−28 and <10−9 respectively) (Figure 5a). These chemokine-encoding genes are expressed by B cells in response to BCR activation [26] and result in the migration of CCR5-expressing cells such as T cells and monocytes [27] to the tissue microenvironment. This suggests that BCR cross-linking and signaling may be more prominent in the proximal colon. Cecal plasma B cells were also enriched for CXCR4 (log fold change of 0.44, adjusted P < 10−10) (Figure 5a), which encodes a chemokine receptor highly expressed by germinal center (GC) B cells and important for the movement of plasmablasts to the GC-T zone interface post-GC responses [28].

Figure 5

B cells are more abundant, clonally expanded and mutated in the sigmoid colon.

a) Mean expression of key differentially expressed genes by IgA+ plasma cells in cecum, transverse colon and sigmoid colon (n=5 donors). Point size shows the fraction of cells with nonzero expression. b) Proportion of CD27+ B cells of total B cells from cecum, transverse colon and sigmoid colon determined by flow cytometry (n=4 donors). Bar represents the mean and connected points represents values of each donor. Analysis is a two-tailed paired t test. c) UMAP visualization of B cells for which matched single-cell VDJ libraries were derived using 10x Genomics 5’ scRNA-seq (n=2 donors) colored by cell type annotation (left), tissue (middle) and antibody isotype (right). Bar plot of antibody isotype frequencies per annotated cell type (far right). d) UMAP visualization of somatic hypermutation frequencies of IgH sequences as in (c). e) Quantitation of somatic hypermutation frequencies of IgH sequences from B cell types and gut regions as in (c). f) Estimated clonal abundances per donor for members of expanded B cell clones in B cell types and gut regions as in (c). g) Binary count of co-occurrence of expanded B cell clones identified by single-cell VDJ analysis shared across gut regions as in (c). h) Co-occurrence of expanded B cell clones identified by bulk B cell receptor sequencing across gut regions (n=3 donors). Statistics in (e) and (f) are calculated with two-sided Wilcoxon signed-rank tests. Rows and columns in (g) and (H) are ordered by hierarchical clustering. * P <0.05; ** P <0.01; *** P <0.001; **** P <0.0001

Among the genes more highly expressed by plasma cells in the sigmoid colon was CD27 (Figure 5a; log fold change of 0.24; adjusted P <10−40), encoding a member of the TNF receptor family that is expressed by memory B cells and even more highly by plasma cells [29]. We confirmed differential expression of CD27 protein and additionally observed a proximal-to-distal gradient of increasing CD27 expression by plasma cells in the colon (Figure 5b). Targeted homing of B cells from their site of activation in lymphoid tissues to the colonic lamina propria relies on signaling through CCR10 and its cognate ligand, CCL28, and integrin α4β7 [30]. CCR10 (log fold change of 0.24; adjusted P < 10−16), ITGA4 (log fold change of 0.58; adjusted P < 10−24) and ITGB4 (log fold change 0.57; adjusted P < 10−108) were also more highly expressed by sigmoid colon IgA+ plasma cells (Figure 5a).

To determine whether the B cell clonal repertoire changes across the colon, we took advantage of the paired single-cell VDJ-sequencing data available from two donors for which scRNA-seq libraries were generated using 10x Genomics 5’ chemistry. We confirmed the expression of IgM and IgD isotypes by follicular B cells, IgG1 and IgG2 by IgG+ plasma cells, IgA1 and IgM expression by memory cells in the mLN and predominantly IgA2 expression by plasma cells of the colon (Figure 5c). The mutation frequency of the heavy chain variable region was greatest in the plasma cells followed by memory B cells, indicating more somatic hypermutation by these cell types compared to the naive follicular B cells (Figure 5d). Additionally, while mutational frequency was consistent across colon regions and mLN for memory and follicular B cells, it was significantly increased in IgA+ plasma cell of the sigmoid colon compared to the other colon regions (Figure 5e). IgG+ plasma cells also showed a trend towards increased mutational frequency in the sigmoid colon, however their numbers were limiting (Figure 5e).

We then identified clonally-related B cells to explore clonal expansion dynamics of different cellular populations throughout the gut. Clonal expansion was evident in memory B cells, IgA+ plasma cells and IgG+ plasma cells (Figure 5f). While the relative abundance of clonal groups did not differ across the colon regions of memory B cells and IgG+ plasma cells, again this was greatest for IgA+ plasma cells in the sigmoid colon (Figure 5f). This was supported by bootstrapped VDJ sequence diversity analysis of clonally-related IgA+ plasma cells, which showed that the diversity of BCR sequences was consistent between donors and that there was a trend for decreased diversity (consistent with higher rates of clonal expansion) of IgA+ plasma cells in the sigmoid colon (Supplementary Figure 5a,b). Although some clones were shared between B cell types (i.e. memory and IgA+ plasma), indicating that alternate B cell fates can derive from a single precursor cell, most expanded clones within the gut were of the same cell type (Supplementary Figure 5c). Finally, we found many examples of B cell clones shared between all three colonic regions, and to a lesser extent the mLN, for both donors (Figure 5g), indicating dissemination of B cells throughout the colon as previously reported [31]. Our observation of B cell dissemination throughout the colon was replicated with bulk BCR sequencing from whole tissue (Figure 5h). Increased clonal variability of sigmoid colon B cells was reflected in a greater spread of BCR variable chains expressed compared with cecum and transverse colon (Supplementary Figure 5d), altogether suggesting a more active response in the distal versus proximal colon.

These data indicate a highly activated state of plasma cells in the distal colon compared with proximal colon plasma cells, characterized by greater accumulation, somatic hypermutation, clonal expansion and stronger homing to the colon mucosa.

IgA is responsive to a richer microbial community in the sigmoid colon

Previous reports have shown that IgA is secreted by plasma cells in response to the presence of specific bacteria species, rather than as a general response to the microbiota [32]. In light of this, we examined whether the increased plasma cell activation we observed in the sigmoid colon was linked to the differences in bacteria species. To this end, we assessed IgA-opsonization of bacteria from the donor microbiota samples (Figure 6a). A greater proportion of bacteria from the sigmoid colon was positive for IgA-binding compared with bacteria of the cecum and transverse colon (Figure 6b). Furthermore, shotgun sequencing of IgA-opsonized bacteria revealed a richer community of species in the sigmoid colon (Figure 6c and Supplementary Table 5). Diversity of IgA-bound bacteria, which considers relative abundance of each species, was lower in the sigmoid colon compared to cecum (Figure 6d).

Figure 6

Increasing number of microbiota species recognized by antibodies in the sigmoid colon.

a) Experimental workflow for assessing Ig-opsonized colon bacterial species. b) Representative histogram of IgA1/2-bound Hoechst+ bacteria and summary plot of bound bacteria as a proportion of total bacteria (n=13 donors). Positive binding is set against an isotype control. c) Richness of bacteria species determined as the number of unique species and d) diversity of species identified from shotgun sequencing of Ig-opsonized bacteria from (b) (n=6 donors). P values were calculated using one-tailed paired t tests. * P <0.05

These data suggest that, compared with the cecum, IgA+ plasma cells of the sigmoid colon respond to a rich and unevenly represented community of bacterial species, likely contributing to their increased activation status, strong homing to the colon and trend towards greater clonal diversity (Supplementary Figure 6).

Discussion

In this study we performed the first simultaneous assessment of the colonic mucosal microbiome and immune cells in human donors at steady-state. This enabled us to compare lymphoid and peripheral tissue immunity and explore how immune cells and their neighboring microbiome change along the colon within the same individuals. In doing so, we highlight previously unappreciated regional differences in both cellular communities. Our unique annotated colon immune single-cell dataset is available at www.gutcellatlas.org, where users can visualize their genes of interest.

We describe a shift in the balance of T helper subsets, with a predominance of TH17 in the cecum and TH1 in the sigmoid colon. Decreasing abundance of TH17 cells has similarly been shown from proximal small intestine to colon of mice [5]. Additionally, simultaneous increase of the genus Bacteroides and TH1 numbers in the sigmoid colon are in line with findings that polysaccharide from Bacteroides fragilis preferentially induces TH1 differentiation in the intestine of germ-free mice [9]. An alternative or complementary explanation for skewed T helper proportions and variation in transcriptional profiles is offered in a study by Harbour et al., which demonstrated that TH17 cells can give rise to a IFNy+ ‘TH1-like’ cell in response to IL-23 production by innate cells [33]. These findings demonstrate the complexity of external signals shaping colonic Th responses leading to regional changes in their numbers and differentiation.

In contrast to conventional T cells, Treg cells are evenly represented across the colon. Treg cell subpopulations within the mLN and colon tissue are analogous to those we have recently described in mouse [22]. We also identify an additional population, termed Treg-like cells, reminiscent of a CD25−FOXP3loPD1hi Treg population in the peripheral blood, although the latter was described to also express Ki67 [34]. This population could represent uncommitted Treg cells experiencing transient loss of FOXP3 while retaining regulatory potential [23] or permanent loss of FOXP3 and adoption of a more pro-inflammatory phenotype after repeated stimulation [35]. Our pseudospace analysis of Treg cells suggests a continuum of activation states from resting cells in the lymphoid tissue through to highly suppressive cells in the periphery. We identify genes underlying this transition including chemokine receptors that are strongly expressed on arrival to the intestine and enable interaction with, and suppression of, TH1 and TH17 cells [36]. Similarly to mouse, pseudospace from lymphoid to peripheral tissue correlates with transcriptional markers of migratory potential and suppression of effector cells.

We find that IgA+ plasma cells are more abundant and have greater expression of colon-specific migration markers (CCR10, ITGA4 and ITGB7) in the sigmoid colon compared to the proximal colon. This adds finer resolution to previous reports describing increasing abundance of IgA+ plasma cells from the small intestine to the colon [37]. Our B cell repertoire analysis demonstrates extensive clonal expansion within each colon region and to a lesser extent between regions, arguing for colonic dissemination of B cells from the same precursor pool, followed by local expansion. Furthermore, more clonal sharing existed between regions of the colon than with mLN, consistent with recent work showing that while mLN clones can also be detected in blood, the intestines are host to a unique B cell clonal network [38]. Sigmoid colon plasma cells, in particular, exhibit greater mutational abundance and clonal expansion. This is consistent with a trend towards reduced plasma cell clonal diversity in the sigmoid colon. While this was not statistically significant, this may have been due to the relatively low number of VDJ sequences obtained from these samples, thus decreasing the power of this analysis. Previous work has shown the mutational frequency of B cells is consistently high between duodenum and colon and are primarily driven by dietary and microbiome antigens respectively [39]. Thus, we suggest that enhanced plasma cell accumulation, mutation and expansion in the sigmoid colon is in response to continued stimulation from the local microbiome. This may happen through increased engagement of sigmoid colon plasma cells in T cell-mediated germinal center reactions in gut-associated lymphoid structures or in T cell-independent somatic hypermutation in local isolated lymphoid follicles [40] followed by local expansion. Yet the exact mechanisms require further study.

Finally, we show greater IgA binding to the microbiota in the sigmoid colon compared to proximal sites. One possible explanation for this observation is the accumulation of upstream secreted IgA and IgA-bound bacteria in the sigmoid colon. However, our simultaneous observation of enhanced plasma cell responses in the sigmoid colon suggests that IgA is locally produced as a result of immune poising. Possible scenarios contributing to a richer immune-reactive microbiome in the sigmoid colon are bacteria derived externally via the rectum. Alternatively, environmental pressures (i.e. lower water and nutrient levels [2]) could restrict outgrowth of dominant gastrointestinal species of the proximal colon, providing space for smaller communities of opportunistic species. The IgA response in the colon is antigen-specific rather than a general response to the presence of bacteria [32]. Thus, the overall increased number of unique species recognized by host IgA antibodies in the sigmoid colon is fitting with the enhanced clonal expansion and mutation of plasma cells at this site.

Together, our simultaneous analyses of microbiome and neighboring immune cells highlight the significance of environmental signals in shaping and maintaining regional adaptive immune cell composition and function in the intestine at steady-state. Dysregulation of T helper cells [41] and plasma cells [42], has been implicated in susceptibility to inflammatory bowel disease. Observations of the linked compartmentalization of these immune cells and microbial species along the colon at steady-state may provide a platform for understanding the mechanisms underpinning the tropism of different intestinal diseases to specific regions of the gut, such as Crohn’s disease and ulcerative colitis.

Methods

Colon and mesenteric lymph node tissue retrieval

Human tissue was obtained from deceased transplant organ donors after ethical approval (reference 15/EE/0152, East of England Cambridge South Research Ethics Committee) and informed consent from the donor family. Fresh mucosal tissue from the cecum, transverse colon and sigmoid colon, and lymph nodes from the intestine mesentery, were excised within 60 minutes of circulatory arrest and colon tissue preserved in University of Wisconsin (UW) organ preservation solution (Belzer UW® Cold Storage Solution, Bridge to Life, USA) and mLN stored in saline at 4°C until processing. Tissue dissociation was conducted within 2 hours of tissue retrieval. Four individuals (287c, 296b, 403c and 411c) had received antibiotics in the two weeks prior to death (Supplementary Table 1).

Tissue dissociation for flow-sorting separation and MAC separation

Tissue pieces from donors 290b, 298c, 302c, 364b and 411c were manually diced and transferred into 5 mM EDTA (Thermo Fisher Scientific)/1 mM DTT (Sigma-Aldrich)/10 mM HEPES (Thermo Fisher Scientific)/2% FBS in RPMI and incubated in a shaker (~200 rpm) for 20 minutes at 37°C. Samples were briefly vortexed before the media renewed and incubation repeated. Tissue pieces were washed with 10 mM HEPES in PBS and transferred into 0.42 mg/ml Liberase TL (Roche)/0.125 KU DNase1 (Sigma)/10 mM HEPES in RPMI and incubation for 30 minutes at 37°C. The digested samples were passed through a 40 µm strainer and washed through with FBS/PBS.

Flow-sorting

Cells from donor 290b, 298c and 302c were pelleted and resuspended in 40% Percoll (GE Healthcare). This was underlayed with 80% Percoll and centrifuged at 600g for 20 minutes with minimal acceleration and break. Cells at the interface were collected and washed with PBS. Cells were stained for fluorescent cytometry using Zombie Aqua Fixable Viability Dye (Biolegend, cat. 423101; diluted 1:200), CD45-BV605 (clone HI30, Bioegend, cat. 304043; dilution 1:100), CD3-FITC (clone OKT3; Biolegend, cat. 317305; dilution 1:100), CD4-BV421 (clone SK3; Biolegend, cat. 344631; dilution 1:100), CD8-PE-Cy7 (clone SK1, Biolegend, cat. 344711; dilution 1:100), CD19-APC-Cy7 (clone HIB19, Biolegend, cat. 302217; dilution 1:100), IgD-PE Dazzle (clone IA6-2, Biolegend, cat. 348207; dilution 1:100), CD27-BV711 (clone M-T271, Biolegend, cat. 356429; dilution 1:100), HLA-DR- BV785 (clone L243, Biolegend, cat. 307641; dilution 1:100), CD14-APC (clone 63D3, Biolegend, cat. 367117; dilution 1:100) and CD11c-PE (clone 3.9; eBioscience, cat. 12-0116-42; dilution 1:100). Non-CD4+ T immune cells were sorted as live singlet CD45+, CD3- and CD4−. CD4+ T cells were sorted as live singlet CD45+, CD3+ and CD4+. Each faction was manually counted using 0.4% Trypan Blue (gibco) and a haemocytometer and diluted to 500 cells/µl in PBS. Sorting was carried out on a BD FACS ARIA Fusion. Analysis of flow-sorting data was conducted with FlowJo Software package (version 10.4).

MACS cell enrichment

Cells from donor 411c were pelleted for 5 minutes at 300g and resuspended in 80 µl of ice-cold MACS buffer (0.5% BSA (Sigma-Aldrich Co. Ltd), 2 mM EDTA (ThermoFisher) in DPBS (Gibco)) and 20 µl of CD45 Microbeads (Miltenyi Biotech). Cells were incubated for 15 minutes at 4°C before being washed with 2 ml of MACS buffer and centrifuged as above. Cells were resuspended in 500 µl of MACS buffer and passed through a pre-wetted MS column on QuadroMACS magnetic cell separator (Miltenyi). The column was washed 4 times with 500 µl of MACS buffer, allowing the full volume of each wash to pass through the column before the next wash. The column was removed from the magnet and the cells eluted with force with 1 ml of MACS buffer into a 15 ml tube. Cells were pelleted as above and cell number and viability were determined using a NucleoCounter NC-200 and Via1-Cassette (ChemoMetec). Cells were resuspended at 500 cells/µl in 0.04% BSA in PBS.

Tissue dissociation for Ficoll separation

Tissues from donor 390c were manually diced and <5.0 grams was added per Miltenyi C tube with 5mL tissue dissociation media (Liberase TL (0.13 U/mL; Roche) and DNase (10 U/mL Benzonase nuclease; Merck)) in 1% FCS and 20 mM HEPES in PBS (Lonza). Samples were dissociated with GentleMACS Octo for 30 minutes homogenizing/37°C cycle. Enzymatic digestion was stopped with the addition of 2 mM EDTA in tissue dissociation media. Digested samples were then passed through a 70 µM smart strainer (Miltenyi) before being washed with PBS and pelleted at 500g for 10 minutes. Cells were resuspended in PBS, layered onto FicollPaque Plus (GE Healthcare) and spun at RT 400g for 25 minutes. Mononuclear cells were retrieved from Ficoll layer and washed with PBS. Cells were filtered through 0.2 µM filter (FLowmi cell strainers, BelArt). Cells were manually counted using a hemocytometer and diluted to a concentration of 1000 cells/µl in 0.04% BSA in PBS.

10x Genomics Chromium GEX library preparation and sequencing

Cells were loaded according to the manufacturer’s protocol for Chromium single cell 3’ kit (version 2) or 5’ gene expression (version 2) in order to attain between 2000-5000 cells/well. Library preparation was carried out according to the manufacturer’s protocol. For samples from donors 290b, 298c and 302c, eight 10x Genomics Chromium 3’ libraries were pooled sequenced on eight lanes of an Illumina Hiseq 4000. For samples from donors 390c and 417c, sixteen 10x Genomics Chromium 5’ libraries were pooled and sequenced on 2 lanes of a S2 flowcell of Illumina Novaseq 6000 with 50 bp paired end reads.

10x Genomics Chromium VDJ library preparation and sequencing

10x Genomics VDJ libraries were generated from the 5’ 10x Genomics Chromium cDNA libraries as detailed in the manufacturer's protocol. BCR libraries for each sample were pooled and sequenced on a single lane of Illumina HiSeq 4000 with 150 bp paired end reads.

Plate-based scRNA-seq

Plate-based scRNA-seq was performed with the NEBNext Single Cell/Low Input RNA Library Prep Kit for Illumina (New England Biolabs Inc, E6420L). Cells from donor 364b were snap frozen in 10% DMSO in 90% BSA. Cells were thawed rapidly in a 37°C water bath and diluted slowly with pre-warmed 2% FBS in D-PBS. Cells were pelleted for 5 minutes at 300g and washed with 500 µl of DPBS and pelleted as before. Cells were resuspended in 200 µl of CD25-PE (clone M-A251, Biolegend, cat. 356102; diluted 1:200), CD127-FITC (clone eBioRDR5, eBioscience, cat. 11-1278-42; diluted 1:200), CD4-BV421 (clone SK3, Biolegend, cat. 344632; diluted 1:200) and TCRα/β-APC (clone 1p26, Biolegend, cat. 306718; diluted 1:200) and Zombie Aqua Fixable Viability Dye (diluted 1:400) and incubated for 30 minutes in the dark at room temperature. Cells were washed twice with 500 µl of 2% FBS in D-PBS before being filtered through a 100 µM filter. Single, live, TCRβ+ cells were FACS sorted into a pre-prepared 384-well plate (Eppendorf, Cat. No. 0030128508) containing 2 µl of 1X NEBNext Cell Lysis Buffer. FACS sorting was performed with a BD Influx sorter with the indexing setting enabled. Plates were sealed and spun at 100g for 1 minute then immediately frozen on dry ice and stored at -80°C.

cDNA generation was then performed in an automated manner on the Agilent Bravo NGS workstation (Agilent Technologies). Briefly, 1.6 µl of Single Cell RT Primer Mix (New England Biolabs Inc) was added to each well and annealed on a PCR machine (MJ Research Peltier Thermal Cycler) at 70°C for 5 minutes. 4.4 µl of Reverse Transcription (RT) mix was added to the mixture and further incubated at 42°C for 90 minutes followed by 70°C for 10 minutes to generate cDNA. 22 µl of cDNA amplification mix containing NEBNext Single Cell cDNA PCR MasterMix and PCR primer was mixed with the cDNA, sealed and spun at 100g for 1 minute. cDNA amplification was then performed on a PCR machine (MJ Research Peltier Thermal Cycler) with 98°C 45 s, 20 cycles of [98°C 10 s, 62°C 15 s, 72°C 3 mins], 72°C 5 mins. The 384-well plate containing the amplified cDNA was purified with an AMPure XP workflow (Beckman Coulter, Cat No. A63880) and quantified with the Accuclear Ultra High Sensitivity dsDNA kit (Biotium, Cat. No. 31028). ~10 ng of cDNA was stamped into a fresh 384-well plate (Eppendorf, Cat. No. 0030128508) for sequencing library preparation.

Sequencing libraries were then generated on the Agilent Bravo NGS workstation (Agilent Technologies). Purified cDNA was fragmented by the addition of 0.8 µl of NEBNext Ultra II FS Enzyme Mix and 2.8 µl of NEBNext Ultra II FS Reaction buffer to each well and incubated on a PCR machine (MJ Research Peltier Thermal Cycler) for 72°C at 15 minutes and 65°C for 30 minutes. A ligation mixture was then prepared containing NEBNext Ultra II Ligation Master Mix, NEBNext Ligation Enhancer and 100 µM Illumina compatible adapters (Integrated DNA Technologies) and 13.4 µl added to each well of the 384-well plate. The ligation reaction was incubated on the Agilent workstation at 20°C for 15 minutes and then purified and size selected with an AMPure XP workflow (Beckman Coulter, Cat No. A63880). 20 µl of KAPA HiFi HS Ready Mix (Kapa Biosystems, Cat. No. 07958927001) was then added to a pre-prepared 384-well plate (Eppendorf, Cat. No. 0030128508) containing 100 µM i5 and i7 indexing primer mix (50 µM each) (Integrated DNA Technologies). The indexing primers pairs were unique to allow multiplexing of up to 384 single cells in one sequencing pool. The plate containing the PCR Master Mix and indexing primers was stamped onto the adapter ligated purified cDNA, sealed and spun at 100g for 1 minute. Amplification was performed on a PCR machine (MJ Research Peltier Thermal Cycler) with 95°C for 5 minutes, 8 cycles of [98°C 30 seconds, 65°C 30 seconds, 72°C 1 minute], 72°C 5 minutes. The PCR products were pooled in equal volume on the Microlab STAR automated liquid handler (Hamilton Robotics) and the pool purified and size selected with an AMPure XP workflow (Beckman Coulter, Cat No. A63880). The purified pool was quantified on an Agilent Bioanalyser (Agilent Technologies) and sequenced on one lane of an Illumina HiSeq 4000 instrument.

Single-cell RNA sequencing data alignment

10x Genomics gene expression raw sequencing data was processed using CellRanger software version 3.0.2 and 10X human transcriptome GRCh38-3.0.0 as the reference. 10x Genomics VDJ immunoglobulin heavy and light chain were processed using cellranger vdj (version 3.1.0) and the reference cellranger-vdj-GRCh38-alts-ensembl-3.1.0 with default settings. NEB sequencing data was processed using STAR 2.5.1b into HTSeq (version 0.11.2) and mapped to 10X human transcriptome GRCh38-1.2.0.

Single-cell RNA sequencing quality control

Single cell read counts from all samples were pooled and filtered considering number of UMIs - keeping genes expressed in minimum of 3 cells, keeping cells where genes detected are in a range 700-6000. Non-immune cells were excluded from the final analysis based on the absence of PTPRC and presence of markers such as EPCAM (epithelial cells) and COL1A1 (fibroblasts).

Cell type annotation

Cells were clustered using Scanpy (version 1.4) processing pipeline [43]. In short, the counts were normalized to 10,000 reads per cell (sc.pp.normalise_per_cell) and log transformed (sc.pp.log1p) to be comparable amongst the cells. The number of UMIs and percentage of mitochondrial genes were regressed out (sc.pp.regress_out) and genes were scaled (sc.pp.scale) to unit variance. The normalized counts were used to detect highly variable genes (sc.pp.highly_variable_genes). Batch correction between the donors was performed using bbknn (version 1.3.6) [44] on 50 PCs and trim parameter set to 100. Clusters were then identified using Leiden graph-based clustering (resolution set to 1). Cell identity was assigned using known markers shown in Figure 2b and the top differentially expressed genes identified using Wilcoxon rank-sum test (sc.tl.rank_genes_groups function, Supplementary Table 4). CD4 T cells, myeloid cells, B cells were subclustered for identification of subsets within each cluster. Treg cells annotated above were further subclustered using the “FindClusters” and “RunUMAP” functions from Seurat (version 3.0.1)[45] (Figure 4). The number of PCs used for Treg cell clustering were estimated by the elbow of a PCA screen plot, in combination to manual exploration of the top genes from each PC. Clustering of mLN Treg cells was performed using 1-14 PCs and resolution of 0.3 and colonic Treg cells with 1-11 PCs and resolution of 0.3.

‘Pseudospace’ analysis

‘Pseudospace’ analysis of Treg cells was performed using Monocle (version 2.10.1) [46]. Data was log normalized and cell ordered based on DDRTree reduction on highly variable genes with donor effect regression. The heatmap in Figure 4c was generated using the “plot_pseudotime_heatmap” function in Monocle (version 2.10.1) [46]. Genes contributing to ‘pseudospace’ were first filtered for exclusion of mitochondria and immunoglobulin genes, expression in at least 15 cells and qval < 0.001. Gene expression was smoothed into 100 bins along pseudospace using a natural spline with 3 degrees of freedom. Matching column annotation for each bin was determined as the most prevalent tissue region origin of the cells within that bin. Genes were grouped into 2 clusters and ordered by hierarchical clustering. Cells were ordered through pseudotime.

Sampling of microbiome

Swabs were taken immediate of the mucosal surface of excised tissue using MWE Transwab Cary Blair (catalogue number: MW168). Swabs were maintained at 4°C. Working in a biosafety cabinet, swabs were washed in 500 µl of anaerobic PBS and mixed with anaerobic 50% glycerol before being snap frozen on dry ice and stored at -80°C until use.

Microbiota profiling and sequencing

The microbiota vials were defrosted on ice. Approximately 100 µl of each was transferred into new Eppendorf tubes. DNA was extracted from microbiome samples using the MP Biomedical FastDNA SPIN Kit for soil (catalogue number 116560200). 16S rRNA gene amplicon libraries were made by PCR amplification of variable regions 1 and 2 of the 16S rRNA gene using the Q5 High-Fidelity Polymerase Kit supplied by New England Biolabs as described in [47]. Primers 27F AATGATACGGCGACCACCGAGATCTACAC (first part, Illumina adaptor) TATGGTAATT (second part, forward primer pad) CC (third part, forward primer linker) AGMGTTYGATYMTGGCTCAG (fourth part, forward primer) and 338R CAAGCAGAAGACGGCATACGAGAT (first part, reverse complement of 3′ Illumina adaptor) ACGAGACTGATT (second part, golay barcode) AGTCAGTCAG (third part, reverse primer pad) AA (fourth part, reverse primer linker) GCTGCCTCCCGTAGGAGT (fifth part, reverse primer) were used. Four PCR amplification reactions per sample were carried out; products were pooled and combined in equimolar amounts for sequencing using the Illumina MiSeq platform, generating 150 bp reads. Analysis of partial 16S rRNA sequences was carried out using SILVA v132 and mothur MiSeq SOP v1.42.3 [48]. The 16S rRNA gene alignments were used to determine a maximum likelihood phylogeny using FastTree (version 2.1.10) [49]. Phylogenetic trees were visualized and edited using iTOL (version 5) [50].

T cell clonal sharing analysis

T cell receptor sequences generated using the Smartseq2 scRNA-seq protocol were reconstructed using the TraCeR software as previously described [51].

Bulk B cell receptor sequencing

Small portions of samples were taken from excised tissues and snap frozen in 1ml of RNAlater (Ambion). RNA extracted from tissue using the QIAshreddar and QIAGEN Mini Kit (50). RNA concentration was measured using a Bioanalyzer. B Cell Receptor (BCR) heavy chain sequences of all B lineage subsets present in the tissue were amplified as previously described [52]. Briefly, RNA was reverse transcribed using a barcoded reverse primer set capturing all antibody (sub)classes. Targeted heavy-chain amplification was performed with a multiplex set of IGHV gene primers to FR1 and a universal reverse primer using HiFi qPCR KAPA Biosystems. After adapter filtering and trimming, BCR sequences were assembled and aligned using MiXCR (version 3.0.1) [53]. It is worth noting that detected BCR sequences are biased towards those included in the reference database and while there is a continuous discovery of novel germline alleles, no database is currently a complete reflection of the human IGH locus diversity. Only in-frame and IGH sequences with at least 3 read counts were kept for the analysis. To calculate the CDR3 nucleotide shared repertoire, the tcR package (version 2.2.1) was used [54].

10x Genomics single-cell VDJ data processing, quality control and annotation

Poor quality VDJ contigs that either did not map to immunoglobulin chains or were assigned incomplete by cellranger vdj were discarded. For additional processing, all IgH sequence contigs per donor were combined together. We further filtered IgH contigs as to whether they had sufficient coverage of constant regions to ensure accurate isotype assignment between closely related subclasses using MaskPrimers.py (pRESTO version 0.5.10) [55]. IgH sequences were then further annotated using IgBlast (version 1.12.0) [56] and reassigned isotype classes using AssignGenes.py (pRESTO) prior to correction of ambiguous V gene assignments using TIgGER (version.03.1) [57]. Clonally-related IgH sequences were identified using DefineClones.py (ChangeO version 0.4.5) [57] with a nearest neighbor distance threshold of 0.2, as determined by visual inspection of the output of distToNearest (Shazam version 0.1.11) [57]. CreateGermlines.py (ChangeO version 0.4.5) was then used to infer germline sequences for each clonal family and observedMutations (Shazam version 0.1.11) was used to calculate somatic hypermutation frequencies for each IgH contig. Estimated clonal abundances and IgH diversity analyses within each donor were performed using estimateAbundance, rarefyDiversity and testDiversity of Alakazam (version 0.2.11) [57] with a bootstrap number of 500. Finally, the number of quality filtered and annotated IgH, IgK or IgL chains were determined per unique cell barcode. If more than one contig per chain was identified, metadata for that cell was ascribed as “Multi”. The subsequent metadata table was then integrated with the single-cell RNA-seq gene expression objects for annotation of IgH contigs with B cell types and downstream analysis. Co-occurrence of expanded clone members between tissues and/or cell types was reported as a binary event for each clone that contained a member within two different tissues and/or cell types.

Quantifying IgA binding of bacteria

Microbiome frozen in 50% glycerol were defrosted on ice before being washed in PBS and pelleted for 3 minutes at 8,000 rpm. Bacteria were then stained on ice for 30 minutes with IgG-PE (Biolegend, clone HP6017, cat. 409304; dilution 1:100), IgA1/2 biotin (BD, clone G20-359, cat. 555884; dilution 1:50), followed by 20 minutes with streptavidin-APC (Biolegend, cat. 405207; dilution 1:100). For isotype controls, mouse IgG1 κ Isotype- Biotin (BD, cat. 550615; dilution 1:50) or mouse IgG2a, κ Isotype- PE (Biolegend, cat. 400212; dilution 1:100) were used. The bacteria were washed before and after DNA was stained with Hoechst33342 (Sigma-Aldrich). The stained bacteria were sorted as Hoechst+ and IgG+IgA+ or IgG−IgA+ into PBS using the BD Influx and then stored at -80°C until DNA extraction.

Genomic DNA Extraction for shotgun sequencing of IgA-bound bacteria

Bacteria in PBS was defrosted on ice before being pelleted at 3900 rpm for 5 minutes at 4°C. The supernatant was removed and the pellet resuspended in 2 ml of 25% sucrose in TE Buffer (10mM Tris pH8 and 1mM EDTA pH8). 50 µl of 100 mg/ml lysozyme in 0.25M Tris was added and incubated at 37°C for 1 hour. 100ul of Proteinase K (18 mg/ml),15 µl of RNAase A (20 mg/ml), 400 µl of 0.5 M EDTA (pH 8), and 250 µl of 10% Sarkosyl were then added. This was left on ice for 2 hours and then incubated at 50°C. The DNA was then mixed with 5 ml of Tris-EDTA buffer and purified using four cycles of addition of 5 ml Phenol: Chlorophorm: Isoamyl Alcohol in 15 ml MaXtract High Density Phase Lock Gel (PLG) tubes (QIAGEN), centrifuging at 2800 rpm for 5 min and retreiving the aqueous phase into a new PLG tube. DNA in the final aqueous phase was precipitated with 10 ml of 100% ethanol at -20 °C overnight, centrifuged at 3900 rpm for 20 min at 4 °C and washed with 10 ml of 70% ethanol, before being centrifuged at 4500 rpm for 10 min at 4 °C and gently dried at 50 °C overnight. DNA was then resuspended in 200 µl of Tris-EDTA buffer and isolated by running through a gel. Samples were quantitated by qbit and pooled. Pooled DNA was sequenced on the Illumina Hiseq2500 platform. Inverse Simpson (diversity) and chao (richness) of IgA-opsonized bacteria was determined by using R package “microbiome”.

Statistical analysis

Sample sizes for each experimental group and statistical tests used are included in the relevant figure legends. Measurements were taken from distinct donors in all experiments.