Whole Transcription Profile of Responders to Anti-TNF Drugs in Pediatric Inflammatory Bowel Disease

Background: Up to 30% of patients with pediatric inflammatory bowel disease (IBD) do not respond to anti-Tumor Necrosis Factor (anti-TNF) therapy. The aim of this study was to identify pharmacogenomic markers that predict early response to anti-TNF drugs in pediatric patients with IBD. Methods: An observational, longitudinal, prospective cohort study was conducted. The study population comprised 38 patients with IBD aged < 18 years who started treatment with infliximab or adalimumab (29 responders and nine non-responders). Whole gene expression profiles from total RNA isolated from whole blood samples of six responders and six non-responders taken before administration of the biologic and after two weeks of therapy were analyzed using next-generation RNA sequencing. The expression of six selected genes was measured for purposes of validation in all of the 38 patients recruited using qPCR. Results: Genes were differentially expressed in non-responders and responders (32 before initiation of treatment and 44 after two weeks, Log2FC (Fold change) >0.6 or <−0.6 and p value < 0.05). After validation, FCGR1A, FCGR1B, and GBP1 were overexpressed in non-responders two weeks after initiation of anti-TNF treatment (Log2FC 1.05, 1.21, and 1.08, respectively, p value < 0.05). Conclusion: Expression of the FCGR1A, FCGR1B, and GBP1 genes is a pharmacogenomic biomarker of early response to anti-TNF agents in pediatric IBD.


Introduction
Inflammatory bowel disease (IBD), which includes ulcerative colitis (UC) and Crohn disease (CD), is a multifactorial autoimmune disorder in which a quarter of patients are diagnosed when aged under 18 years [1,2]. IBD, when diagnosed in children, is linked with more extensive disease and greater complications compared to patients whose disease first appears in adulthood [3]. Since children with pediatric IBD (pIBD) have to take current therapy for longer, treatment must be optimized.
The use of biological therapy, such as anti-Tumor Necrosis Factor (anti-TNF) agents, has dramatically changed the treatment of autoimmune disease, including IBD. The use of these drugs is often linked to more severe symptoms of pIBD [4]. The only anti-TNFs approved for pIBD are infliximab (IFX) and adalimumab (ADL). However, treatment with biological drugs very often fails. Thus, up to 41% of children with moderate to severe CD and treated with IFX do not achieve clinical remission [5].
Mucosal healing is the best outcome in pIBD. However, given that pIBD is a chronic disease whose response cannot be monitored using regular biopsies, the use of non-invasive biomarkers is highly recommended.
Trough serum anti-TNF levels and antidrug antibodies, among other serological biomarkers, are usually measured in IBD to monitor anti-TNF treatment response [6,7]. However, neither can be measured prior to starting or during the first two weeks of anti-TNF treatment. Trough serum anti-TNF levels as soon as six weeks after initiation of treatment were recently reported to predict remission [8]. No earlier biomarkers have been identified to date.
Identification of genomic biomarkers could be useful to identify groups of pIBD patients who are less likely to respond in early stages of treatment or even before initiation. The mRNA levels in some genes have been identified as pharmacogenomic biomarkers of the activity of anti-TNF drugs in the inflamed tissues of adults diagnosed with IBD or other autoimmune disorders [9][10][11][12][13]. However, these biomarkers are identified using invasive techniques that are not suitable for monitoring. In addition, pharmacogenomic biomarkers of response to anti-TNF drugs have not been extensively investigated in pIBD. Identification of biomarkers in blood facilitates monitoring. In a recent comparison with healthy people, several genes were differentially expressed in the blood of children, but not adults, diagnosed with IBD during an active phase of the illness [14]. On the other hand, some studies in adults have revealed the usefulness of biomarkers of gene expression from whole blood in the assessment of response to anti-TNF agents [15,16].
In the present study, we analyzed whole gene expression profiles using next-generation sequencing of RNA in whole blood from children diagnosed with IBD. Differential gene expression before and after two weeks of treatment with IFX or ADL was analyzed with the aim of identifying very early biomarkers of response to IFX or ADL in pIBD.

Patient Samples
This study was prospective and multicentric and recruited 38 IBD patients aged < 18 years between March 2017 and May 2019, (30 with CD and eight with UC, 17 treated with ADL and 21 with IFX) [17]. The groups analyzed were matched for age and sex. Patients with a confirmed diagnosis of inflammatory bowel disease, aged between 1-17 years, and who had started treatment with infliximab (5 mg/kg, 0-2-6 weeks) or adalimumab (160/80 mg in those patients weighing more than 40 kg and 80/40 mg for those weighing 40 kg or less) were included.
Age, sex, type of IBD, anti-TNF drug, and specific disease activity scores, such as Pediatric Crohn Disease Activity Index (PCDAI) and Pediatric Ulcerative Colitis Activity Index (PUCAI), were collected to measure anti-TNF response, which was defined as a decrease of at least 15 points in PCDAI or PUCAI from the start of treatment to weeks 14 (IFX) or 26 (ADL). Study data were collected and managed using Research Electronic Data Capture (REDCap) tools hosted at Hospital General Universitario Gregorio Marañón, Madrid, Spain [18].

Ethics Statement
This study was approved by the Ethics Committee of Hospital General Universitario Gregorio Marañón with the number LAL-TNF-2019-01. Written, informed consent was obtained from the patients and parents or legal guardians.

Extraction of Total RNA from Whole Blood
Blood samples were collected in Paxgene tubes (PreAnalytics, Hombrechtikon, Switzerland) and whole blood RNA was extracted using PAXgene Blood RNA kit (PreAnalytics) at two different points: before the first administration of ADL or IFX (week 0) and after two weeks of the first administration of the drug (week 2) following manufacturer's recommendations. The total RNA concentration was measured by spectrophotometry, and the integrity of RNA was verified by electrophoresis. Only RNA samples with an RNA integrity number > 7 were used.

RNA Sequencing
The quality and integrity of each RNA sample were checked using both a Bioanalyzer and a Nanodrop device before proceeding to the RNA sequencing (RNAseq) protocol. Poly A+ RNA from 100 nanograms of total RNA was reverse transcribed and barcoded RNAseq libraries were constructed using NEBNext Ultra II Directional RNA Library Prep Kit (New England Biolabs) following manufacturer's recommendations. The quality of each RNA sample library was checked using a Bioanalyzer and a Qubit.
Sequencing reads were aligned to the human reference transcriptome (GRCh38 v91) and quantified with RSem v1.3.1 (Li and Dewey 2011). Raw counts were normalized using transcripts per million and the trimmed mean of M values, transformed into log2 expression (log2[rawCount+1]), were compared to calculate fold-change and corrected p value. Only those genes expressed with at least one count in at least 12 samples were taken into account. As there are no gene expression changes with an associated Benjamini and Hochberg-adjusted p value < 0.05, we considered candidates to be confirmed as such by qPCR genes with |log2FC| > 0.6 and a non-adjusted p value < 0.05. The RNAseq data have been deposited with the accession number GSE159034 in the Gene Expression Omnibus database (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE159034) [19].

Quantitative Reverse Transcription-Polymerase Chain Reaction (qRT-PCR)
Total RNA was reverse transcribed and amplified, and relative expression of GBP1, GBP5, IGHG2, GNLY, FCGR1A, FCGR1B, ACTB, and RPL4 was quantified, as described in Salvador-Martín [17]. ACTB and RPL4 were used for normalization and three technical replicates were used for each sample. The oligonucleotide sequences used for gene amplification are shown in Table 1. Primer pair efficiency was used for correction and relative expression calculated using the 2 −∆∆ Ct method.

Statistical Analysis
Individual gene expression analyses at t = 0 and t = 2 weeks were performed using ExpressionSuite v1.1 (Applied Biosystems, Foster City, CA, USA), using the responder sample (D005) as relative quantification 1 in both times of comparison. Comparison of gene expression changes from t = 0 to t = 2 in responder versus non-responders was performed using GraphPad Prism (GraphPad Software, San Diego, CA, USA), using the responder sample (D005) at t = 2 as relative quantification 1. The mean relative quantification on the triplicated samples was used for expression and the unpaired t test applied for analyzing responder versus non-responder groups. P values were corrected using the false discovery rate with a confidence level of 95%. Categorical and numerical variables were compared using the t test and the Fisher exact test, respectively. For all tests, a p value < 0.05 was considered as statistically significant.
The statistical review was performed by a biomedical statistician. The positive predictive value (PPV), negative predictive value (NPV), sensitivity, specificity, and diagnostic odds ratios for relative expression of GBP1, FCGR1A, and FCGR1B were calculated as described elsewhere [20]. The + and − likelihood ratios were calculated with a 95% confidence interval (CI).

Patients' Characteristics
Thirty-eight patients (29 responders and nine non-responders) met the inclusion criteria and were included in the study. The failure rate was 23.7%. The characteristics of both groups of patients are summarized in Table 2. Patients were mainly male (52.6%; median age at diagnosis, 10.5 years), diagnosed with CD (78.9%), and treated with IFX (55.3%). The statistical differences between both groups were in the PCDAI and C-reactive protein (CRP) levels at initiation of treatment The demographic and clinical variables of the six responders and six non-responders selected for RNAseq were more homogeneous than those of the total population (Supplementary Materials Table S1). In these patients only PCDAI was statistically significant between both groups (p = 0.025).

Differential Gene Expression Using RNAseq in the Response of Anti-TNF Agents Prior to Starting Treatment
Twenty genes were overexpressed and two downregulated in non-responders versus responders (Log2FC [Fold change] > 0.6 or <-0.6 and p value < 0.05) immediately prior to the first administration of the anti-TNF agent (Table 3). For this analysis, the relative expression of the whole transcriptome using RNAseq was measured in six responders and six non-responders. Responders were used as the reference group. The most overexpressed gene in non-responders was GNLY (2.8 fold). The most downregulated gene in the same patients was DNAJC13 (2.4 fold).

Differential Gene Expression in Response to Anti-TNF Agents at Week 2 Post-Treatment
Twenty-six genes were overexpressed in non-responders and 16 were downregulated in responders (Log2FC [Fold change] >0.6 or <−0.6 and p value < 0.05) at two weeks post-treatment with anti-TNFs (Table 4). Responders were used as the reference group. The most overexpressed gene in non-responders was GBP5 (2.4 fold). The most downregulated gene in the same patients was IGLV1-44 (2.4 fold).

Functional in Silico Analysis
Functional analysis of differentially expressed genes was performed using Ingenuity Pathways Analysis (IPA, Qiagen, Germany).
Developmental disorder (p value range 1.96 × 10 −2 to 9.23 × 10 −15 ), cell-to-cell signaling interaction (p value range 1.96 × 10 −2 to 4.30 × 10 −9 ), and hematological system development and function (p value range 1.96 × 10 −2 to 4.30 × 10 −9 ) were found to be the most significant disease and biofunctions represented by all of the selected genes prior to initiation of treatment (Supplementary Materials Table S2). Figure 1 shows the main network generated by IPA and based on known interactions between the genes expressed differentially between responders and non-responders prior to initiation of anti-TNF treatment. The main diseases and functions associated with this network were developmental disorder, hereditary disorder, and metabolic disease.  The results for all of the genes selected after two weeks of anti-TNF treatment showed inflammatory response (p value range 1.30 × 10 −2 to 5.74 × 10 −23 ), cellular function and maintenance (p value range 1.19 × 10 −2 to 5.11× 10 −19 ), and humoral immune response (p value range 1.23× 10 −2 to 5.74× 10 −23 ) to be the most significant disease and biofunctions represented by the selected genes after two weeks of treatment (Supplementary Materials  Table S3) The most informative network generated by IPA and based on known interactions between the genes differentially expressed between responders and non-responders after two weeks of anti-TNF treatment is represented in Figure 2. The top disease and functions associated with this network were consistent with those obtained for all of the selected genes, namely, cellular function and maintenance, humoral immune response, and inflammatory response. For this reason, we decided to select most of the genes for validation from among the genes that were differentially expressed after two weeks of treatment. The results for all of the genes selected after two weeks of anti-TNF treatment showed inflammatory response (p value range 1.30 × 10 −2 to 5.74 × 10 −23 ), cellular function and maintenance (p value range 1.19 × 10 −2 to 5.11 × 10 −19 ), and humoral immune response (p value range 1.23 × 10 −2 to 5.74 × 10 −23 ) to be the most significant disease and biofunctions represented by the selected genes after two weeks of treatment (Supplementary  Materials Table S3) The most informative network generated by IPA and based on known interactions between the genes differentially expressed between responders and non-responders after two weeks of anti-TNF treatment is represented in Figure 2. The top disease and functions associated with this network were consistent with those obtained for all of the selected genes, namely, cellular function and maintenance, humoral immune response, and inflammatory response. For this reason, we decided to select most of the genes for validation from among the genes that were differentially expressed after two weeks of treatment.

Validation of Differentially Expressed Genes by qRT-PCR
Eight genes were selected from the RNAseq analyses for validation by real-time PC (one differentially expressed at T0 and seven at T2). Semiquantitative real-time PCR w performed to assess RNA from patients' blood prior to treatment and after two weeks treatment (Figure 3). None of the genes studied was expressed differentially between r sponders and non-responders before initiation of anti-TNF treatment.

Validation of Differentially Expressed Genes by qRT-PCR
Eight genes were selected from the RNAseq analyses for validation by real-time PCR (one differentially expressed at T0 and seven at T2). Semiquantitative real-time PCR was performed to assess RNA from patients' blood prior to treatment and after two weeks of treatment (Figure 3). None of the genes studied was expressed differentially between responders and non-responders before initiation of anti-TNF treatment.  GBP1 was overexpressed in non-responders compared with responders after two weeks of treatment (LogFC = 1.08, p value 0.006, False Discovery Rate (FDR) 0.032). In addition, FCGR1A and FCGR1B were also induced 1.05 fold more in non-responders (p value = 0.006, FDR 0.035) and 1.21 fold more in responders (p value = 0.005; FDR 0.032) (in LogFC). No other statistically significant changes were detected.
A comparison of the results for changes in gene expression using RNAseq and qRT-PCR (Table 5) revealed a good correlation, mainly in the data on genes selected after two weeks of anti-TNF treatment (R 2 = 0.83). The differences were statistically significant using both techniques in three cases.

Prediction of Response to Anti-TNF Therapy Based on Expression of GBP1, FCGR1A, and FCGR1B after Two Weeks of Treatment
Expression of GBP1, FCGR1A, and FCGR1B mRNA at T2 was higher in non-responders than in responders (Figure 4). The PPV, NPV, sensitivity, specificity, diagnostic odds ratio, positive likelihood ratio (+LR), and negative likelihood ratio (-LR) are presented in Table 6. The best diagnostic odds ratio corresponds to FCGR1B expression

Differences in Gene Expression between Responders and Non-Responders during the First Two Weeks of Anti-TNF Therapy
Changes in gene expression from T0 to T2 were measured for the eight selected genes when responders and non-responders were compared (Figure 4). Only FCGR1A changed its expression during the first two weeks of treatment (p value < 0.05). An increase in FCGR1A expression was observed in non-responders, while a decrease was observed in responders. Similar trends were observed for GBP1 and FCGR1B, although the differences were not statistically significant.

Discussion
The use of biologic drugs such as anti-TNF agents has dramatically transformed the treatment of autoimmune diseases, including IBD. However, more than 20% of patients do not respond correctly to this therapy [6]. The identification of the patients in whom therapy is more likely to fail in early stages or even before initiation would enable therapy to be personalized. The benefits of personalization in terms of safety and efficacy are of particular interest in children, who are necessarily treated for longer. Finding specific biomarkers for children is necessary since genetics has a more important role in pIBD than in adult disease. Consequently, several gene polymorphisms have been involved in susceptibility to pIBD or to very-early-onset ulcerative colitis [21,22]. The serum levels of proteins, such as clusterin and ceruloplasmin, also differ between children and adults with IBD [23]. As for response, genetic polymorphisms in genes such as ATG16L1, CDKAL1, ICOSLG, BRWD1, and HLA-DQA1 have been associated with the response to anti-TNF treatment in children [24]. In contrast, several genes associated with expression and response to anti-TNF response have been identified only in adults with IBD [9,10,13]. Our group recently showed expression of SMAD7 in blood to be a biomarker of the response to anti-TNF agents two weeks after initiation of treatment [17]. To our knowledge, ours is the first study to assess whole gene expression profile by RNAseq in pIBD using biomarkers of response to anti-TNF agents. We identified putative gene networks and pathways involved in early response in pIBD and validated GBP1, FCGR1A, and FCGR1B as potential pharmacogenomic markers. Although we are still at a very early stage, the identification of these non-invasive biomarkers in pediatric IBD could revolutionize the selection of biological drug treatment in these patients in the future.
Humoral immune response, inflammatory response, and maintenance of cellular function were associated with the differentially expressed genes in responders and nonresponders after two weeks of anti-TNF treatment. These diseases and biofunctions are clearly associated with IBD and indicate that our strategy could be helpful for identifying biomarkers of response to anti-TNF agents [25,26]. However, since these diseases and biofunctions are too extensive and provide little information in terms of prediction, we decided to focus on the differentially expressed genes verified by qPCR.
Certainly, only three out of eight genes were statistically validated. However, the correlation between log2 ratios (RNAseq) and relative quantification (RT-qPCR) was as high as R 2 = 0.83. This value was lower than those found in other works with ideal conditions for comparison [27]. Nevertheless, it was similar or even higher than the values found in other works with more variable samples, such as tissues from living organisms [28,29].
We found that expression of GBP1 mRNA was upregulated in the peripheral blood cells of non-responders after two weeks of anti-TNF treatment. GBP1 is an interferonstimulated, guanylate-binding protein involved in defense against pathogens and inflammation; its levels are elevated in the mucosa of patients with active disease [30]. The fact that this gene was overexpressed in non-responders after two weeks of treatment suggests a higher likelihood of inflammation and a higher risk of treatment failure.
Similarly, GBP1 was identified as differentially expressed in colitis-susceptible mice and in colitis-resistant mice [31].
Our results also revealed greater expression of the FCGR1A and FCGR1B genes in non-responders after two weeks of anti-TNF treatment. FCRG1A, also known as CD64, is upregulated in adults and children diagnosed with clinically active IBD [32] and has been related to calprotectin level [33]. FCRG1A and FCGR1B are expressed on the surface of neutrophils and have been suggested to be potential therapeutic targets [34]. Furthermore, both ADL and IFX are less effective in the peripheral blood mononuclear cells of adult IBD patients who express elevated levels of CD64 [35]. Here, we demonstrated a similar usefulness of these genes as biomarkers of response to anti-TNF drugs in children with IBD.
There are several limitations of this study. First, the sample size is small for a study of these characteristics [12,16]. Second, our results do not distinguish between the two drugs administered to patients, infliximab and adalimumab. Thus, although the mechanism of action of both drugs is similar, we cannot rule out a differential effect. Third, it is necessary to define the role of the PCDAI index in the causality of the response to anti-TNFs. Finally, the comparison with other works is complicated by differences in the study population, as well as in the evaluation criteria of the response [9,13].
Future research will require more studies involving larger populations to confirm our findings. Single-cell RNAseq might be useful to rule out the effect of mixed-cell populations. In addition, a larger sample size could help to differentiate biomarkers by anti-TNF drug and by type of pIBD. In spite of these limitations, the identification of GBP1, FCGR1A, and FCGR1B as blood biomarkers of response to anti-TNF agents shows great potential for the personalization of therapy in pIBD.

Conclusions
We identified the expression levels of GBP1, FCGR1A1, and FCGR1B1 genes as potential biomarkers of response to treatment with IFX and ADL in children with IBD.
Supplementary Materials: The following are available online at https://www.mdpi.com/1999-492 3/13/1/77/s1, Table S1: Characteristics of patients selected for RNAseq. Table S2: Top diseases and biofunctions of genes expressed differentially between responders and non-responders before initiation of anti-TNF treatment. Table S3: Top diseases and biofunctions of genes expressed differentially between responders and non-responders two weeks after initiation of anti-TNF treatment. Data Availability Statement: The RNAseq data have been deposited with the accession number GSE159034 in the Gene Expression Omnibus database (https://www.ncbi.nlm.nih.gov/geo/query/ acc.cgi?acc=GSE159034) [19].