Characterization of the Regulatory Landscape in Crohn’s Disease Reveals microRNA-Associated Alterations that Shape Anti-TNF Response

Abstract

Background

MicroRNAs (miRNAs) play a key role in regulating gene expression in Crohn’s disease (CD). Although several studies have identified miRNAs with biomarker potential, an exhaustive characterization of the miRNAome in CD is still lacking. We performed the largest miRNA profiling effort to date to analyze miRNA variability across intestinal tissues, disease activity status, and infliximab treatment in CD.

Methods

We generated 119 transcriptomic profiles from the terminal ileum and left colon biopsies of 30 individuals (10 with active CD, 10 with quiescent CD, and 10 healthy controls). Half of the samples were cultured ex vivo with infliximab, and the remaining half with basal medium. Using variance analyses and linear mixed differential expression models, we explored the determinants of miRNAome variability in CD. We also generated infliximab response signatures to identify candidates and examine interactions between miRNAs and the coding transcriptome.

Results

Tissue location, and patient-specific effects, were the main factors in miRNA variability in CD, with some differentially expressed miRNAs involved in epithelial–mesenchymal transition (miR-200s, miR-429). We identified 9 miRNAs with treatment-responsive behaviour, particularly to the terminal ileum of active CD cases. Although the changes observed in active CD cases suggest that many alterations are not offset by infliximab incubation, we described 13 miRNAs–mRNA pairs with potential involvement in the anti-tumor necrosis factor (TNF) treatment, 7 of which have been already validated.

Conclusions

A comprehensive miRNA profiling revealed significant intestinal tissue-specific variability and identified key alterations in the miRNA–mRNA interactome that might be involved in therapeutic response to anti-TNF in CD.

Graphical Abstract

Lay Summary

This article provides an extensive characterization of the miRNA variability in Crohn’s Disease biopsies incubated ex vivo with infliximab, revealing tissue-specific regulatory programs and identifying miRNA interactions associated with anti-TNF response, with potential therapeutic applications.

anti-TNF drug levels, transcriptomics, microRNA, Crohn’s disease, gene regulation

Key Messages

• What is already known: The number of works addressing the full characterization of the miRNA regulatory mechanisms in Crohn’s Disease, including response to anti-TNF, remains limited.

• Findings: The comprehensive characterization of the miRNA profiles determined that tissue effects were predominant, more than disease subtype effects while identifying 9 infliximab-responsive miRNAs and 7 validated interacting pairs of miRNAs and transcript targets with a potential therapeutic role.

• Implications for patient care: This study provides new venues for the development of therapeutic targets, which could improve disease management and health outcomes for Crohn’s Disease patients.

Introduction

Crohn’s disease (CD) is a chronic immune-mediated inflammatory disease characterized by the affection of the gastrointestinal tract in relapse-remission cycles.¹ The development of high-throughput molecular profiling technologies has revolutionized our ability to characterize the pathogenic processes involved in CD. For instance, transcriptomic profiling provides a readout of the cell regulatory activity to identify gene expression changes associated with disease processes and treatment responses.^2,3 However, the translation of these findings into clinical actionability is difficult. To date, fecal calprotectin and C-reactive protein are the only biomarkers routinely used for disease monitoring in the clinic.⁴

MiRNAs are non-coding RNA molecules, typically 19-25 nucleotides long, that repress protein expression by direct inhibition or cleavage of the targets. The primary miRNA transcripts (pri-miRNA) are transcribed in the nucleus and then cleaved by the microprocessor complex into the precursor miRNA (pre-miRNA) in the canonical biogenesis pathway.^5,6 They are then exported to the cytosol, where mature miRNAs are formed by Dicer-mediated cleavage and loaded to an AGO protein to form the RNA-induced silencing complex.

The study of the miRNAome provides insights into the gene expression alterations observed in CD patients. Over the past decade, several studies have attempted to identify miRNAs with biomarker potential by profiling both tissues directly involved with CD activity as well as from less invasive samples, such as circulating blood and stools.⁷ These efforts have crystallized into individual miRNAs that could help to track processes associated with CD, including inflammation, and disease progression.^7,8 Besides the potential as biomarkers, the miRNAome can provide deep insights into the mechanistic molecular aspects that shape transcriptomic regulation in this disease (eg, regulation of intestinal permeability and autophagy).^9,10 While many studies have focused on individual miRNAs, an exhaustive characterization of the miRNAome and its determinants in CD is still missing.

The drug portfolio approved for CD has improved constantly in the last decade. Besides corticosteroids and immunosuppressors, an increasing number of biological treatments, such as anti-cytokine agents like anti-tumor necrosis factor (TNF) (infliximab, adalimumab, certolizumab, golimumab), and ustekinumab; and anti-integrin molecules (vedolizumab) are available for gastroenterologists.^11,12 However, regardless of the treatment of choice, a relatively high proportion of patients fail to achieve long-term remission.¹³ Understanding the molecular aspects that determine drug response is a particular pressing need in CD.¹⁴ Several reports have focused on identifying miRNAs that can predict therapeutic response.^15,16 However, most of these studies have followed a “candidate-gene” approach, which focuses on short lists of miRNAs selected according to previous knowledge. This strategy misses the discovery opportunities provided by a hypothesis-free approach based on profiling the complete miRNAome. Moreover, these studies usually focus on specific patient settings (eg, a single tissue type). Characterizing the broader context and variability across diverse tissues and patient settings is hence a must to dissect the regulatory effects of the miRNAome that are most relevant for drug response.

To fill these gaps, in our study, we carried out a comprehensive evaluation of the miRNAome to characterize the sources of variation that are most relevant in CD, with a particular focus on detecting differences associated with intestinal tissue location, disease activity, and treatment with infliximab.

Materials and Methods

Experimental Design and Patient Recruitment

We collected gut biopsies from 3 groups of patients: 10 individuals with active CD (aCD), 10 with quiescent CD (qCD), and 10 who were healthy controls (HC). For each patient, we collected 4 samples: 2 ascertained from the terminal ileum and 2 from the left colon, for a total of 120 specimens in the overall study (119 analyzed). With these biopsies, we carried out an ex vivo assay by splitting them into 2 groups. Half of them were incubated with infliximab, while the other half was cultivated with basal medium only (see “Ex Vivo Incubation with Infliximab”). The experimental design aimed to balance across the 3 variables of interest, namely tissue, disease activity, and treatment, hence achieving 10 samples per each of the 12 possible combinations. Figure 1 is a diagram summarizing the experimental design. Sample collection was carried out at the Gastroenterology Department of Hospital Universitario de La Princesa (Madrid, Spain). This ex vivo experimental setting with the same biopsy samples was already harnessed to characterize the regulatory role of long non-coding RNAs (lncRNAs) in a previous work from our group.¹⁷

Figure 1.

Diagram outlining the experimental design.

Crohn’s disease patients were classified as active or quiescent according to the Simple Endoscopic Score for Crohn’s Disease (SES-CD), with this score < 3 for qCD, and ≥3 for aCD. Eligible CD patients were detected through a colonoscopy, having inflamed endoscopic and mucosal states, while HC were described as macroscopically and histologically normal. The latter group is referred from colorectal cancer screenings, or changes in bowel transit studies, without having a diagnosis for any inflammatory or autoimmune disease. Exclusion factors that were considered for patient enrollment were age (<18), presence of significant pathologies, changed medication in the last 3 months, alcohol or drug abuse, and pregnancy and lactation status. Crohn’s disease patients were also required not to have received any biological or drug altering the immune system significantly. The main clinical and demographic records of patients are summarized in Table 1.

Table 1.

Demographic and clinical characteristics of the study groups. In brackets, percentage of samples per category.

Characteristics	HC (n = 10)	qCD (n = 10)	aCD (n = 10)
Female gender, n (%)	8 (80)	6 (60)	5 (50)
European ethnicity, n (%)	10 (100)	10 (100)	10 (100)
Non-smokers at diagnosis, n (%)	8 (80)	8 (80)	4 (40)
Median age at diagnosis, years (IQR)	N/A	29 (22-37)	27 (22-49)
Median age at enrollment, years (IQR)	51 (45-55)	44 (34-48)	28 (24-57)
CD characteristics
Location
L1: ileal	N/A	4 (40)	3 (30)
L2: colic	N/A	4 (40)	3 (30)
L3: ileocolic	N/A	2 (20)	3 (30)
L4: upper CD	N/A	0 (0)	1 (10)*
Behavior
B1: inflammatory	N/A	9 (90)	9 (90)
B2: stricturing	N/A	1 (10)	1 (10)
B3: penetrating	N/A	0 (0)	0 (0)
Perianal disease	N/A	5 (50)	2 (20)
Treatment at sample collection
None, n (%)	N/A	1 (10)	4 (40)
Aminosalicylates, n (%)	N/A	2 (20)	1 (10)
Immunomodulators, n (%)	N/A	6 (60)	4 (40)
Biologics, n (%)	N/A	4 (40)	4 (40)
Anti-TNF-α agents	N/A	4 (40)	1 (10)
1 anti-TNF-α	N/A	4 (40)	1 (10)
2 anti-TNF-α	N/A	0 (0)	0 (0)
Anti-integrin	N/A	0 (0)	2 (20)
Anti-IL-12/23	N/A	0 (0)	1 (10)

HC: healthy control; qCD: quiescent Crohn’s disease; aCD: active Crohn’s disease; CD: Crohn’s disease; N/A: not applicable; TNF: tumor necrosis factor; *Ileocolic and upper CD; IQR, interquartile range.

Ex Vivo Incubation with Infliximab

Sample collection was carried out in sterile vials and stored in ice-cold complete medium (composed of RPMI 1640 from Sigma Aldrich, 100 μg/mL penicillin/streptomycin, 2 mM L-glutamine, 50 μg/mL gentamicin, and 10% fetal bovine serum). Biopsies were then washed in HBSS buffer and cultured in complete medium at 37 °C for 18 h, either in the presence of infliximab (4 mL complete medium + 4 µL infliximab stock at 10 mg/mL), or just in a control medium (4 mL complete medium). After culturing, all samples were introduced in 250 μL RNAlater solution and stored at –80 °C.

Sample Preparation and RNA Sequencing

For RNA isolation, intestinal tissues were suspended in TRIZOL Reagent (Tri Reagent; Invitrogen, Waltham, MA, USA) and chloroform (Panreac). The mixture was then homogenized and processed using the RNeasy Mini Kit 250, following the manufacturer’s instructions. To evaluate the quantity and quality of miRNAs and mRNAs, we used Qubit RNA HS Assay Kit (Thermo Fisher Scientific, Cat. #Q32855) and Agilent RNA 6000 Nano Chips (Agilent Technologies, Cat. #5067-1511), respectively. Sequencing libraries for miRNAs were prepared following the protocol included with the kit “NEXTflexTM Small RNA-Seq Kit v3” (©Bioo Scientific Corp. Catalog #5132-05). Briefly, between 20 and 100 ng of total RNA was incubated for 2 min at 70 °C, then 3’ 4N-adenylated adapter and ligase enzyme were added, and ligation was conducted by incubation of this mix overnight at 20 °C. After excess 3’ adapter removal, 5’ adapter was added alongside with ligase enzyme, and the mix was incubated at 20 °C for 1 h. The ligation product was used for the reverse transcription with the M-MuLV Reverse Transcriptase in a thermocycler for 30 min at 42 °C and 10 min at 90 °C. Next, enrichment of the cDNA was performed using PCR cycling: 2 min at 95 °C; 19-22 cycles of 20 s at 95 °C, 30 s at 60 °C, and 15 s at 72 °C; a final elongation of 2 min at 72° and pause at 4 °C. PCR products were resolved on 8% Novex TBE PAGE gels (Cat. #EC6265BOX, Thermo Fisher Scientific), and fragments of mature miRNA ~150 bp cut from the gel. After extraction of small RNAs, and a cleaning up procedure, libraries were visualized on an Agilent 2100 Bioanalyzer using Agilent High Sensitivity DNA Kit (Agilent Technologies, Cat. #5067-4626) and quantified using Qubit dsDNA HS DNA Kit (Thermo Fisher Scientific, Cat. #Q32854).

On the other hand, total RNA (mRNA plus lncRNA) sequencing libraries were prepared using “Illumina Stranded Total RNA Prep, Ligation with Ribo-Zero Plus” kit (Illumina, Inc. Part #20040525) following the “Illumina Stranded total RNA Prep Ligation with Ribo-Zero Plus Reference Guide.” Starting from 58-300 ng of total RNA, rRNA was depleted, and remaining RNA was purified, fragmented, and primed for cDNA synthesis. cDNA first strand synthesis was carried out for 10 min at 25 °C, 15 min at 42 °C, 15 min at 70 °C and pause at 4 °C and cDNA second strand was synthesized at 16 °C for 1 h. Following A-tailing, pre-index anchors were ligated to the ends of the double-stranded cDNA fragments to prepare them for dual indexing. A subsequent PCR amplification step to add the index adapter sequences (30 s at 98 °C; 12-15 cycles of 10 s at 98 °C, 30 s at 60 °C, 30 s at 72 °C; 5 min at 72 °C and pause at 4 °C) was performed. After a final library clean-up, libraries were quantified using Qubit ds DNA HS DNA Kit (Thermo Fisher Scientific) and visualized on an Agilent 2100 Bioanalyzer using Agilent High Sensitivity DNA kit (Agilent Technologies).

Exploratory Data Analyses

To generate the miRNA-enriched libraries, FASTQs were trimmed following the recommendations of the NEXTflex Small RNA-Seq Kit manufacturers. We used Bowtie¹⁸ to align the reads against the human genome GRCh38, allowing no mismatches. We chose miRBase v22 to quantify the mature miRNAs, harnessing the Partek Flow application (version 7.0). To generate mRNA libraries, raw reads were trimmed with skewer followed by quality control analysis via FastQC. Trimmed reads were aligned to the reference CRCh38 genome via STAR v2.7.1 aligner indexed to GENCODE v26 GTF library using -quantMode GeneCounts for feature counting. We obtained 2 expression matrices with the counts quantified for each feature. In total, we obtained 761 miRNAs in the first matrix and 58 278 features in the whole RNA dataset, among which 18 260 are protein-coding (mRNA) genes. We applied a threshold of 20 counts across all samples per feature, leading to the inclusion of 603 miRNAs and 16 555 protein-coding mRNAs in the final study.

The structure of expression and the role of the relevant clinical variables were explored using Principal Component Analysis (PCA), with the “prcomp” R function, using rlog-transformed counts. We also performed Principal Variance Component Analysis (PVCA) to obtain the proportion of variance explained by each variable. Expression data were normalized according to “voom” method (normalized by log2 of counts per million).

Variance Analysis and Differential Expression Models

All statistical analyses were carried out using the R programming language (version 4.2.3). Differential Expression Analysis was performed with linear mixed models implemented in the “dream” package. This approach permits to delve into the transcriptional changes while accounting for the effects of individuals as random effects.¹⁹ We created a general additive model, and another one with interactions including specific effects of tissue, disease activity, and treatment, to detect differentially expressed miRNAs (DE-miRNAs) and mRNAs (DEG). The False Discovery Rate (FDR) threshold for significance was set up at 5%, using an absolute fold change of 2. For treatment response, we used a more liberal FDR threshold of 10% to increase the number of features detected. We plotted heatmaps to represent the voom-normalized counts of the top 5% miRNAs ranked according to Variance Partition (VP) and to the Differential Expression model (DE-miRNAs with the lowest adjusted P-value). Clustering of samples was performed using the “ward.D2” method.

For pathway enrichment analysis, we functionally annotated the selected miRNAs with miRPathDB.²⁰ To incorporate biologically meaningful miRNAs, we used a compilation of summary statistics of pathway enrichment results of many human miRNAs, including only miRNAs whose function has been validated in miRTarBase.²¹ Over-representation analysis of enriched miRNAs with the enricher function available in the “clusterProfiler” package,²² while Gene Set Enrichment Analysis (GSEA) results were represented with ridgeplots.

Compilation of miRNA–mRNA Pairs to Study the Regulatory Interactome in CD

To characterize the functional impact of the miRNAs available in our study, we applied a computational pipeline to obtain bona fide miRNA–mRNA pairs. MiRNA and mRNA sequences were downloaded from miRbase and the ENSEMBL site, respectively. We used the aligner Bowtie to obtain all possible pairwise alignments between the miRNAs and mRNAs available in the normalized expression datasets.¹⁸ No mismatches were allowed (-v 0). From a total of 1527 candidates, we selected 517 miRNA–mRNA unique pairs. To gain insights into the changes in the miRNA–mRNA interactome, we calculated Pearson correlation coefficients across the full dataset and for each combination under study (119 and 10 samples each, respectively). The summary statistics for all miRNA–mRNA pairs in each combination are available in Table S1. To overcome the low statistical power for each combination, we explored the changes in the sign of the correlation coefficients according to disease activity and treatment status. Statistical significance across these analyses was obtained through Fisher exact tests.

Infliximab Response Signatures to Identify miRNAs Involved in Drug Response

To focus on miRNAs with potential role in infliximab response, we created 5 different candidate lists or signatures including dysregulated miRNAs in the presence of infliximab.

Four out of 5 signatures were constructed by harnessing the study by Arijs et al. as a reference dataset (GSE16879).²³ This longitudinal study of patients treated with infliximab included responders and non-responders, with 2 sampling timepoints (before and after treatment), and 2 tissues (ileum and colon). Although Ulcerative Colitis (UC) cases and HC were available, we only used the CD samples in this study. The first signature was created by identifying miRNAs listed in the miRbase database that were differentially expressed using “dream” linear mixed models, for a total of 678 miRNAs (Table S2). For the second signature, we integrated pathways from the Reactome database that are enriched in genes matched with miRNAs on a sequence-based basis, indicating potential regulation by miRNAs. We used the “hgu133plus2probe” package to retrieve affy probe sequences from the GSE16879 reference panel. After obtaining the cDNA sequences from the gene site and translating them with the “biomaRt” R tool, we aligned them to miRbase miRNA sequences to map miRNAs to all genes with the previous configuration of Bowtie, for a total of 464 pathways (Table S3). For the third signature, we computed a VP model to identify the top 10% of miRNAs from the reference study that explain the highest proportion of variance for the following variables: tissue (ileum vs. colon), week (after vs. before treatment), response (response vs. non-response), and subtype (CD vs. control), for a total of 93 miRNAs (Table S4). The fourth signature was constructed using a list of pathways enriched in miRTarBase-validated miRNAs according to their role in TNF regulation, which we deemed as potential hallmarks of infliximab response, for a total of 45 pathways (Table S5).

Finally, to compile the fifth signature we leveraged a recent set of abstracts from over 25 000 articles focused on inflammatory bowel disease (IBD).²⁴ We used the keywords “microRNA,” “miRNA,” and “miR” to perform a broad search in the abstract section of each article and select 8 articles dealing with miRNAs (Table S6). Finally, for each molecule in each of the 5 signatures, we calculated the pairwise correlation coefficient in each of the 12 combinations and ascertained candidates as explained in the Results section.

Results

To characterize the miRNAome in CD and the interplay with the changes observed in the coding transcriptome, we generated 119 RNA-Seq profiles from gut biopsies cultivated ex vivo (one sample excluded). We prepared 2 libraries for each sample, 1 enriched for regulatory miRNAs and the other for coding mRNAs. Our analyses focused on understanding the influence of clinical and demographic factors on the variability of the transcriptomic profiles, with particular emphasis on tissue location, disease activity, and the impact of incubation with infliximab. The cohort comprises 10 CD patients with active endoscopy, 10 qCD, and 10 HCs. For each individual, we generated 4 profiles to ensure a balanced representation across tissues (left colon and terminal ileum) and according to incubation with infliximab (experimental design diagram in Figure 1). After data curation, the profiles for each sample included 603 miRNAs and 16 555 mRNAs.

Intestinal Tissue Location is the Main Determinant Driving Expression Variability

We first performed a PCA on each of the 2 transcriptomic datasets. PC1 accounted for 20.69% and 42.64% of variance in the miRNA and mRNA datasets, respectively. In both cases, the samples were distributed across PC1 according to tissue (Figure 2A and B). PC2 explained a smaller fraction of variance (9.6% for miRNAs and 16.93% for mRNAs), without a clear correlation according to tissue. None of the other 2 main variables of interest, namely disease activity and infliximab status, showed distinct separation along PC1 or PC2, although the most active CD cases were distributed along the negative component of PC2 for both miRNA and mRNA datasets (Figure S1).

Figure 2.

Variance exploratory analyses of the miRNA and mRNA expression datasets. (Left) The principal component analysis spaces of the miRNA (A) and mRNA (B) expression matrices are represented on the left. In both cases, tissue is the main factor driving expression variability along PC1 (20.69% and 42.64%, respectively). (Right) Principal variance component analysis estimation bar chart for miRNA and mRNA datasets (C). The tissue effect carries most of the average variance in both models (66%, and 7.4%, respectively), although the residual effect is particularly large in the miRNA dataset (74.3% of the variance). Interactions between the 3 variables of interest are also included. L: left colon; resid: residual values, unaccounted by tissue, disease activity, treatment status and interaction effects.; T: terminal ileum.

We next conducted PVCA to estimate the overall contributions of the 3 abovementioned variables. For mRNA variability, tissue arose as the main determinant followed by the tissue-by-patient interaction term (66% and 20.3%, respectively, Figure 2C). In contrast, the contribution of tissue was notably lower for the variability in the miRNAome (7.4%), which showed an overwhelmingly dominant residual term (74.3%). This suggests that the 3 variables of interest are not the main drivers of miRNAome variability.

Remarkably, for miRNAs the interaction terms for tissue-by-patient and treatment-by-patient accounted for 6.8% and 3.5% of the variance, respectively (Figure 2C). This hints at variability that stratifies independently for each patient. In both datasets, disease activity played a minor role (2.4% and 1.7% for mRNA and miRNA, respectively). The main effect of incubation with infliximab was particularly small, implying that stratification of the transcriptome according to this variable is negligible and inconsistent across conditions.

Variable miRNAs Reveal Heterogeneity in Co-regulation Across Tissues

We used a VP model to identify subsets of miRNAs (n = 30, top ~5%) whose variability is predominantly influenced by each of the 3 studied variables. Tissue consistently explained an average of 33.2% of the variability among the top regulators selected for this variable (Figure 3A). The figures for the top 30 miRNAs selected based on disease activity and treatment were particularly lower (9.5% and 4%, respectively, Figure S2). The corresponding PCAs for the 3 sets of miRNAs confirmed these patterns (Figure S3). Notably, the 3 groups of miRNAs showed minimal overlap, with at most 2 miRNAs selected in more than 1 set.

Figure 3.

Differences in the miRNA regulatory landscape between terminal ileum and left colon. (A) Fraction of variance explained by each clinical factor for the top 30 (5%) miRNAs for tissue. (B) Heatmap of normalized counts of the top 30 miRNAs for tissue. (C) Heatmap with Spearman’s pairwise correlation coefficients of the top 30 miRNAs detected in the left colon. (D) Heatmap with Spearman’s pairwise correlation coefficients of the top 30 miRNAs detected in the terminal ileum. aCD: active CD, HC: healthy controls, IFX: infliximab; L: left colon; qCD: quiescent CD; T: terminal ileum.

A heatmap based on the top 30 miRNAs selected for tissue showed strong consistency across samples (Figure 3B). An enrichment analysis using miRPathDB Reactome pathways determined that the 15 miRNAs up-regulated in the left colon are associated with cellular senescence, stress, and external stimuli response, while the 15 up-regulated miRNAs in the terminal ileum are linked to transcription and replication processes (Table S7).

To explore the landscapes of co-regulation, we calculated pairwise Spearman correlation values among all selected miRNAs in each of the 2 tissues. In the left colon, we observed pervasive high positive correlation within each cluster (Figure 3C). Despite some miRNAs belonging to the same precursor family (eg, miRNA-200s), some of the most strongly correlated pairs do not belong to the same family, including miR-383-5p with miR-4525 (ρ = 0.99) and miR-212-3p with miR-4525 and miR-383-5p (ρ = 0.84, and ρ = 0.82, respectively). This suggests shared regulation across up-regulated miRNAs in the left colon. Conversely, in the terminal ileum, we observed less clustering (Figure 3D). The weaker correlation among members of the same family, or functionally similar families, hints at multiple regulatory programs shaping the variability of the miRNAome in this tissue. MiR-328-3p and miR-615-3p consistently exhibited negative correlations with the rest of the miRNAs, suggesting potential antagonistic repressor functions that counteract other miRNA-controlled processes.

MiRNAs are Involved in Different Biological Functions Across Tissues

To delve deeper into the biological characterization of miRNA variability across tissues, we performed a differential expression analysis, fitting an additive linear mixed model to control for the impact of patients as a random effect. We identified 109 DE-miRNAs between the 2 locations (Figure 4A). The majority (88.1%) of the miRNAs were up-regulated in the terminal ileum, including miRNAs associated with epithelial to mesenchymal transition (eg, 196a, 196b, 378a, 200b, 00c, and 203a).^25,26 The large overlap with the top miRNAs selected in the VP model (23 out of 30, 76.7%) validates the robustness of both methodologies.

Figure 4.

Biological annotation of the differentially expressed miRNAs (DE-miRNA) detected in the terminal ileum vs. left colon comparison. (A) Heatmap of the normalized counts of the top 50 DE-miRNAs ranked according the P-value. (B) Bar plot of Reactome pathways available in miRPathDB enriched among the 96 detected DE-miRNAs. (C) Dotplot depicting the gene ratios for the same pathways. (D) Ridgeplot depicting the log of the fold-changes for the miRNAs included in pathway analysis based on gene set enrichment analysis. Positive values (logFC > 0) indicate miRNAs that are up-regulated in the terminal ileum. aCD: active CD, HC: healthy controls, IFX: infliximab; L: left colon; qCD: quiescent CD; T: terminal ileum.

To further characterize their biological role, we performed a general enrichment using the miRTarBase-validated Reactome pathways available in miRPathDB. Several pathways, including “Cellular responses to external stimuli,” “Cellular responses to stress,” and “Developmental Biology,” showed enrichment in DE-miRNAs (Figure 4B and C). A GSEA, that considers the entire distribution of expression changes for each pathway, also revealed translation-associated pathways and DNA repair processes up-regulated in the terminal ileum, while signaling activity was down-regulated (Figure 4D).

Association of miRNA Expression Patterns with CD Severity

The evidence for expression variability apportioning according to disease activity was lower. We observed 154 (25.5%) and 203 (33.7%) miRNAs with intermediate expression in qCD in the terminal ileum and left colon, respectively. The corresponding values for the other comparisons, namely HC and aCD with intermediate expression averages, ranked from 25.2% to 49.3%. These figures suggest that the average expression of miRNAs do not correlate clearly according to the stages of disease progression.

In both tissue locations, 22 miRNAs showed increasing expression toward aCD patients. In contrast, 16 miRNAs showed progressively decreasing expression levels from aCD to HC (Figure S4). Enrichment analyses revealed pathways related to the cell cycle, interleukins, and signaling, suggesting that these miRNAs play distinct roles in the immune response in CD patients depending on the modulation by disease severity or inflammatory activity.

Incubation with Infliximab Alters miRNA Levels in the Terminal Ileum of Active CD Cases

To unravel the changes in the miRNAome in response to incubation with infliximab, we employed a linear mixed model with tissue, disease activity, and treatment status as fixed effects, along with their interactions. This approach enabled us to identify differences between treated and non-treated samples following 18 h of infliximab incubation for a specific tissue and disease activity category. Notably, we only identified significant changes in terminal ileum samples from active CD cases, underscoring the high specificity of treatment effects on miRNA-mediated regulation.

The differential expression model identified 9 miRNAs significantly up-regulated (10% FDR) upon infliximab incubation (Figure 5A). This included 378a-3p, 194-5p, and 200b-5p, which are involved in cellular responses to stimuli, stress, and general signaling regulation in the miRPathDB pathways analysis (Figure 5B). However, as discussed above, the variance explained by treatment status across the miRNAome was low, leading to unclear clustering analysis (Figure S2B), and several enriched pathways associated with signaling and gene expression that include a low number of significant features (not shown).

Figure 5.

MiRNAs detected as differentially expressed in response to infliximab in terminal ileum biopsies from active CD cases. (A) Heatmap of normalized counts of the 9 miRNAs detected. (B) List of miRPathDB pathways for which the 9 DE-miRNAs are annotated. aCD: active CD, IFX: infliximab; T: terminal ileum.

Infliximab Fails to Restore the miRNA–mRNA Interactome in the Terminal Ileum of Active CD

To further characterize the impact of infliximab incubation in the regulation of the miRNAome and its overall regulatory impact upon the coding transcriptome, we capitalized on the availability of mRNA profiles generated for the same specimens in each patient. In homeostatic conditions, miRNAs typically negatively regulate the levels of target mRNAs. We hypothesized that the chronic inflammation associated with CD might disrupt the canonical negative correlation between miRNA and mRNA levels. In turn, following incubation with infliximab, this relationship could normalize towards enhanced correlation of their levels of expression.

To evaluate these hypotheses, we first obtained a list of genome-wide pairs of miRNA–mRNA that can interact according to their sequence. Then, to focus on pairs potentially involved in infliximab response, we leveraged a longitudinal study of IBD patients treated with infliximab and a public repository of IBD-associated omic literature. In total, we selected 517 pairs of interest for which miRNA and mRNA counts were available in our normalized expression datasets (see “Methods” for details on the overall strategy). Among these pairs, 83 (16.1%) showed a nominally significant negative correlation across the 119 specimens, with a further 17 (3.3%) reaching the Bonferroni threshold of significance (Table S1). In contrast, when considering the whole universe of all possible miRNA–mRNA pairs (266 772 pairs), only 11.3% and 0.9% showed a nominally significant negative correlation and met the Bonferroni threshold, respectively. Besides confirming that miRNAs and mRNAs are globally coregulated, these proportions imply that a notable fraction of the selected miRNAs impact their predicted mRNA targets across the different combinations in our dataset.

By exploring the 517 pairs in each combination (summary statistics available in Table S1), we gained several insights into the miRNA interactome in CD. First, we compared across the basal combinations according to disease subtype. The number of miRNA–mRNA pairs with negative correlation was larger in the samples from HC than in those from CD (averages of 47.7% and 42%, respectively, P = .0031, Fisher exact test). Notably, we did not observe differences between aCD and qCD (41.7% and 42.4%, respectively, P = .79, Fisher exact test).

Next, we explored the effect of infliximab incubation onto these interactions. For each pair, we calculated the correlation between miRNA and mRNA levels in each of the 12 combinations stemming from our experimental design (Figure S5). An interesting pattern emerged when classifying combinations according to treatment. Across the 6 combinations of variables involving incubation with infliximab, we observed 1649 and 1453 positively and negatively correlated pairs, respectively. The same figures for the 6 corresponding basal combinations were 1740 and 1362, respectively. Hence, incubation with infliximab was associated with an increase in the number of miRNA–RNA with negative correlation estimates (43.9%-46.8%, ∆2.9%, P = .021, Fisher’s exact test), implying recovery of the canonical expectation.

However, these analyses uncovered an unexpected pattern that appeared restricted to the samples from the terminal ileum from active CD cases. Rather than increasing like for the other combinations, the number of miRNA–mRNA pairs moving towards negative correlation estimates decreased by 3.8% (220 vs. 200 out of 517 pairs). This divergent behavior became more evident when looking at miRNA–mRNA pairs with largest differences between the estimates of correlation in the basal and incubated with infliximab samples (ie, |r_INFL – r_BASAL| > 0.3). For the samples from the terminal ileum from active CD cases, we observed 265 pairs with a large difference in correlation estimates of the basal vs. infliximab samples. Of those, 159 (60%) corresponded to pairs moving strongly towards a positive correlation in the abundance of miRNA and mRNA in the incubated group. In contrast, the opposite trend was observed for all 5 remaining pairwise comparisons between basal and infliximab samples. Specifically, more miRNA–mRNA pairs with large changes would move towards negative correlation (an average of 141.2, 55.4%) than towards positive correlation (an average of 113.8, 44.6%). This includes the left colon samples from active CD cases, that did not show the divergent behavior observed for their terminal ileum counterparts.

Identification of miRNA Signatures as Potential Biomarkers for Infliximab Response

We developed a strategy to identify miRNAs with relevant participation in the miRNA–mRNA interactome and that can be associated with infliximab response. Briefly, we used available IBD literature to compile 5 lists of miRNAs involved in treatment response (see “Methods” and Tables S2–S6).

From these lists, we identified 74 miRNA–mRNA pairs (Table S8) that showed negative pairwise correlation (P < .05) in any of the 12 combinations in our study. Out of 35 unique pairs, we filtered out 21 in which the aligned mRNA genes corresponded to the actual host gene that generates the miRNA of interest, and a further pair in which the mRNA corresponds to a novel transcript without known function (ENSG00000275180). Table S9 lists the final 13 miRNA–mRNA pairs selected for their potential involvement in drug response, 7 of which have been described as functionally validated in the MiRPathDB (2.0) database. Considering that several of the 13 miRNAs identified in our study have been previously associated with infliximab response status (eg, miR-7704, miR-24, and miR-147b families),^23,27,28 future attempts to characterize further the interaction with their mRNA targets can be a venue to clarify regulatory mechanisms associated with drug response in CD.

Discussion

Thanks to the developments in high-throughput genomic technologies, we can now generate comprehensive profiles for different molecular layers relevant for CD patients. This includes the miRNAome, a key component of the post-transcriptional regulation of gene expression. Several works have attempted to identify miRNAs that can serve as biomarkers to track disease activity and treatment response in CD.^29–31 However, the lack of an exhaustive characterization of the miRNAome and its variability across demographic and clinical factors relevant for CD, including tissue, disease subtype, and the effect of treatment, remains a critical gap. Towards this purpose, we generated 238 miRNA and mRNA profiles from specimens obtained from a comprehensive cohort of CD patients and HC. We focused particularly on exploring the changes associated with anti-TNF treatment and its impact on the miRNA–mRNA interactome.

The different analyses of variance confirmed that intestinal tissue is the clear primary driver of miRNAome variability in CD. The patterns of correlation among pairs of miRNAs differed between the terminal ileum and the left colon, suggesting that regulatory programs are specific to each tissue. Through pathway annotation, we identified programs associated with cellular response to external stimuli, stress, and developmental biology having a role in tissue variability. We also discovered that translation-associated pathways are up-regulated in the terminal ileum, while interleukin and inflammatory signaling are more promoted in the left colon.

A fraction of miRNAs associated with tissue has previously been described as crucial in IBD development. For instance, miR-200b and miR-200c are down-regulated in CD cases compared to normal mucosa, which agrees with the down-regulation of miR-200s in the terminal ileum observed in our study. In turn, we observed up-regulation of the repressors Snail and Slug, belonging to the epithelial–mesenchymal transition (EMT) mechanisms, and involved in fibrosis development, which is more common in CD.^32,33 This also fits with the observed high correlation with same directionality of miR-200s with miR-429, a known player in the abovementioned EMT-associated regulatory programs.³⁴ These miRNAs have been associated with key transcription factors involved in suppression of EMT progression.²⁵ Indeed, the whole miR-200 family plays an important role in cancer and is increased in dysplastic lesions in UC cases.²⁶

We also observed ileal down-regulation of miR-196s, that can suppress autophagic activity via down-regulation of IRGM and LC3-II,^10,31 and which were identified as up-regulated in active and quiescent CD tissues³¹ but in turn with a decreased expression in UC.³⁵ Our observation can hence be linked with the expected tissue-specific changes since UC involves mostly colonic tissue. However, other studies have observed decreased ileal autophagy, suggesting that other factors can override the impact of miR-196s.³⁶ These findings underscore the importance of incorporating spatial differences in miRNAs to characterize the activity of different molecular pathways underlying IBD.

In addition to spatial differences, we identified expression patterns correlated with disease activity, as shown in the PCA space, where PC2 captured part of such activity. This component may reflect the inflammatory status of the patient, considering that samples with negative values in PC2 represent basal conditions with no influence of treatment. The different analyses of variance also rendered a larger residual for the miRNAome than for the coding transcriptome. Hence, other factors not included in our study might play an important role in miRNA variability. In this regard, patient interaction effects emerged as particularly relevant explanatory elements. Although random variation might be a prevalent feature of miRNA profiling, these results highlight the presence of important individual-specific determinants. This observation implies that the miRNAome can be particularly useful to study the dynamic mechanisms that underlie the molecular heterogeneity observed across CD patients. Future omic profiling efforts with larger sample sizes should aim to elucidate this proposition.

Moving beyond variability, we explored the impact of incubation with infliximab, arguably the most important biological therapeutic agent used in CD. The main result of these analyses is the identification of 9 miRNAs that are altered by this drug only in the terminal ileum of active CD cases. The annotation of these molecules suggests the involvement of these miRNAs in the modulation of the immune response and inflammatory processes. MiR200s-associated regulation of EMT and the intestinal barrier function is mediated by ZEB1 and ZEB2, which repress miR200s in a negative feedback loop.³⁴ We hypothesize that infliximab may trigger a general decrease in cytokine production, diminishing ZEB production and up-regulating miR200s levels in active CD.^34,37,38 On the other hand, miR-192, miR-142, miR-320, and miR-122 are involved in the down-regulation of NOD2 and NF-κB activity, as well as being related to autophagy. In particular, miR-192 has been implicated in IBD pathogenesis via the down-regulation of IL-8 and CXCL3, and by abating NF-κB phosphorylation.¹⁰ Our results show an up-regulation of miR-192 upon infliximab incubation (log2FC = 2.37), which can be linked with the expected reduction in TNF expression that induces overexpression of miR-192.^10,36,39 We conclude that this miRNA might be a potential target for therapy development.

Notably, we did not detect similar effects for the other combinations, especially for specimens stemming from HC. This suggests that anti-TNF treatments have strongest impact on the miRNAome alterations that are associated with disease onset and disease activity. Hence, a renewed focus on the functional implications associated with miRNA alterations that are specific to active CD cases might be a venue to improve our understanding of the variability in anti-TNF response across patients.⁴⁰

We also attempted to characterize the determinants of the miRNA–mRNA interactome in CD. To uncover regulatory networks associated with the miRNAome, we compiled a new database of miRNA–mRNA pairs obtained through various approaches, including pathway enrichment, and VP models, from several publications and 1 reference study.^23,24 Aiming at the effects of infliximab, we geared this effort towards miRNA–mRNA pairs with highest chance of involvement in infliximab response signatures. Through these analyses, we tested the proposition that the different factors relevant to CD also lead to alterations in the miRNA–mRNA interactome. We confirmed that CD onset and disease activity are associated with a loss of the homeostatic expectation for miRNA–mRNA interactions. In contrast, incubation with infliximab was associated with a recovery of the expected negative correlation in the levels of miRNA and mRNA. Despite the limitations due to the low sample size available for each combination, namely 10 individuals per group, these observations confirm that infliximab exerts a general positive impact on the relationship between miRNAs and the coding transcriptome.

Although several treatment-responsive miRNAs have been proposed, the validation status of most candidates is unclear. Moreover, genetic studies have not described miRNA-altering variants associated with anti-TNF response.¹⁴ Rather than focusing on proposing single miRNAs as biomarkers for therapy response, our results suggest that global signatures associated with positive changes in the miRNA–mRNA interactome can help to interpret studies that stratify between responders and non-responders. This aspect is reinforced by the lack of a positive impact of infliximab incubation upon the interactome of the terminal ileum of active CD cases. Again, this finding hints that alterations specific to this context are not easily amenable to anti-TNF treatment.

On a positive note, we shortlisted 13 miRNA–mRNA pairs that can be informative about the molecular changes associated with drug response. For this exercise, we compiled a diverse set of signatures of interest. We focused particularly on those that can recapitulate the effects of the miRNAome in components of the coding transcriptome that have been shown to be relevant for drug response, especially on those whose interaction has already been predicted in functional databases such as MiRPathDB. In most of these 13 miRNA–mRNA pairs, only 1 of the 2 participating molecules has been robustly associated with drug response. As discussed above, rather than attempting to find biomarkers based on exploration of a single omic layer, future works integrating miRNA and mRNA in larger cohorts will help to narrow down specific mechanistic aspects that underlie anti-TNF response in CD.

Our study has 3 relevant strengths. First, the experimental design allowed us to provide the most comprehensive characterization of the landscape of miRNA variability in the context of CD and treatment with infliximab to date. Second, we applied deep modeling to consider all the factors at the same time, hence avoiding losses in statistical power. Indeed, the overlap in detected signatures across the different methods confirms the reliability of our approach. Third, our study uniquely incorporated mRNA profiles obtained from the same samples, allowing for the characterization of the regulatory effects of miRNAs on the coding transcriptome.

Nonetheless, this work has several noteworthy weaknesses that illustrate challenges relevant for the design of future miRNAome studies in IBD. First, there is an inherent limiting aspect associated with the ex vivo design. Even if providing a controlled setting to assess the regulatory changes induced by infliximab, this design cannot fully replicate the complexity of in vivo conditions. The shorter incubation times, and localized drug exposure in a medium enriched with infliximab, necessarily differs from the prolonged and systemic effects observed with intravenous administration in patients. Considering that in vivo settings involve longer drug regimens, it is possible that the ex vivo setting leads to weaker and less sustained anti-inflammatory responses, hampering our ability to detect molecular effects associated with treatment. In addition, the ex vivo design provides only a snapshot of regulatory changes and lacks the longitudinal dimension needed to evaluate whether these alterations correlate with the long-term outcomes seen in patients, such as remission or relapse. Incorporation of longitudinal sampling schemes is a must to validate the observations and clinical implications observed in this study.

Second, due to the complex experimental design aimed at dissecting several relevant factors in the same cohort, we achieved a sample size of only 10 individuals for each combination under study. This design limited the statistical power, particularly for our aim of characterizing the dynamics of the miRNAome rather than focusing on a limited set of miRNA candidates. More importantly, CD presents ample heterogeneity across patients, and complex patient-specific factors, such as lifestyle, diet, and genetics, may have influenced the biological variability observed in the dataset. We accounted for patient-specific effects using linear mixed models, but other unaccounted factors may have affected the generalizability of the trends described here. As mentioned above, this modest sample size was necessary to ensure the complex experimental design. Nonetheless, the experience gained from this study suggests that future research, rather than just profiling larger cohorts, will also maximize their discovery power by adopting similarly balanced and multifactorial frameworks.

Finally, the experimental design based on focusing on 2 timepoints, namely before and after drug incubation, necessarily provides only a partial snapshot of miRNA changes associated with infliximab response. While revealing immediate molecular effects, this approach may not capture the dynamic changes that occur over longer periods. The importance of longitudinal profiling is increasingly recognized, especially given the growing understanding that each patient possesses unique omic “personalities” that influence disease progression and therapeutic outcomes. As is the case across all omic layers, incorporating temporal analyses into future studies will be crucial for capturing these complexities.

The ultimate goal of this study is to generate the foundation for incorporating miRNA–mRNA interactions into clinically actionable insights for CD. By identifying specific miRNAs and their regulatory networks, we highlight potential biomarkers that could guide therapeutic stratification, fueling the development of minimally invasive diagnostic tools. After nearly 2 decades of omic-based efforts, molecular profiling is being slowly incorporated into well-designed clinical trials. Although our observations suggest that miRNA profiling have potential to be geared towards such goals, this omic has received less attention and it will require sustained efforts through basic research.

In this regard, through our study, we identified a general problem due to the incomplete knowledge of the regulatory networks associated with the miRNAome. Current knowledge is biased towards cancer studies, that have a longer tradition of evaluating the miRNAome for biomarker discovery purposes. For this reason, it is often difficult to interpret the observations gathered from annotation databases for specific candidates. Besides the benefits associated with larger sample sizes, as a community we need to develop catalogs that incorporate the functional implications of the miRNAome and the miRNA–mRNA interactome that are most relevant to inflammatory settings like CD.

Conclusions

Our study based on intestinal biopsies represents a significant step towards characterizing the miRNAome in the context of CD and infliximab treatment. Along with tissue and patient-specific effects, we conclude that the effects of anti-TNF treatment in the miRNA–mRNA interactome are particularly altered in active CD cases. Further research, potentially with larger sample sizes and with improved annotation databases, is essential for unraveling the intricate miRNA regulatory networks that play a role in treatment response in CD. Additionally, the lack of available datasets for validation underscores the need for future unbiased studied profiling the miRNAome in CD and other IBDs, particularly those designed to evaluate longitudinal treatment effects. Expanding efforts to validate the candidates identified here in independent cohorts, along with the global integration of miRNAome profiles with other omic layers, will be critical for advancing their potential as biomarkers of drug response in CD.

Author Contributions

Conceptualization, M.C., J.P.G., and U.M.M.; Formal analysis, L.C.-S. and M.B.-M., Investigation, S.F.-T. and L.O.M., Resources, J.J.L., A.M.A., M.C., and J.P.G.; Writing—original draft, L.C.-S. and U.M.M.; Writing—review & editing, M.B.-M. and J.P.G. All authors have read and agreed to the last version of the manuscript.

Funding

This work was supported by a fellowship from “la Caixa” Foundation (ID 100010434) to LC-S and UMM, with fellowship codes LCF/BQ/DR20/11790008 and LCF/BQ/LI18/11630001, respectively. U.M.M. also acknowledges the Spanish Ministry of Science for grants (PID2019-108244RA-I00, RYC2020-030632-I, PID2023-151379OB-I00). M.B.-M. is supported by a Postdoctoral Contract from CIBEREHD. S.F.-T. was funded by Sara Borrell Postdoctoral fellowship CD17/00014 from the Instituto de Salud Carlos III (ISCIII). L.O.M. was funded by the Community of Madrid and Universidad Autónoma de Madrid (Ayudas Atracción de Talento modalidad 2 BMD-5800). We also thank the financial support of grants from ISCIII (PI15/00560, PI19/01034), CIBEREHD and Grupo Español de Trabajo en Enfermedad de Crohn y Colitis Ulcerosa (GETECCU). CIC bioGUNE support was provided from The Department of Industry, Tourism and Trade of the Government of the Autonomous Community of the Basque Country (Emaitek and Elkartek Research Programs 2010–2023), the Innovation Technology Department of the Bizkaia County and MINECO for the Severo Ochoa Excellence Accreditation (SEV-2016-0644).

Conflicts of Interest

Dr. Gisbert has served as speaker, consultant, and advisory member for or has received research funding from MSD, Abbvie, Pfizer, Kern Pharma, Biogen, Mylan, Takeda, Janssen, Roche, Sandoz, Celgene/Bristol Myers, Gilead/Galapagos, Lilly, Ferring, Faes Farma, Shire Pharmaceuticals, Dr. Falk Pharma, Tillotts Pharma, Chiesi, Casen Fleet, Gebro Pharma, Otsuka Pharmaceutical, Norgine and Vifor Pharma. Dr. Chaparro has served as speaker, consultant or research or education funding from MSD, Abbvie, Hospira, Pfizer, Takeda, Janssen, Ferring, Shire Pharmaceuticals, Dr. Falk Pharma, Tillotts Pharma, Biogen, Gilead, and Lilly. The remaining authors declare no financial or non-financial competing interests.

Ethical Considerations

Informed consent was obtained from all subjects involved in the study. The study was conducted in accordance with the Declaration of Helsinki and approved by the Institutional Ethics Committee of Hospital Universitario de La Princesa (Protocol GIS-INH-2015; date of approval: October 5, 2017).

Data Availability

The RNA-seq data for the 30 individuals in this study, including miRNA and mRNA profiles for 119 specimens, have been deposited in the Gene Expression Omnibus (GEO) and are accessible through the GEO Subseries GSE270663.

Author notes

Luis Cervera-Seco and Montse Baldán-Martín share co-first authorship.

Javier P Gisbert and Urko M Marigorta led the study and share the last authorship.