Authors: Leonard Nettey (1Program in Immunology, Harvard Medical School, Boston, MA 02115, USA; 2Institute for Medical Engineering & Science, Massachusetts Institute of Technology, Cambridge, MA 02142, USA; 3Ragon Institute of MGH, MIT, and Harvard, Cambridge, MA 02139, USA; 4Broad Institute of MIT and Harvard, Cambridge, MA 02142, USA; 5Department of Chemistry, Massachusetts Institute of Technology, Cambridge, MA 02139, USA), Hengqi Betty Zheng (6Division of Gastroenterology and Hepatology, Seattle Children’s Hospital and University of Washington, Seattle, WA 98105, USA), Joseph C. Devlin (7Regeneron Pharmaceuticals Inc., Tarrytown, NY 10591, USA), Julie E. Horowitz (7Regeneron Pharmaceuticals Inc., Tarrytown, NY 10591, USA), Sara Rosenbaum (8Division of Gastroenterology, Hepatology, and Nutrition, Boston Children’s Hospital, Boston, MA 02115, USA), Nathalie Fiaschi (7Regeneron Pharmaceuticals Inc., Tarrytown, NY 10591, USA), Wei Keat Lim (7Regeneron Pharmaceuticals Inc., Tarrytown, NY 10591, USA), Kyle Kimler (2Institute for Medical Engineering & Science, Massachusetts Institute of Technology, Cambridge, MA 02142, USA; 3Ragon Institute of MGH, MIT, and Harvard, Cambridge, MA 02139, USA; 4Broad Institute of MIT and Harvard, Cambridge, MA 02142, USA; 5Department of Chemistry, Massachusetts Institute of Technology, Cambridge, MA 02139, USA; 8Division of Gastroenterology, Hepatology, and Nutrition, Boston Children’s Hospital, Boston, MA 02115, USA; 9Division of Pediatric Hematology/Oncology, Boston Children’s Hospital, Boston, MA 02115, USA), Christina Adler (7Regeneron Pharmaceuticals Inc., Tarrytown, NY 10591, USA), Annabelle Lee (7Regeneron Pharmaceuticals Inc., Tarrytown, NY 10591, USA), Min Ni (7Regeneron Pharmaceuticals Inc., Tarrytown, NY 10591, USA), Yi Wei (7Regeneron Pharmaceuticals Inc., Tarrytown, NY 10591, USA), Peter J. Ehmann (7Regeneron Pharmaceuticals Inc., Tarrytown, NY 10591, USA), Shane E. McCarthy (7Regeneron Pharmaceuticals Inc., Tarrytown, NY 10591, USA), Manuel A.R. Ferreira (7Regeneron Pharmaceuticals Inc., Tarrytown, NY 10591, USA), William Galbavy (7Regeneron Pharmaceuticals Inc., Tarrytown, NY 10591, USA), Eli Stahl (7Regeneron Pharmaceuticals Inc., Tarrytown, NY 10591, USA), Sokol Haxhinasto (7Regeneron Pharmaceuticals Inc., Tarrytown, NY 10591, USA), Sandra Coetzee (7Regeneron Pharmaceuticals Inc., Tarrytown, NY 10591, USA), Yoko Yabe (7Regeneron Pharmaceuticals Inc., Tarrytown, NY 10591, USA), Michael Dobosz (7Regeneron Pharmaceuticals Inc., Tarrytown, NY 10591, USA), Paula Keskula (9Division of Pediatric Hematology/Oncology, Boston Children’s Hospital, Boston, MA 02115, USA), Jillian Zavistaski (9Division of Pediatric Hematology/Oncology, Boston Children’s Hospital, Boston, MA 02115, USA), Alexandre Albanese (9Division of Pediatric Hematology/Oncology, Boston Children’s Hospital, Boston, MA 02115, USA), Zoë Steier (2Institute for Medical Engineering & Science, Massachusetts Institute of Technology, Cambridge, MA 02142, USA; 3Ragon Institute of MGH, MIT, and Harvard, Cambridge, MA 02139, USA; 4Broad Institute of MIT and Harvard, Cambridge, MA 02142, USA; 5Department of Chemistry, Massachusetts Institute of Technology, Cambridge, MA 02139, USA), Joshua de Sousa Casal (1Program in Immunology, Harvard Medical School, Boston, MA 02115, USA; 4Broad Institute of MIT and Harvard, Cambridge, MA 02142, USA; 8Division of Gastroenterology, Hepatology, and Nutrition, Boston Children’s Hospital, Boston, MA 02115, USA), Veronika Niederlova (3Ragon Institute of MGH, MIT, and Harvard, Cambridge, MA 02139, USA; 4Broad Institute of MIT and Harvard, Cambridge, MA 02142, USA; 8Division of Gastroenterology, Hepatology, and Nutrition, Boston Children’s Hospital, Boston, MA 02115, USA; 10Laboratory of Adaptive Immunity, Institute of Molecular Genetics of the Czech Academy of Sciences, Prague, Czech Republic), Lusine Ambartsumyan (6Division of Gastroenterology and Hepatology, Seattle Children’s Hospital and University of Washington, Seattle, WA 98105, USA), Ghassan Wahbeh (6Division of Gastroenterology and Hepatology, Seattle Children’s Hospital and University of Washington, Seattle, WA 98105, USA), David L. Suskind (6Division of Gastroenterology and Hepatology, Seattle Children’s Hospital and University of Washington, Seattle, WA 98105, USA), Vanessa Mitsialis (8Division of Gastroenterology, Hepatology, and Nutrition, Boston Children’s Hospital, Boston, MA 02115, USA; 11Brigham and Women’s Hospital Department of Medicine, Division of Gastroenterology, Hepatology, and Endoscopy, Boston, MA 02115, USA), Sara Hamon (7Regeneron Pharmaceuticals Inc., Tarrytown, NY 10591, USA), Stephen M. Carpenter (12Division of Infectious Diseases and HIV Medicine, Department of Medicine, University Hospitals, Case Western Reserve University School of Medicine, Cleveland, OH 44106, USA), Jennifer D. Hamilton (7Regeneron Pharmaceuticals Inc., Tarrytown, NY 10591, USA), Alex K. Shalek (1Program in Immunology, Harvard Medical School, Boston, MA 02115, USA; 2Institute for Medical Engineering & Science, Massachusetts Institute of Technology, Cambridge, MA 02142, USA; 3Ragon Institute of MGH, MIT, and Harvard, Cambridge, MA 02139, USA; 4Broad Institute of MIT and Harvard, Cambridge, MA 02142, USA; 5Department of Chemistry, Massachusetts Institute of Technology, Cambridge, MA 02139, USA; 13Harvard Medical School, Boston, MA 02115, USA; 14Koch Institute for Integrative Cancer Research, Massachusetts Institute of Technology, Cambridge, MA 02139, USA), Scott B. Snapper (8Division of Gastroenterology, Hepatology, and Nutrition, Boston Children’s Hospital, Boston, MA 02115, USA; 11Brigham and Women’s Hospital Department of Medicine, Division of Gastroenterology, Hepatology, and Endoscopy, Boston, MA 02115, USA; 13Harvard Medical School, Boston, MA 02115, USA), Matthew A. Sleeman (7Regeneron Pharmaceuticals Inc., Tarrytown, NY 10591, USA), George D. Kalliolias (7Regeneron Pharmaceuticals Inc., Tarrytown, NY 10591, USA), Andrea T. Hooper (7Regeneron Pharmaceuticals Inc., Tarrytown, NY 10591, USA), Jose Ordovas-Montanes (1Program in Immunology, Harvard Medical School, Boston, MA 02115, USA; 3Ragon Institute of MGH, MIT, and Harvard, Cambridge, MA 02139, USA; 4Broad Institute of MIT and Harvard, Cambridge, MA 02142, USA; 8Division of Gastroenterology, Hepatology, and Nutrition, Boston Children’s Hospital, Boston, MA 02115, USA; 15Harvard Stem Cell Institute, Cambridge, MA 02138, USA), Leslie S. Kean (9Division of Pediatric Hematology/Oncology, Boston Children’s Hospital, Boston, MA 02115, USA; 16Dana Farber Cancer Institute, Boston, MA 02215, USA)
Categories: Article
Source: medRxiv
Authors: Leonard Nettey, Hengqi Betty Zheng, Joseph C. Devlin, Julie E. Horowitz, Sara Rosenbaum, Nathalie Fiaschi, Wei Keat Lim, Kyle Kimler, Christina Adler, Annabelle Lee, Min Ni, Yi Wei, Peter J. Ehmann, Shane E. McCarthy, Manuel A.R. Ferreira, William Galbavy, Eli Stahl, Sokol Haxhinasto, Sandra Coetzee, Yoko Yabe, Michael Dobosz, Paula Keskula, Jillian Zavistaski, Alexandre Albanese, Zoë Steier, Joshua de Sousa Casal, Veronika Niederlova, Lusine Ambartsumyan, Ghassan Wahbeh, David L. Suskind, Vanessa Mitsialis, Sara Hamon, Stephen M. Carpenter, Jennifer D. Hamilton, Alex K. Shalek, Scott B. Snapper, Matthew A. Sleeman, George D. Kalliolias, Andrea T. Hooper, Jose Ordovas-Montanes, Leslie S. Kean
Pediatric Inflammatory Bowel Disease (IBD) remains challenging to treat and difficult to prognosticate. Although multiple immune cell types coordinate pathology in both Crohn’s disease (CD) and ulcerative colitis (UC), specifying which cell types and cell states portend better or worse response to major IBD treatment strategies, including anti-TNF therapies (the only FDA-approved therapy for pediatric IBD), remains challenging. Here, we present the results of the PREDICT study, which enrolled 79 treatment-naïve pediatric patients at the time of diagnostic endoscopy, enabling a comprehensive transcriptomic, histologic, and serologic analysis of 40 CD patients and 16 UC patients, as well as 23 patients with functional gastrointestinal disorders (FGID) who served as pediatric non-inflamed controls. Leveraging these data, we performed a comprehensive analysis of colonic immunology in each of these clinical conditions. Our results indicate that, within the complex landscape of immune pathology in pediatric CD and UC, there is a coordinated shift across the TH1-to-TH17 immune activation continuum among T cells that is pertinent to anti-TNF response. For CD, this shift defines partial response to anti-TNF treatment. For UC, the landscape is more complex, with both TH17 and TFH biology defining disease, and pre-treatment TH17 biology contributing to anti-TNF treatment resistance. Related to this, polyreactive TCR phenotypes within UC TFH cells are correlated with both germinal center activity and the frequency of IgG1 plasma cells, yet opposed to the TH17 signatures associated with anti-TNF nonresponse. The elucidation of these distinct mechanisms of T cell-dependent disease pathology, treatment response, and TCR polyreactivity suggests a model in which sustained TH17 signaling in CD and baseline TH17 signaling in UC are associated with disease pathogenesis, forming a basis for a generalized understanding of IBD pathogenesis and underscoring the need for endotype-specific approaches to IBD therapy.
Inflammatory bowel disease (IBD), for which Crohn’s disease (CD) and ulcerative colitis (UC) are the most common subtypes, encompasses a group of debilitating inflammatory gastrointestinal conditions with increasing incidence^1^. Although pediatric and adult disease have many clinical and biologic commonalities, pediatric IBD can be particularly devastating, given that this disease is often treatment refractory, more predisposed to rapid progression^2^, and requires effective treatment for many decades. Recent efforts have catalogued immunologic associations with treatment resistance in biologic-, but not treatment-naïve, adults; however, no study has yet defined the comparative landscape of these diseases at diagnosis (naïve to all treatments), nor provided a multiomic map of these diseases in pediatric patients^3^. Understanding the cellular and molecular immunopathologic networks causing pediatric disease, elucidating the mechanisms of treatment failure in these young patients, and identifying novel avenues for therapeutic intervention^4^ remain major unmet clinical needs.
Acknowledging the central role that T cells play in both CD and UC, multiple therapies targeting T cell trafficking, IL-12/IL-23, and lymphocyte activation are in use. However, anti-TNF monoclonal antibodies remain first-line biologic treatment for both CD and UC^5,6^, and the only FDA-approved biologic therapies for pediatric IBD. The success of anti-TNF agents is highly variable, and there remains no diagnostic paradigm for prioritizing other agents due to limited understanding of TNF-resistant cell states in CD or UC. Such functional knowledge is urgently required, given that differing mechanisms of IBD pathology may be better treated by biologic therapies targeted towards disease endotypes.
To address these knowledge gaps, we created a prospective natural history study, the Precision Diagnostics in Inflammatory Bowel Disease, Cellular Therapy, and Transplantation study (“PREDICT”, NCT03369353), which enrolled and evaluated a cohort of treatment-naïve pediatric patients (≥ 6 years old) with GI distress requiring endoscopy (n = 79), and who were diagnosed with one of four non-inflammatory ‘Functional GI Disease’ (FGID), ileal-only CD (iCD), colon-involving CD (cCD), or UC. Key to the design of this study was the enrollment of patients prior to diagnostic endoscopy and any treatment, such that the distinctive immunobiology of these conditions could be determined, unaffected by concomitant immunosuppressive treatments. PREDICT combines detailed clinical evaluation and follow-up with multiomic biologic evaluation including 5’-scRNA-seq/scTCR-seq, multiplexed immunohistochemistry (IHC), and human proteome array-based autoantibody detection. This study enabled us to uncover distinct cellular disease mechanisms that drive pathology and treatment resistance within pediatric CD and UC, providing an evidence-based foundation for future evaluations of molecular disease-targeted therapies in these patients.
PREDICT IBD study participants were enrolled from November 2017 to January 2022, during which time 23 FGID, 7 iCD, 33 cCD (1 with perianal-only disease, and 2 without successfully processed diagnostic scRNA-seq biopsies), and 16 UC patients were enrolled (Fig. 1a). In addition to the diagnostic biopsies, 4 iCD, 10 cCD, and 2 UC patients contributed follow-up biopsies that were processed for downstream analysis. At diagnosis, enrolled patients underwent clinical and laboratory evaluation, and both upper and lower endoscopy to evaluate disease and obtain biopsies. Following endoscopy, standard clinical and histologic features were used by the treating physicians to classify participants as FGID, iCD, cCD, or UC, and all treatment decisions were made by the treating physicans^7^. After their initial diagnosis, patients with IBD were followed clinically for up to 2 years. Patients with FGID were followed clinically only as needed. Of the FGID, iCD, cCD, and UC patients analyzed, the median age at diagnosis was 15.7, 17.6, 11.4, and 12.8, respectively (p = 0.013, ANOVA, Supplemental Table 1). Males represented 56.5, 71.4, 54.5, and 56.3% of patients with FGID, iCD, cCD, and UC, respectively (p = 0.88, chi-squared test, Supplemental Table 1).
The majority (70%) of total IBD patients were treated with anti-TNF, with both infliximab and adalimumab used, based on treating physician decision-making (Fig. 1a, Supplemental Table 1). Among patients with successfully processed diagnostic biopsies, 5/7 iCD, 25/31 cCD, and 8/16 UC patients received anti-TNF. For these patients, we defined anti-TNF response status at 2 years post-diagnostic endoscopy. Anti-TNF treatment response categorization [Full Response (FR), Partial Response (PR), Nonresponse (NR)] was adjudicated by two pediatric gastroenterologists (HBZ and SR) based on the following FR was defined as clinical symptom control and biochemical response (normalization of C-reactive protein (CRP) and erythrocyte sedimentation rate (ESR)) on maintenance anti-TNF therapy with dose adjustments only due to insufficient blood drug levels. PR was defined as improvement in clinical symptoms or CRP/ESR without full clinical symptom control or CRP/ESR normalization, with documented escalation of anti-TNF therapy or additional therapeutic agents added. NR was defined as lack of clinical and biochemical response to anti-TNF, using the criteria defined above for FR and PR. In PREDICT, 4/5 iCD, 15/25 cCD, and 0/8 UC patients achieved FR. 1/5 iCD, 10/25 cCD, and 4/8 UC patients achieved PR. 0/5 iCD, 0/25 cCD, and 4/8 UC patients were designated as NR (Fig. 1a, Supplemental Table 1). In addition, for each disease there were patients for whom the treating physician chose not to treat with anti-TNF agents (2/7 iCD, 6/31 cCD, and 8/16 UC). These patients are designated as “not on anti-TNF” (“NOA”), and were treated with a variety of other immunomodulators and antimetabolic agents, based on their treating physician’s choice.
We performed cellular profiling on fresh (non-frozen) diagnostic and follow-up colon biopsies through linked 5’ scRNA-seq and scTCR-seq (Fig. 1a). During sample processing, biopsies were separated into lamina propria and epithelial fractions, and samples from each fraction were separately prepared for scRNA-seq^8,9^. The use of 5’ chemistry enabled us to simultaneously interrogate the gene expression profiles and TCR identities of the profiled T cells using scTCR-seq. Following quality control and cellular annotation, 650,501 high-quality gene expression transcriptomes from diagnostic and follow-up biopsies were further analyzed (Supplemental Table 1).
To annotate cells, we first identified major “cell types” through low-resolution clustering (see Methods)^10^. These cell types were analyzed agnostic to GI layer of origin and included colonic epithelium; stromal and endothelial cells; myeloid cells; plasma cells; B cells; and T cells, NK cells & innate lymphoid cells (Fig. 1b–g, Supplemental Fig. 1a–g, 2a–f, Supplemental Table 2). Each cell type was separately analyzed to identify additional levels of detail, which we termed “cell subtypes” and “cell states” (Supplemental Fig. 1h). We manually annotated a total of 82 cell states within 31 cell subtypes across the 6 cell types (Supplemental Table 3)^8,11,12^. Each final annotated cell state was represented by multiple samples, regardless of the number of cells comprising the cell state (Supplemental Fig. 1a, Supplemental Fig. 2a–f). In addition, IHC was performed on biopsies from 44 PREDICT participants (8 FGID, 6 iCD, 15 UC, 15 cCD). Two immunophenotyping panels were used, identifying B cells (CD20), T cells (CD3), myeloid cells (CD68), Treg cells (FOXP3), CD8 T cells (CD8), plasma cells (BCMA), eosinophils (ECP), neutrophils (MPO), dendritic cells and monocytes (CD11c), NK cells (NCR1), and epithelial cells (panCK) (Fig. 1h, i).
To statistically characterize major differences in the cell state frequencies between FGID controls and IBD at diagnostic endoscopy, we used Milo, a tool that generates differential abundance statistics for scRNA-seq data by assigning cells to partially overlapping neighborhoods on a k-nearest neighbor graph^13^. Milo neighborhoods were created within each cell type to determine which cell states contributed most to the variation observed amongst FGID, cCD, and UC colon biopsies. When comparing IBD samples with colonic inflammation (both cCD and UC) to FGID controls, we found enrichment in multiple activated cell types, including inflamed epithelial cells, inflammatory fibroblasts, inflammatory macrophages, mononuclear phagocytes (MNPs), plasmacytoid dendritic cells (pDCs), neutrophils, and multiple activated and inflammatory T cell populations (Fig. 2a, Supplemental Table 4). These findings demonstrate a diverse mucosal inflammatory response in both cCD and UC. IHC analysis corroborated the overarching scRNA-seq findings, identifying increased densities of CD11c+ cells, MPO+ cells, total CD3+ cells, as well as CD3+CD8-FOXP3+ cells (Tregs) in IBD samples (Fig. 2b, Supplemental Table 5).
While pediatric cCD and UC can be challenging to distinguish endoscopically and histologically^14^, Milo differential abundance identified multiple cell states that could clearly distinguish these two IBD entities (Fig. 2c). Compared to UC, cCD biopsies were enriched in epithelial cell states suggestive of retained epithelial function, including non-inflamed absorptive colonocytes and goblet cells (Fig. 2c). In contrast, most epithelial cell states were depleted in diagnostic samples from UC patients, and the cells that remained in UC were skewed towards LGR5+ epithelial stem cells, CCL24-expressing goblet cells, and an epithelial population with expression of antimicrobial and TH17-promoting genes (PI3 and SAA1, respectively), previously demonstrated to be induced by proinflammatory cytokines, including TNF, IL-1β, and IL-6 (Fig. 2c, Supplemental Fig. 2a)^15–19^. These cell states are consistent with characteristic UC histopathology findings that reveal regenerating stem cells in the setting of inflamed epithelium lacking normally functioning absorptive and secretory cells^20^.
Consistent with the epithelial damage in UC, the fibroblasts that support differentiating epithelium (annotated as ABCA8+ and VSTM2A+ fibroblasts) were generally absent from UC biopsies, but present in cCD biopsies (Fig. 2c). We also identified inflammatory fibroblast populations that were nearly exclusive to either cCD (“TFPI2Hi fibroblasts”) or UC (“PDPNHi fibroblasts”) (Fig. 2c). TFPI2Hi fibroblasts exhibited greater expression of chemokines CCL13 and CXCL10, both of which are induced by IFNγ signaling (Supplemental Fig. 2b)^21–23^. In contrast, PDPNHi fibroblasts were enriched for genes involved in extracellular matrix organization and vascular remodeling, as have been described in adult UC (Supplemental Fig. 2b)^12^. In line with increased vascular remodeling, UC biopsies demonstrated an enrichment in the venule cell state and pericytes involved in capillary constriction and vascular remodeling (Fig. 2c)^11,24^.
While enrichments in inflammatory myeloid populations were often shared between both cCD and UC biopsies compared to FGID (with notable elevations in MNPs, inflammatory macrophages, and pDCs, Fig. 2a), UC biopsies were more skewed towards these monocyte, macrophage, and dendritic cell inflammatory states compared to cCD (Fig. 2c). IHC analysis also revealed a marked increase in MPO+ cells in UC, indicating neutrophilia that was not found in cCD samples (Fig. 2d). Plasma cells were significantly shifted towards IgG1 and IgG3 in UC (Fig. 2c), as has been previously reported^25–27^. We also observed that, concomitant with the increase in IgG, there was a coordinated decrease in both IgM and IgA plasma cells among UC biopsies, with the latter dominated by a reduction in IgA2, rather than all IgA cells (Fig. 2c). Potentially connected to the shifts in the immunoglobulin balance, UC could be distinguished from cCD by the presence of proliferating B cells and germinal center-like B cells, which were both increased in UC versus cCD (Fig. 2c). UC samples were also enriched for PD1Hi TFH cells, which are likely to be germinal center TFH cells^28,29^ (Fig. 2c). The combination of increases in these cell state frequencies suggested elevated germinal center activity in colonic mucosa from UC patients compared to both FGID and cCD, which likely contributed to the observed plasma cell shifts. These scRNA-seq findings were corroborated by IHC analysis which demonstrated increased CD20+ (B cells) and BCMA+ (plasma cells) cell density in UC relative to cCD (Fig. 2d).
T cells play a major role in IBD, coordinating interactions with immune, stromal, and parenchymal cells that contribute to inflammation, with known instructive roles for TH1 and TH17 cells in CD and UC^3,8,9,30–32^. We found that both TH1 and TH17 cell states were enriched in cCD and UC relative to FGID (Fig. 2a). Moreover, testing for differences between cCD and UC highlighted a notable shift towards TH1 cells in cCD and towards TH17 cells in UC (Fig. 2c). Milo-based differential abundance analysis was corroborated with per-biopsy comparisons of FGID, iCD, cCD, and UC cell state frequencies (Fig. 2e–n, Supplemental Table 5). Given the well-established connections between inflammatory T cell signaling and T cell proliferation^33^, we explicitly interrogated proliferating T cell states in CD versus UC. Because proliferating T cell states did not comprise enough cellular neighborhoods for Milo analysis (see Methods), this analysis was limited to per-biopsy comparisons. Our findings were consistent with the observation described above, of more PD1Hi TFH in UC versus CD (Fig. 2c), with a trend towards more proliferating TFH cells also observed in UC (Fig. 2o). Because proliferating T cells often express transcriptomic hallmarks of both TH1 and TH17 cells^34,35^, we defined a proliferating T cell state with features of both (including CXCR3, IFNG, IL17A, and GNLY, termed “Prlf. TH17.1”). Interrogation of Prlf. TH17.1 cells revealed that they were enriched in both cCD and UC compared to FGID, and of the two disease states, most enriched in patients with UC (Fig. 2p, Supplemental Fig. 2f).
Having identified increased TH1 frequency in cCD and increased TH17 frequency in UC, we investigated whether detailed analysis of these T cell populations could further elucidate distinct signatures in each IBD subtype. Given the plasticity of helper T cells and their divergent association with cCD and UC, we performed differential gene expression comparing UC and cCD within the “T helper” cell subtype, which includes TH1, TH17, and Prlf. TH17.1 cells. This analysis identified genes enriched in UC versus cCD related to T cell activation, co-stimulation, or exhaustion, including TOX, TNFRSF4 (OX40), TNFRSF18 (GITR), LAG3, and CD247 (CD3ζ) (Supplemental Fig. 3a, Supplemental Table 6). We further delineated these DE genes within the TH1 and TH17 cell states. We found that these genes, as well as HAVCR2 (TIM-3), were increased in both UC TH17 and UC TH1 cells, confirming that their enrichment in the TH cell subtype was not solely due to differences in cell state abundance (Supplemental Fig. 3b, Supplemental Table 7).
Next, we investigated the ambiguity of the Prlf. TH17.1 cell state by analyzing the expression of TH1- and TH17-identifying genes within the Prlf. TH17.1 cell state, as well as within TH1, and TH17 cell states. Even amongst cells annotated as TH1, we found increased expression of TH1-associated genes (IFNG, TNF, ANXA1, CXCR3, GZMK) in cCD compared to UC (Supplemental Fig. 3c). Importantly, TNF, ANXA1, CXCR3, and GZMK were also increased in cCD Prlf. TH17.1 cells, suggesting that TH1 cells may be the primary cell state of origin in cCD Prlf. TH17.1 cells (Supplemental Fig. 3c). Conversely, genes IL17A, IL17F, IL26, and GNLY were increased in UC Prlf. TH17.1 cells relative to cCD, suggesting an association with the TH17 cell state (Supplemental Fig. 3d), and distinguishing the origin of these cells in the two IBD endotypes. Notably, we observed increased IL17A expression in UC TH1 cells. In the context of increased T cell activation and exhaustion in UC, these data suggest that UC TH1 cells may actually be ex-TH17 cells, known to take on a TH1-like phenotype under prolonged inflammatory conditions^36^.
We took advantage of the 5’ scRNA-seq chemistry to simultaneously track linked TCR sequences, to determine the extent to which TH1 or TH17 clones shared TCRs with proliferating TH17.1 cells in cCD or UC (Supplemental Fig. 3e–g, Supplemental Table 8). We found that Prlf. TH17.1 clones came from both TH1 and TH17 cell states in cCD, with a trend towards increased representation of TH1 clones (Supplemental Fig. 3h). In contrast, nearly all clonal overlap in UC was between Prlf. TH17.1 cells and TH17 cells, reinforcing the likelihood that actively proliferating TH cells in UC are primarily of TH17 origin. We also noted CD8A expression in the CD4+ TH17 cell state (Supplemental Fig. 3i). Given our previous investigation demonstrating expansion of CD8+ TH17 cells in adult UC^8^, we sought to determine their activity in pediatric UC. We found that dual CD4+CD8A+ TH17 cells were significantly enriched for clonal overlap with proliferating TH17.1 cells, suggesting active clonal expansion of this cell substate in UC (Supplemental Fig. 3j).
To explore the extent to which other immune pathways interacted with TH1, TH17, and PD1Hi TFH populations, we performed a cell state abundance correlation analysis between these three CD4 subpopulations and each of the other annotated cell states (Fig. 3a, Supplemental Table 9). Our analysis revealed several cell states across all cell types that were significantly associated with each of these CD4 subpopulations. TH1 cells were associated with NK cells, inflammatory macrophages, and TFPI2-high fibroblasts, consistent with a type 1 immune response^37,38^. Amongst TH17 cells, significant correlations were observed for a variety of cell states, with the strongest correlations between MNPs and venules, as well as between TH17 cells and the PI3+ inflamed epithelium or PDPN+ remodeling fibroblasts. Importantly, while both TH17 and PD1Hi TFH cells were enriched in UC patients, their correlation networks were largely distinct, with PD1Hi TFH being primarily associated with IgG1 plasma cells, GC-like and proliferating B cells, and plasmablasts, consistent with their role in T cell-mediated B cell immune activation.
We calculated the per-biopsy geometric mean of the cell states involved in each correlation network to determine the extent to which each network aligned with each disease state. We found that cCD was distinct from FGID, UC, and iCD in its enrichment for TH1-associated cell states (Fig. 3b). Because IFNγ is a predominant TH1 cytokine, we next determined whether IFNγ may be associated with the cell states in the TH1 network. We first examined STAT1 expression in all cell states, given that STAT1 both mediates and is induced by IFNγ signaling^39,40^. We found that 4 of the 6 top STAT1-expressing cell states (Inflm. Mac, TH1, MNP, TFPI2Hi Fibro.) were present in the TH1 correlation network (Supplemental Fig. 4a). We then compared IFNG expression of IFNG-expressing cell subtypes (including TH, NK, CTL, and InvariantT; Supplemental Fig. 4b–c) to STAT1 expression in TH1-associated cell states within the same biopsy, finding significant correlations between TH
IFNG expression and STAT1 expression in inflammatory macrophages, TH1 cells, Tregs, and XCL+ NK cells (Fig. 3c, Supplemental Table 10). The specificity of this correlation was underscored by the fact that IFNG from other IFNG-expressing cell subtypes was not associated with STAT1 in any TH1-associated cell states. (Supplemental Fig. 4d). These correlations were reproduced when we restricted the dataset from all enrolled patients to focus only on cCD and UC biopsies, suggesting that the TH IFNγ association with STAT1 expression is not solely due to the overall status of tissue inflammation (Supplemental Fig. 4e).
Although cCD was significantly enriched for TH17-associated states compared to both FGID and iCD, this network was most enriched in UC patients (Fig. 3d). To further interrogate TH17 signaling, we next focused on IL-17A, the key TH17 cytokine. IL-17A is known to have multiple effects, including enhanced recruitment of peripheral leukocytes and increased antimicrobial defense^36,41–43^. We therefore created gene set scores for IL-17-induced antimicrobial function, chemokine production, cytokine production, and tissue remodeling using genes downstream of IL-17 signaling in the KEGG IL-17 signaling pathway (Supplemental Fig. 5a)^44–46^. These gene sets were active across multiple cell types (Supplemental Fig. 5b, Supplemental Table 11). Given IL-17 effects on barrier function and the presence of TH17 cells within the epithelial fraction of our biopsies (Fig. 3e, Supplemental Fig. 1a), we focused our analyses on IL-17 effects in colonic epithelium. The antimicrobial and chemokine gene sets were active in colonic epithelium (Fig. 3f). We found strong associations between IL17A in TH cells and activity of IL-17 gene sets in colonic epithelium, suggesting either a direct influence of IL17A on this epithelium, or a common process contributing to expression of these genes in distinct cell types (Fig. 3g). Strong correlations remained when we restricted analysis to cCD and UC biopsies only, but the correlations were not observed between IL17A and active IL-17 gene sets in other cell types (Supplemental Fig. 5c). The antimicrobial and chemokine gene sets, and associated antimicrobial and chemokine genes, were active in similar colonic epithelial cells, were significantly enriched in cCD compared to FGID, and were further enriched in UC (Fig. 3h, Supplemental Fig. 5d–g). Because many of the chemokines in the IL-17-induced chemokine gene set are known neutrophil chemoattractants, we contextualized the scRNA-seq results with multiplexed IHC and found profound neutrophilia in UC epithelium (Supplemental Fig. 5i). A similar increase in potentially inflammatory conventional CD4 T cells (CD3+CD8-FOXP3-), was not observed histologically in UC epithelium relative to cCD (Supplemental Fig. 5j).
In addition to their enrichment in TH17 networks, UC patients also demonstrated significant enrichment in the PD1Hi TFH network (Fig. 3i). This group of differentially abundant cell states was highly suggestive of increased T cell-dependent B cell activation and tertiary lymphoid structure activity. The frequencies of GC-like B cells and proliferating B cells were highly correlated (R = 0.8) (Fig. 3a, Supplemental Table 9). We considered that these cell states may be part of tertiary lymphoid structures (TLS). We thus combined both cell states into a “TLS-like B Cell” annotation and observed this grouping of cells to be increased in UC biopsies relative to both FGID and cCD (Fig. 3j). We then examined expression of the germinal center-organizing cytokine CXCL13. Its expression was increased in UC TFH cells and was greatly enriched in UC follicular dendritic cells (FDCs) (Fig. 3k). Notably, CXCL13 expression among FDCs correlated with TLS-like B cell frequency, but only in UC (Fig 3l). Histologically, we identified lymphoid aggregates of B and T cells, which contained a center of CD20+ cells surrounded by a ring of CD3+ cells (morphologically resembling a germinal center), in biopsies across all disease groups (Fig. 3m–n). Importantly, amongst biopsies with any identified TLS, we observed an increase in the number of TLS per biopsy in UC, as well as increased CD20+ cell density within these structures in UC relative to cCD (Fig. 3o–p).
Because anti-TNF agents are the first line (and only FDA-approved) biologic therapy in pediatric IBD, we next investigated whether we could determine differential expression signatures associated with response to these agents from data in diagnostic biopsies (obtained prior to any treatment) from IBD patients. We first performed per-cell type differential expression (DE) analysis between cCD patients who would be treated with anti-TNF and those not treated with anti-TNF in the 2-year follow-up period (the ‘NOA’ patients). This analysis identified a paucity of DE differences, with the only DE gene identified being TNFAIP6, for which the enrichment was largely driven by a single outlier sample (Supplemental Fig. 6a–b, Supplemental Table 6). We next performed a similar DE analysis to compare cCD FR to PR patients at the time of diagnostic biopsy. This analysis identified only the plasma cell gene PHGDH (Supplemental Fig. 6c–d, Supplemental Table 6) differentially expressed between the two groups. This enzyme, critical for serine synthesis, was enriched in cCD PR vs. FR^47^. These results underscore the challenge in creating a predictive signature at diagnosis for anti-TNF response in cCD^48,49^.
Unlike the lack of differences in diagnostic cCD biopsies, an immunologic signature did emerge when post-treatment cCD biopsies were compared between FR and PR patients. To determine the cell types most likely to reflect an immunologic impact from anti-TNF therapies, we first scored each cell type for pre-treatment TNF pathway activity (Fig. 4a). Both myeloid cells and TNKILCs demonstrated elevated TNF pathway scores and were thus chosen for more detailed analyses. We found that TNF pathway scores were decreased for both myeloid cells and TNKILCs in post-treatment cCD relative to pre-treatment for both response groups, consistent with a biologic impact of anti-TNF on all treated patients (regardless of anti-TNF response) (Fig. 4b–c)^50,51^. To more deeply interrogate the molecular impact of anti-TNF on myeloid and TNKILCs, we conducted a DE analysis in these cell types, comparing pre- and post-treatment cCD, without regard to treatment response (Fig. 4d). Several immune genes were found to be upregulated in pre-treatment cCD myeloid cells and TNKILCs relative to post-treatment, including STAT1, cytolytic genes (GNLY, PRF1), calprotectin (S100A8, S100A9), and IFNγ-induced genes (CXCL9, CXCL10, GBP1, GBP2, GBP5)^52^. Moreover, “Interferon gamma signaling” was the top Reactome term enriched in pre-treatment cCD differential expression, consistent with IFNγ-induced inflammatory states suggested in the TH1 network (Fig. 3c, 4e, Supplemental Table 6). Based on these findings, we compared pre- and post-treatment cell state frequencies of TH1-associated cell states among cCD biopsies, finding decreases in each cell state (Fig. 3a, 4f). Additionally, the two other top STAT1-expressing cell states, MNPs and TH17 cells, were also decreased in post-treatment cCD, when all treatment response groups were analyzed together (Fig. 4f). Overall, these findings suggest successful abrogation of IFNγ signaling through anti-TNF therapy in cCD, regardless of clinical treatment response.
With overall differences in pre- and post-treatment cCD delineated, we investigated whether post-treatment biopsies would illuminate differences in cCD FR and cCD PR. This analysis revealed decreases in the TH1-associated inflammatory macrophage cell state in both FR and PR (Fig. 4g), and that TH1 cells themselves were significantly decreased in FR with a trend towards reduction in PR (Fig. 4h). In contrast, the TH17-associated MNP cell state (shown in Fig.3a), as well as TH17 cells themselves, were decreased only in cCD FR, while being maintained in PR (Fig. 4i–j). To more fully evaluate canonical TH1 and TH17 cytokine gene expression, we analyzed the expression of key TH1 cytokine genes (IFNG, TNF) and TH17 cytokine genes (IL17A, IL17F, IL21, IL22, IL26) within the TH cell subtype. Consistent with the findings above, a gene signature created from these cytokines was significantly reduced in the TH cells from cCD FR patients post-treatment relative to pre-treatment, but significantly increased in cCD PR patients post-treatment TH cells relative to pre-treatment (Fig. 4k)^50^. Single gene analysis of pre- and post-treatment TH cells revealed trends toward reduction in TNF, IFNG, IL22, and IL26 in cCD FR (Fig. 4l). In contrast, cCD PR TH cells demonstrated unchanged TNF and IFNG, along with trends towards increases in IL17A, IL17F, IL21, and IL22 after anti-TNF treatment (Fig. 4l). Together, these data point to persistent TH17 immunopathology associated with anti-TNF PR in cCD.
In the PREDICT cohort, the number of UC patients with repeat biopsies (2/16) was too few to rigorously define the post-treatment impact of anti-TNF therapies on immune cells and interactions. However, unlike for CD, UC diagnostic biopsies were discriminative of anti-TNF response. To discover this, we first performed an analysis comparing diagnostic biopsies from UC patients later treated with anti-TNF to NOA patients, utilizing the per-cell type DE analysis paradigm described above. Unlike in cCD patients, where informative gene expression changes were not observed, in UC, several notable genes were transcriptomically distinct in patients for whom the treating physician initiated anti-TNF therapy versus the NOA patients. Among colonic epithelial cells in pre-treatment biopsies, genes enriched in UC NOA patients were notable for CA1, CA2, and SLC26A3, representing normal colonic absorptive function (Supplemental Fig. 7a, Supplemental Table 12). Genes enriched in the colonic epithelium from UC patients who would eventually be treated with anti-TNF therapies included cell cycle-associated genes, such as RRM2, ZWINT, and CENPM (Supplemental Fig. 7a, Supplemental Table 12). Reactome pathways enriched for the DE results included GTPase cell signaling pathways for UC NOA patients, and DNA synthesis and replication pathways for UC anti-TNF patients (Supplemental Fig. 7b–c). Additional notable DE genes in UC patients treated with anti-TNF therapies included the neutrophil chemoattractant CXCL6 in stromal cells, and the inflammatory cytokine-induced ferritin transport genes FTL and FTH1 in myeloid cells^53,54^ (Supplemental Fig. 7a, Supplemental Table 6).
We next evaluated diagnostic biopsies from UC patients who were eventually treated with anti-TNF for per-cell type DE that could differentiate PR versus NR. Of note, PR was the best response achieved in the PREDICT UC cohort (although other studies have documented UC patients achieving FR to anti-TNF)^55–57^. In patients that partially responded to anti-TNF versus NR patients, pre-treatment UC colonic epithelium and TNKILCs were enriched for genes in cellular metabolism and maintenance pathways (Fig. 5a–b). MIF, which promotes TNF induction and whose release is induced by microbial proteins, was enriched in UC PR versus NR myeloid cells, B cells, and TNKILCs (Fig. 5a)^58–60^. Additional notable genes enriched in UC PR versus NR included LYZ (an antimicrobial enzyme) in myeloid cells and interferon-stimulated genes IFI6, IFITM1, and ISG15 in plasma cells (Fig. 5a)^61^.
In contrast, UC NR versus PR colonic epithelium demonstrated enrichment in pathways governing extracellular degradation and remodeling (Fig. 5a, c). Notably, FCER2 (encoding CD23) and JAK3 (induced upon CD40 ligation or bacterial stimulation, necessary for lymphocyte proliferation and differentiation, and for which there is an FDA-approved therapeutic in adult UC) were enriched in UC NR B cells^62–65^. Many of the genes enriched in UC NR TNKILC are involved in differentiation, activation, and migration of CD4 T cells, including RORA (TH17 cytokine production), CCR6 (TH17 and Treg migration), FOXP3 (Treg differentiation), PRDM1 (encoding BLIMP-1, which mediates chronic activation/exhaustion), IL12RB1 (subunit of IL-12 and IL-23 receptors), and JAK3 (Fig. 5a)^66–71^. We also examined per-biopsy expression of several TNKILC DE genes in the TNKILC cell type, which recapitulated the DE findings described above (Fig. 5d). Of note, in addition to therapeutics directed at JAK3, IL12RB1 is involved in pathways targeted by anti-IL-12p40/anti-IL23p19 therapies, underscoring the potential availability of alternatives to anti-TNF in UC patients bearing anti-TNF nonresponse markers^64,72^.
We further evaluated the role of TH17 cytokine signaling in defining UC treatment response. We found that, compared to pre-treatment UC PR, UC NR TH cells were enriched in TH17 cytokines, with all UC NR biopsies having greater TH
IL17A and IL21 than UC PR biopsies (Fig. 5e). Thus, IL17A was identified as a marker of anti-TNF therapy failure in both cCD and UC.
In addition to TH17 cells, adult UC has previously been associated with TPH cells^26^. Thus, we explored whether TPH cells were present in the pediatric IBD patients profiled in PREDICT. We found that neither the PD1Lo TFH nor the PD1Hi TFH population clearly expressed TPH markers (Supplemental Fig. 8a–b, Supplemental Table 7)^73^. Additionally, no clear subpopulation of cells expressed a combination of canonical TPH markers^73^, assessed through UMAP visualization and through gene expression across a pseudotime diffusion vector spanning CD4 T cells (Supplemental Fig. 8c–e, Supplemental Table 7). This may represent a phenotypic difference in adult and pediatric UC.
While no evidence for an association between TPH and UC patients was observed, we did identify a subset of TFH cells that were altered in PREDICT’s UC We found that although total TFH cells were present in similar relative frequencies in all PREDICT disease groups, UC could be distinguished by the fact that TFH cells in these patients were more shifted towards the PD1Hi TFH cell state (Fig. 6a–b). Moreover, amongst TFH genes involved in B cell help^73^, we found CXCL13 to be increased in UC PD1Hi TFH relative to cCD (Fig. 6c, Supplemental Fig. 8a). During cell state annotation, we also observed that B cell help genes were present in TH17 cells (Supplemental Fig. 2f), with MAF, PDCD1, IL21, and CXCL13 being increased in UC TH17 cells compared to cCD TH17 cells (Fig. 6d).
Given these data, we considered a model wherein UC TH17 cells, but not cCD TH17 cells, exhibit TFH-like plasticity. To test this model, we assessed which set of paired CDR3αβ sequences were detectable between non-proliferating TFH cells and other TH cells within a biopsy (Fig. 6e–f). Our analysis demonstrated that, in 6/14 UC diagnostic biopsies, there was evidence for TFH-TH clonal overlap, while only 3/30 cCD diagnostic biopsies demonstrated TFH-TH clonal overlap (Fig. 6g–h, Supplemental Table 8). Clonal sharing was primarily between TH17 and PD1Hi TFH, with only one clonal pair involving TH1 or PD1Lo TFH cells (Supplemental Table 8). These data indicate that, although a small absolute number of shared clonotypes was identified, UC biopsies were enriched for identifiable TFH-TH clonal overlap (Fig. 6g–h).
To determine whether the coordination of these phenotypic shifts would be distinct in cCD versus UC, we compared the relative frequencies of TH17 cells among all TH cells, and PD1Hi TFH cells among all TFH cells. Notably, we found opposite patterns in the two IBD endotypes, with TH17 and PD1Hi TFH cells being positively correlated in cCD, but inversely correlated in UC (Fig. 6i, Supplemental Table 10). Considering the implications of differences in T cell phenotypes on B cell help and the development of antibody-secreting cells, we explored whether PD1Hi TFH cells or TH17 cells were associated with the shift towards IgG1 plasma cells observed in UC (Fig. 2i). We found that in both cCD and UC, PD1Hi TFH cell proportion and CXCL13 expression were positively correlated with IgG1 plasma cells (Fig. 6j–k). However, in TH17 cells, we found a strong inverse correlation of TH17 cell state frequency or IL17A expression and IgG1 plasma cells in UC, despite both TH17 and IgG1 plasma cells both being increased in UC compared to CD (Fig. 2i, 2n, 6j–k). These data suggest an antagonism between TH17-associated and TFH-associated signaling in UC.
Because of the prevalence of T-B interactions observed in UC, and known clinical associations with autoreactivity among IBD patients^74^, we hypothesized that UC TCRs may bear features that distinguish them from cCD TCRs, and that TCR phenotypes could play a role in the plasticity observed among UC T cells.
Previous studies have suggested that increased CDR3β hydrophobicity may predispose TCRs to polyreactivity or autoreactivity through non-specific interactions^75–77^. Pointing towards poly- and/or autoreactivity in UC T cells, we observed a relative increase in usage of large, aliphatic, and hydrophobic amino acid residues, with a corresponding decrease in polar amino acids, in the polymorphic middle region of CD4 T cell CDR3β sequences in patients with UC (Fig. 7a–b, Supplemental Table 8). To relate these changes to functional differences in TCRs, we applied the previously-defined middle-TCR intrinsic regulatory potential (mTiRP) score^78^. mTiRP is a TCRβ-based metric created from the association of CD4 TCRs with the Treg phenotype, given that Tregs derived from the thymus are trained to recognize self and commensal antigens, and are thus more likely to be self-reactive^78,79^.
Our analysis demonstrated that UC biopsies exhibited increased mTiRP scores within CD4 T cells overall, and within TH cells, Treg cells, and TFH cells (Fig. 7c–d). Moreover, mTiRP scores were further increased in UC PD1Hi TFH cells relative to PD1Lo TFH cells, and in clonally expanded UC PD1Hi TFH cells relative to singlets, suggesting that more active TFH populations in UC bear increased features of polyreactivity (Fig. 7e–f). PD1Hi TFH mTiRP scores were significantly associated with TLS-like B cell (R^2^ = 0.53, p = 0.011) and IgG1 plasma cell (R^2^ = 0.65, p = 0.003) frequencies in UC, but not in FGID or cCD (Fig. 7g–h), potentially suggesting this TCR phenotype helps drive germinal center reactions in UC.
Considering the overlapping and opposing TH17 and PD1Hi TFH phenotypes observed in UC, we examined whether PD1Hi TFH mTiRP would have any relationship to TH17 gene expression. Among canonical TH17 cytokines, we identified a strong inverse association between UC PD1Hi TFH mTiRP score and UC TH17 IL17A expression (R^2^ = 0.87, p < 0.001) (Fig. 7i). We note that UC NR pre-treatment (and post-treatment) biopsies bore the lowest average PD1Hi TFH mTiRP score among UC participants, further suggesting orthogonal mechanisms of pathology (Fig. 7j).
Taken together, these data demonstrate that UC TFH cells display a distinct phenotype from TFH cells in both inflammatory (cCD) and non-inflammatory (FGID) study participants. These phenotypes are associated with polyreactivity and correlate strongly with UC-specific increases in the development of humoral immune activation. However, they are inversely correlated with the TH17 features associated with anti-TNF nonresponse, suggesting non-overlapping mechanisms of immune pathology in UC.
In previous studies, in addition to germinal center activation, adult patients with UC have demonstrated enhanced autoantibody production along with other signs of humoral immune pathogenesis^26,27^. Because of this, we performed a proteome-wide screen of plasma from PREDICT participants to determine the burden of autoantibodies on a per-patient basis. To accomplish this, we used the HuProt microarray (CDI Labs)^80–83^, which assesses immunoglobulins for binding to full-length recombinant human proteins produced in yeast, with coverage against approximately 80% of the human proteome^80–83^. For comparison, we included samples (N=5) from patients diagnosed with autoimmune polyendocrinopathy-candidiasis-ectodermal dystrophy (APECED), a rare condition caused by mutations in AIRE that is characterized by a broad range of autoantibodies, which can be considered positive controls^84^. As shown in Supplemental Fig. 9a and Supplemental Table 13, compared to APECED controls, the most striking finding in PREDICT patients was the high intra-group variability that was evident when comparing the median number of autoantigen hits per patient in each patients were present in all PREDICT cohorts with low and high median numbers of autoantigen hits. Overall, the median (interquartile range [IQR]) number of IgG hits was 185.5 (135.75–284.75) for FGID, 244 (124.5–361.5) for iCD, 301.5 (167–391) for cCD, and 215 (127–322) for UC, with each cohort containing patients demonstrating similar numbers of hits compared to the APECED controls (379 [364–437]). Similar variability was observed among IgA autoantigens (Supplemental Table 13). While many factors may contribute to the high patient-to-patient variability observed in PREDICT, some of this variability was likely explained by differences in antigen presentation driven by specific HLA:TCR interactions. An example of this, wherein 2 UC patients expressing the HLA-DRB1*15:01 allele developed a distinct pattern of autoantigens, is described below.
Given the association of TCR characteristics and UC phenotypes, together with the pronounced inter-patient variability in autoantibody burden we discovered with HuProt analysis, we investigated whether HLA antigen presentation may influence the TFH/TH17 dichotomy, as well as the presence of patient-specific auto-antibody patterns, that we observed in UC. To accomplish this, we analyzed genomic data from PREDICT UC participants to map their HLA alleles, an analysis which detected several previously-identified HLA-DRB1 alleles associated with UC risk (Supplemental Fig. 9b, Supplemental Tables 1, 12)^85^. HLA-DRB115:01, classically associated with multiple sclerosis, but also identified as a risk allele in UC, was present in two UC participants (participant IDs 001–036 and 001–038, Supplemental Fig. 9b)^85,86^. Evaluation of TCR phenotypes revealed that biopsies from these two individuals had the highest PD1Hi TFH mTiRP scores (Supplemental Fig. 9c). We also found substantially decreased frequency of TH17 cells and decreased expression of TH17 IL17A in these patients, relative to UC participants without the HLA-DRB115:01 allele (Supplemental Fig. 9d–e), consistent with the dichotomous relationship between TH17- and TFH biology we previously documented in UC (Fig. 6i–k, 7i). As predicted, PD1Hi TFH frequency and CXCL13 expression in both TH17 cells and PD1Hi TFH cells were increased in these individuals (Supplemental Fig. 9d, f). IgG1 plasma cells were also increased in these individuals, while IgA2 plasma cells were decreased (Supplemental Fig. 9g). Interestingly, the “PI3+ Epi.” cell state, an inflamed epithelial cell state likely reflecting microbial interaction, was also increased in the individuals with the HLA-DRB1*15:01 allele (Supplemental Fig. 9h).
Further analysis of the HuProt microarray data demonstrated that these participants (IDs 001–036 and 001–038) were jointly reactive to multiple antigens, some of which were uniquely present in this pair. These include autoantibodies directed against ENO1, an autoantigen previously identified in multiple autoimmune disorders^87–91^, as well as two novel targets, YWHAQ^92^, MED18^93^ (Supplemental Fig. 9i–j). Together, these data are consistent with a model wherein the combination of HLA type, antigen presentation, and T cell biology contribute to the unique autoantibody profiles observed in UC patients and may further explain the intra-group autoantibody heterogeneity observed across PREDICT.
The combined scRNA-Seq, scTCR-Seq, IHC, and autoantibody analysis have led to a concerted model of pediatric IBD that encompasses three axes of immunopathology, with prominent features associated with TH1, TH17, and TFH cells (Supplemental Fig. 10a, Supplemental Table 14) differentially represented across IBD endotypes. While cCD was primarily enriched in TH1 features with a moderate increase in TH17 features compared to FGID, UC was strongly enriched in TH17 and TFH features compared to both patient groups (Supplemental Fig. 10b). Although we observed no baseline differences in cCD patients who were eventually treated with anti-TNF that would predict their anti-TNF response status, analysis of post-treatment biopsies identified persistent TH17 immunopathology as a signature of less-responsive cCD patients (Supplemental Fig. 10c–d). TH17 cells were further implicated in anti-TNF response when we examined pre-treatment UC biopsies, which demonstrated an increase in TH17 features among NR compared to PR patients (there were no FR to anti-TNF in the PREDICT UC cohort) (Supplemental Fig. 10e). Paradoxically, although both TFH and TH17 features were enriched in UC, we found orthogonal skewing of these patients with greater enrichment for TFH features had decreased association with TH17 features (Supplemental Fig. 10f)^87–91^. Taken together, these data provide evidence for multiple disease endotypes within classically defined IBD that can be distinguished by their guiding immunopathologic signatures.
Here we present a multiomic examination of pediatric IBD from patients enrolled in the PREDICT study, identifying the distinguishing cellular and molecular features of cCD and UC, and their response to treatment with anti-TNF therapies. Our in-depth profiling provides a comprehensive cellular and molecular landscape of parenchymal, mesenchymal, and immune cells in the pediatric IBD colon prior to biologic or non-biologic treatment, with critical contextualization provided by non-colon-involved CD, and by pediatric FGID patients, who have GI symptoms, but lack IBD-associated inflammation. Most cellular profiling of adult or pediatric IBD has focused on either CD or UC, without direct comparison between the two disease states ^8,9,30,32,94–97^. Although these studies have been crucial for CD and UC characterization, the lack of direct comparison within the same study prevents the elucidation of critical nuances both between and within the two diseases. The PREDICT data underscore the importance of these direct comparisons, given that, while there were multiple immune states that were dysregulated in both CD and UC compared to FGID, close examination revealed distinct patterns of pathology between these two subtypes, revealing cell states and immune pathways specifically enriched in CD or UC. Previously published studies that have directly compared single-cell transcriptomic profiles of CD and UC have been of much smaller scale (often comparing 6 or fewer patients in one or both groups), lack contextualization with non-inflamed contemporaneously-enrolled controls, and lack a fully treatment-naïve group^31,98–100^, which would limit the disease-specific inferences that could be drawn.
This work is most complementary to a recent study of transcriptomic differences in biologic-naïve (but pre-treated) adults with CD and UC, which utilized 3’ scRNA-Seq to describe 109 cell states with similar phenotypes associated with either CD (including expansion of TH1 cells and robust interferon stimulation signatures), or UC (including IgG plasma cells, CXCL13+ TFH/TPH cells, and TH17 cells)^3^. The similarity of these overarching features suggests the pathogenic mechanisms governing the generation of these cell states are preserved between pediatric and adult populations. Our study has additionally enabled an in-depth analysis of T cell plasticity (providing evidence for ex-TH17 cells comprising the UC TH1 cell state, as well as clonal proliferation of CD8+ TH17 cells in UC), TCR phenotypes (demonstrating increased polyreactivity in UC TFH cells), and TFH vs. TH17 dynamics (identifying an inverse correlation between TFH features and TH17 features in UC), highlighting the central role T cells play in delineating IBD subtypes and treatment response paradigms.
PREDICT demonstrated a profound reduction of TH1 cells and TH1/IFNG-associated cell states upon anti-TNF treatment, consistent with the hypothesis that T cell abrogation is a primary mechanism of anti-TNF efficacy^101–103^. In vitro attempts to validate this hypothesis have further suggested that anti-TNF’s ability to prevent TNFRII pro-survival signaling, likely from myeloid cells, is a primary mechanism of drug efficacy^104,105^. Our results are consistent with a model whereby anti-TNF limits TH1 differentiation, proliferation, or survival, reducing pathologic IFNγ in the inflamed tissue and contributing to a shift away from inflammatory T cell, NK cell, and macrophage populations. Whereas anti-TNF successfully abrogates TH1-mediated inflammation, our results suggest that it is less efficacious at correcting TH17-associated pathology^105^. In PREDICT, IL17A gene expression was associated with lack of response to anti-TNF, either prior to treatment (in UC), or due to a shift in T cell phenotype following treatment (in cCD). Previous studies have identified a transition to increased TH17 features, including increased IL17A expression, in treatment-refractory adult CD^105,106^.
Given these results, it is, perhaps, counterintuitive that selected case studies have reported new onset IBD following anti-IL-17A treatment of other rheumatological diseases, and that the only IL-17A cytokine blockade trial in IBD was halted due to adverse events in CD patients receiving the anti-IL17A treatment^107–112^. Moreover, due to IL-17’s other roles in host defense and tissue remodeling, it is likely that this cytokine contributes to both persistence of inflammation and resolution of injury^36,113^. Rather than targeting IL-17 directly, our data suggest that preventing TH17 differentiation, for example via JAK inhibition or IL-12/IL-23 blockade, may be more effective in treatment-naïve IBD patients who exhibit greater TH17 features, inclusive of IL17A expression. This is further reinforced by the emerging use of JAK and IL-12/IL-23 inhibition in IBD treatment, including in those who have failed anti-TNF therapy^72,114–118^. Additionally, neutrophils have also been previously linked to anti-TNF nonresponse^12^. In PREDICT, IHC data demonstrated that neutrophil density was significantly increased in UC. Preventing neutrophil recruitment or activation, for example through IL-1 inhibition, may also be a therapeutic strategy for these patients^12^.
Our analysis provides a framework for combining gene expression and TCR analysis to identify distinct, cell state-influenced TCR phenotypes in autoimmune diseases. Although UC CD4 T cells overall demonstrated an increased disposition towards polyreactive TCR features (described primarily through mTiRP^78^), this signature was particularly enriched in UC PD1Hi TFH cells. This could represent increased reactivity to self-antigens, or potentially to microbial antigens, such as those of the commensal microbiome^78,79^. Additionally, the association of these TCR phenotypes with TLS-like B cells and IgG1 plasma cell frequencies may suggest that lymphoid aggregate formation in UC is due to UC-specific antigenic targets in the tissue. Previous studies have suggested that antibodies in UC primarily target commensal antigens, although some targeting of pathogenic bacteria has been observed as well^27,119,120^.
Our evaluation of polyreactive TCR features in TFH highlighted the opposing TH17 and TFH biology within UC, which was consistent with increased germinal center activity in some UC patients. To further evaluate this, we performed a proteome-wide autoantibody screen, to determine if the increased T cell polyreactivity in UC was mirrored by increased antibody-mediated autoreactivity. Our results demonstrated a complex autoantibody landscape in all PRECICT patients, with substantial inter-patient variability being the most prominent feature of the data, and with examples of patients with high autoantibodies in FGID, CD and UC. These data are consistent with multiple factors likely influencing the largely patient-specific autoantibody profile that was measured in PREDICT, including HLA genotype, variations in dietary and microbial exposures, and robustness of humoral immune response. The influence of HLA genotype on the autoantibody profile was highlighted by an evaluation of two UC patients who both expressed the HLA-DRB1*15:01 allele, which has been associated with T cell autoreactivity and cross-reactivity in multiple autoimmune diseases^86,121–124^. These two UC patients were both characterized by TFH features rather than TH17 features, and both developed autoantibodies to several antigens not present in any other PREDICT patient, including to ENO1 (a known autoantigen in autoimmune diseases including Hashimoto’s encephalopathy, asthma, juvenile idiopathic arthritis, and Behçet’s disease)^87–91^. These patients highlight the (often unmeasured) influence that HLA genetics can have on the propensity for autoantibody production, underscoring the complex integration of multiple factors in IBD patients that ultimately produce each patient’s unique clinical presentation.
Intriguingly, an experimental murine model of periodontal disease, which is dominated by TH17 immunopathology, demonstrated that TH17-to-TFH plasticity was critical for generating IgG, restricting growth of oral commensals, and decreasing pathologic tissue damage^125^. Given the similarities to the findings in the PREDICT study, we hypothesize that polyreactive TCR phenotypes may predispose T cells to TH17-to-TFH plasticity that could limit TH17-mediated colonic inflammation and prevent nonresponse to anti-TNF therapy. While further work is required to address the relationship between TH17 and TFH cell biology, it is interesting to note that, in contrast to recent adult UC cohorts which revealed a dominant TPH cell signature together with plasma cell infiltration in UC, pediatric UC harbors a co-existence of TH17 and TFH cells of related clonotypes. Thus, this careful study of pediatric participants may be helping to reveal an earlier “transition state” in this disease prior to the formation of more organized lymphoid aggregates.
In summary, our results provide novel insights into the overlapping and distinct biology of pediatric cCD and UC, identifying three prominent axes of immunopathology (TH1, TH17, and TFH) whose interactions and relative magnitude define these two IBD subtypes. They delineate which immune pathways portend suboptimal response to anti-TNF therapy and nominate data-driven approaches to IBD treatment.
All raw and processed gene expression data will be made available upon publication of the peer-reviewed version of this manuscript.
Pediatric patients were enrolled in the Precision Diagnostics in Inflammatory Bowel Disease, Cellular Therapies, and Transplantation (PREDICT) trial (clinicaltrials.gov #NCT03369353), as described previously^9^. Inclusion criteria required participants to be 6 years or older, weigh 10 kg or greater, and to be undergoing evaluation for new diagnosis of inflammatory bowel disease (IBD) or functional gastrointestinal disorder (FGID). Enrollment took place in accordance with Fred Hutchinson Cancer Center Institutional Review Board (Protocol #9730, ethical approval given) and Boston Children’s Hospital Institutional Review Board (Protocol IRB-P00030890, ethical approval given) approved protocols, with written informed consent and assent when applicable. Enrollment for the trial began on May 1, 2017. For this analysis, only enrollees who underwent diagnostic endoscopy prior to February 1, 2022 were included. Diagnoses of FGID, Crohn’s disease (CD), and ulcerative colitis (UC) were made by the patient’s clinician based on endoscopic, histopathologic, and clinical features. Patients were designated as ileal-only Crohn’s disease (iCD) if inflammation was present in the terminal ileum and not in any section of the colon. Designation of colon-involved Crohn’s disease (cCD) required inflammation anywhere in the large intestine or in perianal area. Patients were excluded from the analysis if diagnosed with indeterminate IBD (2 cases), inflammation due to infectious etiology (1 case), or a non-FGID disorder (1 case revised to neurogenic bowel dysfunction). Clinical variables were recorded from enrollment through 2 years post-diagnostic endoscopy. Clinical variables recorded included age, sex, race, disease location and phenotype (using Montreal Criteria), clinical disease severity (calculated using Pediatric Crohn’s Disease Activity Index (PCDAI) or Pediatric Ulcerative Colitis Activity Index (PUCAI) for CD and UC, respectively), anti-TNF use, and anti-TNF response category.
For patients administered anti-TNF (5 iCD patients, 27 cCD patients, 8 UC patients), therapy was initiated within 90 days of diagnostic endoscopy. IBD patients who were not administered anti-TNF within this window were denoted as not on anti-TNF (NOA). Anti-TNF response status was determined at 2 years after diagnostic endoscopy. Full response (FR) was defined as clinical symptom control and biochemical response with a PCDAI score less than 12.5 or PUCAI score less than 10 on maintenance anti-TNF therapy. Partial response (PR) was defined as a lack of complete clinical symptom control and biochemical response with documented dose escalation of anti-TNF therapy and a possible change to a new therapeutic. Nonresponse (NR) was defined as complete lack of clinical and biochemical response to anti-TNF therapy.
Colon and whole blood biopsies were obtained during the diagnostic endoscopy visit. Whole blood was drawn into EDTA anticoagulant and frozen at −80° C until required for genomic DNA experiments. Per protocol, up to 2 pinch biopsies were taken from non-necrotic areas of the colon, with preference for areas with marked erythema or edema. If no inflammation was observed, biopsies were obtained from left colon by default. Pinch biopsies were immediately placed into RPMI 1640 (ThermoFisher, 21870–076) on ice before subsequent processing.
Pinch biopsy processing was performed using a modified version of a previously published protocol^9,126^. While intact, biopsy pinches were handled using a P1000 pipette applying gentle suction, and all centrifugation steps done in a temperature-controlled 4°C centrifuge. Biopsy pinches were first rinsed in 25 mL PBS (ThermoFisher 14190–144) and allowed to settle. Each individual pinch was then transferred to 10 mL epithelial cell solution (ECS) (HBSS Ca/Mg-Free [ThermoFisher 14175–103], 100 U/mL penicillin [ThermoFisher 15140–122], 100 μg/mL streptomycin [ThermoFisher 15140–122], 10 mM HEPES [ThermoFisher 15630–080], and 2% FBS [ThermoFisher SH3007103HI]), freshly supplemented with 10 mM EDTA [ThermoFisher AM9261]. Separation of the epithelial layer (EPI) from the underlying lamina propria (LP) was performed for 15 minutes at 37°C with rotation at 700 RPM. The tube was then removed and placed on ice immediately for 10 minutes before shaking vigorously 20 times and undergoing 10 seconds on high vortex. Visual macroscopic inspection of the tube at this point yielded visible epithelial sheets.
The LP tissue fraction was carefully removed and placed into a large volume (45 mL) of ice-cold PBS to rinse before transferring to 10 mL of enzymatic digestion mix (Base: RPMI1640 [ThermoFisher 21870–076], 100 U/ml penicillin [ThermoFisher 15140–122], 100 μg/mL streptomycin [ThermoFisher 15140–122], 10 mM HEPES [ThermoFisher 15630–080], 2% FBS [ThermoFisher SH3007103HI], and 50 μg/mL gentamicin [ThermoFisher 15750–060]), freshly supplemented immediately prior with 100 μg/mL of Liberase TM [Roche 5401127001] and 100 μg/mL of DNase I [Sigma D5025–150KU]), at 37°C with rotation at 700 rpm for 30 minutes. LP enzymatic dissociation was quenched by addition 80 μL of 0.5M EDTA and placed on ice for five minutes. Samples were typically fully dissociated at this step and after gentle trituration with a P1000 pipette filtered through a 40 μm cell strainer (VWR 21008–949) into a new 50 mL conical tube and rinsed with ice-cold PBS to 35 mL total volume. This tube was spun down at 500g for 10 minutes and resuspended in 1.2 mL ECS. The tube was then spun down at 800g for 2 minutes and resuspended in 500 μL of ACK lysis buffer (Gibco A10492–01) and placed on ice for 3 minutes. 500 μL 0.2% FBS (ThermoFisher SH3007103HI) in DPBS (ThermoFisher 14-190-144) was added to post-lysis solution before centrifugation at 800g for 2 minutes and resuspension in 700 μL ice-cold PBS with 0.4% BSA (ThermoFisher AM2616). LP cells were filtered through a 40 μm strainer then rinsed with an additional 700 μL ice-cold PBS with 0.4% BSA. LP cells were centrifuged at 500g for 5 minutes, then resuspended in 60 μL PBS + 0.4% BSA.
During LP processing, ECS was added to EPI-containing tube to 35 reach mL before centrifugation at 600g for 10 minutes and resuspension in 1 mL ECS. Cells were transferred to a 1.5mL Eppendorf tube pre-coated with 0.4% BSA in PBS in order to minimize time spent centrifuging and obtain a more concentrated cell pellet. Cells were spun down at 800g for 3 minutes, resuspended with 1 mL ECS, and centrifuged again at 800g for 3 minutes. Pellets were resuspended in ACK lysis buffer [ThermoFisher A1049201] for 2 minutes on ice to remove red blood cells, even if no RBC contamination was visibly observed, in order to maintain consistency across samples. 500 μL ECS without FBS was added to EPI cells before centrifugation at 800g for 5 minutes. Cells were washed (resuspended in 1 mL ECS without FBS, centrifuged at 800g for 3 minutes) and resuspended in 1 mL 37° C TrypLE express enzyme [ThermoFisher 12604–013]. After adequate mixing, another 400 μL TrypLE solution was added before additional mixing. Cells were incubated in TrypLE for 4 minutes in a 37°C bath followed by gentle trituration with a P1000 pipette. Cells were incubated for another 3 minutes in a 37°C bath followed by gentle trituration, and an additional 3 minute incubation at 37° C. Cells were centrifuged at 800g for 3 minutes then resuspended in 700 μL PBS + 0.4% BSA. Resuspended cells were filtered through 40 μm strainer and rinsed with 700 μL PBS + 0.4% BSA before being transferred to a new BSA-coated tube. Cells were centrifuged at 500g for 5 minutes and resuspended in 50 μL PBS + 0.4% BSA.
Cells from both EPI and LP fractions were counted and prepared as a single-cell suspension for scRNA-seq.
Single cells were loaded onto 5’ library chips as per the manufacturer’s protocol for Chromium Single Cell 5’ Library and Gel Bead Kit (10X Genomics). Biopsies were processed with 5’ v1 kits November 2020 and prior (Chromium Single Cell 5ʹ Library and Gel Bead Kit, PN-1000006; Chromium Single Cell V(D)J Enrichment Kit, Human T Cell, PN-1000005), or 5’ v2 kits (Chromium Next GEM Single Cell 5ʹ v2 Library and Gel Bead Kit, PN-1000263; Chromium Single Cell V(D)J Enrichment Kit, Human T Cell, PN-1000252) since November 2020. The EPI and LP fractions were processed in separate channels of the 10X Chromium Single Cell Platform. An input of 20,000 single cells was added to each channel. Briefly, single cells were portioned into Gel Beads in Emulsion (GEMs) in the Chromium controller with cell lysis and barcoded reverse transcription of RNA, followed by cDNA amplification, enzymatic fragmentation and 5’ adaptor and sample index attachment.
For single cell libraries prepared using the 5’ v2 assay from 10x Genomics, paired-end sequencing was performed on Illumina NovaSeq 6000 for RNA-seq libraries (Read 1 26-bp for UMI and cell barcode, Read 2 80-bp for transcript read, with 10-bp i7 and 10-bp i5 reads) and for V(D)J libraries (Read 1 150-bp, 10-bp i7, 10-bp i5, Read 2 150-bp). For single cell libraries prepared using the 5’ v1 assay from 10x Genomics, paired-end sequencing was performed on Illumina NovaSeq 6000 for RNA-seq libraries (Read 1 26-bp for UMI and cell barcode, Read 2 80-bp for transcript read, with 8-bp i7 and 0-bp i5 reads) and for V(D)J libraries (Read 1 150-bp, 8-bp i7, 0-bp i5, Read 2 150-bp). Quality-filtered base calls were converted to demultiplexed FASTQ files. For all RNA-seq libraries, Cell Ranger Single Cell Software Suite (10X Genomics, 7.2.0) was used to perform transcriptome alignment, filtering, and UMI counting. The human GRCh38–1.2.0 genome assembly was used for the alignment. Intronic reads were not included in alignment. For all V(D)J libraries, Cell Ranger Single Cell Software Suite (10X Genomics, 7.2.0) was used to perform de novo assembly of read pairs into contigs followed by alignment and annotations of contigs against the germline segment V(D)J reference sequences from the GRCh38-alts-ensembl-2.0.0 genome assembly.
For all aligned RNA-seq libraries, CellBender “remove-background” was used to reduce contaminating signal due to ambient mRNA^127,128^. Raw Cell Ranger-aligned counts matrix were used as input. Default remove-background settings were used for all samples, with exception of learning rate (set to 3×10^−5^ for all samples, instead of default 10×10^−5^)
Cell x gene matrix outputs from CellBender were aggregated into a single object using Seurat (version 5.0)^129^. Samples were converted to AnnData format using zellkonverter^130^. The AnnData object was processed with scanpy log-normalization, variable gene selection (flavor = seurat_v3, batch = 10x version, n_top_genes = 3000), PCA dimensional reduction, and coarse leiden clustering (n_neighbors = 20, n_pcs = 20, leiden resolution = 0.1)^10^. Coarse clusters were grouped into one of 6 Colonic Epithelium (ColEpi); Stromal and Endothelial Cells (Stroma); Myeloid Cells (Myeloid); Plasma Cells (Plasma); B Cells (BCell); and T Cells, NK Cells, and Innate Lymphoid Cells (TNKILC). Clusters that were low-quality based on mitochondrial gene expression, number of features detected, number of UMIs detected, or a combination of these metrics, were not included in any of the above groups and discarded from further processing. Sample-specific clusters were also not included.
After partitioning into coarse cell types, each cell type was log-normalized and underwent 3 rounds of scVI model training, clustering, differential expression, annotation, and low-quality cell removal. For each round, an scVI model was trained and the scVI latent space was used for dimensional reduction, and between-cluster differential expression performed within the scVI framework^131^. Settings for scVI: n_layers = 2, n_latent = 30, n_hidden = 128, dropout_rate = 0.1, gene_likelihood = “nb”, max_epochs = 250 (with early stopping), batch_key = “version_10x”. AnnData objects were reduced to top n variable genes prior to scVI model training. Number of variable gene selected differed by cellgroup (number of genes used for scVI input in final round of annotation and QC: ColEpi = 2000; Stroma = 5000; Myeloid = 5000; Plasma = 750; BCell = 3000; TNKILC = 5000). For Plasma and BCell types, additional gene filtering was performed to exclude variable immunoglobulin genes (pattern matching for “ÎG[HKL][VJC]”). This step was taken to reduce the ability of variable immune receptor genes to dominate model fitting and latent space generation in plasma cells/plasmablasts (constant genes were retained). Additional genes matching the following patterns were removed from scVI training as they appeared to dominate clustering without clear biological “ĈH17’”, “^RP11”, “ÎGLL”, “ÂC[0–9]{5}”. Conversely, genes matching the following pattern were explicitly included as input (if not already an highly variable gene) for plasma cell scVI model training, to ensure constant genes would be “IGH[ADEGM][1–4]*$”
For cluster annotation, annotation was based on marker genes identified through literature review. Annotations were defined at additional levels within each cell type (level 1): cell subtype (level 2) and cell state (level 3). If appropriate and supported by literature, clusters were recombined into overarching cell types to obtain “level 2 annotations.” Conversely, if clusters were observed to be heterogeneous through differential expression analysis and visualization, subclustering was attempted, along with differential expression and annotation to obtain level 3 annotations (“cell state”). Annotations from final round of scVI processing were used for this study.
After Cell Ranger alignment to the V(D)J reference, the filtered contigs files from each sample were combined into a single data frame. In some cases, more than one alpha chain was detected per barcode. In very rare cases, more than one beta chain was detected per barcode. To simplify data analysis, a single alpha or beta chain was selected for downstream processing. First, barcodes with only a CDR3 listed were excluded (these chains lacked values for CDRs 1 and 2 as well as framework regions 1–4). The selection preference order (1) which chain had the most UMIs, (2) which chain had the most reads, (3) which chain is listed first (this condition was rarely reached in the dataset). Cell annotations from the RNA-seq processing step were merged with TCR data through paired barcode sequences.
A fully automated multiplex immunohistochemistry assay was performed on the Ventana Discovery ULTRA platform (Ventana Medical Systems, Tucson, AZ) as previously described^9,132^. Optimal concentrations of each antibody were determined, and they were applied in the following sequence and detected with the indicated fluorophore.
Following staining, the tissue was counter-stained and cover slipped with Invitrogen ProLong Gold Antifade Mountant with NucBlue. Whole slide imaging was performed on the Zeiss Axioscan which was equipped with a Colibri light source.
For analyses of common variants, we used Illumina GSA array genotyping data and performed imputation with the TOPMed imputation server and reference panel^133^. Exome sequencing was performed at the Regeneron Genetics Center using a custom automated sample preparation approach. Samples were captured with IDT xGen v1 probes and sequenced using Illumina HiSeq 2500-v4 or Illumina NovaSeq instruments, with 75-bp paired-end reads and two index reads. The GRCh38 human genome reference sequence and Ensembl version 100 gene definitions were used for variant calling using DeepVariant and for annotation^134^.
We imputed HLA alleles using the Michigan Imputation Server^133^, HLA reference panel multiethnic-hla-panel-4digit-v2@1.0.0^135^. Prior to imputation, in order to maximize the number of informative variants for HLA imputation, the PREDICT IBD sample was merged with n=1,000 random individuals from a Regeneron Genetics Center multi-ethnic cohort genotyped on the same array version. Variants in the MHC region 28–34Mb and present in the HLA reference panel were extracted from TOPMED imputed genotype probabilities and hard-called with certainty threshold 0.9. The resulting dataset was QC-filtered for 12,101 variants with less than 0.1% missing genotype calls and HWE Pval > 1e-10. (All samples had less than 1% missing genotype calls.)
The HuProt proteome microarray (CDI Labs, USA) was used to assay plasma for autoreactivity against the human proteome^80–83^. HuProt arrays and plasma samples were diluted with CDI proprietary CDISampleBuffer (CDI Labs, USA) at 1000 and incubated with CDIArrayBlock (CDI Labs, USA) at room temperature with gentle shaking for 1 hour. After the blocking period, each sample was probed onto a HuProt microarray at room temperature for 1 hour with gentle shaking. The arrays were then washed with TBS-Tween (1x TBS / 0.1% Tween-20) 3 times for 10 minutes each. Following washes, the arrays were probed with fluorescent anti-human IgA or IgG secondary for 1 hour in the dark with gentle shaking. 3 washes were performed with TBS-Tween for 10 minutes each, followed by 3 rinses with ddH20. The arrays were dried with compressed CO2. GenePix 4000B scanner was used to scan arrays. Scanned images were aligned to the array layout. Raw fluorescence intensity values were extracted using GenePix software for further data processing.
Gene set enrichment analysis was performed using fgsea^136^. The DESeq2 “stat” value (Log2FC/SElog2FC) was used as gene-level statistical input for fgsea. 1000 permutations were used for fgsea calculation of enrichment scores. Reactome gene sets were used as input (version 2023.2, obtained through MSigDB C2 curated gene set)^137,138^. Gene sets were included if they had at least 15 genes and no more than 100 genes.
Per-sample mean cell frequencies were calculated using custom code. Cells were grouped first by their sample or biopsy of origin (relevant grouping listed in text and/or figure captions; “sample” will be used throughout methods), then by the annotation of interest (e.g. CD4 Naïve T cells). The number of cells per sample and annotation was tallied. Unless otherwise stated, the number of cells was divided by a relative grouping, e.g. a cell subtype (e.g. Naïve T cells), a cell type (e.g. T cells), or the entire independent sample (e.g. “Relative to Sample”). For per-biopsy group comparisons, biopsies were only included for analysis if the relative grouping had at least 10 cells. This was an empirical measure taken to prevent statistical conclusions driven by low overall cell detection. Statistical testing was performed with ggpubr’s “stat_compare_means()” function^139^.
Per-sample quantitative feature (e.g. gene expression, module score, mTiRP, etc.) values were calculated using custom code. Data were aggregated based on sample of origin, then by cellular annotation(s) of interest. Mean values were calculated per combination of sample and cell annotation. For gene expression, log-normalized values were used as input for mean calculation. Samples were required to meet a minimum number of cells within sample and cell annotation combination to be included for analysis. 10 cells were required for calculation of gene-based metrics, and 50 cells were required for TCR-based metrics.
Only baseline, non-iCD samples were used as input for analysis with MiloR^13^. Milo analysis was performed on a per cell type basis (ColEpi, Stroma, Myeloid, Plasma, BCell, TNKILC). For each cell type, the data object was downsampled to the number of cells in the smallest disease category between FGID, cCD, and UC. This step was performed to reduce the effect of differing total number of cells skewing differential abundance results. Downsampling was performed on a per disease category basis, irrespective of sample of origin within cell types. After downsampling, any remaining biopsies with fewer than 5 cells were removed before proceeding with analysis. Following biopsy removal step, the Milo object was created and kNN graph generated with the cell group scVI latent space used as the reduced dimension (all 30 scVI dimensions used as input). Because there were many samples in this analysis (~70), we set k = 200 to ensure the average neighborhood would have cells numbering greater than 5 x number of samples (reducing the likelihood of a single sample comprising a neighborhood), per package maintainer recommendations^13^. We verified that each cell group fulfilled this heuristic, then calculated neighborhood distances and built the neighborhood graph, using Milo functions with default inputs and latent space inputs as described above. Following graph construction, we tested for differential abundance with the scVI latent space used as the reduced dimension and a design formula specifying 10x version as a batch “~version_10x + disease_group”. Alpha value was set to 0.1 for identifying differentially abundant neighborhoods. Neighborhoods with annotation fraction < 0.7 were excluded from visualization. 12 cell states did not meet the annotation fraction threshold in any neighborhood and were not Lymph. Endo., IgG4 PC, Prlf. IgG PC, Prlf. IgM PC, IgD PC, IgE PC, Plasmablast, Prlf. Th17.1, Prlf. Tfh, Prlf. Treg, Vg9d2 TC, and NKT.
Cell state correlation analysis was performed using per-cell type frequencies of annotated cell states. Cell states were excluded if they did not have at least 100 cells (2 cell states IgD PC and IgE PC). Only baseline biopsies were included for correlation analysis. The previously described cell frequency calculation method was used to generate per cell type frequencies for each cell state in each biopsy. Pearson correlation was used to compare biopsy frequencies of TH1, TH17, and PD1hi TFH cells against all cell states excluding IgD PC and IgE PC. Test statistics were generated using “cor.test” function from the R stats package (version 4.3.1). Cell state correlations with a p value < 0.001 were selected as part of TH1, TH17, or PD1Hi TFH network. Cell states in network were visualized with ggraph (version 2.2.1). Only edges with a correlation coefficient > 0.35 were visualized. Graph layout algorithm set to “kk”, maximum iterations = 150. Node location manually adjusted when necessary for visualization of all cell states.
To provide a quantitative measure of the usage of cell states, we calculated the geometric mean of cell state frequencies, as has been performed previously^30^. Briefly, we took the logarithm value of the associated cell state frequencies in each biopsy, calculated the mean of all cell state frequencies in logarithm space, then exponentiated the average exp(∏i=1klog(xi))
Where x is the cell state proportion, i is the index of the cell state in the correlation network, and k is the number of total cell states in the correlation network. Because calculation of the geometric mean requires a logarithmic transformation, we added a “pseudoproportion” of 1/1000 to any cell state which was not identified in a biopsy.
For differential expression (DE) analysis, a single-cell object containing aggregated, annotated cells from each cell type was filtered to include only the biopsies of interest for that comparison. Raw counts from each cell type within each biopsy were aggregated into pseudobulk format using the Seurat “AggregateExpression()” function^140^. Gene selection was performed by filtering for genes with ≥ 250 counts in at least as many samples as were in the smaller of the two groups being compared (e.g. if comparing 15 cCD FR and 10 cCD PR samples, genes were included if there were ≥ 250 counts in at least 10 samples overall in the cell type of interest). For the combined myeloid and TNKILC pre- vs. post-treatment analysis, a 500 count minimum was used. All input genes are included in supplemental DE tables. DESeq2 was run with default settings, using UMI data from the pseudobulk object^141^.
Decoupler python package (version 1.5.0) was used for all gene set scoring^50^. For TNF scoring, the PROGENy database was used to select genes involved in TNF pathway activity^51^. Top 500 genes were obtained initially. The genes were then filtered to only include those with p-value < 1×10^−10^, to increase confidence that the genes included in the network play a role in TNF signaling. This resulted in 225 input genes. These genes were used as input for decoupler’s multivariate linear model (mlm), with weights provided by PROGENy database entries^50,51^. Decoupler was performed on individual cell types (e.g. Myeloid, TNKILC) and mlm estimate values were used for data analysis. For IL-17 pathway scoring, genes that were listed as cell extrinsic and downstream of IL-17RA/IL-17RC signaling in Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway “IL-17 SIGNALING PATHWAY” were investigated for gene expression in this dataset^46^. Genes with robustly detectable expression were retained for further analysis and grouped into KEGG-provided gene Chemokines (CXCL1, CXCL2, CXCL5, CXCL8, CXCL10, CCL2, CCL7, CCL20), Cytokines (IL6, TNF, PTGS2, CSF2, CSF3), Anti-microbial (DEFA5, DEFA6, DEFB1, DEFB4A, DEFB4B, DEFB124, MUC5B, S100A7, S100A8, S100A9, LCN2), Tissue remodeling (MMP1, MMP3, MMP9, MMP13). An unweighted network of these genes was created with the gene group as the source and the gene as the target. This network was used as input with each cell type for decoupler’s weighted sum algorithm^50^. However, equal weight was given to all genes in the network. Normalized wsum values were used for data analysis. The same analytical steps as IL-17 gene set scoring were used for TH1 and TH17 cytokine (IFNG, TNF, IL17A, IL17F, IL21, IL22, IL26) scoring.
Comparison of multiple metrics through regression analysis was performed using a combination of custom code and ggpmisc (version 0.5.5)^142^. Cell frequency and feature quantification was performed as described in “Cell frequency analysis” and “Quantitative feature analysis.” The regression analysis was performed using “stat_poly_eq” (for statistics) and “stat_poly_line” (for visualization) functions. The method was set to “lm” for all regressions. A 95% confidence interval used for shading.
A T cell was defined as clonally expanded if its paired CDR3α and CDR3β amino acid sequences were observed in at least one other cell in the same biopsy. The number of cells in a clonotype (paired CDR3α+CDR3β) per biopsy was quantified and reported as the level of clonal 1 (singlet), 2, 3, 4, or 5+. Cells that did not have a CDR3α or CDR3β were not considered for clonal expansion analysis.
For TH-TFH cell clonal overlap analysis, the dataset was reduced to TH1, TH17, PD1Hi TFH, and PD1Lo TFH cell states. Additional filtering was performed to ensure paired CDR3α+CDR3β was present in each cell. Cells were grouped by paired CDR3 per biopsy and the number of cell subtypes within each biopsy+CDR3 pair was determined. Clonotypes spanning more than one cell subtype were considered to have clonal overlap. Visualization of clonal overlap was performed with ggraph (version 2.2.1).
Clonal overlap between Prlf. TH17.1 and TH1 or TH17 was performed in an identical manner, except that cell states were used to define overlap rather than cell subtypes.
Quantitative image analysis was performed using HALO Indica Labs Hyperplex module (IndicaLabs, Albuquerque, NM). For each sample, images from the pan-immune panels and the pan-CK marker were fused to generate a single image. A classifier was first applied to detect the tissue in each section and an automatic annotation was generated as “whole section.” A second classifier was applied using panCK and DAPI channels to define the epithelium layer (panCK+) and the lamina propria (panCK-). Automated annotations were generated as “epithelium” and as “lamina propria”. Tertiary lymphoid structures (TLS) were annotated by a pathologist. Their number, frequency, size, and cell composition were measured. Immune cells were detected in each area of interest (whole section, epithelium, and lamina propria for pan-immune panel 1; whole section, epithelium, lamina propria, and tertiary lymphoid structures for pan-immune panel 2). T cells were defined as CD3+, B cells as CD20+, myeloid cells as CD68+. For each T cell subset, CD8 T cells were defined as CD3+CD8+FOXP3-, CD4 as CD3+CD8-, Tregs as CD3+CD8-FOXP3+ and CD4 conventional (non-Treg CD4) as CD3+CD8-FOXP3-. Plasma cells were defined as BCMA+. Eosinophils were defined as ECP+. Dendritic cells were defined as CD11c+. Neutrophils were defined as MPO+. NK cells were defined as NCR1+. Numbers of positive cells for each immune subset were counted and their density measured.
Diffusion pseudotime was calculated using scanpy^10^. For pseudotime analysis, a subset of the TNKILC cell group was created from CD4 T cell states (NaiveCD4, MemoryCD4, Int. TH, TH1, TH17, Prlf. TH17.1, PD1Lo TFH, PD1Hi TFH, Prlf. TFH, Treg, LEF1+ Treg, LAG3+ Treg, Prlf. Treg, CD4 CTL). Nearest neighbors were recalculated on this subset, with n_neighbors = 20 and all 30 scVI latent dimensions used. The starting position for the pseudotime analysis was selected as the cell with the most positive value in the UMAP2 dimension, which corresponded to a cell in the NaiveCD4 cell state. Diffusion pseudotime was calculated using scanpy’s “tl.dpt” function, with n_branchings = 0 and n_dcs = 10. Gene expression of cells ordered by pseudotime was visualized using scanpy’s “pl.paga_path” function, with “normalize_to_zero_one = True”.
CDR3 length was quantified as the number of amino acids in the CDR3, as determined by Cell Ranger alignment and TCR classification. Previous studies have highlighted that the middle residues of the CDR3β, determined by N and P nucleotides between regions belonging to the V gene and J gene, are most likely to directly participate in contacts between TCRs and antigenic peptides^75^. Thus, we focused our subsequent biophysical analysis on the middle region of the CDR3β. The CDR3β middle region was calculated only for TCRs with CDR3β lengths between 12–17 amino acids, in accordance with definition for T cell regulatory potential (TiRP) score^78^. The Cell Ranger definition of CDR3 includes junctional amino acids at IMGT positions 104 and 118^143^. The CDR3β middle region was defined as all amino acids from the 5^th^ amino acid from the beginning to the 6^th^ from the end, corresponding to P108 and P112 in IMGT numbering^143^. mTiRP score was calculated using publicly available code from Lagattuta et al. (https://github.com/immunogenomics/TiRP), with Cell Ranger-determined V gene and CDR3β used as input^78^. CDR3β middle-region biophysical properties, including hydrophobicity (“grand average of hydrophobicity”, or “gravy”), size, aliphatic propensity, and polarity, were calculated with “alakazam” R package (version 1.3.0). CDR3β middle-region amino acids were used as input to alakazam aminoAcidProperties() function with nucleotides set to FALSE, all other settings default. Normalized z-scores were generated through transformation to standard score within the z=X-μσ where X is the vector of all metric values, μ is the metric’s mean within the dataset, and σ is the metric’s standard deviation within the dataset. Given the inherent heterogeneity in TCR features, we limited sample-comparison analysis to biopsies with at least 50 CDR3β sequences detected in the cell grouping of interest.
Data from HuProt proteome microarray was processed for autoantigen reactivity among PREDICT participants^80–83^. Briefly, median fluorescence intensity (MFI) of each replicate spot pair on the protein microarray was averaged. Replicates with values differing by two-fold or greater were discarded from further analysis. log2(MFI) was calculated to stabilize variance across cohorts. Microarray spots with readings 4-fold above secondary-only median were considered false-positives and removed. Spots that were not detected were also removed. Z-scores were calculated within each sample array for all spots passing QC. Autoantigen detection (binary variable) was defined as spots with readings 4-fold above secondary-only median for the relevant spot, and with a Z-score > +3 within that sample array. Arithmetic mean of values was used for visual representation when an antigen was present more than once on the proteome microarray, and thus had multiple independent readings. An antigen was considered a hit if the antigen met Z-score-based hit criteria in 1 or more positions on the microarray. ComplexHeatmap (version 2.24.0) was used to generate antigen heatmaps and any accompanying dendrogram using default settings^144,145^.
Features that were associated with TH1 (TH1 cell state frequency, geometric mean of TH1-associated cell state frequencies, TH
IFNG expression, myeloid and TNKILC cell type STAT1 expression, and myeloid and TNKILC TNF pathway score), TH17 (TH17 cell state frequency, geometric mean of TH17-associated cell state frequencies, TH
IL17A expression, ColEpi antimicrobial gene set score, ColEpi chemokine gene set score), and TFH (PD1Hi TFH cell state frequency, geometric mean of PD1Hi TFH-associated cell state frequencies, TFH
CXCL13 expression, IgG1 plasma cell state frequency, and TFH mTiRP score) were used as inputs to generate per-biopsy summaries of overall TH1, TH17, and TFH tendencies. Per-biopsy values for each feature were generated as described for the relevant feature. Each feature was scaled from 0 to 1, corresponding to the minimum and maximum values of that feature within the dataset. The arithmetic mean of each feature within each feature set (TH1, TH17, or TFH) was then averaged per biopsy. The average value of each scaled feature was used to represent biopsies in the per-biopsy radar plot. For group-based radar plots created with ggradar (version 0.2), the arithmetic mean of the scaled feature set among all biopsies within the group was used to represent the group value on the radar plot.