Mucosal Gene Transcript Signatures in Treatment Naïve Inflammatory Bowel Disease: A Comparative Analysis of Disease to Symptomatic and Healthy Controls in the European IBD-Character Cohort

Introduction

Inflammatory Bowel Disease (IBD) is a chronic inflammatory disease of the gut with heterogeneous manifestations, and a clinical presentation as Crohn’s disease (CD), ulcerative colitis (UC), or IBD unclassified (IBD-U).^1–3

The disease course of IBD shows large individual variation in terms of severity, behavior, progression, and response to treatment. Both clinical variables^4–6 and molecular signatures^7–11 have been associated with different aspects of the disease course, but only a few have been shown to be robust and reliable enough to be used for decision-making in clinical practice with regard to therapy.^12–14

A transcriptome-wide analysis is a hypothesis-free strategy to provide insights into the differences in gene regulation in biological samples. Various studies have associated transcriptional changes with IBD.^15–19 These studies have paved the way towards the first map of mRNA expression changes characterizing the disease, and may be used for the identification of novel therapeutic targets. Recently, efforts have been put into the exploration of global gene expression signatures to facilitate personalized IBD care. In particular, responses to anti-Tumor Necrosis Factor (anti-TNF) agents, have been explored, demonstrating that gene expression signatures hold predictive biomarker potential.^20–25

In order to advance our understanding of IBD and to increase the diagnostic and prognostic precision of the disease, the EU project “Inflammatory Bowel Disease CHARACTERization by a multi-modal integrated biomarker study” (IBD-Character)²⁶ was launched in 2013. IBD-Character is a collaborative effort, involving six European university clinics (Akershus, Norway; Ørebro, Sweden; Linkøping, Sweden; Edinburgh, United Kingdom; Zaragoza, Spain; Maastricht, Netherlands). We performed mucosal transcriptome analyses of an inception cohort of newly diagnosed, mainly treatment naïve (92%) IBD patients, to create a transcriptional snapshot of IBD in its early manifestation. The RNA transcriptional status as well as demographic and clinical data were compared between patients and control groups, and between inflamed and non-inflamed tissue in IBD patients. In addition, we examined how gene expression may be associated with the clinical outcome defined as the need for treatment escalation and response to anti-TNF agents.

Methods

Patient Recruitment, Sample Collection and Follow-Up

Patients referred to the hospital with suspected IBD were recruited prospectively at six gastrointestinal clinics across Europe during 2012–2015, as part of the EU FP7 IBD-Character study.²⁶

At inclusion, clinical data were recorded and blood samples were obtained. Mucosal biopsies were collected during ileocolonoscopy according to a study protocol and stored in Allprotect (Qiagen) stabilizing reagent at −80°C. For CD and UC patients, biopsies from both inflamed and non-inflamed tissue were collected.

In total, 323 subjects were recruited and grouped into five groups; CD (N = 75), UC (N = 87), IBD-U (N = 3), symptomatic controls (SC, N = 131), and healthy controls (HC, N = 27). Table 1 demonstrates demographics and clinical characteristics.

Table 1 Demographics, Clinical Characteristics and Sample Material

The diagnosis of IBD was based on internationally accepted diagnostic criteria, following clinical, endoscopic, histologic, and radiologic examinations.^1,3 The Montreal classification was used to characterize the clinical phenotype.^27,28 The patients were mainly diagnosed at the date of biopsy sampling. Mean time from diagnosis to inclusion was 14 days, while the median was 0 days. The patients were mainly treatment naïve (92%), and they were further followed by their local gastroenterologist according to national and European guidelines.^1,3,29 Follow-up data on disease course were available in 154/165 (93%) patients, 72 CD and 82 UC. Treatment escalation was defined as introduction of anti-TNF agents, cyclosporine or surgery introduced for a disease flare or disease refractory to primary treatment.

Response to anti-TNF therapy during follow-up was evaluated 14 weeks after introduction of anti-TNF. C-reactive protein (S-CRP) was analyzed in serum samples, and fecal calprotectin (F-calprotectin) was assessed according to the manufacturer’s protocol (BÜHLMANN fCAL™, Schönenbuch, Switzerland),³⁰ at the date of introduction of anti-TNF treatment, and in week 14, respectively.

For CD, response was defined as no use of concomitant steroids, a decrease in the Harvey-Bradshaw Index (HBI) of at least 3 points, or to ≤4, and at least one of the following two criteria: (a) a reduction in S-CRP by at least 50%, or to ≤3 mg/l from baseline, or (b) a reduction F-calprotectin by at least 30%, or to <250mg/kg. Remission was defined as HBI ≤3 points and S-CRP ≤3 mg/l, and no concomitant steroids.

For UC, response was defined as no use of concomitant steroids, a reduction from the baseline partial Mayo score of at least 3 points, and at least 30%, and a reduction in the rectal bleeding subscore of at least 1 point or an absolute rectal bleeding score of 0 or 1. Additionally, at least one of the following three criteria was achieved: (a) a reduction in S-CRP fall by at least 50%, or to ≤3 mg/l from baseline, or (b) a reduction in F-calprotectin by at least 30%, or to <250mg/kg, or (c) a reduction of at least 1 in the Endoscopic Mayo score. Remission was defined as partial Mayo score 0–1, and at least one out of three criteria achieved: F-Calprotectin <250mg/kg, S-CRP ≤3 mg/l, and/or Endoscopic Mayo score 0–1.

Individuals with gastrointestinal symptoms, but with no endoscopic or histologic signs of IBD inflammation at inclusion, and no discernible evidence of IBD or other inflammation during follow-up, were classified as symptomatic controls. Seventy-five percent of these patients had diarrhea with 4–15 bowel movements per day. The remaining 25% had abdominal pain, bloating and some had constipation. In symptomatic controls with persisting symptoms 1–2 months after inclusion, a follow-up F-calprotectin was assessed. MRI and capsule endoscopy were performed with normal findings on selected cases with high F-calprotectin. Celiac serology, b-hemoglobin, and s-ferritin were assessed in all patients. Testing for fecal occult blood (FOB) was done on clinical indication. Upper endoscopy was performed in selected cases with either/or/and: Upper GI-symptoms, anemia, positive FOB, or those with elevated anti-TTG/anti deaminated-gliadin IgG, suggestive of celiac disease. If there were findings in upper GI that could be related to symptoms, the patient was excluded.

Seventy-eight percent of the symptomatic controls ended up with a diagnosis of irritable bowel syndrome (IBS), whereas in 22% of the cases no certain diagnosis could be set, as they did not meet the IBS Rome IV criteria, and other disease was also ruled out. Patients with a possible, but no definite IBD diagnosis, were excluded from the study. Patients with ischemic colitis, microscopic colitis and infectious enterocolitis were also excluded.

Healthy controls had no gastrointestinal symptoms or other disease, and they went through rectoscopy with mucosal pinch biopsies. Additionally, fecal samples and blood samples were collected from healthy controls.

Sample Processing and Gene Expression Profiling

RNA Isolation

Total RNA was isolated using QIAsymphony RNA Kit (Qiagen, Hilden, Germany) and the protocol RNA CT 400. Prior to loading into the QIAsymphony, the tissue biopsies were homogenized in Buffer RLT Plus using a 5 mm steal bead in the FastPrep-24™ Classic Instrument (MP Biomedicals, Irvine, CA, USA) at 4 m/s for 30 sec, centrifuged and the lysate transferred to new tubes. Further processing was per the manufacturer’s protocol. RNA concentration was measured using the NanoDrop ND-1000 spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA) and using the OD260 for calculation. The RNA quality, indicated by the RNA integrity number (RIN), was assessed with an Agilent 2100 Bioanalyzer, Agilent 2100 Expert software and Agilent RNA 6000 Nano Kit (Agilent Technologies Inc., Santa Clara, CA, USA). Exclusion criteria for RNA samples were RIN number below 7.5 and RNA concentration lower than 66.7 ng/μL. Four hundred and ninety-eight RNA samples passed these criteria.

Microarray Analysis

mRNA amplification, labelling, and hybridization of the RNA samples were performed following the manufacturer’s instructions (Agilent One-Color Microarray-Based Gene Expression Analysis, Low Input Quick Amp Labeling; version 6.9). In short, 100 ng RNA was amplified and labelled with Cy3 using the Low Input Quick Amp Labeling Kit, one-color (product number 5190–2305), and the labelled cRNA was purified using the Qiagen RNeasy Mini Kit.

Amplification and labelling efficiency were controlled using NanoDrop ND-1000 spectrophotometer. The cRNA quality was assessed by measuring the RNA absorbance ratio 260nm/280nm. The specific activity of labeled cRNA (pmol Cy3 per ug cRNA) was defined to be 6 or above according to the Agilent recommendations. Four hundred and eighty-one cDNA samples passed these criteria and were eligible for microarray analyses.

Six hundred ng cRNA per sample was fragmented and applied to Agilent SurePrint G3 human Gene Expression 8x60k v2 Microarrays (G4851B). The cRNA samples were randomized on the slides while attempting to balance the potential confounding variables clinical center, diagnosis, gender, age and smoking status. After hybridization for 17 h at 65°C the microarray slides were washed, scanned with the Agilent Microarray Scanner and the data extracted using Agilent Feature Extraction Software.

Statistical Analysis

Background correction (using the normexp method with offset 16) and quantile normalizing between arrays were done using the limma R package.³¹ After removing samples with sex chromosome mismatch, unclear diagnosis and missing gut location, gene expression profiles from 437 cRNA samples remained for statistical analyses. The distribution of samples regarding inflammation status and location for the diagnostic groups is demonstrated in Table 1.

A linear regression analysis using limma was used to compare the different diagnosis groups and inflammation statuses, while adjusting for paired samples from the same subject and different gut sampling location. Transcripts with the 90% quantile below 5 were not analyzed as they were likely to be below the limit of detection. Differentially expressed genes (DEGs) were identified at a false discovery rate of 5% (p = 0.05) and log2FC >±1.

The period from the date of IBD diagnosis, when the first-line treatment was initiated, until treatment escalation, was analyzed using Cox proportional hazards models. The analyses were performed on all patients with IBD and by the inflammation status and UC and CD separately.

Associations between gene expression levels and response to anti-TNF therapy were analyzed using logistic regression analysis. As in differential abundance analyses, low abundance transcripts were not analyzed.

We used the fgsea R package¹¹ to test for enriched Gene Ontology (GO) terms in ranked lists based on the differential expression test.

Curated gene sets representing pathways known, or predicted, to be of pathophysiological relevance in IBD in general, or for the effect of the therapeutic options used for IBD patients in this study, were downloaded from the MsigDB database,³² and gene expression signatures predicting activity of these pathways were transformed into PathTracer deregulation scores.³³

Ethics

Regional committees for medical and health research ethics (REC) in South-East Norway, and Norwegian Centre for Research Data accepted the ethical aspects of the research project. All clinical centers were granted local ethics approval for this study, and all patients gave written and informed consent prior to participating in this study. This study was conducted in accordance with the Declaration of Helsinki.

Data Transparency Statement

Clinical and descriptive data of the cohort, analytic methods and microarray raw data can be shared on request to the corresponding author ([email protected]) with permission from the IBD-Character data access committee.

Author Approval

All authors had access to the study data and reviewed and approved the final manuscript.

Results

In total, 323 subjects were included: Patients with treatment naïve IBD consisting of CD (N = 75), UC (N = 87) and IBD-U (N = 3). Additionally, symptomatic controls (N = 131) and healthy controls (N = 27) were recruited. Demographics and clinical characteristics, as Montreal Classification, S-CRP, F-calprotectin are shown in Table 1.

Principal component analysis (PCA) was used to assess to which extent inflammation status and gastrointestinal segment influenced expression variability. The PCA analysis confirmed that the diagnosis groups form one combined cluster (Figure 1A). Importantly, PCA analyses showed that the main determinant of the overall expression profile was gut segment, primarily separating terminal ileum vs colorectum, although different segments of the colorectum (proximal colon vs rectum) could also be observed (Figure 1B). Inflammation status also separated the data to some extent in PCA (Figure 1C). Neither inflammation status, segment, nor diagnosis was, however, able to provide definite separation of the data. Moreover, there was no evident subgrouping in the normalized gene expression data among neither UC (Supplementary Figure 1A) nor CD (Supplementary Figure 1B), and healthy controls as judged from hierarchical cluster analyses, although for UC patients, there was a tendency that biopsies from inflamed tissue clustered together. In CD and UC, biopsies from terminal ileum clustered together, while biopsies from different segments of the colorectum were distributed across the entire cluster.

Figure 1 Principal Component Analysis plots, demonstrating separation of samples based on (A) diagnostic group, (B) gut segments and (C) inflammation status separated by color (inflamed samples in red, vs non-inflamed samples in blue) and gut segment separated by shape (samples from colon demonstrated with circles, vs samples from ileum demonstrated with triangles. Abbreviations: CD, Crohn’s disease; UC, ulcerative colitis; IBD-U, IBD-unclassified.

Figure 2 Continued.

Figure 2 Continued.

Figure 2 Crohn's disease. (A) Hierarchical clustering of differentially regulated genes. (B) Volcano plot comparing inflamed (CDi) and non-inflamed (CDni) samples in CD. (C) Volcano plot of inflamed biopsies from CD (CDi) compared to healthy controls. (D) Volcano plot of non-inflamed biopsies from CD (CDni) compared to healthy controls. (E) Gene set enrichment analyses given as Normalised Enrichment Scores (NES) in CDi (top panel) and CDni (bottom panel) versus healthy controls. “Orange” represents overrepresented GO terms, with positive NES. “Green” represents underrepresented GO terms, with negative NES. Abbreviations: CDi, inflamed CD; CDni, non-inflamed CD; HC, healthy controls; NES, Normalised Enrichment Score.

The significantly differentially expressed genes (DEGs) at 5% False Discovery Rate (FDR) and the correlation between the comparisons are presented in Table 2. A high number of DEGs were observed in inflamed samples, regardless of IBD phenotype. This picture was strikingly different in non-inflamed IBD: In non-inflamed CD 401 transcripts were significantly differentially expressed compared to symptomatic controls. In non-inflamed UC, there were no differentially expressed genes that reached statistical significance when compared to symptomatic controls. Compared to healthy controls, however, non-inflamed CD and UC samples revealed 3216 and 3431 DEGs, respectively. The CD and UC results were highly correlated (r = 0.95).

Table 2 Number of Differentially Expressed Genes at a 5% False Discovery Rate and the Correlation Between the Comparisons. Up and Down Refers to Whether the Genes are Lower or Higher Expressed in the First Group Compared with the Second

These data suggest that gene expression profiles of symptomatic individuals resemble that of IBD patients, especially UC, and diverge from that of healthy controls. Hence, healthy controls were used as a reference dataset in the downstream analyses in the two following sections. Since biopsies from healthy controls were obtained from rectum only, all downstream analyses included correction for location effects.

Crohn’s Disease

Hierarchical clustering of DEGs showed extensive heterogeneity in the gene expression signatures from patients with CD. Some tendency of clustering of the inflamed samples was observed. Although the gut segment origin explained much of the variation in the data based on PCA analyses, no systematic clustering of samples from the different locations was found (Figure 2A).

Comparing gene expression status of the inflamed CD tissues with the non-inflamed CD tissues revealed 7669 DEGs (Table 2). Volcano plots show that there was a skewed distribution towards upregulated genes (Figure 2B). Gene Set Enrichment Analyses (GSEA) of the highly differentially regulated genes (log2 FC >±1) revealed that biological processes related to inflammation were overrepresented whereas mitochondrial respiration was underrepresented in inflamed CD samples compared to non-inflamed CD samples (Supplementary Figure 2A). Thus, inflammation had a strong effect on gene expression, consistent with the PCA analyses (Figure 1C).

To test whether the strong influence on global gene expression by active inflammation could mask underlying differences among the diagnosis groups, we compared gene expression between inflamed, and non-inflamed CD, respectably, with healthy controls. The comparison of inflamed CD and healthy controls showed 418 highly dysregulated genes with log2 FC >1 and 197 genes with log2 FC <-1 (Figure 2C). The highest regulated genes (up and down) primarily represented genes involved in transmembrane transport and ion transport. GSEA showed that enrichment of processes related to proliferation in addition to immune responses. Metabolism and vesicle transport were the most significantly underrepresented gene ontology (GO) terms in inflamed CD (Figure 2E).

In non-inflamed CD, few genes were highly dysregulated compared to healthy controls, with only 67 genes with log2 FC >1 and 40 genes with log2 FC <-1 (Figure 2D). Nevertheless, GSEA revealed dysregulation of several processes related to basic cellular functions between non-inflamed CD and healthy controls possibly suggesting underlying differences in gene regulation (Figure 2E). Mitophagy (autophagy of mitochondria), as well as cellular response to vitamin D, were among the most negatively enriched GO terms. Conversely, mitochondrial translation elongation and termination and mitochondrial respiratory chain complex I assembly, were among the most positively enriched GO terms.

Ulcerative Colitis

Hierarchical clustering of genes differentially regulated in UC gave a reasonable separation of the inflamed and non-inflamed biopsies (Figure 3A). Differential expressed analysis between inflamed UC and non-inflamed UC, resulted in 17,704 DEGs (Table 2). Volcano plots showed a skewing of the data as 907 genes were strongly upregulated, while only 357 were strongly downregulated (log2 FC >±1) (Figure 3B). GSEA applied to the highly upregulated genes (log2 FC >±1), identified overrepresentation of processes related to inflammation and underrepresentation of biological processes related to mitochondrial respiration in inflamed UC compared to non-inflamed UC, largely overlapping the GSEA results in CD (Supplementary Figure 2B).

The comparison of inflamed UC and healthy controls showed 1088 genes with log2 FC >1 and 625 genes with log2 FC <-1 (Figure 3C). GSEA showed strong enrichment of processes related to proliferation and immune responses in inflamed UC. The most underrepresented GO terms in inflamed UC were dominated by vesicle transport and mitochondrial functions (Figure 3E).

Comparing non-inflamed UC and healthy controls, 111 highly dysregulated genes were found with 58 genes with log2FC >1 and 53 genes with Log2 FC <-1. Volcano plots showed that there was a slight skewing towards upregulated genes (Figure 3D). GSEA revealed dysregulation of several processes related to basic cellular functions between non-inflamed UC and healthy controls. Mitochondrial translation elongation and termination were among the most positively enriched GO terms (Figure 3E).

Figure 3 Continued.

Figure 3 Continued.

Figure 3 Ulcerative colitis. (A) Hierarchical clustering of differentially regulated genes. (B) Volcano plot comparing inflamed (UCi) and non-inflamed (UCni) samples in UC. (C) Volcano plot of inflamed biopsies from UC (UCi) compared to healthy controls. (D) Volcano plot of non-inflamed biopsies from UC (UCni) compared to healthy controls. (E) Gene set enrichment analyses given as Normalised Enrichment Scores (NES) in UCi (top panel) and UCni (bottom panel) versus healthy controls. “Orange” represents overrepresented GO terms, with positive NES. “Green” represents underrepresented GO terms, with negative NES. Abbreviations: UCi, inflamed UC; UCni, non-inflamed UC; HC, healthy controls; NES, Normalised Enrichment Score.

Ulcerative Colitis versus Crohn’s Disease

Comparing inflamed UC with inflamed CD yielded 204 highly (log2FC >±1) differentially expressed transcripts (146 upregulated and 58 downregulated, data not shown). The two gene expression signatures were highly correlated (r = 0.93), supporting that inflammation might mask underlying biological differences among the diagnosis groups. On a biological pathway level, the comparisons involving inflamed biopsies from UC and CD were remarkably similar with respect to normalized enrichment scores, regardless of diagnosis or whether the healthy or symptomatic controls were used in the comparison (Figure 4). In contrast, mitochondria-associated pathways had negative normalized enrichment score (NES) in inflamed UC compared to inflamed CD.

Figure 4 Gene set enrichment analyses. Ranking of the biological processes that are among the top 30 enriched in at least two datasets based on Normalized Enrichment Scores (NES). NES scores span from −2 to 3, where negative scores define underrepresented GO terms (green) and positive scores define overrepresented GO terms (Orange). Color intensity correlates with NES score. Significantly regulated biological processes are labelled with an asterisk (*). Abbreviations: CDi, inflamed CD; CDni, non-inflamed CD; UCi, inflamed UC; UCni, non-inflamed UC; SC, symptomatic controls; HC, healthy controls; NES, Normalised Enrichment Score.

Comparisons of non-inflamed tissue from UC and CD with healthy controls seemed to better reveal underlying biological differences between UC and CD (Figure 4). Non-inflamed UC tissue showed negative NES, or no enrichment, for several processes affecting immune response regulation. However, both showed enrichment of TNFα - and NF-κB signaling. Furthermore, the comparisons of non-inflamed UC or CD with healthy controls revealed shared positive normalized enrichment scores for processes related to cell proliferation, mRNA splicing/translation and mitochondrial respiration, and negative normalized enrichment scores for cell adhesion, extracellular matrix organization and collagen fibril organization.

Taken together, comparing IBD with healthy controls, revealed qualitative and quantitative differences in gene expression signatures. While mostly quantitative differences relating to immune cell activation could be seen in the inflamed biopsies, differences in gene expression signatures in non-inflamed samples suggested underlying biological differences between the diagnosis groups.

Symptomatic Controls

The unique signatures of non-inflamed CD, and non-inflamed CD compared to healthy controls, were not found in the comparison with symptomatic controls (Figure 4). When comparing non-inflamed CD with symptomatic controls, we found an underrepresentation of mitochondrial pathway, which is opposite to the overrepresentation of the same GO terms when compared with healthy controls. In non-inflamed UC, the comparison with symptomatic controls was less discriminative. Nevertheless, some GO terms were significantly different also in this comparison (Figure 4), despite no differentially expressed genes at a 5% false discovery rate (Table 2).

Treatment Escalation and Anti-TNF Response

The median follow-up period was 752 days. Forty-five IBD patients (30%) experienced treatment escalation to anti-TNF agents and/or surgery during the follow-up period. Characteristics and measurements of relevant biomarkers at inclusion to the study according to whether the patient escalated or not, are shown in Supplementary Table 1.

Single-Gene Based Prediction

Transcriptome-wide analysis at 0.05 FDR did not reveal any association between transcription levels of any single gene with need for treatment escalation, when controlled for country, gut segment, or both (data not shown).

Of the genes identified by Lee et al,⁸ as predictive of treatment escalation in CD, four (FOXO3, HLA-B, IGFBP1  and IGFBP3) were available on the microarray in this study (Supplementary Table 2). Both inflamed and non-inflamed samples from UC and CD were explored. IGFBP1 expression was found to be nominally significant in inflamed UC samples (p = 0.045), possibly predicting treatment escalation. No significant difference in expression level was found for the candidate transcripts FOXO3, HLA-B and IGFBP3 when comparing samples from escalators and non-escalators.

Forty-two IBD patients (30 CD, 12 UC) underwent treatment escalation to anti-TNF agents. We could not find any significant association between single transcripts and anti-TNF response among these patients (data not shown).

Prediction of Treatment Escalation by Gene Set Signatures

The analyses above illustrate the limited power of gene expression changes of single genes to predict functional impact, which is reasonable as compensation and balancing of single-gene expression changes occur at several levels in biological systems. Thus, we explored whether the distinct signatures revealed by hierarchical clustering of the DEGs in CD and UC (Figures 2A and 3A, respectively) were correlated with the need of treatment escalation. As the above analyses suggested that the inflammation might mask underlying differences in gene regulation between the diagnosis groups, we performed analyses of the non-inflamed and inflamed samples falling in the same clusters for CD and UC, separately.

In non-inflamed CD, none of the patients in cluster 2 experienced treatment escalation. In contrast, 16 out of 38, and 7 out of 19 patients escalated in cluster 1, and cluster 3, respectively (Supplementary Figure 3A). However, the difference in clinical outcome between clusters was not statistically significant (p = 0.175).

In non-inflamed UC, there were also no statistical significance (p = 0.832) when comparing gene expression signatures and their relation to treatment escalation (Supplementary Figure 3A). Further, the inflamed biopsies did not predict need for or time to treatment escalation, neither in CD nor in UC (Supplementary Figure 3B).

Prediction of Anti-TNF Response by Gene Set Signatures

We explored whether the distinct signatures revealed by hierarchical clustering of the DEGs in UC and CD (Figures 2A and 3A, respectively) were correlated with anti-TNF response. When analyzing non-inflamed and inflamed samples separately for UC and CD, no signatures were found to predict anti-TNF response or anti-TNF remission (data not shown).

Prediction of Treatment Escalation and Anti-TNF Response by Pathway Signatures

It has been proposed that dysregulation of signaling pathways, might better predict the changes likely to have functional impact. Thus, to identify novel differentially regulated pathways in the IBD patients, we used PathTracer which has been shown to reveal subtle changes in pathway activity in cancer data.³³ A PathScore heatmap (Supplementary Figure 3B) showed six clusters with different patterns of high and low dysregulation scores. An expected enrichment of pathways with high dysregulation scores of inflammatory responses, kRAS signaling, IFN-γ response and TNF-α signaling in inflamed samples of UC, CD and IBD-U from different locations were found (Supplementary Figure 4A). PathTracer scores did not separate UC and CD, but separated IBD from symptomatic and healthy controls.

On group level, the PathTracer clusters that separated high and low inflammatory response dysregulation scores (Supplementary Figure 4B) did not predict time for treatment escalation (Supplementary Figure 4C), although dysregulation scores of hypoxia and unfolded protein response pathways predicted need for treatment escalation in inflamed mucosa in UC (HR 1.2, p = 0.02, and HR 1.3, p = 0.01, respectively).

Cluster analyses based on pathway signature analyses revealed no association between the mucosal transcriptome and anti-TNF response, as the transcripts of patients exposed to anti-TNF therapy were scattered and there was no systematic clustering separating responders from non-responders (Supplementary Figure 4A).

Discussion

In accordance with earlier studies, location,³⁴ and inflammation status,^35,36 were important drivers of the transcriptome signatures separating the samples. In inflamed tissue, the immune responses appear to mask potentially important biological differences. Our study adds to the present knowledge by suggesting that the non-inflamed areas of the intestine, close to the inflammation, could have an even higher biomarker potential than inflamed tissue alone, supporting earlier studies demonstrating distinct gene expression profiles in non-inflamed tissue.^15,35,37,38

In the present study, the non-inflamed samples from both CD and UC revealed dysregulation of a large spectrum of biological pathways. Potent inflammatory pathways, as TNF-α signaling, were dysregulated. In non-inflamed CD, cellular response to vitamin D was among the top hits of negatively enriched GO terms, strengthening the knowledge base of vitamin D as a central factor in the pathogenesis of CD.³⁹ Of note, major disturbances in mitochondria-associated pathways were also found in non-inflamed biopsies from patients with CD and UC, compared to healthy controls, suggesting a widespread disturbance of intestinal energy regulation.

The molecular mechanistic background for mitochondrial dysfunction in the pathogenesis of IBD has been demonstrated both in CD and UC. Sünderhauf et al⁴⁰ recently showed how mitochondrial functions are essential for normal differentiation of mucin-producing Goblet cells, and how this mechanism was dysregulated in UC patients in remission, supporting the concept of UC as a barrier disease, and underscoring the value of investigating non-inflamed samples to avoid blurring of patterns caused by general inflammatory drive. The link between mitochondrial dysfunction and a weakened mucosal barrier has also been investigated in non-inflamed CD by Khaloian et al.⁴¹ Khaloian et al revealed mitochondrial dysfunctions to impact epithelial cell differentiation, leading to the development of aberrant Paneth cells with reduced barrier function properties.

In a publication by Haberman et al, suppression of mucosal mitochondrial functions was demonstrated in active UC.³⁸ In the present study, we present data that support dysregulation of mitochondrial functions also in non-inflamed UC samples. Furthermore, underrepresentation of GO terms related to mitochondrial function can also be seen in CD when comparing inflamed to non-inflamed CD samples and when comparing inflamed CD samples to healthy controls. Conversely, the non-inflamed samples in CD and UC showed enrichment of biological processes related to mitochondrial function compared to healthy controls.

Furthermore, when non-inflamed CD and UC samples were compared to healthy controls, negative enrichment scores for cell adhesion, extracellular matrix organization, and collagen fibril organization were found for both CD and UC, demonstrating dysregulation of pathways central for maintenance of robust gut barrier functions. Taken together, the present study supports a role of mitochondrial dysfunction and barrier dysfunction in the pathogenesis of both CD and UC.

In this study, we have used two well-defined control groups, both symptomatic and healthy controls. The symptomatic controls represented a large control group and were mainly (78%) represented by patients with irritable bowel syndrome (IBS). Our data demonstrated that the tissue responses among symptomatic controls resemble IBD on a transcript level. Furthermore, the comparisons between IBD and symptomatic controls presented pathway signatures that were less influenced by inflammatory status of the IBD-biopsy, than the comparisons with healthy controls. With respect to GO terms, we found differences in mitochondria-related pathways in the comparisons of IBD-samples with healthy controls and symptomatic controls, respectively. This finding might indicate similarities between IBS and IBD in context of a dysfunctional gut barrier. Interestingly, there have also been demonstrated major similarities in fecal microbiota composition between IBD and symptomatic controls in the IBD-character cohort, while the healthy controls separated to a greater extent.⁴² The diagnostic approach of IBS are symptom-driven, based on the Rome IV criteria.⁴³ IBS has been reported to affect approximately 10–20% of the population globally,^44,45 and severe morbidity, high economic costs and reduced quality of life have been reported among these patients.^46,47 The therapeutic armamentarium available for IBD patients has steadily expanded the last years. Unfortunately, this has not been the case for IBS. The unravelling of molecular patterns in symptomatic controls, could be a source to mine for biomarkers and novel therapeutic targets also for patients with irritable bowel syndrome (IBS). This is the first study with a large number of well-defined symptomatic controls, demonstrating the large mucosal transcriptional overlap with IBD in this group, compared to healthy individuals. These findings emphasize the importance of having strict criteria when defining a control group, securing homogeneity.

Development of personalized treatment strategies, based on molecular markers, is a central aim in IBD research.^14,48 However, at present, there are no reliable biomarkers to guide individually tailored treatment options. In the present paper, we presented transcriptomic profiling of newly diagnosed IBD patients with an aim of finding marker genes or marker clusters predicting treatment escalation, as well as response to anti-TNF agents.

The four genes identified by Lee et al,⁸ as predictive of prognosis in CD, were available on the microarray in this study (Supplementary Table 2). Among them, IGFBP1 expression was found to be associated with treatment escalation in inflamed UC samples (p = 0.045). IGFBP1 has been demonstrated to prolong the half-lives of insulin-like growth factors I and II; proteins involved in processes including immunity.⁴⁹ A previous study also demonstrated a connection between IGFBP1 loci and prognosis in rheumatoid arthritis, another immune-mediated disease.⁵⁰ Taken together, these results are suggestive of a common involvement of the IGFBP1 locus in IBD and rheumatoid arthritis, with a substantial impact on prognosis.

We further explored if gene expression signatures could predict treatment escalation without making significant findings. Nevertheless, in non-inflamed CD, the hierarchical cluster analyses might be suggestive of a potential to predict clinical outcome as none of the patients in gene set cluster 2 experienced treatment escalation. Noteworthy, the most discriminative transcripts, generating the clustering were DEFB4A, encoding antimicrobial peptides, and PCK1, encoding the metabolic enzyme phosphoenolpyruvate carboxykinase, an integral in energy metabolism. Studies have already revealed PCK1 to be involved in inflammation. Ko et al⁵¹ found that a deletion of the PCK1 gene gave a more proinflammatory phenotype when they explored metabolism of macrophages. We find it interesting that expression of transcripts involved in mitochondrial functions and barrier functions are dominating the clustering, even though the comparison of the signatures regarding their inclination to treatment escalation did not reach statistical significance.

The role of the transcriptome as a predictor of anti-TNF response has been extensively investigated. Some studies find signatures,^21,23 others do not.²⁰ In a multi-omic study of treatment naïve pediatric IBD patients, distinct transcriptome signatures at diagnosis were found to predict treatment escalation to anti-TNF-therapy, and anti-TNF response.³⁷ Kugathasan et al, also studying a pediatric CD cohort, found gene expression signatures to be related to clinical outcome in an integrative model.¹⁰ In our study, no single gene association could be found to predict treatment response to anti-TNF agents. Furthermore, we did not identify gene set signatures or pathway signatures predicting anti-TNF response in CD or UC.

A strength regarding the anti-TNF response evaluation in this study, is the use of objective markers as f-calprotectin and CRP to assess response, more robust than symptom-based scoring alone (clinical activity indices as HBI and partial Mayo score). Nevertheless, there might be a potential to include more stringent criteria, defining remission, rather than response, to identify predictive gene expression signatures, as implemented in the PANTS study by Kennedy et al.⁵²

The gut transcriptome is highly influenced by inflammation status,^35,36 bowel segment³⁴ and medical treatment.^21,25,53,54 A strength of this study is, therefore, the use of a standardized protocol for biopsy sampling during ileocolonoscopy in this treatment naïve cohort. Strong influence on gene expression by active inflammation could mask underlying differences among the diagnosis groups, reinforcing the importance of a standardized protocol for biopsy sampling.^55,56 The sampling of both inflamed and non-inflamed samples made it possible to look for precursor molecular changes in the ileocolonic mucosa close to segments with active inflammation. The patients were newly diagnosed, mainly treatment naïve (>92%), and the impact of a longstanding inflammation on the transcriptome, could possibly be less dominant in this cohort.

A limitation in the present study is the lack of application of endoscopic severity scores such as Simple endoscopic score for CD and endoscopic Mayo score, as well as histopathological evaluation, to grade inflammation in the mucosal biopsies, as earlier studies have demonstrated these factors influence on the gene expression profiles.^38,57 This could restrict the opportunity to classify molecular phenotypes and relate these molecular phenotypes to endoscopic and histological phenotype, probably playing an important role when performing treatment response analyses.

In conclusion, gene expression signatures of the non-inflamed tissue in treatment naïve CD and UC reveal dysregulation of biological pathways that might potentially be central in disease initiation and prognosis. In particular, major disturbances in mitochondria-associated pathways and barrier dysfunction were confirmed. The gene expression signatures of the control groups were related to if they were symptomatic or not, which may have important implications for future study designs. Larger studies of treatment-naïve IBD with standardized sampling and prospective follow-up of clinical outcome are relevant to identify gene expression signatures applicable for supporting clinical decision-making.