Next Article in Journal
Insight into Factors Influencing Wound Healing Using Phosphorylated Cellulose-Filled-Chitosan Nanocomposite Films
Previous Article in Journal
The Biological Role and Therapeutic Potential of NK Cells in Hematological and Solid Tumors
Previous Article in Special Issue
Dietary Nanoparticles Interact with Gluten Peptides and Alter the Intestinal Homeostasis Increasing the Risk of Celiac Disease
 
 
Search for Articles:
Title / Keyword
Author / Affiliation / Email
Journal
Article Type
 
 
Advanced Search
 
Section
Special Issue
Volume
Issue
Number
Page
 
Journals
IJMS
Volume 22
Issue 21
10.3390/ijms222111382
ijms-logo
Submit to this Journal Review for this Journal Propose a Special Issue
► Article Menu

Article Menu

share Share announcement Help format_quote Cite question_answer Discuss in SciProfiles thumb_up
...
Endorse textsms
...
Comment
first_page
settings
Order Article Reprints
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

A Combined mRNA- and miRNA-Sequencing Approach Reveals miRNAs as Potential Regulators of the Small Intestinal Transcriptome in Celiac Disease

by
Ineke Luise Tan
1,2,
Donatella Barisani
3,
Roberto Panceri
4,
Rutger Modderman
1,
Marijn Visschedijk
2,
Rinse K. Weersma
2,
Cisca Wijmenga
1,
Iris Jonkers
1,†,
Rodrigo Coutinho de Almeida
5,† and
Sebo Withoff
1,*
1
Department of Genetics, University Medical Center Groningen, University of Groningen, 9700 RB Groningen, The Netherlands
2
Department of Gastroenterology and Hepatology, University Medical Center Groningen, University of Groningen, 9700 RB Groningen, The Netherlands
3
School of Medicine and Surgery, University of Milano-Bicocca, 20900 Monza, Italy
4
Fondazione MBBM, S. Gerardo Hospital, 20900 Monza, Italy
5
Section Molecular Epidemiology, Department of Biomedical Data Sciences, Leiden University Medical Center, 2300 RC Leiden, The Netherlands
*
Author to whom correspondence should be addressed.
†
These authors contributed equally.
Int. J. Mol. Sci. 2021, 22(21), 11382; https://doi.org/10.3390/ijms222111382
Submission received: 26 September 2021 / Revised: 14 October 2021 / Accepted: 17 October 2021 / Published: 21 October 2021
(This article belongs to the Special Issue Molecular Research on Celiac Disease)
Browse Figures
Versions Notes

Abstract

:
Celiac disease (CeD) is triggered by gluten and results in inflammation and villous atrophy of the small intestine. We aimed to explore the role of miRNA-mediated deregulation of the transcriptome in CeD. Duodenal biopsies of CeD patients (n = 33) and control subjects (n = 10) were available for miRNA-sequencing, with RNA-sequencing also available for controls (n = 5) and CeD (n = 6). Differential expression analysis was performed to select CeD-associated miRNAs and genes. MiRNA‒target transcript pairs selected from public databases that also displayed a strong negative expression correlation in the current dataset (R < −0.7) were used to construct a CeD miRNA‒target transcript interaction network. The network includes 2030 miRNA‒target transcript interactions, including 423 experimentally validated pairs. Pathway analysis found that interactions are involved in immune-related pathways (e.g., interferon signaling) and metabolic pathways (e.g., lipid metabolism). The network includes 13 genes previously prioritized to be causally deregulated by CeD-associated genomic variants, including STAT1. CeD-associated miRNAs might play a role in promoting inflammation and decreasing lipid metabolism in the small intestine, thereby contributing unbalanced cell turnover in the intestinal crypt. Some CeD-associated miRNAs deregulate genes that are also affected by genomic CeD-risk variants, adding an additional layer of complexity to the deregulated transcriptome in CeD.

1. Introduction

In patients with celiac disease (CeD), immune-mediated destruction of the small intestinal villous structure takes place as a response to the presence of dietary gluten. CeD occurs in 1–2% of the Caucasian population in genetically predisposed individuals [1,2]. Nearly all CeD patients carry the genetic risk human leukocyte antigen (HLA) haplotypes HLA-DQ2 or HLA-DQ8. Besides the HLA region, more than 40 non-HLA CeD risk loci have been identified by genome-wide association studies (GWAS) [3,4].
Upon gluten intake, antigen-presenting cells present deamidated gliadin peptides to CD4+ cells in the context of the predisposing HLA-types (HLA-DQ2/8). The activated CD4+ gluten-specific T-cells subsequently stimulate CD8+ T-cells, which eventually move to the intraepithelial compartment and destroy the villous structure of the small intestine [1,5].
Previous studies have shown that the transcriptome of the small intestine is affected in CeD [6,7,8]. Affected pathways include immune-related pathways for B- and T-cells and metabolic pathways (e.g., fatty-acid metabolism). A recent study showed that CeD-associated genomic variants (single nucleotide polymorphisms (SNPs)) influence gene expression and compiled a list of 118 prioritized genes likely to play a role in CeD pathophysiology [9]. This list includes genes involved in immune cell migration, activation, and differentiation. In addition to SNPs, gene expression can be regulated by multiple other factors, both on the transcriptional and post-transcriptional level, including microRNAs (miRNAs) [10]. In the current study, we examined the role of small intestinal miRNAs in CeD pathophysiology.
miRNAs are small non-coding RNAs consisting of 19–24 nucleotides. The physiological role of cellular miRNAs is to fine-tune gene expression by directing the RNA-induced silencing complex to target messenger RNAs (mRNAs) [11]. miRNA target recognition takes place via pairing of the seed region of the miRNA (nucleotides 2–7) with the 3′ untranslated region of the target mRNA, which eventually results in inhibition of translation or promotion of decay of the target transcript [11].
Recently, extracellular miRNAs in the circulation have been proposed as biomarkers for CeD [12]. The small intestine could be the source of these circulating miRNAs, as it has been shown that the small intestinal miRNA profiles are different between CeD and controls [13,14,15,16,17]. However, the functional roles of both circulating and intracellular miRNAs in CeD pathophysiology are not yet fully understood. To uncover the functions and pathways in which CeD-associated miRNAs play a role, the transcripts targeted by the deregulated miRNAs need to be identified. One expects that when miRNAs are overexpressed, target transcripts are downregulated [14,17]. It is also important to investigate these interdependencies in a tissue-specific context, as miRNA expression is highly cell type‒specific [10,18,19]. In the context of CeD, only a few miRNA‒transcript pairs have been investigated thus far as no high-throughput studies have searched for such miRNA‒target transcript pairs in CeD [14,17].
To get insight into how miRNAs affect target genes and pathways in the small intestine of CeD patients, we performed miRNA- and RNA-sequencing on small intestinal biopsies from CeD patients and controls. Using these samples, we investigated which miRNAs were inversely correlated with predicted and experimentally validated target transcripts and used these data to construct a miRNA‒transcript interaction network for CeD. Subsequent pathway-analyses on the target genes revealed that CeD-associated miRNAs are associated with increased inflammatory processes and unbalanced cell turnover in the intestinal crypt. Several miRNA target genes revealed by our analysis have previously been shown to be affected by CeD-risk SNPs, adding another level of complexity to the deregulated gene expression in CeD.

2. Results

2.1. Differential Expression Analysis in CeD Identifies Immune-Related Genes and miRNAs

To start, we used principal component analysis (PCA) to investigate whether miRNA profile correlates with CeD status. No outliers were observed in the PCA. The PCA showed a clear significant association with disease status (Figure 1A), independent of sex, age, and RNA isolation method.
Differential expression analysis identified 111 miRNAs to be differentially expressed between CeD and controls (FDR < 0.05): 52 that were decreased and 59 that were increased in CeD (Supplementary Table S2). These included miRNAs that had been previously identified to be decreased in CeD: miR-338-5p (Log2 Fold Change (LFC) = −2, FDR = 1.4 × 10−11), miR-192-5p (LFC = −1.2, FDR = 1.1 × 10−8), miR-194-5p (LFC = –0.8, FDR = 1.2 × 10−4), miR-31-3p (LFC = −1.0, FDR = 7.0 × 10−7) (Supplemental Table S2) [14,15,20], and ‘immune miRNA’ miR-155-5p (LFC = 1.2, FDR = 1.3 × 10−5) [11,21,22,23,24].
No outliers were observed in the mRNA expression profile in a PCA. The mRNA expression profile of CeD biopsies differed significantly from that of control biopsy according to the PCA (Figure 1B), independent of sex, age, and RNA isolation method.
In total, 3869 genes were differentially expressed between CeD and controls (FDR < 0.05), of which 2209 were downregulated and 1660 were upregulated (Supplementary Table S3). We compared this differentially expressed gene list to the list of Loberman-Nachum et al., who previously compiled genes consistently reported to be involved in CeD by multiple studies (Supplementary Table S4) [8]. The majority of genes in the Loberman-Nachum et al. consensus set (n = 403) [8] were also concordantly and significantly differentially expressed in our dataset (n = 334) (334/403 = 83%). These genes included the five genes that showed the highest discriminatory value between CeD and controls in their dataset: LPL (encoding lipoprotein lipase, plays an important role in lipid clearance, utilization, and storage), BIRC3 (C-IAP-2/baculoviral IAP repeat containing 3, an anti-apoptotic protein binding to TRAF-1 and TRAF-2), UGT2B7 (UDP glucuronosyltransferase family 2 member B7, involved in conjugation and subsequent elimination of potentially toxic compounds), THSD4 (thrombospondin type 1 domain containing 4, a peptidase involved in matrix homeostasis), and HMGCS2 (3-hydroxy-3-methylglutaryl-coA synthase 2, involved in lipid metabolism) [8,25]. The consensus set included genes that are associated with genetic risk loci for CeD, including STAT1 (LFC 1.6, FDR = 2.7 × 10−13), which encodes a transcription factor important in the response to interferon signaling [8,9].
To identify pathways related to CeD, we performed pathway-enrichment analysis on the differentially expressed genes in our study. GO terms associated to the down- and up-regulated genes in our study are displayed in Supplementary Figure S1. The transcripts upregulated in CeD are associated with several immune pathways, including B- and T-cell activation (e.g., ‘regulation of lymphocyte activation’ (FDR = 4.6 × 10−13) and ‘interferon-gamma pathways’ (FDR = 4.5 × 10−8)), and with pathways related to the cell cycle (e.g., ‘positive regulation of cell cycle process’ (FDR = 3.4 × 10−11)). Downregulated genes were associated to several metabolic pathways, including lipid metabolism (‘Fatty acid metabolic process’ (FDR = 2.7 × 10−27) and ‘digestion’ (FDR = 1.94 × 10−5)). These and similar pathways have previously been linked to CeD [8].

2.2. miRNA‒Target Transcript Interaction Network

In order to identify miRNAs and their target transcripts, we focused on those miRNAs that were differentially expressed between CeD and controls. We selected miRNA‒target pairs from two prediction tools (TargetScan and microTCDS) and from databases that list experimentally validated miRNA‒target pairs (TarBase and miRTarbase) (Figure 1). Afterwards, we performed a Pearson correlation between the selected miRNAs and genes. This identified 2030 negatively correlated miRNA‒target transcript pairs (R < −0.7 and p-value < 0.05). The resulting miRNA‒target transcript interaction network consists of 31 miRNAs connected to 1344 genes (Figure 2, supporting data displayed in Supplementary Table S5). As mentioned before, we identified 3869 genes to be differentially expressed between CeD and controls, 2209 downregulated genes and 1660 upregulated genes. In total, 39% of all downregulated genes (866/2209) and 29% of all upregulated genes (478/1660) are targeted by at least one differentially expressed miRNA. Of the 2030 miRNA‒transcript pairs, only 423 miRNA‒transcript interactions have been experimentally validated. Transcripts that were targeted by more than one miRNA showed a slight, but significantly stronger, negative correlation with miRNA levels compared to genes targeted by only one miRNA (respectively, mean correlation median −0.78 interquartile range (IQR) [−0.75; −0.81] versus correlation −0.76 IQR [−0.73; −0.82], MWU p-value = 0.002).
MiRNA families consist of miRNAs that share homology in their seed sequence, resulting in target transcript similarities [11]. In our miRNA‒target transcript network, multiple members of the same miRNA family could be identified amongst both upregulated miRNAs (miR17 family: miR-18a-3p and miR-17-5p; miR15 family: miR-15a/b-5p and miR-16-5p) and downregulated miRNAs (miR30 family: miR-30a/e-3p; miR28 family: miR-151b and miR-28-5p).

2.3. Pathway Enrichment Analyses

To assess which pathways might be affected by CeD-associated miRNAs, we performed pathway enrichment analyses based on the transcripts present in the miRNA‒target transcript network. This identified 360 pathways that were significantly associated with the miRNA‒target transcript network (Supplementary Table S6). Figure 3 shows the top 30 pathways associated with the network, split into up- and downregulated transcripts and sorted based on how many different miRNA target transcripts were present in the pathway. Downregulated genes in the miRNA‒target transcript network were associated with pathways related to metabolism, such as lipid metabolism, whereas upregulated genes were associated with immune pathways related to T-cells, response to interferon-gamma, and cell-cycle (Figure 3 and Table S6).
We prioritized miRNA target transcripts based on the criteria that they are targeted by multiple miRNAs and that there is supporting experimental evidence for the miRNA‒target transcript interaction. This analysis identified 492 transcripts that were targeted by at least two miRNAs, and there is experimental evidence for at least one of the interacting miRNAs for 34% (168/492) of these target transcripts. The downregulated targets (118/168) are significantly enriched with transcripts involved in metabolic pathways (phospholipid metabolic process, FDR = 0.03) and epithelial cell maturation (FDR = 0.046), whereas the upregulated transcripts targets are involved in cell-cycle pathways (e.g., mitotic nuclear division, FDR = 0.003), positive regulation of type I interferon (FDR = 0.03), and the NOTCH pathway (FDR = 0.023).
To investigate whether the individual miRNAs in the network regulate specific pathways, we also performed gene set enrichments per miRNA (Figure 4). This revealed that, in our miRNA‒target transcript network, multiple individual miRNAs target similar pathways even though their target transcripts vary. For example, the 44 target transcripts of miR-31-3p (see Figure 4) are significantly enriched in transcripts involved in lymphocyte differentiation, and this was also the case for the target transcripts of miR-22-5p and miR-30a-3p.

2.4. Cell Type‒Enrichment Analysis

As small and bulk RNA sequencing results can be influenced by tissue composition, we also assessed the association between the tissue composition of the biopsy, the disease status of the patient, and the miRNA expression.
Using the xCell package, we calculated cell type‒enrichment scores based on the RNA-seq data for multiple cell types relevant in CeD pathophysiology: epithelial cells, CD4- and CD8-positive T-cells, and B-cells (Figure 5A,B). Although the differences between CeD (n = 6) and controls (n = 5) did not reach significance for any of these cell types, immune cells did show a higher enrichment in CeD biopsies (B-cells: p = 0.03; CD8+ T-cells p = 0.076; immune cell score p = 0.094), and epithelial cells were significantly depleted in CeD compared to controls (p = 0.03) (Figure 5A). To investigate in which cell types the miRNAs in the miRNA‒target transcript network are expressed, we correlated the expression levels of the miRNAs with the cell type–enrichment scores. Figure 5B shows that miRNAs downregulated in CeD were correlated with enrichment for epithelial cells and that miRNAs upregulated in CeD were correlated with enrichment for immune cells.
We then explored whether the CeD miRNAs are expressed in a cell type‒specific manner. A previous overview by De Rie et al. of the enrichment of specific miRNAs in purified human cell types had shown that miRNA expression is highly cell type‒specific [26]. We used the De Rie et al. data to explore whether the 111 up- and downregulated miRNAs in CeD are likely expressed by immune cells or intestinal epithelial cells. Here, we observed a significant correlation between the fold change of miRNAs between CeD and controls in our study and the enrichment scores for CD3+ T-cells (Pearson’s R = 0.27 p = 0.02) and intestinal epithelial cells (Pearson’s R = −0.29, p = 0.01) reported by De Rie et al. Enrichment scores of the other cell subtypes that were analyzed did not reach significance: CD4+ T-cells (p = 0.97), CD8+ T-cells (p = 0.59), and B-cells (p = 0.894). Because not all the CeD miRNAs we identified were included in De Rie et al., we calculated a miRNA-based enrichment score for each miRNA profile by calculating enrichment for the top 10 CD3+ miRNAs and the top 10 intestinal epithelial cell miRNAs from De Rie et al. Our CeD miRNA profiles showed an enrichment for De Rie et al.’s CD3+ miRNAs and a depletion of De Rie et al.’s intestinal epithelial cell miRNAs (Figure 5C). Moreover, correlation of the miRNA-based enrichment scores with the expression of the miRNAs that were included in our miRNA‒target transcript network led to similar results (Figure 5D). miRNAs downregulated in CeD were enriched for miRNAs expressed in intestinal epithelial cells, whereas miRNAs upregulated in CeD were enriched for the immune cell‒specific miRNAs.
These results suggest that cell type–composition can partially explain the differences in miRNA expression between CeD and controls. Therefore, we again tested the differences between CeD and controls (Supplementary Table S2), including the miRNA-based cell type–enrichment scores in the statistical model to correct for cell type‒composition. Of the 31 miRNAs in the miRNA‒target transcript network, nine showed differences between CeD and controls (FDR < 0.1) that were independent of enrichment of CD3+ or intestinal epithelial cell types.

2.5. Linking the miRNA‒Target Transcript Network to Genes in CeD-Associated Genetic Risk Regions

A previous study integrated multiple in-silico approaches to prioritize genes in CeD-associated genomic risk regions and identified 118 genes as likely causal in CeD [9]. Our data show a clear enrichment of the transcripts of these prioritized genes in CeD biopsies (Figure 6A). Of the 118 genes, 102 were expressed in our biopsy transcriptome dataset, and 13 of the 102 are also found in our miRNA‒target transcript network (Figure 6B). This suggests that, in addition to being affected by genetic factors that predispose to CeD, these genes might also be regulated post-transcriptionally by CeD-associated miRNAs. This network includes the previously mentioned STAT1 [9]. In the miRNA‒target transcript network, STAT1 is associated with miR-22-5p (R = −0.86; p = 6.6 × 10−4) and miR-30a-3p (R = −0.88; p = 3.9 × 10−4), and there is other previous experimental evidence (TarBase: HITS-CLIP) to support the interaction with miR-30a-3p (Supplementary Table S5).

3. Discussion

To our knowledge, this is the first study to use next-generation sequencing to generate a CeD-specific miRNA‒target transcript interaction network, thereby providing the first unbiased analysis of miRNAs and their targets in the context of CeD. For our analyses, we used public databases of predicted miRNA‒target transcript pairs or experimentally validated miRNA‒target transcript pairs that showed a strong negative correlation in our expression dataset. This resulted in a network of 2030 miRNA‒target transcript pairs that provides a starting point for understanding the complex relations between miRNAs and target transcript regulation in CeD pathophysiology.
Our analyses show that miRNA target transcripts are involved in many pathways important in CeD pathogenesis. Downregulated transcripts appear to be involved in, for instance, (lipid) metabolism pathways [12], and upregulated (de-repressed) transcripts appear to play a role in cell cycle pathways and immune-pathways (e.g., T- and B-cell-pathways, interferon pathways) [27,28,29]. The network also includes several sets of miRNAs belonging to the same miRNA families (families miR17, miR15, miR30, and miR28). MiRNA families are defined by homology in their seed sequence, resulting in the targeting of the same transcript [11]. Additionally, we identified different miRNAs that target different transcripts in the same pathways, which suggests that miRNAs from different families can cooperate to enhance repression of particular pathway, as previously suggested by others [19,30].
Our results suggest that CeD-associated miRNAs are involved in regulating barrier homeostasis, a process that is affected in CeD [12]. Some miRNAs appear to do this by affecting lipid metabolism in the small intestine. It has been shown, for instance, that mice with a small intestinal DICER-knock out display abnormal absorption and processing of lipids [31,32]. Lipid metabolism is also important in the maintenance of the regenerative capacity of the small intestinal crypt [33,34]. One of the lipid metabolism transcripts targeted in our small-intestinal miRNA‒target transcript network is the transcription factor PRDM16, which is targeted by miR-500a-3p and miR-361-3p. In a murine model, Stine et al. showed that loss of Prdm16 inhibits transcription of many fatty-acid oxidation genes, resulting in villous atrophy of the small intestine, which is a hallmark of CeD [35]. We observed that a number of fatty-acid oxidation genes that are regulated by PRDM16 (CPT1A, ECI1, CD36, ACSL1, ACAA2, and HADH) are also targeted by CeD-associated miRNAs (miR-361-3p, miR-155-5p, miR-15a-5p, miR-18a-3p, miR-425-5p, and miR-138-5p). This shows that deregulated small intestinal miRNAs in CeD patients contribute to villous atrophy by regulating genes related to fatty-acid metabolism. Another miRNA prioritized in our study, miR-31-3p, is downregulated in CeD and has been previously associated with cell cycle and immune pathways [14,15,20]. Tian et al. showed that miR-31-3p is highly expressed in the intestinal epithelial crypt and that it is an important restorative factor in the intestine through maintenance of the homeostasis of cell turnover from the intestinal crypt to the villous tip [36]. MiR-31-3p knockout mice display more severe intestinal inflammation as a response to chemically induced colitis (DSS), and this response can be dampened by administration of miR-31-3p [37]. Other target transcripts that play a role in intestinal epithelial maintenance are RXRA (targeted by the upregulated 155-5p, miR-1260b, miR-132-3p, miR-425-5p, miR-18a-3p, and miR-425-5p; has a key role in retinol signaling in the differentiation of mature enterocytes [38]), VAV2 (targeted by 155-5p, miR-15a/b-5p, and miR-17-5p; has a role in wound repair in the intestine and in differentiation and migration of mature enterocytes along the crypt-villous axis via RAC1 [39,40]), CUX1 (miR-132-3p; transcription factor targeting VAV2 [40]), and PACSIN2 (miR-155-5p, miR-1260b, miR-138-5p, and miR-361-3p; controls morphology of the microvilli [41]).
MiRNAs also appear to deregulate the immune response in CeD, for instance by affecting interferon signaling, which is key in CeD pathophysiology [9,28,29]. Above, we discussed how miR-31-3p is involved in intestinal barrier homeostasis, but downregulation of miR-31-3p has also been shown to enhance the response of CD8+ cells to viral triggers, leading to a higher pro-inflammatory interferon response [42]. The target genes of miR-31-3p include HMGB2 and PRKDC, which are both nucleic acid sensors that are important in eliciting an interferon response after sensing of cytosolic DNA [43,44]. Interestingly, a paralog of HMGB2, HMGB1, has been proposed as a biomarker for CeD [45,46,47]. Another prioritized and upregulated miRNA is miR-155-5p, a well-described immune miRNA that also has a role in enhancing the interferon response [21,22,23,24]. UBXN1 is a target gene of miR-155-5p in our miRNA‒target transcript network. A previous study has shown that UBXN1 inhibits pro-inflammatory NF-κB signaling and the interferon response to viral stimuli [48,49]. Furthermore, knockdown of the miR-155-5p target gene JUNB in regulatory T-cells has been shown to lead to increased production of pro-inflammatory cytokines, such as IFN-gamma, in colonic tissue [50,51]. Regulatory T-cell function has previously been shown to be affected in individuals with CeD [28]. In addition to miR-155-5p and miR-31-3p, other miRNAs also play a role in interferon signaling. Our network connects six additional downregulated miRNAs (miR-192-3p, miR-30a-3p, miR-653-5p, miR-22-5p, miR-28-5p, and miR-151b) to transcripts involved in interferon-gamma response. Two of these miRNAs, miR-22-5p and miR-30a-3p, were previously described to enhance interferon signaling [52,53]. Altogether, our analysis shows that the deregulated interferon response in CeD is partially regulated by a CeD-specific miRNA profile.
Interestingly, 13 of the CeD-associated target transcripts identified by our network have previously been described to be potentially causally deregulated by CeD-associated SNPs [9]. Two transcripts have also been experimentally validated to be the target of the CeD-associated miRNAs in our network: STAT1 has been shown to be regulated by miR-30a-3p and ERRFI1 by miR-138-5p [9]. MiR-22-5p and miR-30a-3p, which we described above as regulating genes associated with interferon signaling, also affect the transcripts of genes in genetic risk loci for CeD such as STAT1 (miR-22-5p and miR-30a-3p) and TRAFD1 (miR-30a-3p) [9]. Several CeD-associated SNPs can cause altered binding sites for miRNAs, thereby affecting the binding between miRNA and gene, depending on genotype [54]. The current study could not detect these kinds of differences due to the lack of genotype information and, more importantly, limited sample size. However, taken together, these results suggest that miRNAs and genetic risk SNPs associated with CeD cooperate in deregulating the expression of transcripts in CeD pathophysiology.
A limitation in miRNA research is that most miRNA‒target transcript interactions have been predicted by target-prediction algorithms but have not been validated experimentally. In the current network, 423 out of 2030 interactions have experimental evidence supporting these interactions. The best way to get insight into physiological miRNA‒target transcript interactions would be via crosslinking-based methods such as HITS-CLIP, in which miRNAs are crosslinked to target transcripts and subsequently sequenced [11,55]. Both miRNAs and mRNA are expressed in a cell type‒specific manner, and therefore expression in a biopsy is expected to depend on cell type differences between CeD and controls (e.g., expansion of lymphocytes). Indeed, we observed that the predicted cell type composition explained part of the miRNA differences that we observed between CeD and controls. However, even though miRNAs may not be produced in the same cells as the target transcripts, these miRNAs might still be able to affect transcripts with different target cells [56,57]. Because of these cell type‒dependent expression levels, miRNA‒target transcript interactions should ideally be performed at single cell‒level [58,59]. While the currently available single cell techniques are not yet capable of capturing the full spectrum of miRNA‒target transcript interactions, they do hold potential for the future [60].
Interestingly, our biopsy-focused approach identified three miRNAs—miR-21-3p, miR-500a-3p, and miR-15b-5p—that also can be detected in circulation as biomarkers for CeD [13,15,16,61] (Tan et al. submitted manuscript). We previously observed that these three miRNAs were upregulated in the circulation of CeD patients compared to controls up to two years the rise in current serological antibodies (anti-tissue transglutaminase) could be detected (Tan et al. submitted manuscript).
Taken together, the results of our paired miRNA‒target transcript sequencing study of small intestinal biopsies of CeD patients versus controls have revealed that the miRNAs deregulated in CeD could play a role in metabolic, cell cycle, and immune pathways that are deregulated in the small intestine in CeD. Our study is an exploratory, hypothesis generating study to investigate the potential role of miRNAs in CeD. Future functional studies should be performed to validate and confirm the role of miRNAs candidates in the pathogenesis of CeD, preferably in a cell-type specific manner. Moreover, to get more insights in the specificity of these regulatory miRNA-gene interaction pathways, it would be of value to also include other types of intestinal inflammation (e.g., Crohn’s disease) in future studies. A better understanding of the role of these miRNAs in CeD pathogenesis could aid in the search for biomarkers relevant to disease processes and the identification of novel therapeutic options for CeD.

4. Materials and Methods

4.1. Sample Collection

Pediatric patients and controls were included at the San Gerardo Hospital, Monza, Italy. The parents of all participants provided informed consent for the study. Duodenal biopsies were collected from untreated CeD patients at time of diagnosis (n = 33). CeD diagnosis was established based on serology (anti-transglutaminase antibodies) and histopathological examination (villous atrophy and influx of intraepithelial lymphocytes). Biopsies were also collected from control individuals (n = 10) who underwent upper endoscopies for other indications and did not show signs of CeD in the histopathological examination of small-intestinal biopsies. Clinical characteristics are described in Supplementary Table S1. The study was conducted after approval of the San Gerardo Hospital Ethical Committee, Monza, Italy.

4.2. RNA Isolation, Small RNA, and mRNA-Sequencing

RNA was isolated from the small-intestinal biopsies using either the miRVana isolation Kit (Ambion, Carlsbad, CA, USA) or qiazol lysis reagent (Qiagen, Hilden, Germany, 79306). The proportion of CeD and controls did not differ between both isolation methods (χ2 P = 0.43). RNA quality was assessed using the Caliper GX bioanalyzer (Agilent, Santa Clara, CA, USA). Small RNA-libraries were generated starting from 500 ng total RNA using the TruSeq Small RNA Sample Prep kit (Illumina, San Diego, CA, USA), implementing 15 amplification cycles. RNA libraries were prepared as previously described in Zorro et al. [9,62], using the Illumina TruSeq stranded total RNA library kit with a riboZero rRNA depletion step. After measurement of the cDNA concentration, libraries were pooled equimolarly per lane on a HiSeq 2500 (Illumina San Diego, CA, USA). Raw reads were aligned to miRbase 22 (small-RNA-seq) and human_g1k_v37 ensemble Release 75 (RNA-seq), as described previously (Tan et al. submitted manuscript; [62]). For all 43 samples, small-RNA libraries were generated and sequenced. Bulk RNA-sequencing was performed for 5/10 of the controls and 6/33 of the CeD patients (GEO accession number GSE146190) [9,62].

4.3. Prioritizing Genes Targeted by miRNAs in the Small Intestinal Biopsies

4.3.1. Differential Expression Analyses in miRNA Sequencing and RNA Sequencing

To find miRNAs and genes that are deregulated in CeD in small intestinal biopsies, we performed differential expression analysis between CeD patients and controls for the miRNA profiles and RNA profiles separately (R-package DESeq2, version 1.26.0 [63]), while correcting for the covariates age and sex. p-values were adjusted for multiple testing using the Benjamini-Hochberg correction for False Discovery Rate (FDR) [64]. Regularized log-normalized miRNA and RNA counts were used in all downstream analyses.

4.3.2. miRNA‒Target Transcript Network

To build a miRNA‒target transcript network (see Figure 1), we first identified miRNA‒target transcript pairs by combining data provided by two prediction databases (TargetScan version 7.2 [65] and microTCDS version 7.0 [66]) and by two experimentally-validated miRNA target databases (TarBase version 7.0 [67] and miRTarbase version 7.0 [68]). We then calculated Pearson’s correlations for all miRNA‒target transcript pairs, and only those pairs with negative Pearson’s correlations (R < −0.7) and a p-value < 0.05 were subsequently integrated into the miRNA‒target transcript network and visualized with the RedeR package [69].

4.4. Enrichment Analyses

We performed pathway analyses to identify which pathways were associated to (1) transcripts that were differentially expressed in CeD versus controls and (2) all the transcripts combined that are targeted by miRNAs in the miRNA‒target transcript network. To zoom in on specific functions of the individual miRNAs in the miRNA‒target transcript network, we also performed pathway analyses using separate transcript lists per miRNA. Pathway analyses were performed using the R-package clusterProfiler (version 3.14.3) using Gene Ontology (GO) terms (Biological Process) [70]. The online tool REVIGO was used to reduce the number of redundant GO terms from the long lists of significantly associated GO terms (settings: small 0.5) [71].

4.5. Cell Type‒Enrichment Analyses

To gain insight into the cell types that might contribute to the miRNA and bulk transcript differences between CeD and controls, we calculated enrichment scores for cell types using two approaches. In the first approach, we used xCell, an in-silico method that can be used to reliably calculate enrichment of certain cell types in bulk RNA-seq samples [72,73,74]. In the second approach, we used cell type–specific miRNA information to calculate cell type–enrichment scores from the miRNA-seq data. Here, we extracted the top 10 most-enriched miRNAs specific for immune cells and intestinal epithelial cells from a publicly accessible miRNA expression atlas [26]. This list was then used as input for the R-package GSVA (version 1.34.0) to perform single sample enrichment analyses to calculate: (1) the miRNA-based cell type enrichment score and (2) the enrichment for previously prioritized CeD genes [9].

Supplementary Materials

The following are available online at https://www.mdpi.com/article/10.3390/ijms222111382/s1.

Author Contributions

Conceptualization and study design: I.L.T., D.B., R.P., R.M., M.V., R.K.W., C.W., I.J., R.C.d.A., S.W.; sample collection, D.B. and R.P.; sample processing, R.C.d.A. and R.M.; data analysis and visualization, I.L.T. and R.C.d.A.; data analysis supervision, R.C.d.A., I.J. and S.W.; writing—original draft preparation, I.L.T. and R.C.d.A.; writing—review and editing: I.L.T.; D.B., R.P., R.M., M.V., R.K.W., C.W., I.J., R.C.d.A., S.W.; supervision: I.J., R.C.d.A. and S.W. All authors read and agreed to the published version of the manuscript.

Funding

I.T. is supported by a MD/PhD scholarship from the Junior Scientific Masterclass (Graduate School of Medical Sciences, University Medical Center Groningen and University of Groningen) and a Royal Netherlands Academy of Arts and Sciences (KNAW) Ter Meulen grant. DB is supported by 2016-ATE-0312; I.J. by a Rosalind Franklin Fellowship from the University of Groningen and a Netherlands Organization for Scientific Research (NWO) VIDI grant (no. 016.171.047). S.W. and C.W. were supported by The Netherlands Organ-on-Chip Initiative, an NWO Gravitation project (024.003.001) funded by the Ministry of Education, Culture, and Science of the government of The Netherlands; and European Research Council advanced grant (FP7/2007-2013/ERC Advanced Grant Agreement 2012-322698).

Institutional Review Board Statement

The study was conducted after approval of the San Gerardo Hospital, Monza, Italy [9,62].

Informed Consent Statement

The parents of all participants provided informed consent for the study.

Data Availability Statement

GEO accession number GSE146190) [9,62].

Acknowledgments

We would like to thank all the participants and their families for the donation of biomaterials and phenotypes. We would like to acknowledge the research group of M. Swertz for the high-performance computing infrastructure and data storage. We also thank Kate Mc Intyre for editing the manuscript.

Conflicts of Interest

The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

References

  1. Lindfors, K.; Ciacci, C.; Kurppa, K.; Lundin, K.E.A.; Makharia, G.K.; Mearin, M.L.; Murray, J.A.; Verdu, E.F.; Kaukinen, K. Coeliac disease. Nat. Rev. Dis. Primers 2019, 5, 1–18. [Google Scholar] [CrossRef]
  2. Singh, P.; Arora, A.; Strand, T.A.; Leffler, D.A.; Catassi, C.; Green, P.H.; Kelly, C.P.; Ahuja, V.; Makharia, G.K. Global Prevalence of Celiac Disease: Systematic Review and Meta-analysis. Clin. Gastroenterol. Hepatol. 2018, 16, 823–836.e2. [Google Scholar] [CrossRef] [Green Version]
  3. Ricaño-Ponce, I.; Gutierrez-Achury, J.; Costa, A.F.; Deelen, P.; Kurilshikov, A.; Zorro, M.M.; Platteel, M.; van der Graaf, A.; Sugai, E.; Moreno, M.L.; et al. Immunochip meta-analysis in European and Argentinian populations identifies two novel genetic loci associated with celiac disease. Eur. J. Hum. Genet. 2020, 28, 313–323. [Google Scholar] [CrossRef]
  4. Trynka, G.; Hunt, K.A.; Bockett, N.A.; Romanos, J.; Mistry, V.; Szperl, A.; Bakker, S.F.; Bardella, M.T.; Bhaw-Rosun, L.; Castillejo, G.; et al. Dense genotyping identifies and localizes multiple common and rare variant association signals in celiac disease. Nat. Genet. 2011, 43, 1193–1201. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  5. Jabri, B.; Abadie, V. IL-15 functions as a danger signal to regulate tissue-resident T cells and tissue destruction. Nat. Rev. Immunol. 2015, 15, 771–783. [Google Scholar] [CrossRef] [PubMed]
  6. Bragde, H.; Jansson, U.; Jarlsfelt, I.; Söderman, J. Gene expression profiling of duodenal biopsies discriminates celiac disease mucosa from normal mucosa. Pediatr. Res. 2011, 69, 530–537. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  7. Bragde, H.; Jansson, U.; Fredrikson, M.; Grodzinsky, E.; Söderman, J. Celiac disease biomarkers identified by transcriptome analysis of small intestinal biopsies. Cell. Mol. Life Sci. 2018, 75, 4385–4401. [Google Scholar] [CrossRef] [Green Version]
  8. Loberman-Nachum, N.; Sosnovski, K.; Di Segni, A.; Efroni, G.; Braun, T.; BenShoshan, M.; Anafi, L.; Avivi, C.; Barshack, I.; Shouval, D.S.; et al. Defining the Celiac Disease Transcriptome using Clinical Pathology Specimens Reveals Biologic Pathways and Supports Diagnosis. Sci. Rep. 2019, 9, 16163. [Google Scholar] [CrossRef] [Green Version]
  9. van der Graaf, A.; Zorro, M.; Claringbould, A.; Vosa, U.; Aguirre-Gamboa, R.; Li, C.; Mooiweer, J.; Ricano-Ponce, I.; Borek, Z.; Koning, F.; et al. Systematic prioritization of candidate genes in disease loci identifies TRAFD1 as a master regulator of IFNγ signalling in celiac disease. Front. Genet. 2020, 11, 1–16. [Google Scholar] [CrossRef] [Green Version]
  10. Chen, P.; Liao, K.; Xiao, C. MicroRNA says no to mass production. Nat. Immunol. 2018, 19, 1040–1042. [Google Scholar] [CrossRef] [PubMed]
  11. Bartel, D.P. Metazoan MicroRNAs. Cell 2018, 173, 20–51. [Google Scholar] [CrossRef] [Green Version]
  12. Ramírez-sánchez, A.D.; Tan, I.L.; Gonera-De Jong, B.C.; Visschedijk, M.C.; Jonkers, I.; Withoff, S. Molecular biomarkers for celiac disease: Past, present and future. Int. J. Mol. Sci. 2020, 21, 8528. [Google Scholar] [CrossRef]
  13. Capuano, M.; Iaffaldano, L.; Tinto, N.; Montanaro, D.; Capobianco, V.; Izzo, V.; Tucci, F.; Troncone, G.; Greco, L.; Sacchetti, L. MicroRNA-449a overexpression, reduced NOTCH1 signals and scarce goblet cells characterize the small intestine of celiac patients. PLoS ONE 2011, 6, e29094. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  14. Magni, S.; Comani, G.B.; Elli, L.; Vanessi, S.; Ballarini, E.; Nicolini, G.; Rusconi, M.; Castoldi, M.; Meneveri, R.; Muckenthaler, M.U.; et al. MIRNAs affect the expression of innate and adaptive immunity proteins in celiac disease. Am. J. Gastroenterol. 2014, 109, 1662. [Google Scholar] [CrossRef] [PubMed]
  15. Buoli Comani, G.; Panceri, R.; Dinelli, M.; Biondi, A.; Mancuso, C.; Meneveri, R.; Barisani, D. miRNA-regulated gene expression differs in celiac disease patients according to the age of presentation. Genes Nutr. 2015, 10, 482. [Google Scholar] [CrossRef] [Green Version]
  16. Bascuñán, K.A.; Pérez-Bravo, F.; Gaudioso, G.; Vaira, V.; Roncoroni, L.; Elli, L.; Monguzzi, E.; Araya, M. A miRNA-Based Blood and Mucosal Approach for Detecting and Monitoring Celiac Disease. Dig. Dis. Sci. 2019, 65, 1982–1991. [Google Scholar] [CrossRef] [PubMed]
  17. Comincini, S.; Manai, F.; Meazza, C.; Pagani, S.; Martinelli, C.; Pasqua, N.; Pelizzo, G.; Biggiogera, M.; Bozzola, M. Identification of autophagy-related genes and their regulatory miRNAs associated with celiac disease in children. Int. J. Mol. Sci. 2017, 18, 391. [Google Scholar] [CrossRef] [PubMed]
  18. Hsin, J.P.; Lu, Y.; Loeb, G.B.; Leslie, C.S.; Rudensky, A.Y. The effect of cellular context on miR-155-mediated gene regulation in four major immune cell types. Nat. Immunol. 2018, 19, 1137–1145. [Google Scholar] [CrossRef]
  19. Bracken, C.P.; Scott, H.S.; Goodall, G.J. A network-biology perspective of microRNA function and dysfunction in cancer. Nat. Rev. Genet. 2016, 17, 719–732. [Google Scholar] [CrossRef]
  20. Vaira, V.; Roncoroni, L.; Barisani, D.; Gaudioso, G.; Bosari, S.; Bulfamante, G.; Doneda, L.; Conte, D.; Tomba, C.; Bardella, M.T.; et al. microRNA profiles in coeliac patients distinguish different clinical phenotypes and are modulated by gliadin peptides in primary duodenal fibroblasts. Clin. Sci. 2014, 126, 417–423. [Google Scholar] [CrossRef]
  21. Huffaker, T.B.; Hu, R.; Runtsch, M.C.; Bake, E.; Chen, X.; Zhao, J.; Round, J.L.; Baltimore, D.; O’Connell, R.M. Epistasis between MicroRNAs 155 and 146a during T Cell-Mediated Antitumor Immunity. Cell Rep. 2012, 2, 1697–1709. [Google Scholar] [CrossRef] [Green Version]
  22. Torri, A.; Carpi, D.; Bulgheroni, E.; Crosti, M.C.; Moro, M.; Gruarin, P.; Rossi, R.L.; Rossetti, G.; Vizio, D.D.; Hoxha, M.; et al. Extracellular MicroRNA signature of human helper T cell subsets in health and autoimmunity. J. Biol. Chem. 2017, 292, 2903–2915. [Google Scholar] [CrossRef] [Green Version]
  23. Gracias, D.T.; Stelekati, E.; Hope, J.L.; Boesteanu, A.C.; Doering, T.A.; Norton, J.; Mueller, Y.M.; Fraietta, J.A.; Wherry, E.J.; Turner, M.; et al. The microRNA miR-155 controls CD8+ T cell responses by regulating interferon signaling. Nat. Immunol. 2013, 14, 593–602. [Google Scholar] [CrossRef] [Green Version]
  24. Singh, U.P.; Murphy, A.E.; Enos, R.T.; Shamran, H.A.; Singh, N.P.; Guan, H.; Hegde, V.L.; Fan, D.; Price, R.L.; Taub, D.D.; et al. miR-155 deficiency protects mice from experimental colitis by reducing T helper type 1/type 17 responses. Immunology 2014, 143, 478–489. [Google Scholar] [CrossRef]
  25. Stelzer, G.; Rosen, N.; Plaschkes, I.; Zimmerman, S.; Twik, M.; Fishilevich, S.; Iny Stein, T.; Nudel, R.; Lieder, I.; Mazor, Y.; et al. The GeneCards suite: From gene data mining to disease genome sequence analyses. Curr. Protoc. Bioinform. 2016, 2016, 1.30.1–1.30.33. [Google Scholar] [CrossRef] [PubMed]
  26. De Rie, D.; Abugessaisa, I.; Alam, T.; Arner, E.; Arner, P.; Ashoor, H.; Åström, G.; Babina, M.; Bertin, N.; Burroughs, A.M.; et al. An integrated expression atlas of miRNAs and their promoters in human and mouse. Nat. Biotechnol. 2017, 35, 872–878. [Google Scholar] [CrossRef] [PubMed]
  27. Iversen, R.; Roy, B.; Stamnaes, J.; Høydahl, L.S.; Hnida, K.; Neumann, R.S.; Korponay-Szabó, I.R.; Lundin, K.E.A.; Sollid, L.M. Efficient T cell–B cell collaboration guides autoantibody epitope bias and onset of celiac disease. Proc. Natl. Acad. Sci. USA 2019, 116, 15134–15139. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  28. Jabri, B.; Sollid, L.M. T Cells in Celiac Disease. J. Immunol. 2017, 198, 3005–3014. [Google Scholar] [CrossRef]
  29. Abadie, V.; Kim, S.M.; Lejeune, T.; Palanski, B.A.; Ernest, J.D.; Tastet, O.; Voisine, J.; Discepolo, V.; Marietta, E.V.; Hawash, M.B.F.; et al. IL-15, gluten and HLA-DQ8 drive tissue destruction in coeliac disease. Nature 2020, 578, 600–604. [Google Scholar] [CrossRef]
  30. Briskin, D.; Wang, P.Y.; Bartel, D.P. The biochemical basis for the cooperative action of microRNAs. Proc. Natl. Acad. Sci. USA 2020, 117, 17764–17774. [Google Scholar] [CrossRef]
  31. Gil-Zamorano, J.; Tomé-Carneiro, J.; Lopez de las Hazas, M.C.; del Pozo-Acebo, L.; Crespo, M.C.; Gómez-Coronado, D.; Chapado, L.A.; Herrera, E.; Latasa, M.J.; Ruiz-Roso, M.B.; et al. Intestinal miRNAs regulated in response to dietary lipids. Sci. Rep. 2020, 10, 18921. [Google Scholar] [CrossRef]
  32. McKenna, L.B.; Schug, J.; Vourekas, A.; McKenna, J.B.; Bramswig, N.C.; Friedman, J.R.; Kaestner, K.H. MicroRNAs control intestinal epithelial differentiation, architecture, and barrier function. Gastroenterology 2010, 139, 1654–1664.e1. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  33. Beumer, J.; Clevers, H. Cell fate specification and differentiation in the adult mammalian intestine. Nat. Rev. Mol. Cell Biol. 2021, 22, 39–53. [Google Scholar] [CrossRef]
  34. Wang, B.; Rong, X.; Palladino, E.N.D.; Wang, J.; Fogelman, A.M.; Martín, M.G.; Alrefai, W.A.; Ford, D.A.; Tontonoz, P. Phospholipid Remodeling and Cholesterol Availability Regulate Intestinal Stemness and Tumorigenesis. Cell Stem Cell 2018, 22, 206–220.e4. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  35. Stine, R.R.; Sakers, A.P.; TeSlaa, T.; Kissig, M.; Stine, Z.E.; Kwon, C.W.; Cheng, L.; Lim, H.W.; Kaestner, K.H.; Rabinowitz, J.D.; et al. PRDM16 Maintains Homeostasis of the Intestinal Epithelium by Controlling Region-Specific Metabolism. Cell Stem Cell 2019, 25, 830–845.e8. [Google Scholar] [CrossRef] [PubMed]
  36. Tian, Y.; Ma, X.; Lv, C.; Sheng, X.; Li, X.; Zhao, R.; Song, Y.; Andl, T.; Plikus, M.V.; Sun, J.; et al. Stress responsive miR-31 is a major modulator of mouse intestinal stem cells during regeneration and tumorigenesis. Elife 2017, 6, e29538. [Google Scholar] [CrossRef] [PubMed]
  37. Tian, Y.; Xu, J.; Li, Y.; Zhao, R.; Du, S.; Lv, C.; Wu, W.; Liu, R.; Sheng, X.; Song, Y.; et al. MicroRNA-31 Reduces Inflammatory Signaling and Promotes Regeneration in Colon Epithelium, and Delivery of Mimics in Microspheres Reduces Colitis in Mice. Gastroenterology 2019, 156, 2281–2296.e6. [Google Scholar] [CrossRef]
  38. Lukonin, I.; Serra, D.; Challet Meylan, L.; Volkmann, K.; Baaten, J.; Zhao, R.; Meeusen, S.; Colman, K.; Maurer, F.; Stadler, M.B.; et al. Phenotypic landscape of intestinal organoid regeneration. Nature 2020, 586, 275–280. [Google Scholar] [CrossRef]
  39. Liu, J.Y.; Seno, H.; Miletic, A.V.; Mills, J.C.; Swat, W.; Stappenbeck, T.S. Vav proteins are necessary for correct differentiation of mouse cecal and colonic enterocytes. J. Cell Sci. 2009, 122, 324–334. [Google Scholar] [CrossRef] [Green Version]
  40. Latreille, R.; Servant, R.; Darsigny, M.; Marcoux, S.; Jones, C.; Perreault, N.; Boudreau, F. Transcription factor CUX1 is required for intestinal epithelial wound healing and targets the VAV2-RAC1 Signalling complex. Biochim. Biophys. Acta Mol. Cell Res. 2017, 1864, 2347–2355. [Google Scholar] [CrossRef]
  41. Postema, M.M.; Grega-Larson, N.E.; Meenderink, L.M.; Tyska, M.J. PACSIN2-dependent apical endocytosis regulates the morphology of epithelial microvilli. Mol. Biol. Cell 2019, 30, 2515–2526. [Google Scholar] [CrossRef] [PubMed]
  42. Moffett, H.F.; Cartwright, A.N.R.; Kim, H.J.; Godec, J.; Pyrdol, J.; Äijö, T.; Martinez, G.J.; Rao, A.; Lu, J.; Golub, T.R.; et al. The microRNA miR-31 inhibits CD8+ T cell function in chronic viral infection. Nat. Immunol. 2017, 18, 791–799. [Google Scholar] [CrossRef] [Green Version]
  43. Ferguson, B.J.; Mansur, D.S.; Peters, N.E.; Ren, H.; Smith, G.L. DNA-PK is a DNA sensor for IRF-3-dependent innate immunity. eLife 2012, 1, e00047. [Google Scholar] [CrossRef] [PubMed]
  44. Yanai, H.; Ban, T.; Wang, Z.; Choi, M.K.; Kawamura, T.; Negishi, H.; Nakasato, M.; Lu, Y.; Hangai, S.; Koshiba, R.; et al. HMGB proteins function as universal sentinels for nucleic-acid-mediated innate immune responses. Nature 2009, 462, 99–103. [Google Scholar] [CrossRef] [PubMed]
  45. Palone, F.; Vitali, R.; Trovato, C.M.; Montuori, M.; Negroni, A.; Mallardo, S.; Stronati, L. Faecal high mobility group box 1 in children with celiac disease: A pilot study. Dig. Liver Dis. 2018, 50, 916–919. [Google Scholar] [CrossRef] [PubMed]
  46. Manti, S.; Cuppari, C.; Tardino, L.; Parisi, G.; Spina, M.; Salpietro, C.; Leonardi, S. HMGB1 as a new biomarker of celiac disease in children: A multicenter study. Nutrition 2017, 37, 18–21. [Google Scholar] [CrossRef] [PubMed]
  47. Boyapati, R.K.; Rossi, A.G.; Satsangi, J.; Ho, G.T. Gut mucosal DAMPs in IBD: From mechanisms to therapeutic implications. Mucosal Immunol. 2016, 9, 567–582. [Google Scholar] [CrossRef] [Green Version]
  48. Wang, Y.B.; Tan, B.; Mu, R.; Chang, Y.; Wu, M.; Tu, H.Q.; Zhang, Y.C.; Guo, S.S.; Qin, X.H.; Li, T.; et al. Ubiquitin-associated domain-containing ubiquitin regulatory X (UBX) protein UBXN1 is a negative regulator of nuclear factor κB (NF-κB) signaling. J. Biol. Chem. 2015, 290, 10395–10405. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  49. Wang, P.; Yang, L.; Cheng, G.; Yang, G.; Xu, Z.; You, F.; Sun, Q.; Lin, R.; Fikrig, E.; Sutton, R.E. UBXN1 Interferes with Rig-I-like Receptor-Mediated Antiviral Immune Response by Targeting MAVS. Cell Rep. 2013, 3, 1057–1070. [Google Scholar] [CrossRef] [Green Version]
  50. Koizumi, S.I.; Sasaki, D.; Hsieh, T.H.; Taira, N.; Arakaki, N.; Yamasaki, S.; Wang, K.; Sarkar, S.; Shirahata, H.; Miyagi, M.; et al. JunB regulates homeostasis and suppressive functions of effector regulatory T cells. Nat. Commun. 2018, 9, 5344. [Google Scholar] [CrossRef]
  51. Wheaton, J.D.; Ciofani, M. JunB Controls Intestinal Effector Programs in Regulatory T Cells. Front. Immunol. 2020, 11, 44. [Google Scholar] [CrossRef] [Green Version]
  52. Polioudakis, D.; Bhinge, A.A.; Killion, P.J.; Lee, B.K.; Abell, N.S.; Iyer, V.R. A Myc-microRNA network promotes exit from quiescence by suppressing the interferon response and cell-cycle arrest genes. Nucleic Acids Res. 2013, 41, 2239–2254. [Google Scholar] [CrossRef]
  53. Alsaleh, G.; François, A.; Philippe, L.; Gong, Y.Z.; Bahram, S.; Cetin, S.; Pfeffer, S.; Gottenberg, J.E.; Wachsmann, D.; Georgel, P.; et al. MiR-30a-3p negatively regulates BAFF synthesis in systemic sclerosis and rheumatoid arthritis fibroblasts. PLoS ONE 2014, 9, e111266. [Google Scholar] [CrossRef]
  54. de Almeida, R.C.; Chagas, V.S.; Castro, M.A.A.; Petzl-Erler, M.L. Integrative analysis identifies genetic variants associated with autoimmune diseases affecting putative microRNA binding sites. Front. Genet. 2018, 9, 139. [Google Scholar] [CrossRef] [Green Version]
  55. Kern, F.; Backes, C.; Hirsch, P.; Fehlmann, T.; Hart, M.; Meese, E.; Keller, A. What’s the target: Understanding two decades of in silico microRNA-target prediction. Brief. Bioinform. 2020, 21, 1999–2010. [Google Scholar] [CrossRef]
  56. Mittelbrunn, M.; Gutiérrez-Vázquez, C.; Villarroya-Beltri, C.; González, S.; Sánchez-Cabo, F.; González, M.Á.; Bernad, A.; Sánchez-Madrid, F. Unidirectional transfer of microRNA-loaded exosomes from T cells to antigen-presenting cells. Nat. Commun. 2011, 2, 282. [Google Scholar] [CrossRef] [Green Version]
  57. Montecalvo, A.; Larregina, A.T.; Shufesky, W.J.; Stolz, D.B.; Sullivan, M.L.G.; Karlsson, J.M.; Baty, C.J.; Gibson, G.A.; Erdos, G.; Wang, Z.; et al. Mechanism of transfer of functional microRNAs between mouse dendritic cells via exosomes. Blood 2012, 119, 756–766. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  58. Moerkens, R.; Mooiweer, J.; Withoff, S.; Wijmenga, C. Celiac disease-on-chip: Modeling a multifactorial disease in vitro. United Eur. Gastroenterol. J. 2019, 7, 467–476. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  59. Gambardella, G.; Carissimo, A.; Chen, A.; Cutillo, L.; Nowakowski, T.J.; Di Bernardo, D.; Blelloch, R. The impact of microRNAs on transcriptional heterogeneity and gene co-expression across single embryonic stem cells. Nat. Commun. 2017, 8, 14126. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  60. Wang, N.; Zheng, J.; Chen, Z.; Liu, Y.; Dura, B.; Kwak, M.; Xavier-Ferrucio, J.; Lu, Y.C.; Zhang, M.; Roden, C.; et al. Single-cell microRNA-mRNA co-sequencing reveals non-genetic heterogeneity and mechanisms of microRNA regulation. Nat. Commun. 2019, 10, 95. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  61. Amr, K.S.; Bayoumi, F.S.; Eissa, E.; Abu-Zekry, M. Circulating microRNAs as potential non-invasive biomarkers in pediatric patients with celiac disease. Eur. Ann. Allergy Clin. Immunol. 2019, 51, 159–164. [Google Scholar] [CrossRef] [Green Version]
  62. Zorro, M.M.; Aguirre-Gamboa, R.; Mayassi, T.; Ciszewski, C.; Barisani, D.; Hu, S.; Weersma, R.K.; Withoff, S.; Li, Y.; Wijmenga, C.; et al. Tissue alarmins and adaptive cytokine induce dynamic and distinct transcriptional responses in tissue-resident intraepithelial cytotoxic T lymphocytes. J. Autoimmun. 2020, 108, 102422. [Google Scholar] [CrossRef]
  63. Love, M.I.; Huber, W.; Anders, S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014, 15, 550. [Google Scholar] [CrossRef] [Green Version]
  64. Benjamini, Y.; Hochberg, Y. Controlling the false discovery rate: A practical and powerful approach to multiple testing. J. R. Stat. Soc. Ser. B 1995, 57, 289–300. [Google Scholar] [CrossRef]
  65. Agarwal, V.; Bell, G.W.; Nam, J.W.; Bartel, D.P. Predicting effective microRNA target sites in mammalian mRNAs. eLife 2015, 4, e05005. [Google Scholar] [CrossRef] [PubMed]
  66. Paraskevopoulou, M.D.; Georgakilas, G.; Kostoulas, N.; Vlachos, I.S.; Vergoulis, T.; Reczko, M.; Filippidis, C.; Dalamagas, T.; Hatzigeorgiou, A.G. DIANA-microT web server v5.0: Service integration into miRNA functional analysis workflows. Nucleic Acids Res. 2013, 41, W169–W173. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  67. Karagkouni, D.; Paraskevopoulou, M.D.; Chatzopoulos, S.; Vlachos, I.S.; Tastsoglou, S.; Kanellos, I.; Papadimitriou, D.; Kavakiotis, I.; Maniou, S.; Skoufos, G.; et al. DIANA-TarBase v8: A decade-long collection of experimentally supported miRNA-gene interactions. Nucleic Acids Res. 2018, 46, D239–D245. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  68. Chou, C.H.; Chang, N.W.; Shrestha, S.; Hsu, S.D.; Lin, Y.L.; Lee, W.H.; Yang, C.D.; Hong, H.C.; Wei, T.Y.; Tu, S.J.; et al. miRTarBase 2016: Updates to the experimentally validated miRNA-target interactions database. Nucleic Acids Res. 2016, 44, D239–D247. [Google Scholar] [CrossRef] [PubMed]
  69. Castro, M.A.A.; Wang, X.; Fletcher, M.N.C.; Meyer, K.B.; Markowetz, F. RedeR: R/Bioconductor package for representing modular structures, nested networks and multiple levels of hierarchical associations. Genome Biol. 2012, 13, R29. [Google Scholar] [CrossRef] [Green Version]
  70. Yu, G.; Wang, L.G.; Han, Y.; He, Q.Y. ClusterProfiler: An R package for comparing biological themes among gene clusters. Omi. J. Integr. Biol. 2012, 16, 284–287. [Google Scholar] [CrossRef]
  71. Supek, F.; Bošnjak, M.; Škunca, N.; Šmuc, T. Revigo summarizes and visualizes long lists of gene ontology terms. PLoS ONE 2011, 6, e21800. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  72. Aran, D. Cell-type enrichment analysis of bulk transcriptomes using xcell. In Methods in Molecular Biology; Humana: New York, NY, USA, 2020; Volume 2120, pp. 263–276. [Google Scholar]
  73. Kim-Hellmuth, S.; Aguet, F.; Oliva, M.; Muñoz-Aguirre, M.; Kasela, S.; Wucher, V.; Castel, S.E.; Hamel, A.R.; Viñuela, A.; Roberts, A.L.; et al. Cell type–specific genetic regulation of gene expression across human tissues. Science 2020, 369, eaaz8528. [Google Scholar] [CrossRef] [PubMed]
  74. Marderstein, A.R.; Uppal, M.; Verma, A.; Bhinder, B.; Tayyebi, Z.; Mezey, J.; Clark, A.G.; Elemento, O. Demographic and genetic factors influence the abundance of infiltrating immune cells in human tissues. Nat. Commun. 2020, 11, 2213–2214. [Google Scholar] [CrossRef] [PubMed]
Ijms 22 11382 g001 550
Figure 1. miRNA and gene expression profiles were generated for small-intestinal biopsies using a next generation sequencing approach. Principal component analysis using (A) the miRNA-seq profile and (B) the mRNA-seq analysis shows clear separation of controls (black) and active celiac disease (CeD) patients (grey). Differential expression analysis was performed between CeD and controls (FDR = p-value adjusted for multiple testing). For all the differentially expressed miRNAs and genes, we extracted previously described miRNA-target transcript interactions. In total, 2030 miRNA–transcript target pairs remained with a negative Pearson’s correlation R < −0.7.
Figure 1. miRNA and gene expression profiles were generated for small-intestinal biopsies using a next generation sequencing approach. Principal component analysis using (A) the miRNA-seq profile and (B) the mRNA-seq analysis shows clear separation of controls (black) and active celiac disease (CeD) patients (grey). Differential expression analysis was performed between CeD and controls (FDR = p-value adjusted for multiple testing). For all the differentially expressed miRNAs and genes, we extracted previously described miRNA-target transcript interactions. In total, 2030 miRNA–transcript target pairs remained with a negative Pearson’s correlation R < −0.7.
Ijms 22 11382 g001
Ijms 22 11382 g002 550
Figure 2. miRNA–target transcript network in small intestinal biopsies. This network shows the 2030 miRNA–target transcript pairs that were selected by the workflow in Figure 1. Only the gene names of experimentally validated target transcripts are shown. Supporting data for the miRNA-mRNA pairs are shown in Supplementary Table S1.
Figure 2. miRNA–target transcript network in small intestinal biopsies. This network shows the 2030 miRNA–target transcript pairs that were selected by the workflow in Figure 1. Only the gene names of experimentally validated target transcripts are shown. Supporting data for the miRNA-mRNA pairs are shown in Supplementary Table S1.
Ijms 22 11382 g002
Ijms 22 11382 g003 550
Figure 3. Top 30 enriched GO pathways associated with the down- or upregulated genes in CeD that are targeted by at least one miRNA. Sorted by how many miRNAs targeted the transcripts in the pathway. The colorscale indicates significance of pathway enrichment as indicated with p.adjust.
Figure 3. Top 30 enriched GO pathways associated with the down- or upregulated genes in CeD that are targeted by at least one miRNA. Sorted by how many miRNAs targeted the transcripts in the pathway. The colorscale indicates significance of pathway enrichment as indicated with p.adjust.
Ijms 22 11382 g003
Ijms 22 11382 g004 550
Figure 4. Enrichment of the target genes associated with the individual miRNAs in the miRNA–target transcript network. Y-axis shows representative GO terms that summarize the groups of GO terms selected by REVIGO (input: all the GO terms of the miRNAs combined). miRNAs are sorted from lowest fold change (miR-31) to highest fold change (miR-1260b) in the CeD versus controls comparison. The significance of the enrichment is indicated in FDR by the colorscale and the GeneRatio is indicative of the ratio of genes targeted by the miRNA and the total number of genes in the pathway.
Figure 4. Enrichment of the target genes associated with the individual miRNAs in the miRNA–target transcript network. Y-axis shows representative GO terms that summarize the groups of GO terms selected by REVIGO (input: all the GO terms of the miRNAs combined). miRNAs are sorted from lowest fold change (miR-31) to highest fold change (miR-1260b) in the CeD versus controls comparison. The significance of the enrichment is indicated in FDR by the colorscale and the GeneRatio is indicative of the ratio of genes targeted by the miRNA and the total number of genes in the pathway.
Ijms 22 11382 g004
Ijms 22 11382 g005 550
Figure 5. Cell type–enrichment based on Xcell enrichment scores calculated from the RNA-seq data (A,B) or based on the miRNA based cell type–enrichment score (C,D). (A) Immune cell enrichment scores show an increased trend in CeD. Epithelial cells show a lower enrichment in CeD. (B) Heatmap showing the (unscaled) Pearson’s R correlation coefficient between miRNA levels and the cell type–enrichment score. miRNAs that are downregulated in CeD (blue) show a high positive correlation with epithelial cells, whereas upregulated miRNAs show the opposite trend. (C) Enrichment for CD3+ associated miRNAs is increased in CeD. Enrichment for epithelial-specific miRNAs is lower in CeD than in controls. (D) Heatmap showing the (unscaled) Pearson’s R correlation coefficient between miRNA levels and the enrichment for CD3+ or intestinal epithelial cells. Mann-Whitney U p-value: n.s. >0.05, * 0.01–0.05, **** <0.0001.
Figure 5. Cell type–enrichment based on Xcell enrichment scores calculated from the RNA-seq data (A,B) or based on the miRNA based cell type–enrichment score (C,D). (A) Immune cell enrichment scores show an increased trend in CeD. Epithelial cells show a lower enrichment in CeD. (B) Heatmap showing the (unscaled) Pearson’s R correlation coefficient between miRNA levels and the cell type–enrichment score. miRNAs that are downregulated in CeD (blue) show a high positive correlation with epithelial cells, whereas upregulated miRNAs show the opposite trend. (C) Enrichment for CD3+ associated miRNAs is increased in CeD. Enrichment for epithelial-specific miRNAs is lower in CeD than in controls. (D) Heatmap showing the (unscaled) Pearson’s R correlation coefficient between miRNA levels and the enrichment for CD3+ or intestinal epithelial cells. Mann-Whitney U p-value: n.s. >0.05, * 0.01–0.05, **** <0.0001.
Ijms 22 11382 g005
Ijms 22 11382 g006 550
Figure 6. Overlap between genes of which the expression is affected by both celiac disease risk–SNPs and by miRNAs. (A) Single sample enrichment analysis reveals that the genes prioritized by van der Graaf et al. are enriched in active CeD compared to controls [9]. This analysis included the 102 genes (out of 118) with sufficient expression in the biopsy transcriptome dataset. Mann-Whitney U p-value: * 0.01–0.05. (B) Subset of the miRNA-target transcript interaction network. Only genes that overlap with the prioritized CeD genes are displayed. Thick connecting lines indicate that the interaction was previously experimentally validated. After correction for cell type composition (miRNA-based enrichment score for intestinal epithelial cells and CD3+ cells), the miRNAs indicated with asterisk still show a trend between CeD and controls (FDR: * 0.05–0.1, ** <0.05).
Figure 6. Overlap between genes of which the expression is affected by both celiac disease risk–SNPs and by miRNAs. (A) Single sample enrichment analysis reveals that the genes prioritized by van der Graaf et al. are enriched in active CeD compared to controls [9]. This analysis included the 102 genes (out of 118) with sufficient expression in the biopsy transcriptome dataset. Mann-Whitney U p-value: * 0.01–0.05. (B) Subset of the miRNA-target transcript interaction network. Only genes that overlap with the prioritized CeD genes are displayed. Thick connecting lines indicate that the interaction was previously experimentally validated. After correction for cell type composition (miRNA-based enrichment score for intestinal epithelial cells and CD3+ cells), the miRNAs indicated with asterisk still show a trend between CeD and controls (FDR: * 0.05–0.1, ** <0.05).
Ijms 22 11382 g006
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Tan, I.L.; Barisani, D.; Panceri, R.; Modderman, R.; Visschedijk, M.; Weersma, R.K.; Wijmenga, C.; Jonkers, I.; Coutinho de Almeida, R.; Withoff, S. A Combined mRNA- and miRNA-Sequencing Approach Reveals miRNAs as Potential Regulators of the Small Intestinal Transcriptome in Celiac Disease. Int. J. Mol. Sci. 2021, 22, 11382. https://doi.org/10.3390/ijms222111382

AMA Style

Tan IL, Barisani D, Panceri R, Modderman R, Visschedijk M, Weersma RK, Wijmenga C, Jonkers I, Coutinho de Almeida R, Withoff S. A Combined mRNA- and miRNA-Sequencing Approach Reveals miRNAs as Potential Regulators of the Small Intestinal Transcriptome in Celiac Disease. International Journal of Molecular Sciences. 2021; 22(21):11382. https://doi.org/10.3390/ijms222111382

Chicago/Turabian Style

Tan, Ineke Luise, Donatella Barisani, Roberto Panceri, Rutger Modderman, Marijn Visschedijk, Rinse K. Weersma, Cisca Wijmenga, Iris Jonkers, Rodrigo Coutinho de Almeida, and Sebo Withoff. 2021. "A Combined mRNA- and miRNA-Sequencing Approach Reveals miRNAs as Potential Regulators of the Small Intestinal Transcriptome in Celiac Disease" International Journal of Molecular Sciences 22, no. 21: 11382. https://doi.org/10.3390/ijms222111382

APA Style

Tan, I. L., Barisani, D., Panceri, R., Modderman, R., Visschedijk, M., Weersma, R. K., Wijmenga, C., Jonkers, I., Coutinho de Almeida, R., & Withoff, S. (2021). A Combined mRNA- and miRNA-Sequencing Approach Reveals miRNAs as Potential Regulators of the Small Intestinal Transcriptome in Celiac Disease. International Journal of Molecular Sciences, 22(21), 11382. https://doi.org/10.3390/ijms222111382

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop