Blood gene expression biomarkers of response to anti-TNF drugs in pediatric inflammatory bowel diseases before initiation of treatment

Biomedicine & Pharmacotherapy 173 (2024) 116299 Available online 23 February 2024 0753-3322/© 2024 The Authors. Published by Elsevier Masson SAS. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/). Blood gene expression biomarkers of response to anti-TNF drugs in pediatric inflammatory bowel diseases before initiation of treatment Sara Salvador-Martín a , 1 , Gianluca Rubbini a , 1 , Perceval Vellosillo a , b , 1 , Paula Zapata-Cobo a , Marta Velasco c , Laura M. Palomino c , Susana Clemente d , Oscar Segarra d , Ana Moreno-´ Alvarez e , Ana Fern´ andez-Lorenzo e , Bego˜ na P´ erez-Moneo f , Montserrat Montraveta g , Cesar S´ anchez a , Mar Tolín a , In´ es Loverdos h , María Jos´ e Fobelo i , Victor Manuel Navas-L´ opez j , Lorena Magallares k , Ruth García-Romero l , Ricardo Torres-Peral m , Alejandro Rodríguez n , Ferr´ an Bossacoma o , Vicente Merino-Boh´ orquez p , Enrique Salcedo q , Rebeca ´ Alvarez r , Ana Dopazo r , María Sanjurjo-S´ aez a , Luis A. L´ opez-Fern´ andez a , * a Instituto de Investigaci´ on Sanitaria Gregorio Mara˜ n´ on, Hospital General Universitario Gregorio Mara˜ n´ on, Madrid, Spain b Unidad de Investigaci´ on Materno Infantil Fundaci´ on Familia Alonso (UDIMIFFA), Spain c Hospital Universitario Infantil Ni˜ no Jesús, Madrid, Spain d Hospital Vall d ′ Hebr´ on, Barcelona, Spain e Hospital Universitario A Coru˜ na, A Coru˜ na, Spain f Hospital Universitario Infanta Leonor, Madrid, Spain g Hospital Universitario Germans Trias i Pujol, Badalona, Spain h Hospital Universitario Parc Taulí, Sabadell, Spain i Hospital Universitario Virgen de Valme, Universidad de Sevilla, Sevilla, Spain j Hospital Regional Universitario de M´ alaga, Spain k Hospital Universitario la Paz, Madrid, Spain l Hospital Universitario Miguel Servet, Zaragoza, Spain m Complejo Asistencial Universitario de Salamanca, Salamanca, Spain n Hospital Universitario Virgen del Rocio, Sevilla, Spain o Servicio de Gastroenterología, Hepatología y Nutrici´ on Pedi´ atrica, Hospital Sant Joan de D` eu, Barcelona, Spain p Hospital Universitario Virgen de la Macarena, Universidad de Sevilla, Sevilla, Spain q Hospital Universitario Doce de Octubre, Madrid, Spain r Genomics Unit, Spanish National Center for Cardiovascular Disease (CNIC), Madrid 28029, Spain ARTICLE INFO Keywords: RNA-seq Pharmacogenomics Inflammatory bowel disease Anti-TNF drugs Pediatrics ABSTRACT Background/aims: Changes in gene expression profiles among individuals with inflammatory bowel diseases (IBDs) could potentially influence the responsiveness to anti-TNF treatment. The aim of this study was to identify genes that could serve as predictors of early response to anti-TNF therapies in pediatric IBD patients prior to the initiation of treatment. * Corresponding author. E-mail addresses: [email protected] (S. Salvador-Martín), [email protected] (G. Rubbini), [email protected] (P. Vellosillo), paula. [email protected] (P. Zapata-Cobo), [email protected] (M. Velasco), [email protected] (L.M. Palomino), susana.clemente@ vallhebron.cat (S. Clemente), [email protected] (O. Segarra), [email protected] (A. Moreno-´ Alvarez), ana.estefania.fernandez.lorenzo@ sergas.es (A. Fern´ andez-Lorenzo), [email protected] (B. P´ erez-Moneo), [email protected] (M. Montraveta), cesar.sanchez. [email protected] (C. S´ anchez), [email protected] (M. Tolín), [email protected] (I. Loverdos), [email protected] (M.J. Fobelo), [email protected] (V.M. Navas-L´ opez), [email protected] (L. Magallares), [email protected] (R. García-Romero), [email protected] (R. Torres-Peral), [email protected] (A. Rodríguez), [email protected] (F. Bossacoma), [email protected] (V. Merino-Boh´ orquez), [email protected] (E. Salcedo), [email protected] (R. ´ Alvarez), [email protected] (A. Dopazo), [email protected] (M. Sanjurjo-S´ aez), [email protected] (L.A. L´ opez-Fern´ andez). 1 These authors contributed equally to this work. Contents lists available at ScienceDirect Biomedicine & Pharmacotherapy journal homepage: www.elsevier.com/locate/biopha https://doi.org/10.1016/j.biopha.2024.116299 Received 11 November 2023; Received in revised form 6 February 2024; Accepted 17 February 2024 Biomedicine & Pharmacotherapy 173 (2024) 116299 2 Methods: We conducted a prospective, longitudinal, and multicenter study, enrolling 24 pediatric IBD patients aged less than 18 years who were initiating treatment with either infliximab or adalimumab. RNA-seq from blood samples was analyzed using the DESeq2 library by comparing responders and non-responders to anti-TNF drugs. Results: Bioinformatic analyses unveiled 102 differentially expressed genes, with 99 genes exhibiting higher expression in responders compared to non-responders prior to the initiation of anti-TNF therapy. Functional enrichment analyses highlighted defense response to Gram-negative bacteria (FDR =2.3 ×10–7) as the most significant biological processes, and hemoglobin binding (FDR =0.002), as the most significant molecular function. Gene Set Enrichment Analysis (GSEA) revealed notable enrichment in transcriptional misregulation in cancer (FDR =0.016). Notably, 13 genes (CEACAM8, CEACAM6, CILP2, COL17A1, OLFM4, INHBA, LCN2, LTF, MMP8, DEFA4, PRTN3, AZU1, and ELANE) were selected for validation, and a consistent trend of increased expression in responders prior to drug administration was observed for most of these genes, with findings for 4 of them being statistically significant (CEACAM8, LCN2, LTF2, and PRTN3). Conclusions: We identified 102 differentially expressed genes involved in the response to anti-TNF drugs in children with IBDs and validated CEACAM8, LCN2, LTF2, and PRTN3. Genes participating in defense response to Gram-negative bacterium, serine-type endopeptidase activity, and transcriptional misregulation in cancer are good candidates for anticipating the response to anti-TNF drugs in children with IBDs. 1. Introduction Inflammatory bowel diseases (IBDs) comprise a family of inflammatory, chronic, and autoimmune diseases of unknown origin that affect the gastrointestinal tract and are clinically differentiated into Crohn’s disease (CD), ulcerative colitis (UC), and IBD unclassified (IBD-U). The incidence and prevalence of IBDs increased dramatically in the second half of the last century [1]. In 25% of cases, IBDs develop during childhood or adolescence [2]. Genetic factors play an important role in the pathogenesis of IBDs, with 5%−10% of patients reporting a positive family history [3,4]. Analysis of the genes and genetic loci identified in IBDs indicates that intestinal homeostasis depends on a series of pathways, including epithelial barrier function, innate and adaptive mucosal defense, cell migration, autophagy, adaptive immunity and metabolic pathways associated with cellular homeostasis [5,6]. In children, the manifestations of IBDs differ from those seen in adults. Consequently, genetic alterations and effects involved in normal development and growth require specific investigation and a different therapeutic approach [7]. Anti–tumor necrosis factor (TNF) drugs, such as infliximab and adalimumab, function in IBDs by inhibiting TNF, thereby reducing inflammation in the gastrointestinal tract [8]. In pediatrics, anti-TNFs are also used for conditions such as juvenile idiopathic arthritis, which can co-occur with certain types of IBD. In some cases, these conditions may also be associated with uveitis, necessitating a similar biological treatment approach [9]. Younger people respond better to anti-TNF therapy [10]. However, since children have to live with this chronic disease for longer than adults, it is of paramount importance to identify strategies that can optimize pharmacological treatments with the aim of inducing remission of acute phases and maintaining an acceptable state. Therapy with anti-TNFs is recommended to induce and maintain remission in children with IBDs. In 10%−40% of cases, these therapies are characterized by a clinically irrelevant response or primary failure to infliximab or adalimumab [11,12]. These percentages depend on the age, sex, and anti-TNF drug and increase with time [13–15]. The expression of specific genes in RNA in blood can impact the way patients respond to anti-TNF drugs [16–19]. We previously used RNA-seq to analyze gene expression in a small cohort of children with IBDs receiving anti-TNF drugs and identified differentially expressed genes in children who are more likely to respond to anti-TNF therapy and in those who are less likely to benefit after 2 weeks of treatment [18]. On completion of the study, we continued to recruit patients, making it possible to perform a more complete and precise study. The aim of this work was to identify new pharmacogenomic biomarkers that can predict early response in a larger cohort of pediatric IBD patients treated with anti-TNFs, with emphasis on the period prior to initiation of anti-TNF treatment. By analyzing the expression of RNA in children with IBDs, healthcare providers might be able to identify those who are more likely to respond well to anti-TNF therapy and those who are less likely to benefit. 2. Material and methods 2.1. Study design We performed a prospective, longitudinal, and multicenter study. 2.2. Patients Participants were recruited from the pediatric gastroenterology units of 17 hospitals (Suppl. Table S1) in Spain between May 9, 2017 and January 29, 2022. The inclusion criteria were age <18 years, diagnosis of IBD (CD, UC, IBD-U), and initiation of treatment with infliximab or adalimumab. The exclusion criterion was not meeting any of the inclusion criteria. The results of RNA-seq experiments for 12 patients (6 responders and 6 non-responders) included in our cohort have been published elsewhere and are freely available (https://www.ncbi.nlm.nih.gov/geo/ query/acc.cgi?acc=GSE159034). We collected these data and added 12 new patients to the RNA-seq analysis to increase the statistical power of the comparison. These 12 new patients were selected from the cohort to ensure that responder and non-responder data were similar, after adjusting for sex, age, type of anti-TNF drug, and type of IBD. A further 2 non-responders and 16 new responders were selected for validation for reverse transcription real-time PCR (RT-PCR). 2.3. Variables and data collection Demographic and clinical variables (sex, age at diagnosis and at start of treatment, type of IBD, anti-TNF drug, first or second anti-TNF, months until treatment, concomitant immunomodulator) were collected from the clinical records. At the start of treatment, we recorded the Pediatric Crohn’s Disease Activity Index (PCDAI) for CD and the Pediatric Ulcerative Colitis Activity Index (PUCAI) for UC, as well as C- reactive protein (CRP) and fecal calprotectin (FC). The response to the anti-TNF drugs was defined as a decrease of at least 15 points in the PCDAI or PUCAI from the start of treatment to weeks 14 (infliximab) and 26 (adalimumab). Data were collected using the data capture tool (REDCap) [20]. 2.4. Ethical considerations The study was conducted in accordance with the Declaration of Helsinki and approved by the Ethics Committee of Hospital General Universitario Gregorio Mara˜ n´ on (protocol code LAL-TNF-2019–01/FG- S. Salvador-Martín et al. Biomedicine & Pharmacotherapy 173 (2024) 116299 3 2019–01, 23 September 2019). The study was validated at all the participating centers and approved by their respective ethics committees with the same code between October 2019 and December 2020. Written informed consent to perform this research was obtained from the patients and/or their legal guardians. Participants aged >12 years signed an assent form. The approved study authorized the use of samples and clinical data from a repository (C.0003459, Instituto de Salud Carlos III) from a previous project (protocol code LAL-TNF-2016–01, approved on 17 March 2016 by the same ethics committee). In the informed consent, the legal representatives agreed to the transfer of samples and data for future research in digestive diseases. 2.5. Sample collection and RNA isolation Before initiation of anti-TNF therapy, blood samples were collected and preserved in PAXgene tubes (BD Biosciences, Franklin Lakes, NJ, USA). The tubes were received at the processing laboratory and maintained at room temperature no later than 24 hours after extraction. All the tubes were then preserved at –80ºC until processing for RNA isolation using the PAXgene blood RNA kit (Qiagen, Venlo, Netherlands) following the manufacturer’s instructions. Ratio A 260 /A 280 of isolated RNAs were measured using the Quawell-5000 spectrophotometer (Quawell, Quebec, CA, USA). The quality of RNAs was measured using the Agilent RNA6000 kit in a 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA). Only samples with an RNA integrity number (RIN) >8 were used. 2.6. RNA sequencing RNA-seq processing of the 12 new samples was performed as with the previous 12 samples [18] (GSM481800, GSM481802, GSM481804, GSM481806, GSM481808, GSM481810, GSM481812, GSM481814, GSM481816, GSM481818, GSM481820, and GSM481822). Sample selection was based on similarity between responders and non-responders in terms of sex, treatment, and age. Complete transcriptome sequencing was performed at the Spanish Cardiovascular Disease Research Center based in Madrid (Centro Nacional de Investigaciones Cardiovasculares [CNIC]). The RNA-seq data of the new 12 samples have been deposited with the accession number GSE247146 in the Gene Expression Omnibus database (https://www.ncbi.nlm.nih.gob/geo/query/acc.cgi?acc=GSE247146) [21]. 2.7. Differential expression (DE) analysis from RNA-seq data Raw data per-gene counts were analyzed using the DESeq2 library [22] combined with a non-parametrical procedure to identify differential gene-expression from RNA-seq data. After estimating size factors and gene dispersions by local fit type [23], a negative binomial generalize linear model (GLM, negative binomial Wald test [22]) was used to obtain the log2 fold changes (log2FC) for each gene by comparing 2 states between samples: responders and non-responders to anti-TNF drugs. The differential gene-expression log2FC data were obtained by considering non-responder patients as a basal gene-expression reference performed using DESeq2 [22]. This protocol was generated to test the effect of the treatment. The p values provided for each gene by the Wald test were adjusted using the Benjamini-Hochberg false discovery rate (FDR) [24]. Genes were considered significant when the FDR was < 0.05. A log2FC greater than 0.6 or less than –0.6 were chosen to consider those genes as over-expressed or under-expressed respectively. 2.8. Functional characterization of significant DE genes from RNA-seq experiments based on the fold change values With the aim of interpreting profile expression from RNA-seq data at the functional level, we initially performed a gene set enrichment analysis (GSEA) [25] to provide the subset of genes that contribute to the enrichment scores based on log2FC values of a list of genes. We used Gene Ontology (GO) [26] and KEGG [27] resources for the analysis. Of note, the complete analysis––performed using GO resources––took into account GO annotations based on experimental evidence codes. Functional enrichment analysis of DE genes with FDR <0.05 was performed by applying the ClusterProfiler R package [28] using GO biological processes [26] and KEGG pathways [27]. p Values were adjusted using the Benjamini-Hochberg FDR [24] and an FDR <0.05 was considered as significant. A simplified method provided by the ClusterProfiler package was included to remove GO functional redundancy with a cutoff value of 0.5. The genes involved in the enriched KEGG pathways were mapped through the pathview R library [29], rendering the differential gene expression as log2FC within each pathway from the pathway database. 2.9. cDNA synthesis and real-time PCR validation Complementary DNAs (cDNAs) were synthetized from 1 µg of total RNA using the high-capacity cDNA archive kit (Applied Biosystems, Waltham, MA, USA) in a volume of 20 µl following the manufacturer’s instructions. For all samples, a 1/20 dilution was used as the template for the real-time PCRs. Pairs of oligonucleotides specific for each of the selected genes were selected using the KiCqStart® web from Sigma Aldrich (https://www. sigmaaldrich.com/ES/es/semi-configurators/kicqstartPrimers?active Link=selectprimers) and synthetized (Merck, NJ, USA) (Table 1). Real-time PCRs quantification was performed for CEACAM8, CEACAM6, CILP2, COL17A1, OLFM4, INHBA, LCN2, LTF, MMP8, and DEFA4 in a final volume of 10 µl (5 µl qPCRBIO SyGreen Mix 2X [PCRBiosystems, London, UK], 2 µl cDNA, and 1 µM of each primer) using the following amplification parameters: 95ºC for 3 min.; 40 cycles of 95ºC 15 sec., and 60ºC 30 sec. Relative gene expression of PRTN3, AZU1, and ELANE was performed in 10 µl (5 µl FastStart Universal SYBR Green Master [Rox] [Roche Diagnostics, Basel, Switzerland], 2 µl cDNA, 1 µM of each primer) using the following amplification parameters: 95ºC for 10 min., 40 cycles of 95ºC 15 sec., and 60ºC 60 sec. After amplification, the specificity the PCR products was analyzed using melting curves. The ACTB and GAPDH genes were used for normalization. PCR efficiency was calculated from slope. Real-time PCRs were performed in a QuantStudio 3 (Life Technologies, Carlsbad, CA, USA). Relative quantification of gene expression was analyzed using the 2 -ΔΔCt method and QuantStudio and Design analysis software (Life Technologies). Table 1 Oligonucleotide sequences used to quantify gene expression in RT-PCR. Gene Symbol Forward (5’−3’) Reverse (5’−3’) CEACAM8 ATGTGAAATACAGAACCCAG CATTGACAGACCAAGAATACTG CEACAM6 AACTCTTGGTATTACCCTCC TGAGTTTTGTAATTCCAGCC CILP2 ACTCAGAGTCAGACAAGAAC CAAACTCCATTCCATACCAC COL17A1 ATGGACAGAATTGGATGATG CTGTCACAATTTTGGTCTCC INHBA GAGAATTTGCTGAAGAGGAG CAAAAGTTGTTGTGATTGCC LCN2 GGAAAAAGAAGTGTGACTACTG GTAACTCTTAATGTTGCCCAG LTF CTGGCAGAGAACTACAAATC AAGATATCCTTCCACAGGTC MMP8 TGACAGAGACCTCATTTTCC TATCCTTGGGATAACCTTGC OLFM4 CCCAGTTGTTTTCCAATTTC CTTCTGAGAAAGAACATGAGC AZU1 GCTGCTTCCAAAGCCAGA GCTGCTGATGGAAAACGT DEFA4 TGCCTAGCTTGAGGATCTGTC CATCACCTCTTGCCTGGAGT ELANE TTCCTCGCCTGTGTCCTG CTGCAGGGACACCATGAA PRTN3 ACGCGGAGAACAAACTGAAC ACTGGCTGGTCCTGCTGT ACTB AGCCACATCGCTCAGACAC GCCCAATACGACCAAATCC GAPDH CTGTGCTGTGGAAGCTAAGT GATGTCCACGTCACACTTCA S. Salvador-Martín et al. Biomedicine & Pharmacotherapy 173 (2024) 116299 4 2.10. Statistical analysis Qualitative variables were expressed as frequency and percent. Differences between qualitative variables were calculated using the Mann- Whitney test. Quantitative continuous variables were expressed as median and interquartile range (IQR). Differences between quantitative variables were calculated using the chi-squared test or Fisher’s exact test (when n<30), respectively. Normality of distribution was evaluated using Kolmogorov-Smirnov test or Shapiro Wilk (when n<30) and excluded when p value <0.05. The statistical analysis was performed using IBM SPSS Statistics, Version 26.0 (IBM Corp., Armonk, NY, USA). Gene expression was quantified using ExpressionSuite v1.1 (Applied Biosystems, Foster City, CA, USA) with the non-responder sample (GM008) as a reference, and data were reported as the mean relative quantification and maximum and minimum relative quantification; all the assays were run in triplicate. Efficiency was calculated for each primer pair probe and used for correction. The statistical analyses and graphs were performed using GraphPad Prism v5.1 (Dotmatics, Boston, MA, USA). Differences in gene expression between responders and non-responders were tested using the Mann-Whitney test, with a confidence level of 95%. For all tests, a p value <0.05 was considered statistically significant. 3. Results 3.1. Patients A total of 129 children meeting the inclusion criteria were recruited and gave their written informed consent. Complete data for evaluation of the primary response were collected from the clinical records of 124 patients. Characteristics of all recruited patients prior to the initiation of anti-TNF therapy are summarized in Table 2. Twenty-one were classified as non-responders (14 known at the time of the validations experiments) and 103 as responders to anti-TNF drugs. For the RNA-seq study, 12 responders and 12 non-responders were included. The cohort used in this study included up to 12 patients analyzed elsewhere [18]. The demographic characteristics prior to the initiation of anti-TNF therapy of the discovery cohort (n =24) are summarized in Table 3. Most patients were female (58.3%), with a median age at diagnosis of 10.2 years, and had been diagnosed with CD (66.7%) and treated with infliximab (62.5%). The only statistically significant difference between the 2 groups was in months from diagnosis until the start of treatment (1.9 in responders versus 16.3 in non-responders). For the validation study, we used all the non-responders known at that time (n =14) and a cohort of similar size (n =16) of responders adjusted for sex, age, type of IBD, and anti-TNF drug. These patients had not been included in the previous RNA-seq experiments. The demographic characteristics of this validation cohort prior to the initiation of anti-TNF therapy are summarized in Table 4. Statistically significant differences between the groups were recorded for 2 parameters: months from diagnosis until the start of treatment (2.1 in responders versus 16.3 in non-responders) and PCDAI at the start of treatment (36.3 in responders versus 17.5 in non-responders). PCDAI or PUCAI, CRP and FC are shown at the time of response definition (14 weeks infliximab, 26 weeks adalimumab) (Table 5). All of these parameters are increased in non-responders compared to responders. 3.2. Genes associated with response to anti-TNF drugs A total of 102 genes were differentially expressed (FC >0.6 or < –0.6) between responders and non-responders after the bioinformatic analysis (p value adjusted by FDR <0.05). Of these, 99 genes were overexpressed in responders versus non-responders, while only 3 were downregulated (Suppl. Table S2). Among the top 10 we found DEFA3, MMP8, COL17A1, AP002990.1, OLFM4, INHBA, CEACAM6, LCN2, LTF and CILP2 (Table 6), all of which were over-expressed in responders. 3.3. Functional analysis of differentially expressed genes Functional enrichment analyses were performed on the 102 differentially expressed genes (FDR <0.05, Log2FC >0.6, and log2FC < –0.6). The heatmap in Fig. 1 displays, for each gene (x-axis) involved in the respective enriched biological processes from the Gene Ontology database (GO) (y-axis), the resulting differential expression values obtained through the DEseq2 pipeline (defined by Log2 FoldChange values represented by the color gradient on the right side of the figure). The GO biological process enrichment results corresponding to 52 over- and under-expressed genes included 28 entities, while the GO molecular functions included 6 entities (Suppl. Table S3). As seen in the heatmap, some of these genes, namely, AHSA, INHBA, DEFA3, LTF, PRTN3, MMP8, BPI, and ELANE, have high differential expression values (Log2FoldChange values ranging between 2 and 6). Additionally, Suppl. Fig. 1 shows the number of significant DE genes (x-axis) involved in the Table 2 Characteristics of 124 patients. Characteristic Overall (n =124) Responders (n =103) Non- Responders (n =21) p Value Gender Male, n (%) 69 (55.6%) 60 (58.3%) 9 (42.9%) 0.23 Female, n (%) 55 (44.4%) 43 (41.7%) 12 (57.1%) Age (years) At diagnosis, median (IQR, range) 11.9 (3.9, 1–17) 12 (4, 1−17) 11.5 (6.5, 1–15) 0.26 At start of treatment, median (IQR, range) 12.6 (4.5, 1.1–17.3) 12.7 (4.7, 1.4–17.3) 12 (4.6, 1.1–16.3) 0.43 Type of IBD CD, n (%) 97 (78.2%) 83 (80.6%) 14 (66.7%) UC, n (%) Indet, n (%) 24 (19.4%) 3 (2.4%) 18 (17.5%) 2 (1.9%) 6 (28.6%) 1 (4.8%) 0.21 Type of Anti-TNF Infliximab, n (%) 77 (62.1%) 64 (62.1%) 13 (61.9%) 1 Adalimumab, n (%) Line treatment 1st 2nd Months until treatment 47 (37.9) 117 (94.4%) 7 (5.6%) 6.1 (14.4, 0–96.1) 39 (37.9%) 100 (97.1%) 3 (2.9%) 5.3 (12.5, 0–96.1) 8 (38.1%) 17 (81%) 4 (19%) 7.4 (19.4, 0.1–55.5) 0.02* 0.03* PCDAI at start of treatment, median (IQR, range) 32.5 (27.5, 0–75) 35 (27.5, 0–75) 20 (13.7, 8–63) 0.04* PUCAI at start of treatment, median (IQR, range) 47.5 (28.8, 0–75) 50 (23.8, 5–65) 40 (45, 0−75) 0.45 CRP at start of treatment, median (IQR, range) 14.5 (30.4, 0–130) 16.5 (32.6, 0–130) 10.8 (22.1, 1–127) 0.98 FC at start of treatment, median (IQR, range) Concomitant immunomodulator at start of treatment 1800 (2307, 27–9543) 2000 (2248, 27–9543) 1175 (1837.5, 288–7595) 0.30 Azathioprine, n (%) 91 (76.5%) 75 (75.8%) 16 (80%) 0.44 Methotrexate, n (%) 6 (5%) 3 (2.9%) 3 (15%) None, n (%) 22 (18.5%) 21 (21.2%) 1 (5%) IBD, inflammatory bowel disease; CD, Crohn’s disease; UC, ulcerative colitis; IQR, interquartile range; PCDAI, Pediatric Crohn’s Disease Activity Index; PUCAI, Pediatric Ulcerative Colitis Activity Index; CRP, C-reactive protein; FC, fecal calprotectin. * p value <0.05, ** p value <0.01. The available n is not always 124 since data for some variables are not available for all patients S. Salvador-Martín et al. Biomedicine & Pharmacotherapy 173 (2024) 116299 5 enriched GO biological processes obtained (y-axis). The color gradient displays the FDR values for each GO term, with the functions at the bottom of the figure being the most significant (lowest FDR values). Furthermore, the size of the dots represents the ratio of significant DE genes associated with the corresponding function to the total number of significant DE genes obtained from the differential expression analysis performed using the DESeq2 workflow. Suppl. Fig. 1 and Suppl. Table S3 show how the most significant biological processes are associated with a higher number of DE genes. However, other processes, such as positive regulation of cell adhesion, are represented by as many genes as the most significant processes (Suppl. Fig. 1). The most significant biological processes were defense response to Gram-negative bacterium (FDR = 2.3 ×10 −7 ), erythrocyte development (FDR =2.3 ×10 −7 ), myeloid cell differentiation (FDR =2.3 ×10 −7 ), and neuroinflammatory response (FDR =0.0126) (Suppl. Table S3). The same analysis was performed on GO molecular functions. The most significant molecular functions were hemoglobin binding (FDR = 0.002) and serine-type endopeptidase activity (FDR =0.003) (Fig. 2, Suppl. Table S3). The proportion of differentially expressed genes from the total number of genes for each molecular function included the same processes (Suppl. Fig. 2). The KEGG pathways enrichment results for the same differentially expressed genes showed a unique entity, transcriptional misregulation in cancer (FDR 6.83 ×10 −4 ) (Fig. 3). The gradient color represents the expression value (Log2FoldChange values) corresponding to the significant DE genes obtained using the DESeq2 pipeline. We mapped DE genes in some cancer types, such as Ewing’s sarcoma, epithelial cancer, and or T/B lymphoblastic leukemias. In summary, the pathway yielded 8 DE genes, some of which were highly overexpressed, for example, defensins (DEFA1B, DEFA3, and DEFA4), ETS transcription factor (ERG), and neutrophil elastase (ELANE). GSEA results from a pre-ranked list of genes based on log2FC values showed 9 enriched gene sets for specific GO biological processes (Supl. Table S4). The most significant processes were defense response to Gram-positive bacterium (FDR =1.40 ×10 −5 ), antibacterial humoral response (FDR =0.002), and innate immune response in mucosa (FDR = 0.004). GSEA results for enriched gene sets for specific KEGG pathways also showed a significant result for transcriptional misregulation in cancer (FDR =0.016). Supp. Fig. S3 shows statistically significant GSEA results based on a Table 3 Characteristics of patients in RNA-seq. Characteristic Overall (n =24) Responders (n =12) Non- Responders (n =12) p Value Sex Male, n (%) 10 (41.7%) 5 (41.7%) 5 (41.7%) 1 Female, n (%) 14 (58.3%) 7 (58.3%) 7 (58.3%) Age (years) At diagnosis, median (IQR, range) 10.2 (3.33, 1–16) 10.2 (5.4, 8–16) 10.1 (5.4, 1–15) 0.55 At start of treatment, median (IQR, range) 11.4 (4.7, 1.1–16.3) 10.3 (5.2, 7.9–15.9) 12 (4.8, 1.1–16.3) 0.76 Type of IBD CD, n (%) 16 (66.7%) 8 (66.7%) 8 (66.7%) 1 UC, n (%) Unclassified, n (%) 7 (29.2%) 1 (4.2%) 4 (33.3%) 0 3 (25%) 1 (8.3%) Type of Anti-TNF Infliximab, n (%) 15 (62.5%) 7 (58.3%) 8 (66.7%) 1 0.22 0.002** Adalimumab, n (%) Line of treatment 1st anti-TNF, n (%) 2nd anti-TNF, n (%) Months until treatment 9 (37.5%) 21 (87.5%) 3 (12.5%) 5.2 (15.8, 0.1–55.5) 5 (41.7%) 12 (100%) 0 1.9 (5, 0.1–12.2) 4 (33.3%) 9 (75%) 3 (25%) 16.3 (19.8, 0.1–55.5) PCDAI at start of treatment, median (IQR, range) 26.3 (21.3, 0–63) 30 (15, 0−53) 17.5 (15, 8−63) 0.26 PUCAI at start of treatment, median (IQR, range) 55 (20, 40−75) 55 (21.3, 40–65) 55 (na, 45–75) 0.86 CRP at start of treatment, median (IQR, range) 17.4 (19.8, 0–54) 19 (34.3, 0–52) 9.3 (12.3, 5–54) 0.88 FC at start of treatment, median (IQR, range) Concomitant immunomodulator at start of treatment Yes, n (%) No, n (%) Type of immunomodulator 2464 (3831.5, 148–7595) 20 (83.3%) 4 (16.7%) 3136 (3509.9, 148–6000) 8 (66.7%) 4 (33.3%) 2027.2 (3305, 288–7595) 12 (100%) 0 0.35 0.09 Azathioprine, n (%) 18 (75%) 8 (66.7%) 10 (83.3%) 0.19 Methotrexate, n (%) 2 (8.3%) 0 2 (16.7%) None, n (%) 4 (16.7%) 4 (33.3%) 0 IBD, inflammatory bowel disease; CD, Crohn’s disease; UC, ulcerative colitis; Undet, IBD-unclassified; IQR, interquartile range; PCDAI, Pediatric Crohn’s Disease Activity Index; PUCAI, Pediatric Ulcerative Colitis Activity Index; CRP, C-reactive protein; FC, fecal calprotectin. na, not applicable. * p value <0.05, ** p value <0.01. Table 4 Characteristics of patients included in RT-PCR validation. Characteristic Overall (n =30) Responders (n =16) Non- Responders (n =14) p Value Gender Male, n (%) 14 (46.7%) 9 (56.3%) 5 (35.7%) 0.30 Female, n (%) 16 (53.3%) 7 (43.8%) 9 (64.3%) Age (years) At diagnosis, median (IQR, range) 10.5 (5.2, 1–16) 10.6 (5.61, 2–16) 10.14 (5.77, 1–15) 0.42 At start of treatment, median (IQR, range) 11.9 (5, 1.1–16.3) 12.2 (5.5, 2–15.9) 11.9 (51.1–16.3) 0.82 Type of IBD CD, n (%) 23 (76.7%) 13 (81.3%) 10 (71.4%) 0.50 UC, n (%) Undet, n (%) 6 (20%) 1 (3.3%) 3 (18.8%) 0 3 (21.4%) 1 (7.1%) Type of Anti-TNF Infliximab, n (%) 23 (76.7%) 13 (81.3%) 11 (78.6%) 0.68 Adalimumab, n (%) Line treatment 1st 2nd Months until treatment 7 (23.3%) 27 (90%) 3 (10%) 5.7 (20.6, 0.1–55.5) 3 (18.8%) 16 (100%) 0 2.1 (11.5, 0.1–31.2) 7 (21.4%) 11 (78.6%) 3 (28.6%) 16.3 (19.7, 0.1–55.5) 0.09 0.01* PCDAI at start of treatment, median (IQR, range) 25 (31.25, 8–65) 36.3 (25, 18−65) 17.5 (12.5, 8–63) 0.007** PUCAI at start of treatment, median (IQR, range) 55 (27.5, 40–75) 52.5 (na, 40–65) 55 (na, 45–75) 0.80 CRP at start of treatment, median (IQR, range) 11 (12.3, 0–88) 12.5 (15.8, 0–88) 9.3 (12.1, 5–54) 0.92 FC at start of treatment, median (IQR, range) Concomitant immunomodulator at start of treatment Yes No Type of immunomodulator 1594.5 (2048.5, 176–7595) 27 (90%) 3 (10%) 1662.5 (2055.6, 126–6000) 13 (81.3%) 3 (18.8%) 905.5 (2469.5, 288–7595) 14 (100%) 0 0.98 0.14 Azathioprine, n (%) 25 (83.3%) 13 (81.3%) 12 (85.7%) 0.40 Methotrexate, n (%) 2 (6.7%) 0 2 (14.3%) None, n (%) 3 (10%) 3 (18.8%) 0 IBD, inflammatory bowel disease; CD, Crohn’s disease; UC, ulcerative colitis; IQR, interquartile range; PCDAI, Pediatric Crohn’s Disease Activity Index; PUCAI, Pediatric Ulcerative Colitis Activity Index; CRP, C-reactive protein; FC, fecal calprotectin. na, not applicable. * p value <0.05, ** p value <0.01. S. Salvador-Martín et al. Biomedicine & Pharmacotherapy 173 (2024) 116299 6 pre-ranked list of log2FC values for all genes showing enriched gene sets for specific GO biological processes and specific KEGG pathways among significant log2FC values. 3.4. Validation of RNA-seq by semiquantitative real-time PCR The 12 most differentially expressed genes were selected for validation using RT-PCR with KiCqStart® oligonucleotides in a new cohort of responder patients. DEFA3 and AP002990.1 were ruled out because no oligonucleotides were available. In addition, DEFA4, ELANE, and AZU1 were selected, since they were found to be differentially expressed Table 5 Clinical characteristics of patients after 14 weeks (infliximab) or 26 weeks (adalimumab) of treatment. RNA-seq Overall (n ¼24) Responders (n ¼12) Non- Responders (n ¼12) p Value PCDAI at 14 or 26 weeks, median (IQR, range) 10 (16.3, 0–38) 5 (10, 0−15) 18.8 (22.5, 8–38) 0.01* PUCAI at 14 or 26 weeks, median (IQR, range) 10 (30, 0−40) 10 (na, 0–10) 40 (0, 40−40) 0.16 CRP at 14 or 26 weeks, median (IQR, range) 0.9 (8.2, 0–11) 0.4 (0.8, 0–3) 8.7 (4.5, 1–11) 0.003** FC at 14 or 26 weeks, median (IQR, range) 325 (615.6, 32–1838) 264 (583, 32–1178) 326 (na, 324–1838) 0.31 qPCR validation Overall (n =30) Responders (n =16) Non- Responders (n =14) p Value PCDAI at 14 or 26 weeks, median (IQR, range) 8.8 (15, 0−38) 5 (9.4, 0–20) 20 (16.3, 8–38) 0.001** PUCAI at 14 or 26 weeks, median (IQR, range) 10 (na, 0–40) 5 (na, 0–10) 40 (0, 40−40) 0.22 CRP at 14 or 26 weeks, median (IQR, range) 1.2 (7.1, 0–11) 0.7 (1.6, 0–4) 8.7 (4.5, 1–11) 0.003** FC at 14 or 26 weeks, median (IQR, range) 310.6 (456, 37–1838) 171.5 (371, 37−613) 717 (1311, 324–1838) 0.04* All Overall (n =124) Responders (n =103) Non- Responders (n =21) p Value PCDAI at 14 or 26 weeks, median (IQR, range) 2.50 (10, 0−48) 0 (7.5, 0–48) 20 (12.5,8–38) 0.00002** PUCAI at 14 or 26 weeks, median (IQR, range) 5 (20, 0−40) 0 (10, 0−30) 30 (na, 20–40) 0.01* CRP at 14 or 26 weeks, median (IQR, range) 1 (5.4, 0–100) 0.9 (3.7, 0–55) 9.5 (10.3, 1–100) 0.0003** FC at 14 or 26 weeks, median (IQR, range) 293 (522.1, 10–6685) 205 (491, 10–5100) 847 (1514, 85–6685) 0.04* IQR, interquartile range; PCDAI, Pediatric Crohn’s Disease Activity Index; PUCAI, Pediatric Ulcerative Colitis Activity Index; CRP, C-reactive protein; FC, fecal calprotectin. na, not applicable; * p value <0.05; ** p value <0.01. Table 6 Top 35 differentially expressed genes according to response to anti-TNF drugs in children prior to initiation of treatment. Gene Symbol Ensembl Gene Base Mean Log2FC adj p value Validation DEFA3 ENSG00000239839 7793 6.39 5.33 ×10–7 MMP8 ENSG00000118113 824 4.85 1.29 ×10–4 YES COL17A1 ENSG00000065618 41 4.69 4.44 ×10–4 YES AP002990.1 ENSG00000255508 675 4.63 8.43 ×10–4 OLFM4 ENSG00000102837 701 4.57 2.94 ×10–3 YES INHBA ENSG00000122641 13 4.29 2.00 ×10–3 YES CEACAM6 ENSG00000086548 152 4.27 3.44 ×10–3 YES LCN2 ENSG00000148346 1064 4.15 1.30 ×10–4 YES LTF ENSG00000012223 1136 3.96 2.43 ×10–3 YES CILP2 ENSG00000160161 2 3.95 2.55 ×10–2 YES CEACAM8 ENSG00000124469 221 3.91 2.00 ×10–3 YES PRTN3 ENSG00000196415 61 3.84 2.69 ×10–2 YES HTRA3 ENSG00000170801 15 3.71 4.61 ×10–3 LINC02009 ENSG00000283646 45 3.69 8.87 ×10–3 BPI ENSG00000101425 519 3.64 8.44 ×10–4 SCGB1C1 ENSG00000188076 2 3.59 4.31 ×10–2 AC104232.1 ENSG00000254006 7 3.59 4.72 ×10–2 ATP2C2 ENSG00000064270 70 3.56 4.12 ×10–3 ELANE ENSG00000197561 101 3.55 2.72 ×10–3 YES PCOLCE2 ENSG00000163710 11 3.38 3.74 ×10–2 CRISP3 ENSG00000096006 232 3.31 4.56 ×10–3 CHIT1 ENSG00000133063 101 3.27 1.04 ×10–2 DEFA4 ENSG00000164821 189 3.17 1.98 ×10–2 YES DEFA1B ENSG00000240247 8110 3.13 3.28 ×10–2 TARM1 ENSG00000248385 9 2.95 1.04 ×10–2 MS4A3 ENSG00000149516 241 2.92 1.20 ×10–2 RETN ENSG00000104918 199 2.92 4.61 ×10–3 HP ENSG00000257017 717 2.75 9.11 ×10–3 PDGFRA ENSG00000134853 2 -2.73 2.17 ×10–2 AZU1 ENSG00000172232 98 2.69 1.33 ×10–2 YES ABCA13 ENSG00000179869 30 2.64 1.87 ×10–2 AC016526.3 ENSG00000280078 2 2.64 8.87 ×10–3 TCN1 ENSG00000134827 283 2.58 1.22 ×10–2 AC100835.2 ENSG00000275527 37 2.54 5.84 ×10–3 CA1 ENSG00000133742 843 2.50 1.80 ×10–2 S. Salvador-Martín et al. Biomedicine & Pharmacotherapy 173 (2024) 116299 7 in another immune-mediated disease, multiple sclerosis, in a previous work by our group [30]. The genes finally selected for validation were MMP8, COL17A1, OLFM4, INHBA, CEACAM6, CEACAM8, LCN2, LTF, CILP2, PRTN3, ELANE, DEFA4 and AZU1. Expression of CEACAM8 (2.59-fold, p value 0.0236), LCN2 (2.00- fold, p value 0.0126), LTF (2.64-fold, p value 0.0083), and PRNT3 (3.10- fold, p value 0.0235) was significantly higher in responders than in nonresponders before initiation of treatment. However, the results showed higher expression of all genes in responders than in non-responders (Fig. 4). Comparison between fold change values in gene expression between RNAseq and RT-PCR shows that although there is no linearity in correlation, they were all upregulated in both techniques (Suppl. Table S5). 4. Discussion Identifying patients diagnosed with IBDs who will respond to biological therapies is a priority for gastroenterologists. In this study, we analyzed the whole gene expression of 12 responders and 12 nonresponders to anti-TNF drugs in blood from children with IBDs. RNA- seq identified 102 genes that were differentially expressed in these groups before initiation of treatment. A selection of 13 genes that were upregulated in responders to anti-TNF drugs were also more highly expressed in a partially different population and using RT-PCR; statistical differences were found for 4 of them. Defense response to Gramnegative bacterium was the main enriched biological process, while hemoglobin binding and serine-type endopeptidase activity were included among the most significant molecular functions. Functional analysis of DE genes between children who respond and do not respond to anti-TNF drugs reveals a transcriptional misregulation in cancer. In this study, we identified 102 DE genes between children who responded and did not respond to anti-TNF agents in blood from children with IBDs, thus highlighting their potential implication as biomarkers of anti-TNF response. Of these, 99 genes were over-expressed in responders versus non-responders, while only 3 were downregulated. RNA-seq was recently shown to be an efficient tool for identifying gene expression–based markers of response to biological therapy identified in biopsies from adult IBD patients (ACTN1, CXCL6, LAMA4, EMILIN1, Fig. 1. Enriched GO terms for DE genes. Enriched Gene Ontology biological processes (FDR <0.05) in over-expressed and under-expressed genes (log 2 FC >0.6 and log 2 FC <–0.6, respectively). The number of genes involved in the set of enriched biological process is shown within the heatmap and the gradient color defines the log 2 FC for each DE gene. S. Salvador-Martín et al. Biomedicine & Pharmacotherapy 173 (2024) 116299 8 CRIP2, CXCL13, and MAPKAPK2) [31]. A recent study using 5 Gene Expression Omnibus datasets from biopsies of patients with IBDs showed overexpression of SELE, TREM1, AQP9, FPR2, and HCAR3 in those who responded to infliximab [32]. Another study in UC patients revealed a different signature (S100A8, S100A9, TREM1, TLR2, IL1B, and FPR1) [33]. This diversity of results gives us an idea of the difficulty in finding a small set of DE genes of response to anti-TNF drugs. None of these genes was identified in blood from children with IBDs in our study. This fact may reflect strong differences in gene expression biomarkers between tissues, type of IBD, and the patient’s age. In this study, we identified 102 genes that were differentially expressed in blood from children with IBDs who responded or not to anti-TNF therapy before the initiation of treatment. Testing blood rather than tissue constitutes an extraordinary advantage, because it is less invasive. The PANTs study (324 patients) did not identify genes regulated in response to anti-TNF drugs in the blood of patients aged ≥6 years with CD at week 0. The main difference between our study and the PANTs study was the age of the participants, since PANTs also included adults. Similarly, in our previous study, which was based on half the number of the current sample, we did not identify any DE gene with an FDR >0.05 [18]. The 4 genes validated in our cohort, namely, PRTN3, LTF, LCN, and CEACAM8, show the importance of neutrophils in the response to anti- TNF therapy. PRTN3 was one of the 4 genes validated in our study. A relationship has been identified between IBDs and the generation of antineutrophil cytoplasmatic antibodies, such as those against PRTN3 [34,35]. This relationship is even useful for distinguishing UC from CD in children [36,37]. Measurement of proteinase 3 antineutrophil cytoplasmic antibody (PR3-ANCA) is useful not only for diagnosing UC, but also for evaluating disease severity and extension and predicting clinical course. PR3-ANCA ≥3.5 U/mL demonstrated 44.5% sensitivity and 95.6% specificity for the diagnosis of UC and 75% sensitivity and 69.0% specificity for steroid requirements in UC patients [38]. Stool levels of myeloperoxidase, another component of azurophil granules, increased significantly in active IBD patients, correlating with laboratory parameters and endoscopic grade of inflammation [39]. In role in autoimmune diseases, it seems to modify adaptive immune responses and induce vascular permeability [40]. LTF codes for lactoferrin and was another of the validated genes overexpressed in responders to anti-TNF drugs. Fecal lactoferrin (FL) Fig. 2. Enriched GO terms for DE genes. Enriched Gene Ontology molecular functions (FDR <0.05) in over-expressed and under-expressed genes (log 2 FC >0.6 and log 2 FC <–0.6, respectively). The number of genes involved in the set of enriched biological process is shown within the heatmap, and the gradient color defines the log 2 FC for each DE gene. S. Salvador-Martín et al. Biomedicine & Pharmacotherapy 173 (2024) 116299 9 was recently shown to be a useful biomarker of response in IBDs. Patients with FL levels >84.5 µg/mL at week 4 presented better clinical remission, endoscopic improvement, and remission (OR 0.43, 0.46, and 0.40 respectively) [41]. This finding supports the higher expression of LTF in blood we observed in children who responded to anti-TNF therapy before starting treatment and its potential use as a biomarker of response. LCN2 codes for lipocalin-2, a protein contained in neutrophilic granules, and exerts a protective effect against pathogens in intestinal inflammation. It is upregulated in serum and biopsies from colonoscopy in patients diagnosed with IBDs [42]. A 62.1% decrease in urinary levels of LCN2 was observed in CD patients after treatment with infliximab [43]. Our results suggest that individuals with higher LCN2 expression in blood benefit more from anti-TNF therapy. CEACAM8, also known as CD66b, was upregulated in children with IBDs who responded to anti-TNF therapy. The expression of CD66b in neutrophils from mucosa, but not blood, correlated with the severity of CD [44]. The functional analysis revealed defense response to Gram-negative bacterium to be the main biological process enriched in the set of genes evaluated, while hemoglobin binding and serine-type endopeptidase activity were included in the most significant molecular functions. In our study, the genes that were more highly expressed in children with IBDs who responded to anti-TNF drugs than in those who did not and contributed more to defense response to Gram-negative bacterium were LTF, BPI, OPTN, CAMP, DEFA4, AZU1, ELANE, DEFA3, and DEFA1B, with most of them also being serine-type endopeptidases (LTF, MMP8, HTRA3, AZU1, PRTN3, ELANE, and HP). AZU1, PRTN3, ELANE and DEFA4 have been associated with methylprednisolone in relapsingremitting multiple sclerosis [45], and in sepsis [46]. AZU1, ELANE, and PRTN3 are expressed coordinately, and are present, together with DEFA4, in azurophil granules, which are specialized lysosomes of the neutrophil for killing of microorganisms [47]. However, they are also involved in the formation of neutrophil extracellular traps, which, in addition to their role in defense against microorganisms, have also been implicated in cancer and autoimmune processes [47]. An association between neutrophil azurophilic granule exocytosis and TNF-alpha has also been described [48]. According to GO, hemoglobin binding proved to be the most significant molecular function in our study population. A relationship between indicators of iron status and indicators of inflammation has been described in children and adolescents with IBDs [49]. The authors found lower hemoglobin levels and iron-binding capacity in children with IBDs after 2–6 months of treatment with infliximab than in healthy children. Similarly, the presence of anemia is a biomarker of poor response to infliximab [50]. The higher expression of genes involved in hemoglobin binding might explain an improvement in anemia values and contribute to an improvement in IBDs in children receiving anti-TNF therapy. Functional analysis of DE genes between children who respond and do not respond to anti-TNF drugs reveals a transcriptional misregulation in cancer due mainly to the inclusion of defensin genes. Regulation of gene expression plays an essential role in IBD-associated colorectal cancer, as shown in preclinical animal models [51]. The association between chronic diseases and increased risk of cancer is well known [52]. For example, colorectal cancer is 1.5–2 times more common in IBD patients than in the general population and is associated with up to 15% Fig. 3. Mapping of DE genes in transcriptional misregulation in cancer. Overview of differential expression within the transcriptional misregulation in cancer enriched pathway associated with the KEGG database. The DE genes with FDR <0.05 belonging to the pathway are highlighted in red or green depending on their degree of over-expression or under-expression, respectively. S. Salvador-Martín et al.