NGS study of glucocorticoid response genes in inflammatory bowel disease patients.

INTRODUCTION: Despite intensive research and a long history of glucocorticoids being applied in various clinical areas, they still generate a challenge for personalized medicine by causing resistance or dependence in nearly 50% of patients treated. The objective of the present study was to determine the genetic predictors of variable reactions in inflammatory bowel disease patients to glucocorticoid therapy. Therefore, based on the current knowledge on how glucocorticoids act, we have compiled a panel of 21 genes for variant analysis: NR3C1, NLRP1, IPO13, FKBP5, HSPA4, ABCB1, STIP1, HSP90AA1, IL-1A, IL-1B, IL-2, IL-4, CXCL8, IL-10, NFKBIA, JUN, MIF, TNF, MAPK14, CYP3A4, and CYP3A5.
MATERIAL AND METHODS: These genes were analyzed using the amplicon next-generation sequencing method in a group of 139 diagnosed and clinically characterized inflammatory bowel disease patients with a confirmed glucocorticoid response.
RESULTS: Analysis of all the targeted DNA sequences for the whole patient group indicated 121 different functional variants. After association analyses of 31 selected variants, the polymorphism c.1088A>G in the NR3C1 gene was linked with glucocorticoid resistance (p = 0.002), variant c.241+6A>G of the FKBP5 gene with glucocorticoid sensitivity (p = 0.040), and deletion c.306-7delT in the MAPK14 gene with an adverse therapeutic effect (dependency and resistance, p = 0.041) in ulcerative colitis patients. In Crohn's disease, the change c.2685+49T>C of the ABCB1 gene related to glucocorticoid resistance (p = 0.034).
CONCLUSIONS: Among the 21 analyzed genes, four (NR3C1, FKBP5, MAPK14, and ABCB1) revealed a significant impact on the glucocorticoid treatment response, which could result in valuable pharmacogenetic biomarkers after being confirmed in other populations and in functional studies. Copyright:

Introduction

Crohn’s disease (CD) and ulcerative colitis (UC) are inflammatory disorders of the gastrointestinal tract known as inflammatory bowel disease (IBD). In UC the inflammation involves the distal intestine and ileum, and in CD it can affect the whole region of the digestive tract. The development of IBD with advanced symptoms frequently begins at the age of 20–30, and due to the chronic nature of the disorder, patients require long-term treatment, usually for their whole life [1, 2]. Synthetic glucocorticosteroids (GCs) are among the most commonly used first-line conventional drugs, dedicated to remission induction in IBD patients with mild and severe exacerbations [3]. The GCs are widely known for their rapid anti-inflammatory and immunomodulatory effect, which is followed by binding the intracellular GC receptor (GR) and potent inhibition of T cell activation and cytokine secretion. Effective suppression of inflammatory gene expression occurs through several distinct mechanisms [4, 5]. However, not all patients show a positive response to GC therapy. On average, only every second adult patient with IBD is sensitive to GCs, close to 20% of subjects are GC resistant, and approximately 30–40% of patients become dependent on these drugs [6]. Personalized medicine could be very promising in CD and UC therapy outcome, but except for thiopurine drugs with pharmacogenetics dosing guidelines, we do not have genetic biomarkers for IBD treatment [7, 8].

Although the reason for the differential response is multifactorial, including disease severity and complications, as well as environmental factors, it can be expected that the genetic component plays the most crucial role, which is estimated at up to 95% [9, 10]. To date, most pharmacogenetic research on GCs has indicated the NR3C1 gene coding for the GCs receptor (GR) protein and the ABCB1 gene responsible for synthesis of the multidrug resistance protein 1 – a membrane GC efflux transporter whose overexpression results in decreased cytoplasmic GC concentration [11, 12]. Fewer GC response studies also concern the participation of the FKBP5, TNF, and NLRP1 genes coding for co-chaperone 51 kDa FK506-binding protein (FKBP-51), tumor necrosis factor (TNF) and NACHT, LRR and PYD domain-containing protein 1, respectively [13-16]. However, more molecules are engaged in GCs’ mechanism of action, whose polymorphism could relate to inter-individual differences observed in drug efficacy and side effects. These include: 1) genes coding for a GR cellular protein complex consisting of chaperone heat shock protein HSP 90-α and HSP 70 kDa, co-chaperones, including stress-induced-phosphoprotein 1 (STI1) and the immunophilin FKBP-51, which bind, stabilize and organize hormone-free GR in the cytosol; 2) factors involved in the hormonal control of target gene transcription in the nucleus, i.e., transcription factor AP-1, nuclear factor NF-κB p105 subunit, and of histone acetylation, as well as chromatin remodeling, 3) principal effectors of GCs, e.g. cytokines and their regulators: tumor necrosis factor (TNF), interleukin-1β (IL-1β), interleukin-2 (IL-2), interleukin-8 (IL-8) [4, 6, 17]. Polymorphisms in genes involved in the nuclear transport factors (importin IPO13) or in the metabolism of these hormones, i.e. cytochromes P450 (CYP3A4, CYP3A5), have also been suggested as other possible candidates of GC response modulators [11]. On the one hand, science emphasizes the role of genetics in predicting GCs’ clinical reactions and on the other hand, practice indicates the need of personalized medicine for improving the GC efficacy and safety rates [9, 13].

In this study, we aimed to determine the genetic predictors of GCs’ clinical response within a selected group of 21 genes coding for GR and their complex proteins, transcription factors, GC target genes, GC biotransformation enzymes and GC transporters in a group of 139 diagnosed and clinically characterized Polish IBD patients with a confirmed glucocorticoid response. The goal of these investigations is to better understand diverse patient reactions to GCs and risk factors. Our findings, which identify potential new variants associated with GC reactions, provide the basis for further research, including therapy in different diseases, various populations, and functional studies. In this work, we focused on the relationship between altered GC response and functional variants of the following genes: NR3C1, NLRP1, IPO13, FKBP5, HSPA4, ABCB1, STIP1, HSP90AA1, IL-1A, IL-1B, IL-2, IL-4, CXCL8, IL-10, NFKBIA, JUN, MIF, TNF, MAPK14, CYP3A4, and CYP3A5. A novel method, based on next-generation sequencing (NGS) of amplicon libraries, was developed for the purpose of analysis of candidate genes in the current study.

Material and methods

Patients and samples

One hundred and thirty-nine Polish patients with IBD hospitalized in the Department of Gastroenterology, Dietetics and Internal Medicine, Poznan University of Medical Sciences (Poland) were enrolled in this study between May 2013 and February 2018. This group included 77 (55.4%) individuals with Crohn’s disease (CD) and 62 (44.6%) individuals with ulcerative colitis (UC). All of them were of Caucasian ancestry. The characteristics of the patients are shown in Table I [18]. In addition to the response to GC treatment defined in all patients, in over half of the study group (51.1%) side effects were observed, including weakening, osteoporosis/osteopenia, peptic ulcers, stretch marks, acne, difficult wound healing, petechiae, edema, hypertension, thrombosis, infection, headache, glucose metabolism disorder, cataract, glaucoma, mental disorders, and increased weight and appetite.

Table I

Clinical characteristics of patients included in the study

Parameter	Total patients (n = 139)	GC-sensitive (n = 77)#	GC-dependent (n = 38)#	GC-resistant (n = 24)#
Age at study, mean ± SD	36.07 ±14.08	38.2 ±15.63	36.5 ±15.61	33.5 ±10.29
Age at diagnosis, mean ± SD	31.3 ±14.31	32.5 ±14.94	30.5 ±14.02	30.5 ±14.02
Sex, n (%):
Female	58 (41.7)	31 (40.3)	17 (44.7)	10 (41.7)
Male	81 (58.3)	46 (59.7)	21 (55.3)	14 (58.3)
Diagnosis, n (%):
CD	77 (55.4)	40 (52.0)	26 (68.4)	11 (45.8)
UC	62 (44.6)	37 (48.0)	12 (31.6)	13 (54.2)
Disease severity, n (%):
Mild	8 (5.8)	4 (5.2)	1 (2.6)	3 (12.5)
Moderate	72 (51.8)	42 (54.5)	25 (65.8)	5 (20.8)
Severe	59 (42.4)	31 (40.3)	12 (31.6)	16 (66.7)
Intestinal location*, n (%):
L1, ileal	24 (17.3)	13 (16.9)	8 (21.1)	3 (12.5)
L2, colonic	9 (6.5)	4 (5.2)	5 (13.2)	0 (0.0)
L3, ileocolonic	44 (31.7)	24 (31.2)	13 (34.2)	7 (29.2)
E1, proctitis	5 (3.6)	4 (5.2)	1 (2.6)	0 (0.0)
E2, left-side	10 (7.1)	8 (10.4)	2 (5.2)	0 (0.0)
E3, pancolitis	47 (33.8)	24 (31.2)	9 (23.7)	14 (58.3)
Side effects occurrence of GCs, n (%)	71 (51.1)	31 (40.3)	25 (65.8)	15 (62.5)

According to response classification as detailed in ‘Material and Methods, Patients and samples’ description.

Disease location was determined according to the Montreal classification [18]. GC – glucocorticoid, SD – standard deviation, CD – Crohn’s disease, UC – ulcerative colitis, n – number.

The inclusion criteria were as follows: 1) age over 18 years; 2) diagnosis of IBD based on clinical symptoms, endoscopy with histopathological confirmation (Figure 1), radiological images, and laboratory tests (including complete blood count with smear, hematocrit, erythrocyte sedimentation rate, C-reactive protein, iron, ferritin); 3) disease duration over 6 months and 4) at least one episode of treatment with GCs. All individuals gave their written consent to genetic testing and evaluation of biochemical parameters in serum and colonoscopy. The study was approved by the Bioethical Committee of Poznan University of Medical Sciences under Resolution No. 436/13.

Figure 1

Histopathology of colon mucosa in IBD patients. Patient with UC: A – crypt abscess (blue arrow), B – thinned intestinal epithelium (blue arrow). Patient with CD: C, D – granulomas (blue arrows). Stain: hematoxylin and eosin, magnification 200× (A, C, D), 400× (B)

In all patients, the clinical data were obtained with a minimum of 1 year of continuous follow-up evaluation after the start of GC treatment. Some patients were evaluated retrospectively and the rest prospectively. The disease activity was determined according to the Crohn’s Disease Activity Index (CDAI) for patients with CD and according to the Truelove-Witts score for patients with UC. GC therapy consisted mostly of methylprednisolone and hydrocortisone, and in rare cases, budesonide (Table II). The dose of GCs (methylprednisolone and hydrocortisone) was converted to prednisolone using a corticosteroid conversion calculator (https://clincalc.com/corticosteroids/). According to the existing criteria of the European Crohn’s and Colitis Organization (ECCO), the patients were assigned to one of three groups: 1) GC-sensitive, if after the initiation of treatment within a few days they achieved and maintained clinical remission; 2) GC-dependent, if the symptoms of exacerbation recurred after reducing the dose of GCS below 10 mg/24 h prednisolone or within 3 months after GCs were discontinued, 3) GC-resistant, if symptoms of exacerbation persisted despite 4 weeks of GC therapy in doses corresponding to a minimum of 0.75 mg/kg of prednisolone [1, 2]. The GCs were administered as a monotherapy or in combination therapy with 5-aminosalicylic acid preparations.

Table II

GCs – dosages in the study group

Drug	Dose [mg/24 h]		Number of patients, n (%)		Response, n (%)			Sex, n (%)		IBD diagnosis, n (%)
					GC-S	GC-D	GC-R	Female	Male	CD	UC
Methylprednisolone (oral)	Medium	12–24	127 (91.4)	45 (32.4)	31 (22.3)	11 (7.9)	3 (2.2)	21 (15.1)	24 (17.27)	20 (14.4)	25 (18.0)
	High	> 24		82 (59.0)	36 (25.9)	26 (18.7)	2 (1.4)	28 (20.1)	54 (38.9)	35 (25.2)	47 (33.8)
Hydrocortisone (intravenous)	Low	≤ 150	103 (74.1)	4 (2.9)	4 (2.9)	0 (0.0)	0 (0.0)	3 (2.2)	1 (0.7)	1 (0.7)	3 (2.2)
	Medium	> 150–250		21 (15.1)	12 (8.6)	7 (5.0)	2 (1.4)	8 (5.8)	13 (9.4)	7 (5.0)	14 (10.1)
	High	> 250		78 (56.1)	38 (27.4)	22 (15.8)	18 (13.0)	31 (22.3)	47 (33.8)	31 (22.3)	47 (33.8)
Budesonide (oral)	Standard	9	18 (13.0)	18 (13.0)	9 (6.5)	3 (2.2)	6 (4.3)	12 (8.6)	6 (4.3)	18	0 (0.0)
P-value			–	–	0.072			0.469		< 0.001 (0.861*)

GC-S – glucocorticoid sensitivity, GC-D – glucocorticoid dependence, GC-R – glucocorticoid resistance, IBD – inflammatory bowel disease, CD – Crohn’s disease, UC – ulcerative colitis, n – number.

Excluding the budesonide treatment in the statistical analysis, which in the course of this research was available in Poland only for CD patients’ treatment.

Moreover, in total 12 healthy subjects (6 women and 6 men) without drug treatment ranging in age from 21 to 54 years old were used in this study as controls in all molecular, biochemistry and immunohistochemistry methods. Exclusion criteria in the control group included chronic disease or any organ dysfunction, ongoing infections, excessive alcohol consumption, cigarette smoking and receiving any supplementation influencing hepatic and kidney enzymes.

Long-range (LR)-PCR amplification

Genomic DNA from all subjects was isolated from peripheral blood according to standard procedures using the method with guanidine isothiocyanate and stored at 4°C in a TE buffer containing 1.0 mM EDTA and 10 mM Tris-Cl until use. The amplification of the NR3C1, NLRP1, IPO13, FKBP5, HSPA4, ABCB1, STIP1, HSP90AA1, IL-1A, IL-1B, IL-2, IL-4, CXCL8, IL-10, NFKBIA, JUN, MIF, TNF, MAPK14, CYP3A4, and CYP3A5 genes was carried out using primers containing exons, splice junctions and promoters, as well as 5′ and 3′ flanking sequences (Table III). The LR-PCR conditions that were established are presented in Table IV. Primers were designed as part of this work using the NCBI (National Center for Biotechnology Information) Primer-BLAST tool in reference to the human genomic sequence (GRCh37/hg19). The only exceptions were the primers and amplification conditions for the IL-1B, TNF and ABCB1 genes, which originated from our previous studies [19-21]. Due to the limitation of the LR-PCR method in the length of efficiently amplified DNA fragments, the LR-PCR of short genes such as IL-1A (11.48 kb), IL-1B (7.02 kb) and IL-2 (5.25 kb) was performed in one fragment, but the long genes, e.g. NR3C1 (123.76 kb), FKBP5 (115.35 kb) or ABCB1 (210.00 kb), were divided into several fragments (from 2 to 6 parts, Table III). Therefore, the amplification of the DNA sequence of 21 genes generates 54 amplicons. The numbering of the amplified fragments in Tables III and IV is identical (and also corresponds with the lane numbering in Figure 2).

Table III

Characteristics of amplicons and primer sequences for NGS

Frag. no.	Gene	OMIM	Region	Amplicon [bp]	Sequence (5′→3′) primer forward	Sequence (5′→3′) primer reverse
1	NR3C1	138040	5′UTR–prom. 1	2400	ATAAAAAGGCAACAAGTTTCTAGGCG	CACACTTGCCCATCTACAGTTAATTC
2			Exons 1–2, prom. 2–7	6301	AAAGTTACCTCCCCAAACCCTAAAC	TACTTGAGAACTTGCAGGAACATTTG
3			Exons 3–4	5039	AGGAAACATGGGTAGAATTGTAATGA	CACAGATCTCCAAGTAAAAGGTATGC
4			Exons 5–7	6208	CACAGTTTGTAATCCGTTACTCCAAG	TGCTATGTTAACCAATCCCCAATAGT
5			Exons 8–9–3′UTR	2261	TTTGTTATCCTGCATTTTCTCTTGGC	AGTTAGGACTGTTAATTTGCACAACC
6	NLRP1	606636	5′UTR–exon 3	5287	AGTGGATAATTCTGGCTTCTAGACTG	CTCCCTGTGTAAGAGTTAAAAGTGTG
7			Exons 4–5	7575	GACTGAGAGTTGTGTCTGTTTTGTTC	GTCCAAACTATTACCATCTTTGCAGG
8			Exons 6–12	12504	GAGAAAAGTAAAGCAGAAAGGAGCAG	AACATTCTTCGTTAGTTCACTTAGCC
9			Exons 13–15	5023	TAATCACCTCCCATAAAAGCTTCAC	TCCTGTGGGCTATTTCTTCTCATCAA
10			Exon 16–3′UTR	1088	TTAACATCCTCTCATTCCACTTTTC	CATCTCGCATATGATTAAGAACAAAT
11	IPO13	610411	5′UTR – exons 1-2	4357	CTCGAGTCACGTCGTCCTACTAAG	CTTAGCCTTAAGAATCCATGGAGGA
12			Exons 3–14	6294	ATCACTCCACTCTCTACTTTCATCC	CTAGAAAGTACTAAGCCTTGATGCAG
13			Exons 15–20, 3′UTR	4762	AAACCTCAAAAGGAGATACCAAGGAC	CTGCTAAGGGAGAAAAGCAAGAATAG
14	FKBP5	602623	5′UTR–prom. 1	2935	AAAGGAAAAGATTAAAGGCAGCTCAC	GGATTTCCCAAACCGTCCAAAAGG
15			Exon 1–prom. 2	2534	GAGACATGAGTTCTAATCTGGCTATG	TATGCTTTAACTTTTCCCTTGTCCAC
16			Exons 2–3	6516	ATCCTCTGTCAAGAATATGGCTGTAG	GGAAAATACCATGACTACTGGAACAG
17			Exons 4–5	2523	TGGACCCTTTTATGCAATTTATCCTG	CCACACTTGACTAATTCTCTTGATGC
18			Exons 6–7	7518	TGATGTTGTTGTCTTTGGATTTACCA	CTAGCTCAGGTATATGCTGTGTGTG
19			Exons 8–11, 3′UTR	12540	TCTCTTGACTCTGGTCCTTCTAAATC	ACTGTGGGTTCTGTTGGTTAAAAATG
20	HSPA4	601113	5′UTR–exon 1	1290	GACTCTTGAGAGAACCCGACATAAAG	ATTCAACTATTATTAAGCGCCTGCTG
21			Exons 2–7	13489	AATGTTCAGAAGTGGGGTTAAAGTTC	ACAATAACACTCTGGGTAGTTTCTGA
22			Exons 8–13	7051	CTCTTCTTCTTCCCTTTTGTCATTGG	GTATCAACCAGGAGAGAGTTTAAAGG
23			Exons 14–19, 3′UTR	9368	TCCATAGATATTTCATACCCTCCTCC	AATATGAAGAAGCTCAGAGCACACC
24	ABCB1	171050	5′UTR–exon 1	3518	CAGTTCAGAAACAGAACAGAACATGA	AATGGCTTTGGAACTCTTAAGAAAGG
25			Exons 2–5	15827	TTAGGAAGCAGAAAGGTGATACAGA	CTTGTTTTTGCTGCAAGTTTCCAATA
26			Exons 6–10	16809	GATATTTCTGCTGACAGACTACCTCT	TTCTAAAGTCAAGCCAACATTACTGG
27			Exons 11–21	14483	CTGATGGTTTTTCTTCACATTCCTCA	AAGGAGAAAATTAGTTTCATGCTGGG
28			Exons 22–26	17000	TCATTCCTAGTTTGTCAGACTCCTTT	CTCTCTATACCATCTGTCATTGCTGA
29			Exons 27–29, 3′UTR	5967	TAATTGTGCTACATTCAAAGTGTGCT	ACCCCACCTCCTAAAATCTTATATCG
30	STIP1	605063	5′UTR–exon 1	2224	CAGCACTTTCAAGTGAGGGATTAC	CACACACACAAGCATACAAAAAGAGG
31			Exons 2–8	6446	TGCACTTTATATGTTTGAGAAGGGAG	TTGAAATCTTGGCCAATACCAAAAAC
32			Exons 9–14, 3′UTR	5400	AATCTTGCCAGCAATACTCTGAAAAT	TCCATGACCCAACTCCTAAATCTAAC
33	HSP90AA1	140571	5′UTR–exon 1	1404	TTGTAGATATAGAAGTACAGCTCCCG	TTAGAATGTTTCTGTCACTTTAGAAGC
34			Exon 2	1238	TCGTGTCATTGATGCATTTAGCTTAC	TCAAGCCCAAACTCCTCAGCATAG
35			Exons 3–12, 3′UTR	5157	CCGCCCTTCTATTTTCGGTTTACTTA	TACAGCTCAGAACACTTCAATAACAG
36	IL1A	147760	Whole gene	12577	ATTCATTTGCTAAGAGTCTGGTGTTC	CTTCTCCCTTTGACATCACTTTTAGG
37	IL1B	147720	Whole gene	10196	TCTACATCATCATCACTCTTCCACTC	GGATATTTCTAGGTACTTTGCAGCAC
38	IL2	147680	Whole gene	6505	CTGTACTGAGGATGAATGAACACCTA	TTAAAGGGACAAGGGAAGATTACCAG
39	IL4	147780	Whole gene	10003	TGAAACCTCAGAATAGACCTACCTTG	AATCCCTAAACCTTGTGTTCTTGTTC
40	CXCL8	146930	Whole gene	5000	TCTTTCCACTAAGCATACAACTTTCC	AACAATAGCCTATGAGCTAGAAGGTG
41	IL10	124092	Whole gene	6199	ATCCTAATGAAATCGGGGTAAAGGAG	CAGGATAACCCAGGCAATTTGTATTC
42	NFKBIA	164008	Whole gene	5016	CACAATGTAGGTCAGATAGCATAAAG	GGATGAAAACCTCAAGATAGAGCATC
43	JUN	165160	Whole gene	4730	GCAATACAAATCTCTCGGCTTCTAC	TACTGGGTACTAGAACACTGCTAATG
44	MIF	153620	Whole gene	2731	CTTCACTGAGACTGGTATATGGATTG	GGAGCTTCTATTCTCCTTCCTCTTC
45	TNF	191160	Whole gene	8831	TCCATCAGCTCATACCATTAGAGAAG	AACATAATCAAAAGAAGGCACAGAG
46	MAPK14	600289	5′UTR–exon 1	1910	GCATCTCAAAGAGTAACCCATATGTC	GGTATTTTAAGCTCTATGGATCGCTG
47			Exons 2–3	7515	CCTCATACAACCACAGAGTTGAGTAG	GCGATAAAAGAGGAGGGAGTCTATATG
48			Exons 4–8	4900	CATCGTGGTACAGCTACCATTTTAAG	TCCATAAGTAGAACTGATGAGAGTTCC
49			Exons 9–12–3′UTR	13300	CACTGTCATATCACAAGCTCTTTAGG	AACATATTCCTAAAAGGACAGCAGC
50	CYP3A4	124010	5′UTR – exon 3	10626	ATGAAAGTCATGAAAGTGAGACACCTA	ATAACCAGTGTGAATGTTGTTAGAACC
51			Exons 4–13, 3′UTR	15518	TATGAAACATTATGGAGAATGGCATGG	GTCACAGTGGATTAACTTTCACCTATG
52	CYP3A5	605325	5′UTR – exons 4	8829	GAATTTCGTGGGTAAAGATGTGTAGG	ACAATATGCCTACCTTTCCAAAATGG
53			Exons 5–10	7268	CTAGGTCTTGTTCCTGTCTATTCCAC	TCAGTTGGAAATGCAGAAATCATTCG
54			Exons 11–13, 3′UTR	4910	ACATTCATAACTCCTCCACACATCTC	CTGAGTCTACCTATCTGTCACTTACC

Frag. – fragment, bp – base pair, UTR – untranslated region, prom. – promoter.

Table IV

Conditions of LR-PCR reactions

Amplified fragments*	Reaction mixture	PCR program
1, 3–8, 10–13, 16, 17, 20–23, 25–29, 31–41, 45, 47, 48, 52–54	38.4 ng template DNA0.48 μl of 5 μmol/l primers10 μl of GoTaq Master MixWater up to 20 μl
9, 24, 50-51	50 ng template DNA0.36 μl of a 5 μmol/l primers15 μl of GoTaq Master MixWater up to 20 μl
18, 19	20 ng template DNA0.48 μl of a 5 μmol/l primers,4 μl of GC Enhancer10 μl of GoTaq Master MixWater up to 20 μl
43, 44	38.4 ng template DNA0.48 μl of 5 μmol/l primers10 μl of GoTaq Master MixWater up to 20 μl
2, 14, 15, 30, 42, 46, 49	20 ng template DNA0.48 μl of 5 μmol/l primers,4 μl of GC Enhancer [Jena Bioscience]10 μl of GoTaq Master,Water up to 20 μl

The fragment numbers: 1–54 are identical with Table III and Figure 3.

Figure 2

Agarose electrophoresis of all 54 LR-PCR products for one DNA sample. Lane: 1–5 NR3C1 (2400 bp, 6301 bp, 5031 bp, 6208 bp, 2261 bp), 6–10 – NLRP1 (5287 bp, 7575 bp, 12504 bp, 5023 bp, 1088 bp), 11–13 – IPO13 (4357 bp, 6294 bp, 4762 bp), 14–19 – FKBP5 (2935 bp, 2534 bp, 6516 bp, 2523 bp, 7518 bp, 12540 bp), 20–23 – HSPA4 (1290 bp, 13489 bp, 7051 bp, 9368 bp), 24–29 – ABCB1 (3518 bp, 15887 bp, 16809 bp, 14483 bp, 17000 bp, 5967 bp), 30–32 – STIP1 (2224 bp, 6456 bp, 5400 bp), 33–35 – HSP90AA1 (1404 bp, 1238 bp, 5157 bp), 36 – IL1A (12577 bp), 37 – IL1B (10196 bp), 38 – IL2 (6505 bp), 39 – IL4 (10003 bp), 40 – CXCL8 (5000 bp), 41 – IL10 (6199 bp), 42 – NFKBI (5016 bp), 43 – JUN (4730 bp), 44 – MIF (2731 bp), 45 – TNF (8831 bp), 46–49 – MAPK14 (1910 bp, 7515 bp, 4900 bp, 13300 bp), 50–51 – CYP3A4 (10626 bp, 15518 bp), 52–54 – CYP3A5 (8829 bp, 7268 bp, 4910 bp), 1 kb – marker ladder 1000 bp. The gel images were obtained by trimming and color adjusting of the full-length gels in the IrfanView 4.44 program. The lanes numbering corresponds with the numbering of the amplified fragments in Tables III and IV

Figure 3

The pipeline of NGS results analysis and variant selection

SNVs – single nucleotide variants, GQX – genotyping quality, MAF – minor allele frequency.

The amplification effect was controlled by electrophoresis in a 0.6% agarose gel (Sigma). To identify the approximate size of amplified fragments, the 1 kb DNA Ladder (Solis Biodyne) was used.

Library preparation, NGS and Sanger sequencing

Fifty-four amplicons of each patient were pooled in equimolar ratios. According to the manufacturer’s protocol, 1 ng of the pooled DNA fragments was subjected to the Nextera XT procedure (Illumina) using the Nextera XT DNA Sample Preparation Kit (Illumina) and Nextera XT Index Kit (96) (Illumina). Finally, we obtained 139 both-side indexed DNA libraries ready for high-throughput sequencing. Normalization of all libraries was carried out with magnetic beads, according to the Nextera XT procedure. Sequencing on an Illumina MiSeq System was performed as paired-end targeted resequencing using MiSeq Reagent Kit v2 (300 cycles) (Illumina). To verify the variants which were detected, Sanger sequencing for random samples was performed on an Applied Biosystems 3500 Genetic Analyzer using the BigDye Terminator v3.1 Cycle Sequencing Kit. The NGS analyses described were performed for the IBD patient group. However, for comparative purposes of the sequencing data, the minor allele frequency for the European population (MAF EU) from the 1000 Genomes Project was also used.

Bioinformatic and statistical analysis

Secondary bioinformatic analysis of NGS reads generated was performed using MiSeq Reporter analysis software (Illumina), and the resulting variants were characterized and filtered using VariantStudio Software v 3.0 (Illumina), according to the pipeline in Figure 3. In the first step, the NGS results concerning the sequences of 21 genes analyzed were aligned to the human reference sequence (version GRCh37/hg19). Next, the SNVs (single nucleotide variants) were detected and named, and we applied filtering parameters (GQX – genotyping quality ≥ 30, read depth ≥ 10 and heterozygous read ratio ≥ 35%) to select high-quality NGS results. The data obtained were filtered according to the potential consequence at the protein level and the frequency of the minor allele ≥ 2%. Thirty-one polymorphic variants were subjected to statistical analysis. We correlated the genetic results with GC response, including the division into three groups – GC-sensitive (S), -dependent (D) and -resistant (R) – for which the abbreviations S, D, and R, respectively, were used. This analysis aimed to compare all three GC-response groups in terms of genotype distribution for each polymorphic variant. Moreover, regarding GC response, the patients were also divided into two groups in two combinations: 1) correct responsive (GC-sensitive) versus poorly responsive subjects (GC-dependent and GC-resistant) abbreviated as S vs. D + R; as well as 2) GC-resistant versus the rest of the patients (GC-sensitive and GC-dependent) abbreviated as S + D vs. R. Those combinations were necessary to determine the genetic factors responsible for GC sensitivity and related to GC resistance. The calculations were carried out for the whole group of 139 patients, and for diagnosis of CD and UC.

The concordance of genotype distribution was verified using the Hardy-Weinberg equilibrium (HWE). All correlation analyses between GC response and genetic variants (Table V), as well as analysis of clinical data in subgroups (Table II), were verified using the Pearson χ2 and Fisher tests with Bonferroni correction.

Table V

Selected variants with potential GC therapy response significance

No.	Number rs	Gene	Nucleotide change	Protein change	NGS read depth#	HWE	MAF this study (%)	MAF EU (%)	GC response vs. genotype (p-value)
									IBD			CD			CU
									S vs. D vs. R	S vs. D + R	S + D vs. R	S vs. D vs. R	S vs. D + R	S + D vs. R	S vs. D vs. R	S vs. D + R	S + D vs. R
1	41423247	NR3C1	1184+646C>G	Regulatory feature	187	0.58	0.32	0.38	0.543	0.882	0.432	0.555	0.737	0.339	0.139	0.458	0.208
2	6190		68G>A	Arg23Lys	261	0.69	0.32	0.02	0.269	0.163	0.156	0.256	0.109	0.298	0.605	0.800	0.360
3	56149945		1088A>G	Asn363Ser	281	0.25	0.05	0.03	0.111	0.140	0.065	0.148	0.388	0.223	0.019	0.152	0.002
4	6196		2301T>C	Asn767Asn	213	0.37	0.14	0.17	0.276	0.167	0.470	0.212	0.386	0.160	0.695	0.330	0.299
5	12150220	NLRP1	464T>A	Leu155His	181	0.04	0.47	0.47	0.159	0.091	0.288	0.359	0.140	0.540	0.531	0.416	0.467
6	11651595		737C>G	Thr246Ser	136	0.24	0.08	0.05	0.712	0.347	0.743	0.571	0.249	0.894	0.614	0.551	0.771
7	52795654		2345C>G	Thr782Ser	180	0.24	0.08	0.05	0.712	0.347	0.743	0.571	0.249	0.894	0.614	0.551	0.771
8	11657747		2633C>T	Thr878Met	78	0.24	0.08	0.06	0.712	0.347	0.743	0.571	0.249	0.894	0.614	0.551	0.771
9	71368558		3052+197C>T	Unknown	49	0.14	0.08	0.05	0.703	0.436	0.633	0.812	0.466	0.796	0.614	0.551	0.771
10	2301582		3175G>A	Val1059Met	70	0.01	0.36	0.40	0.310	0.127	0.252	0.737	0.438	0.615	0.536	0.224	0.421
11	35596958		3355A>G	Met1119Val	84	0.24	0.08	0.05	0.712	0.347	0.743	0.571	0.249	0.894	0.614	0.551	0.771
12	11651270		3550A>G	Met1184Val	164	0.02	0.47	0.46	0.263	0.147	0.155	0.215	0.120	0.160	0.923	0.741	0.699
13	11653832		3721G>C	Val1241Leu	120	0.24	0.08	0.05	0.712	0.347	0.743	0.571	0.249	0.894	0.614	0.551	0.771
14	2817033	FKBP5	241+6A>G	Splice variant?	318	0.67	0.49	0.48	0.559	0.552	0.749	0.351	0.737	0.140	0.100	0.040	0.089
15	61745470	HSPA4	37T>A	Cys13Ser	443	0.73	0.03	0.02	0.902	0.677	0.713	0.387	0.168	0.559	0.659	0.614	0.785
16	61755724		475G>A	Ala159Thr	54	0.73	0.03	0.02	0.788	0.752	0.713	0.627	0.936	0.402	0.918	0.683	0.838
17	2229109	ABCB1	1199G>A	Ser400Asn	22	0.00	0.03	0.04	0.615	0.525	0.356	0.568	0.540	0.584	0.216	0.080	0.356
18	1128503		1236T>C	Gly412Gly	37	0.05	0.44	0.44	0.168	0.244	0.371	0.281	0.395	0.195	0.276	0.281	0.868
19	2032582		2677T>G	Ser893Ala	30	0.14	0.45	0.44	0.552	0.377	0.486	0.279	0.278	0.372	0.976	0.825	0.942
			2677T>A	Ser893Thr	65	0.72	0.01	0.02
20	2032583		2685+49T>C	Unknown	49	0.39	0.12	0.12	0.117	0.570	0.046	0.086	0.751	0.034	0.690	0.592	0.409
21	1045642		3435T>C	Ile1145Ile	62	0.32	0.45	0.45	0.131	0.107	0.529	0.749	0.670	0.751	0.146	0.082	0.749
22	11231722	STIP1	1245+5C>G	Splice variant?	110	0.04	0.36	0.40	0.730	0.406	0.554	0.966	0.770	0.884	0.687	0.339	0.506
23	8005905	HSP90AA1	211A>T	Met71Leu	405	0.59	0.11	0.11	0.411	0.528	0.526	0.375	0.510	0.163	0.542	0.672	0.595
24	17561	IL1A	340G>T	Ala114Ser	56	0.25	0.32	0.30	0.547	0.673	0.227	0.488	0.352	0.502	0.868	0.762	0.549
25	188378669	CXCL8	91G>T	Glu31Ter	147	0.76	0.03	0.00	0.341	0.143	0.417	0.166	0.139	0.089	0.478	0.777	0.459
26	2233418	NFKB1	548-3C>T	Splice variant?	66	0.80	0.02	0.01	0.520	0.570	0.253	0.797	0.956	0.524	0.567	0.518	0.460
27	2070766	MIF	282-6C>G	Splice variant?	85	0.74	0.11	0.18	0.755	0.957	0.457	0.892	0.576	0.914	0.772	0.366	0.529
28	1800629	TNF	-308G>A	Regulatoryfeature	154	0.00	0.15	0.13	0.270	0.950	0.144	0.213	0.805	0.117	0.537	0.546	0.784
29	361525		-238G>A	Regulatory feature	138	0.47	0.06	0.06	0.297	0.315	0.124	0.095	0.052	0.087	0.376	0.376	0.478
30	61763106	MAPK14	306-7delT	Splice variant?	196	0.05	0.09	0.10	0.590	0.317	0.383	0.066	0.881	0.901	0.155	0.041	0.147
31	776746	CYP3A5	219-213G>A	Splice variant	364	0.80	0.02	0.07	0.520	0.570	0.253	0.750	0.603	0.471	0.605	0.800	0.360

S – GC sensitivity, D – GC dependence, R – GC resistance, HWE – Hardy-Weinberg equilibrium, MAF – minor allele frequency, EU – European population, #mean.

Statistical significance was set at p ≤ 0.05. Bonferroni correction was applied for the number of study groups, but not for the number of analyzed loci. All calculations were performed using Statistica 12.0 software (Stat Soft, 2016).

Results

Clinical outcome

Our study group included 139 IBD patients, of whom 77 (56%) reacted correctly to GC therapy and were classified as GC-sensitive, 24 (17%) demonstrated resistance to treatment and 38 (27%) were GC-dependent (Table I). The most frequently used GC drug was methylprednisolone in a dose above 24 mg/24 h (Table II). Correlation analysis of the treatment and drug dose with the response and gender did not reveal statistically significant differences (p = 0.072, p = 0.469, respectively). There was a relationship between the treatment method and the diagnosis of CD and UC (p < 0.001). However, if budesonide was omitted, the significance disappeared (p = 0.861, Table II).

LR-PCR and NGS results

The target region of NR3C1, NLRP1, IPO13, FKBP5, HSPA4, ABCB1, STIP1, HSP90AA1, IL-1A, IL-1B, IL-2, IL-4, CXCL8, IL-10, NFKBIA, JUN, MIF, TNF, MAPK14, CYP3A4, and CYP3A5 genes was amplified for each DNA sample in 54 fragments. The results of LR-PCR for one example DNA sample are presented in Figure 3. The NGS results obtained indicated on average 430 ±63 SNVs per sample. Among 121 different functional DNA variants in the total group, for correlation analysis we selected 31 polymorphic variants in 14 genes (Table V), which occurred at a minor allele frequency (MAF) above 2%. These have consequences for the protein structure and functionality based on VariantStudio 3.0 Software classification (Figures 1 and 4). Variants in the remaining studied 7 genes that were studied (IL-1B, IL-2, IL-10, IL-4, JUN, CYP3A4, IPO13) did not meet these selection criteria and were not qualified for further analysis. Moreover, in the literature described in response to the GCS context, the loci in the NR3C1 (c.1184+646C>G and p.Asn767Asn), ABCB1 (c.2685+49T>C, p.Gly412Gly and p.Ile1145Ile) and TNF (c.-308G>A and c.-238G>A) genes were also included [11]. The highest average number of NGS reads (above 400×) was observed for variants p.Cys13Ser (rs61745470) and p.Met71Leu (rs8005905) of the HSPA4 and HSP90AA genes, respectively. The lowest coverage occurred for the ABCB1 gene changes (22–65×) (Table V). To verify the genotyping results obtained for this gene, we checked all five loci for the whole group of patients using the second method – Sanger sequencing. These data were overlapping. The SNVs were mostly in HWE. The exception concerned changes p.Val1059Met (rs2301582) and p.Met1184Val (rs11651270) in the NLRP1 gene, and p.Ser400Asn (rs2229109) in the ABCB1 gene (p < 0.05, Table V). The allele frequencies identified and presented as MAF generally did not significantly differ from the European data (EU MAF) originating from the 1000 Genomes Project, except for the change p.Glu31Ter (rs188378669) in the CXCL8 gene, in which in contrast to the European population, the minor allele occurred eight times more often in our patient group (Table V).

Figure 4

Flow diagram of genetic variant selection steps

SNVs – single nucleotide variants, bp – base pair, SD – standard deviation, MAF – minor allele frequency.

Correlation data

In the first step, correlation analysis of obtained genotypes with a response to GC was performed in the whole group of 139 IBD patients (jointly all CD and UC cases). Only for the change c.2685+49T>C (rs2032583) in the ABCB1 gene was an association with GC reaction visible (p = 0.046). The presence of a substitution in position c.2685 correlated with the occurrence of resistance to GCs. Subsequently, the correlation analysis for two patient subgroups CD and UC was carried out, which demonstrated a statistically significant relationship in UC between GC response and the occurrence of the variants c.1088G>A (rs56149945) in NR3C1 (p = 0.019), c.241+6A>G (rs2817033) in the FKBP5 gene and c.306-7delT (rs61763106) in the MAPK14 gene (Table V, Figure 5), while in CD, the change c.2685+49T>C (rs2032583) in the ABCB1 gene was relevant.

Figure 5

Distribution of genotypes for significantly correlated variants: A – NR3C1, c.1088A>G, UC patients, B – FKBP5, c.241+6A>G, UC patients, C – MAPK14, c.306-7delT, UC patients, D – ABCB1, c.2685+49T>C, CD patients

ns – non-significant, *relevant genotypes for the GC-response type, p < 0.05, **p < 0.01, ***p < 0.001.

The allele G in position c.1088 (rs56149945) of the NR3C1 gene was identified in UC in 1 patient with GC sensitivity (allele frequency: 1.35%) compared to 3 patients heterozygous and 1 homozygous (allele frequency: 19.2%) in the GC-resistant group. In patients with GC dependency allele G was not found (Figure 5). Consequently, subjects with genotype AG and homozygous GG showed significantly increased resistance to GCs compared to the rest of the patients with UC (OR = 21.333, 95% CI: 2.130–213.648, χ2 = 11.44, p < 0.001). The splicing variant c.241+6A>G (rs2817033) in the FKBP5 gene was relevant to the GC response when comparing the GC-sensitive group versus the rest of UC patients (p = 0.040, Table V). The frequency of the G allele was 62.2% in the GC-sensitive group, 58.3% in -dependent and 42.3% in -resistant. However, the substitution c.241+6A>G inversely correlates with the response to GCs. The homozygous AA is significantly more frequent than AG + GG genotypes in the GC-sensitive group compared to the rest of the UC patients (OR = 0.188, 95% CI: 0.044–0.799, χ2 = 5.84, p = 0.016, Figure 5).

We also observed that the variant c.306-7delT in the MAPK14 gene (p = 0.040) has an effect on GC sensitivity. The allele with deletion T in position c.306-7 was found in UC with a frequency of 4.05% in the GC-sensitive group, 16.7% in the GC-dependent and 19.2% in the GC-resistant group. Deletion T carriers and homozygotes with a mutation showed significant increased dependency or resistance to therapy compared to GC-sensitive patients with UC (OR = 9.722, 95% CI: 1.922–49.183, χ2 = 9.73, p = 0.002).

In CD patients, only one substitution, c.2685+49T>C in the ABCB1 gene, has an impact on resistance to GCs (p = 0.034), which was also visible in the whole group of 139 IBD patients (p = 0.046). The changed allele C was indicated in 11 (17.3%) GC-sensitive and 9 (13.6%) GC-dependent patients. In patients with GC resistance, this allele was absent. The results show that allele C causes a protective effect on GC resistance (Figure 5). The genotype TC demonstrated significantly lower frequency versus TT in GC-resistant patients compared to rest of the CD patients (OR = 0.099, 95% CI: 0.006–1.755, χ2 = 4.50, p = 0.034).

Discussion

Pharmacogenetics in IBD is the subject of a significant amount of research, which consistently concludes that there is a lack of sufficiently strong biomarkers for GC therapy that are useful in clinical practice, thus emphasizing the need for further research [7, 14, 17]. Although the mechanism of GC action involves dozens of genes, so far intensive pharmacogenetic studies of these drugs have focused mainly on single genes, and have produced contradictory results. An example is the NR3C1 gene, of which the c.1184+646C>G polymorphism shows the effect with the GC response [14, 22], although there is a study which does not confirm this observation [16, 23].

Considering the action, transport and metabolism of GCs, in our research we decided to expand significantly the number of genes studied to include both frequently tested NR3C1 for GR and transporter ABCB1 and cytokine genes IL-1A, IL-1B, IL-2, IL-4, CXCL8, IL-10, TNF, chaperones and co-chaperones HSP90AA1, HSPA4, STIP1, FKBP5, inflammasome activator NLRP1, kinase MAPK14, transcription factor genes JUN and NFKBIA, the macrophage migration inhibitory factor MIF gene, metabolism enzyme genes CYP3A4 and CYP3A5, and the importin gene IPO13. The most important purpose of the present study was to determine whether the consequence of genetic variants of the 21 genes selected affect the GC response by studying a cohort of 139 Polish adult IBD patients treated with methylprednisolone, hydrocortisone and slight budesonide (Table II). It was also examined whether the GC treatment depends on gender and diagnosis of CD or UC and whether it affects the GC response in IBD patients because different treatment strategies, i.e. oral methylprednisolone, intravenous hydrocortisone, and oral budesonide, could create a severe limitation (Table II). It is also worth emphasizing that the bioavailability and hepatic activity of budesonide differ from methylprednisolone and hydrocortisone. At the same time, budesonide is dedicated primarily to a specific form of CD. According to ECCO guidelines, it should be effective in patients with ileocaecal localization of CD [1]. We have found no connection between the method of GC therapy and response to it. The GC drugs used in our studies were similar to those investigated in earlier pharmacogenetic research on patients with IBD conducted by Krupoves et al., who evaluated prednisolone and budesonide [24, 25]. Gender also showed no effect on the response to the GCs, but we did identify a significant difference in treatment between patients with CD and UC. However, oral treatment with budesonide during this study was dedicated in Poland only to patients with CD. Therefore, the UC patients could not receive it, which led to differences in GC treatment in our CD and UC patient groups.

High-throughput NGS technology was very helpful in achieving our molecular genetics aim. Therefore, as an additional goal of this work, we developed a rapid, flexible and economical application of NGS by combining it with LR-PCR as target enrichment for pharmacogenetic studies on GC drugs. A similar solution was used successfully in earlier studies [19, 20]. In bioinformatic analysis for our scientific purposes, in contrast to diagnostic standards of minimum coverage 30× [26], we carefully applied a reads cover filter of 10 to avoid omission of any variants at this stage (Figure 2). For most of the loci, the mean coverage for each position analyzed was higher than the cut-off value, i.e. between 22 and 435 (Table V). When it was below 30×, we verified the genotyping data. The weak point of this method turned out to be only the ABCB1 gene. Therefore, for all changes selected for this gene, we performed a full verification for all patients. In the future, the input of the ABCB1 gene amplicons can be increased in the whole library to improve the results. It must be pointed out, however, that the LR-PCR NGS method we applied is easy to modify and relatively cheap. Our calculations including reagents and materials show that LR-PCR NGS would be at least seven times less expensive and laborious than Sanger sequencing [20].

The results of sequence analysis with a total length of 376 894 kb per patient showed the presence of an average of 430 changes per sample, 31 of which were selected for statistical analysis, and finally 4 significantly associated with the response to GCs (Figure 4). These variants include p.Asn363Ser in NR3C1, c.241+6A>G in the FKBP5, c.306-7delT in the MAPK14 and c.2685+49T>C in the ABCB1 gene. Most of the changes investigated in the context of the GC response in IBD refer to the NR3C1 and ABCB1 genes (Table V). The NR3C1 gene variants have been related to both GC sensitivity and resistance [14, 24, 27, 28], while ABCB1 variants, leading to decreased activity of transporter protein, were associated with GC dependence [25, 29]. Previous studies on pediatric IBD patients showed the significance of substitution c.1184+646C>G, leading to increased sensitivity in GC responders versus dependents [14, 22]. Krupoves et al. demonstrated exonic changes p.Arg23Lys and p.Asn767Asn of the NR3C1 gene in pediatric CD in association with GC resistance [24]. Our investigation does not confirm this outcome and is consistent with studies conducted on a group of 200 Caucasians with IBD, where the researchers did not find any association [16]. In our observation, only variant p.Asn363Ser in the NR3C1 gene, which results in the amino acid change from asparagine to serine in codon 363 in the N-terminal trans-activation domain of the GR protein, revealed a significant difference in genotype distribution in the GC-sensitive, -dependent and -resistant groups (p = 0.019). The genotype AG and homozygous GG correlate with increased resistance to GCs compared to the rest of the UC patients, but not CD subjects (Figure 5). In contrast, independent studies have shown no association between the p.Asn363Ser change and the response to GC therapy [14, 16, 23, 24]. There is only evidence that this substitution influences increased sensitivity for exogenous GCs. Huizenga et al. demonstrated stronger inhibition of lymphocyte proliferation by exogenous dexamethasone in heterozygous AG than in wild-type patients. Carriers of the G allele were characterized by a higher BMI index and higher insulin level in blood after the administration of dexamethasone [30].

In previous research on ABCB1 gene variants, Annese et al. found no association between p.Ser893Ala (rs2032582) and p.Ile1145Ile (rs1045642) polymorphisms and for the response to GC treatment. They investigated a group of 594 patients with IBD treated with GCs and a control group of 450 healthy people [31]. The same observation was made by Cucchiara et al. in their research on groups of 200 children with CD and 186 children with UC [15]. However, in contrast to those data, in a Chinese population among three studied variants c.2677T>G, c.1236C>T (p.Gly412Gly, rs1128503), and c.3435C>T of the ABCB1 gene, the CC genotypes in positions c.1236 and c.3435 were more frequent in GC-dependent CD patients than in those responsive to GC, showing that the CC homozygosity of the two SNPs might be the risk factor for GC dependence in CD [29]. Furthermore, Canadian researchers demonstrated a significant association between the polymorphism 2685+49T>C (rs2032583) and the haplotype C-T-C for changes c.3435C>T, c.117+4196T>C (rs3789243) and 2685+49T>C with the occurrence of GC dependence in 260 juvenile patients with CD [25]. They found that in position c.2685+49 in the ABCB1 gene the allele C and genotype CT are protective for dependency on GC treatment, (p = 0.029 and p = 0.035, respectively) [25]. This observation is partly consistent with our results. After correlation analysis of five SNVs – c.3435C>T, c.2677T>G, c.1236C>T, c.2685+49T>C and c.1199G>A (p.Ser400Asn, rs2229109) – in the ABCB1 gene, we found an inverse association of c.2685+49T>C with GC resistance (p = 0.034). We confirm the protective function of the C allele by GC-sensitive and -dependent patients against GC-resistant patients with CD. For the other four ABCB1 variants, we did not observe any significant changes in the frequency of genotypes in particular GC groups.

In these studies, we also identified two associations in loci which were not described earlier in the context of GC treatment and they are not reported in the ClinVar database (https://www.ncbi.nlm.nih.gov/clinvar/). We found that the splicing variants c.241+6A>G in the FKBP5 and c.306-7delT in the MAPK14 gene were relevant in GCs’ response comparing the GC-sensitive group versus the rest of the UC patients (p = 0.040 and p = 0.041, respectively). Until now, only studies by Italian scientists exist in this area, which demonstrated the influence of a putative promoter variant c.-20+18122T>C (rs4713916) in the FKBP5 gene on GC sensitivity in patients with CD but not in UC [16]. The FKBP-51 protein is critical in the correct binding of GCs to GR and in transporting GR to the cell nucleus. Maltese et al. explain that increased expression of FKBP-51 leads to the decreased ability of GR-α to bind GCs, but also to increased translocation of the nonactive GR-β isoform to the cell nucleus [16, 32]. We found that the wildtype homozygous AA in position c.241+6 is significantly rarer than AG + GG genotypes in GC-sensitive patients compared to the rest of the UC patients. Surprisingly, this would indicate that change promotes sensitivity to treatment, which may be related to the maintenance of optimal proportions in the binding of FKBP-51 protein to the active and non-active isoform of the GR. Therefore, we postulate that a splicing variant of the FKBP5 gene observed in our studies could change the level of functional FKBP-51 protein or modify it, and contribute to the emergence of higher GC sensitivity. However, this issue requires further functional research.

The last significant variant in our studies, c.306-7delT, is located in the intron 3 acceptor region of the MAPK14 gene coding for mitogen-activated protein kinase 14, which is also known as MAP kinase p38 α. We found a correlation between carriers of T deletion and homozygous mutants and increased GC dependency or resistance compared to GC sensitivity in patients with UC (OR = 9.722, p = 0.002). MAP kinase p38, which is activated by IL2 and IL-4, phosphorylates GR, resulting in the inhibition of its transcriptional activity [33]. GCs reduce the MAPK signaling pathway via the induction of MAPK phosphatase 1 (MKP1), which consequently inhibits the expression of different inflammatory genes [34]. Furthermore, Bhavsar et al. have proved in severe asthma patients that GC resistance may be improved by inhibitors of p38 MAPK [35]. Accordingly, we suppose that the intronic variant c.306-7delT of the MAPK14 gene relevant in our study could lead to the creation of an altered protein variant, which may not be efficiently inhibited by MKP1 and consequently GC resistance and dependence occurs [12].

In the literature, DNA variants are described, i.e., p.Leu155His (rs12150220) in the NLRP1 gene [14], and c.-308G>A (rs1800629) in the TNF gene [15, 36], related to the response to GCs. We have not confirmed these observations in our studies. Moreover, our results showed apparent differences in the associations of variants between CD and UC that were analyzed. These immunological disorders, although belonging to one group, IBD, are two different diseases with as yet unclear etiology, a different genetic background and probably different pharmacogenetics [16, 37].

It is worth emphasizing that the NGS results obtained seem to have further potential due to the remaining untranslated regions (UTRs), promoters and intron variants which, according to our assumptions, we essentially filtered out. Thus, in the next step, we are planning to correlate these SNVs with clinical data concerning the GC effect. The especially high importance of non-exonic variants was proved by de Iudicibus et al. in recent research, because of binding sites for microRNA [38]. Moreover, it should be taken into account that for personalized glucocorticoid treatment in IBD, a comprehensive approach integrating both genome and transcriptome research, including coding and non-coding RNA, will be necessary [39].

We are aware of the limited subjects of our study group. Therefore, explaining the functional effect of the variants being analyzed will be of particular importance and the present results might become part of the meta-analyses, which will bring us closer to selecting biomarkers for predicting and individualizing GC treatment.

In conclusion, our results demonstrate the critical role of genetic factors of various genes such as FKBP5 and MAPK14 in the response of IBD patients to GC treatment, which indicates the complex contribution of low penetration genes in the GC response. Moreover, we observed that the genetic background differs between CD and UC. Our results should be confirmed in studies on different populations and further functional investigation.

Conflict of interest

Marcin Gabryel is a recipient of a fellowship for a young researcher from Poznan Medical University (Poland, grant no. 502-14-02223359-41103) and from the Polish National Science Centre (grant no. 2015/17/N/N25/00689). Other authors declare no conflict of interest.