Relation between NOD2 genotype and changes in innate signaling in Crohn’s disease on mRNA and miRNA levels

Crohn’s disease is associated with an altered innate immune response of pathogenic importance. This altered response can be associated to loss-of-function polymorphisms in the NOD2 (nucleotide-binding oligomerization domain-containing protein 2) gene, but also changes in transcriptional and post-transcriptional regulatory layers, including microRNA activity. Here, we characterized the link between NOD2 genotype and inflammatory-mediated changes in innate signaling by studying transcriptional and post-transcriptional activity in response to NOD2-agonist muramyl dipeptide in monocytes from healthy controls, and Crohn’s disease patients with and without NOD2 loss-of-function polymorphisms. We measured the expression of genes and microRNAs in monocytes from these subjects after stimulation with muramyl dipeptide. Gene expression profiles mainly distinguished the actual muramyl dipeptide response, but not the genotype. A hyper-responsive phenotype was found in Crohn’s disease patients without NOD2 mutations, characterized by upregulated cytokine receptors and general downregulation of microRNA expression. Conversely, microRNA expression could identify genotype-specific differences between subject groups but exhibited little change upon muramyl dipeptide treatment. Only two microRNAs showed muramyl dipeptide-induced response, including miR-155, which was found to regulate multiple genes and whose host gene was one of the highest muramyl dipeptide responders. miR-155 was upregulated in Crohn’s disease patients with NOD2 mutations following lipopolysaccharide and Escherichia coli treatment, but the upregulation was substantially reduced upon muramyl dipeptide treatment. While Crohn’s disease patients with NOD2 mutations on average showed a reduced muramyl dipeptide response, the cohort exhibited large individual variance: a small subset had inflammatory responses almost comparable to wild-type patients on both gene and miR-155 regulatory levels.


Isolation and stimulation of monocytes
Peripheral blood mononuclear cells were separated from whole blood by Ficoll-Paque density gradient centrifugation according manufacturer's instructions (GE Healthcare).

RNA purification from stimulated monocytes
For RNA-seq, small RNA-seq and validation analysis in monocytes, RNA were purified using the NucleoSpin®miRNA kit (Macherey-Nagel). The two step purification protocol was followed according to the manufacturer's instructions, with a final elution volume 40 ul.
Small RNA (< 200 bp) and large RNA (> 200 bp) were isolated in two separate fractions.
Purification of total RNA from transfected and stimulated THP-1 cells.
For gene expression validation in THP-1 cells, total RNA was purified using PureLink® RNA Mini Kit (Ambion™, Life Technologies) combined with an on-column DNase treatment according to the manufacturer's instructions using lysis buffer added 1%βmercaptoethanol and a final elution volume of 30 ul.

qPCR analysis of gene expression in monocytes and THP-1 cells
For validation of gene expression in monocytes and THP-1 cells, cDNA was synthesized using SuperScript III reverse transcriptase (Invitrogen) with a total input RNA of 200 ng for monocytes and 436 ng for THP-1 cells in a 20 µl reaction volume. All samples were run in triplicates in a 5 µl reaction volume consisting of 0.22 µl cDNA, 1.5 µl of primer mix (1.67 pmol/µl), 0.78 µl RNase-free water, and 2.5 µl Brilliant SYBR Green QPCR Master Mix (Invitrogen). PCR reactions were quantified using the QuantStudio 6 Flex Real-Time PCR System (Applied Biosystems). The thermal cycler details were as follows: 10 min at 95°C followed by amplification of the cDNA for 40 cycles with melting for 30 s at 95°C, annealing for 1 min at 60°C, and elongation for 45 s at 72°C.

qPCR data analysis of gene expression
The standard deviation (SD) between triplicates in qPCR experiments from both monocytes and THP-1 cells was required to be no more than 0.6. If the SD of triplicates exceeded 0.6, pairwise replicates were examined using the same threshold and the pairs with the minimum SD were kept for further processes. The samples that did not satisfy these criteria in triplicates or in any of the replicate pairs were eliminated from the analyses. For the qPCR validation focusing on MDP responses using paired samples, only the samples satisfy these criteria in both 'before MDP' and 'after MDP' states were shown in the resulting plots. All results were normalized to human RPLP0 (ribosomal protein large p0) by the 2 -ΔCt method 1 . In order to eliminate batch differences between the two rounds of miR-155 experiments in monocytes using MDP stimuli ( Figure 4A, right panel) and using LPS/E.coli stimuli ( Figure 4B), for each unstimulated CD WT and CD NOD2 subjects used in both experiments (5 for CD WT and 4 for CD NOD2 subjects, for details see Supplementary   Table 1), we calculated the qPCR-expression differences between experiments as per subject and the median !"#$! per subject groups was used for multiplying to the qPCR expression values derived from the later LPS/E.coli experiment.

qPCR analysis of miRNA expression in monocytes
For validation of miRNA expression, total RNA (10 ng) was reverse transcribed into cDNA using miRNA-specific primers (Applied Biosystems) and the TaqMan® miRNA reverse transcription kit (Applied Biosystems) according to the manufacturer's instructions. The expression levels of selected miRNAs were measured using commercially available predesigned miRNA-specific TaqMan® miRNA assays (Applied Biosystems) according to the manufacturer's recommendations. PCR reactions were quantified using the LightCycler 480 real-time PCR system (Roche). All samples were run in triplicates in a 5 µl reaction volume and cycles were run as previously described by Coskun et al. 2 . The same SD threshold 0.6 as used in gene qPCR analysis was employed for examining technical triplicates or replicate pairs. Expression levels of each miRNA were normalized to human RNU6B (RNA, U6B small nuclear). The relative expression levels of these miRNAs were determined by the 2 -ΔCt method 2 .

Mapping and Quantification of RNA-seq libraries
After filtering the reads for a minimum sequencing quality of 30 in 50% of the bases with the FASTX-Toolkit (http://hannonlab.cshl.edu/fastx_toolkit/), nucleotides with quality score lower than 20 were subsequently trimmed from 3'end of the RNA-seq reads. The remaining reads were mapped to the hg19 genome assembly using Tophat v2.0.

Mapping and Quantification of Small RNA-Seq
After truncating the 3'-ends adaptors and filtering the reads for a minimum sequencing quality of 30 in 50% of the bases with the FASTX-Toolkit, we used Bowtie v0.12.7 6 (settings: --strata --best -k 10 -3 30) to map small RNA-Seq to the hg19 reference assembly. Because the protocol applied did not distinguish strands, miRNA levels were quantified by the sum of mapped reads falling into the corresponding genomic coordinates of their precursors using UCSC gene annotation, regardless of strand. The lowest 40% expressed miRNAs (including miRNAs with no expression) were removed before running differential analysis in DESeq. The expression of miRNAs was visualized as TPM (tags per million mapped tags).

Association Between miRNAs and Their Potential Target Genes
We only analyzed miRNA-gene interactions supported by the published data in mirTarBase4.3 7 . We required miRNA and gene expression from each interaction to be reversely regulated in any of the differential expression pairs consisting of miRNAs and genes passing the cutoff from differential expression analysis (FDR ≤ 0.05) (arrows in Figure 1).

Gene Ontology Analysis and KEGG pathway analysis
GO analysis was performed using DAVID 8,9 using all expressed genes in all samples as background. GO terms (BP_FAT, CC_FAT and MF_FAT) were sorted by corresponding Bonferroni-adjusted P values. Gene enrichment analysis in KEGG 10 pathways was conducted using Pathview 11 using native KEGG graphs and default settings.

Pseudo counts and outlier truncation
For each analysis involving the use of log transformation when zero values existed in the analyzed dataset, the smallest non-zero value in the analyzed data points was employed as pseudo-count before log transformation. For the ease of visualization, in Figure 4F  Primer pairs used for gene expression validation.

Supplementary Table 3:
List of differentially expressed genes and miRNAs for respective group-wise analysis differentially expressed genes regardless of MDP treatment (defined in Figure 2F). Figure 3):

Supplementary Figure S3 (expanding
Panels A and B show the same analyses as in Figure 3A-B, but use an alternative method for defining differentially expressed gene sets (FDR ≤ 0.05 from Cuffdiff combined with a non-adjusted P-value ≤ 0.05 from DESeq).
A) Enriched significant GO terms (Bonferroni-corrected P values ≤ 0.05) for differential expressed genes between subject groups before MDP treatment. X-axis showslog 2 (Bonferroni-corrected P values) and Y-axis shows GO terms. A) Corresponding small RNA-seq data for the validated miRNA-155 shown in Figure 4A (right panel), expressed as TPM. Error bars indicate standard error of the mean.

B) As
B) Heat maps as in Figure 4C showing genes that are differential expressed upon MDP treatment within a subject group. Figure 4C (right panel).

C) As A, but showing RNA-seq data for the two validated genes shown in
D) The left panel shows the expression of MAVS as measured by qPCR following the conventions in Figure 3C. The right panel shows the corresponding RNA-seq data as in A. E) As D, but showing qPCR data and RNA-seq data for PALD1.