Abstract Background The development of factor VIII (FVIII) inhibitor is a severe complication during replacement therapy for hemophilia A patients. Objectives We investigated the potential risk factors for FVIII inhibitor formation based on genome‐wide RNA‐sequencing and whole‐genome bisulfite sequencing analysis. Methods RNA‐sequencing and whole‐genome bisulfite sequencing analysis were applied on 17 blood samples with F8 intron 22 inversion, including seven with inhibitors and 10 without. Results Altogether, 344 mRNA transcripts and 20 long noncoding RNAs (lncRNA) transcripts were differentially expressed. Among the differentially expressed transcripts, 200 mRNAs and 12 lncRNAs were upregulated, and 144 mRNAs and eight lncRNAs were downregulated. Gene ontology enrichment analysis of differentially expressed mRNAs showed that genes involved in immune stimulation, especially those for T‐cell activation, were upregulated, whereas genes involved in negative immune response regulation were downregulated. Coexpression analysis revealed that the targeted upregulated genes of differentially expressed lncRNA were similarly closely related to immune activation, especially T‐cell activation. Methylation analysis showed inhibitor patients exhibited a slightly lower methylation status in the CpG islands, 5′ untranslated region, and exon regions (p < 0.01). Genes with differentially methylated regions were also related to T‐cell activation. Conclusions There is an upregulation of genes involved in activation of the immune system in hemophilia A patients with inhibitors. The lncRNA and methylation modifications may play important roles in inhibitor production. These findings are potentially to reveal novel therapeutic targets for prevention and treatment of inhibitors. Keywords: DNA methylation, FVIII inhibitor, hemophilia A, long‐noncoding RNA (lncRNA), RNA‐seq __________________________________________________________________ Essentials. * The risk factors for factor VIII (FVIII) inhibitor in hemophilia A have not been fully established. * RNA sequence and DNA methylation information were studied in hemophilia A patients. * In patients with inhibitors, upregulated genes are related to immune response. * Long noncoding RNA and methylation modifications may be involved in FVIII inhibitor formation. 1. INTRODUCTION Hemophilia A is an X chromosome‐linked recessive bleeding disorder characterized by the deficiency or dysfunction of coagulation protein factor VIII (FVIII). The incidence rate of hemophilia A is approximately one in 5000 male births. Recurrent joint and muscle bleeding are the major clinical manifestations. The severity of hemophilia is defined according to FVIII concentration in plasma.[46] ^1 FVIII replacement therapy is the primary choice for hemophilia A treatment. However, neutralizing antibodies against exogenous FVIII may be produced in as many as 20%–30% of patients with severe hemophilia A,[47] ^2 affecting treatment effectiveness, increasing morbidity, and raising the management cost.[48] ^3 FVIII inhibitor emergence therefore represents a serious complication during replacement therapy in hemophilia A. These inhibitors can often be eradicated by immune tolerance induction treatment by repeated FVII infusions (every other day to once or twice daily) over months to years.[49] ^3 However, the treatment burden and high cost limit its application particularly in countries with economic constraints. FVIII inhibitor development is related to several external environmental and intrinsic genetic factors.[50] ^4 , [51]^5 An unequivocal explanation for inhibitor development has not been fully established over the past 2 decades. To improve the efficacy of replacement treatment, it is essential to understand the mechanism and pathways of inhibitor emergence in hemophilia A patients. Considering that the immune system is theoretically remodeled in patients who had developed inhibitors, exploring their gene expression alteration therefore may reveal the risk factors of inhibitor production. Additionally, factors regulating gene expression may also be involved. Epigenetic regulation plays a significant role in gene expression. Furthermore, long noncoding RNAs (lncRNAs) are newly highlighted RNA molecules that regulate gene expression by diverse pathways,[52] ^6 although the functions of most lncRNAs are not well verified. Whether anomalous gene regulation by epigenetic modifications and lncRNAs play a role in FVIII inhibitor production remains to be established. In this study, we used RNA‐sequencing (RNA‐seq) technology in hemophilia A inhibitor and noninhibitor patients, to investigate the systematic changes at transcription levels of peripheral blood mononuclear cells. We also performed whole‐genome bisulfite sequencing and compared the methylation levels between inhibitor and noninhibitor patients. We then performed an integrative analysis through the combined bioinformatic data for deeper understanding of the inhibitor‐producing mechanisms. 2. METHODS 2.1. Ethics statement Ethical approval was granted by the Ethics Committee of Blood Diseases Hospital, Chinese Academy of Medical Sciences (reference no. IIT2017005‐EC‐1). 2.2. Blood sample collection The blood sample of hemophilia A patients visiting the Institute of Hematology and Blood Diseases Hospital were collected from December 2017 to May 2018. The inclusion criteria were as follows: patients with severe hemophilia A (FVIII < 0.01 IU/ml); hemophilia A patients with inhibitors developed within 75 FVIII exposure days; hemophilia A noninhibitor patients with at least 75 FVIII exposure days; no bleeding events within 2 weeks before enrollment; patients carrying F8 gene intron 22 inversion mutation; and patients consenting to participate in our study. 2.3. Definition of positive inhibitor A positive FVIII inhibitor was defined as having an inhibitor titer of ≥0.6 BU/mL as measured by the Bethesda assay, or the Nijmegen modified Bethesda assay.[53] ^7 Low‐titer inhibitor was defined as FVIII inhibitor titer <5 BU/ml. High‐titer inhibitor was defined as inhibitor ≥5 BU/ml. 2.4. RNA‐seq and quantitative real‐time polymerase chain reaction analyses Peripheral blood mononuclear cells of hemophilia A patients were isolated from 3–5 ml EDTA‐K2 anticoagulant venous blood using density gradient centrifugation. Patients were divided into the inhibitor‐positive and inhibitor‐negative control groups. Total RNA was extracted using TRIzol reagent (Life Technologies, Carlsbad, CA, USA) for transcriptomics and quantitative real‐time polymerase chain reaction (qRT‐PCR) analyses (Table [54]S1). 2.5. Gene ontology and Kyoto Encyclopedia of Genes and Genomes pathway analyses The Metascape Gene ontology (GO) enrichment and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway databases were used to analyze all differentially expressed mRNAs as well as lncRNA coexpressing differentially expressed mRNAs to interpret the biological meaning of the transcripts. GO terms and KEGG pathways with corrected p values <0.05 were considered significantly enriched. 2.6. Whole genome bisulfite sequencing Whole genome bisulfite sequencing was conducted for the same 17 patients studied with RNA‐seq. DSS software was used to identify differentially methylated regions (DMRs) between the inhibitor‐positive and inhibitor‐negative groups. Differentially methylated genes were defined as genes in which the gene body regions (from transcriptional start site to transcription end sites) or the promoter regions (upstream 2 kb from the transcriptional start site) had overlaps with the DMRs. GO enrichment analysis was performed for annotated differentially methylated genes. 2.7. Statistical analyses All statistical analyses were performed using SAS software, version 9.4. Wilcoxon rank‐sum test was used to compare age distribution between the two groups. Differences were considered statistically significant when p < 0.05. 3. RESULTS 3.1. General information This study included seven hemophilia A patients with inhibitors and 10 without. Table [55]1 presents the general information of the enrolled patients. TABLE 1. The general information of enrolled patients Patient number[56] ^a FVIII:C[57] ^c (IU/ml) Inhibitor titer (BU)[58] ^e ED (d) Age (y)[59] ^f Age (y)[60] ^g IP1 <0.01 0.66 <75[61] ^d 7 12 IP2 <0.01 1587.20 9 26 28 IP3 <0.01 38.00 42 1 3 1P4 <0.01 1.55 14 16 20 1P5 <0.01 15.20 <75[62] ^d 56 62 1P6 <0.01 200.00 53 1 3 1P7[63] ^b <0.01 26.40 7 49 53 IN1 <0.01 <0.6 10 IN2 <0.01 <0.6 23 IN3 <0.01 <0.6 29 IN4 <0.01 <0.6 28 IN5 <0.01 <0.6 44 IN6 <0.01 <0.6 19 IN7 <0.01 <0.6 51 IN8 <0.01 <0.6 34 IN9 <0.01 <0.6 17 IN10 <0.01 <0.6 14 [64]Open in a new tab Abbreviations: ED, exposure day of factor VII at inhibitor development; IN, inhibitor negative; IP, inhibitor positive. ^^a All patients are ethnic Chinese. ^^b This patient number is listed as IP12 in Figure [65]1B. ^^c The factor VIII:C level before inhibitor formation, for patients with inhibitor. ^^d Patients were unable to recall the specific exposure days. ^^e Inhibitor titer at the time of blood sample taken. ^^f Age at the time of inhibitor formation. ^^g Age at the time of blood sample taken. The median age was 20 years (3–62 years) for inhibitor‐positive patients and 25.5 years (10–51 years) for inhibitor‐negative patients (p > 0.05). The median inhibitor titer for the inhibitor‐positive patients was 38.00 BU (0.66–1587.20 BU). Five of the seven inhibitor‐positive patients had high‐titer inhibitors and two had low‐titer inhibitors. 3.2. RNA‐seq results of patients in the inhibitor‐positive and inhibitor‐negative groups and validation of RNA‐seq results by qRT‐PCR We screened the expression profiles of all transcripts in the peripheral blood mononuclear cells from the two groups of patients. More than 80,000 mRNA transcripts and more than 19,000 corresponding genes were identified. Likewise, nearly 5400 lncRNAs and approximately 3495 corresponding gene symbols were identified. Differentially expressed transcripts between the inhibitor‐positive and inhibitor‐negative groups were first analyzed (Figure [66]1A,B). Altogether, 344 mRNA transcripts and 20 lncRNA differentially expressed transcripts were observed according to the threshold of log2FoldChange ≤−1 or ≥1, fragments per kilobase of transcript sequence per million base pairs sequenced ≥0.5 in more than 70% of patients (p < 0.05). Among these, 200 mRNAs and 12 lncRNAs were upregulated, and 144 mRNAs and 6 lncRNAs were downregulated (Table [67]S2–S4). FIGURE 1. FIGURE 1 [68]Open in a new tab Transcription profile of patients. (A) Volcano plot of gene expression Level of mRNAs based on p‐value (<0.05). (B) Hierarchical clustering analysis of differentially expressed mRNA based on p‐value (<0.05) and log2 fold change (≤−1 or ≥1). (C) Gene ontology (GO) enrichment of upregulated differentially expressed mRNAs. (D) GO enrichment of downregulated differentially expressed mRNAs. (E) Kyoto Encyclopedia of Genes and Genomes enrichment of upregulated differentially expressed mRNAs. IN, inhibitor negative; IP, inhibitor positive; lncRNA, long noncoding RNA GO and KEGG enrichment analyses were performed. The 20 most overrepresented GO terms for the upregulated differentially expressed mRNAs are shown in Figure [69]1C. The upregulated differentially expressed transcripts were mainly enriched in the immune response and cellular biological processes. Six of the 20 GO terms were related to T‐cell activation. The rest were related to the positive regulation of leukocyte or lymphocyte activation and differentiation. Downregulated differentially expressed mRNAs were mainly enriched in covalent chromatin, histone, and protein modifications. Seven GO terms were related to the negative regulation of biological processes (Figure [70]1D). The upregulated and downregulated genes involved in the immune response are shown in Tables [71]2 and [72]3. TABLE 2. Upregulated genes related to T‐cell activation Gene name Description log2 fold change p‐value Reported IL1ß Interleukin 1 beta 3.44 <0.001 Yes[73] ^a EGR1 Early growth response 1 2.60 <0.001 No THBS1 Thrombospondin 1 2.14 0.002 No CD83 CD83 molecule 2.07 0.005 No TNFAIP3 TNF alpha induced protein 3 1.88 0.002 No GPR183 G protein‐coupled receptor 183 1.76 <0.001 No IFNG Interferon gamma 1.72 0.017 No PDE4D Phosphodiesterase 4D 1.43 0.005 No ARG2 Arginase 2 1.41 0.009 No SOX4 SRY‐box transcription factor 4 1.41 0.003 No HLA‐DQB1 Major histocompatibility complex, class II, DQ beta 1 1.37 0.005 Yes[74] ^b DEFA3 Defensin alpha 3 1.33 0.018 No HLA‐DRB5 Major histocompatibility complex, class II, DR beta 5 1.27 0.042 Yes[75] ^b H2BC10 H2B clustered histone 10 1.23 0.036 No DUSP10 Dual specificity phosphatase 10 1.22 0.006 No HLA‐DPA1 Major histocompatibility complex, class II, DP alpha 1 1.11 0.013 No MAFB MAF bZIP transcription factor B 1.11 0.025 No RNASE2 Ribonuclease A family member 2 1.07 0.024 No SAMSN1 SAM domain, SH3 domain and nuclear localization signals 1 1.04 0.005 No ZC3H12A Zinc finger CCCH‐type containing 12A 1.00 0.041 No [76]Open in a new tab ^^a Immune polymorphism of IL‐1ß gene has been reported to be associated with inhibitor development. ^^b Alleles of HLA‐DR and HLA‐DQ genes have been reported to be associated with inhibitor development. TABLE 3. Downregulated genes related to negative regulation of immune response Gene name Description log2 fold change p‐value Reported FPR2 Formyl peptide receptor 2 −4.28 0.006 No SERPING1 Serpin family G member 1 −2.73 0.007 No SEC14L1 SEC14 like lipid binding 1 −2.19 0.021 No PCBP2 Poly(rC) binding protein 2 −1.87 0.009 No XPO1 Exportin 1 −1.74 0.002 No DNMT1 DNA methyltransferase 1 −1.53 0.004 No SETDB2 SET domain bifurcated histone lysine methyltransferase 2 −1.44 0.003 No DDX6 DEAD‐box helicase 6 −1.39 0.039 No PHF8 PHD finger protein 8 −1.39 0.022 No POM121 POM121 transmembrane nucleoporin −1.32 0.008 No EED Embryonic ectoderm development −1.27 0.021 No AGO2 Argonaute RISC catalytic component 2 −1.20 0.022 No SFN Stratifin −2.18 0.001 No AMPD2 Adenosine monophosphate deaminase 2 −1.50 0.028 No TEX9 Testis expressed 9 −1.13 0.007 No POLE2 DNA polymerase epsilon 2, accessory subunit −1.10 0.006 No REC8 REC8 meiotic recombination protein −1.08 0.038 No NEIL3 Nei like DNA glycosylase 3 −1.08 0.017 No [77]Open in a new tab The results of the KEGG enrichment analysis are shown in Figure [78]1E. The four most enriched KEGG pathways were influenza A, the Toll‐like receptor signaling pathway, T‐helper cell 17 (Th17) cell differentiation, and the B‐cell receptor signaling pathway. To validate the reproducibility and repeatability of the transcriptome sequencing results, we randomly selected the following 13 genes for qRT‐PCR analysis (Table [79]S1): Plant homeodomain finger protein 8 (PHF8), cellular nucleic acid‐binding protein, male‐specific lethal complex subunit 3, nucleoporin 98 and 96 precursor, Pinin interacting serine and arginine‐rich protein, SOS Ras/Rac guanine nucleotide exchange factor 1, phosphoribosyl pyrophosphate synthetase associated protein 2, AT‐hook transcription factor, E74 like ETS transcription factor 2, mitogen‐activated protein kinase 4, RAB GTPase activating protein 1 like, SET domain bifurcated histone lysine methyltransferase 2, and Rap guanine nucleotide exchange factor 2 (Figure [80]S1). The results were consistent with those of the RNA‐seq findings. 3.3. Potential targets and function of the dysregulated lncRNAs Target prediction programs were used to investigate whether the differentially expressed lncRNAs regulate genes associated with inhibitor formation. We found that 18 of the dysregulated lncRNAs had predicted target genes (Table [81]S4). On further analysis with the profile of differentially expressed lncRNAs and targeted differentially expressed mRNAs, we found nine upregulated lncRNAs coexpressed with 64 upregulated mRNAs (58 genes)(Figure [82]2A). FIGURE 2. FIGURE 2 [83]Open in a new tab Coexpression analysis of upregulated differentially expressed long noncoding RNAs (lncRNAs) and targeted upregulated differentially expressed mRNAs. (A) Coexpression network of upregulate differentially expressed lncRNAs and upregulated differentially expressed mRNAs. (B) Gene ontology (GO) enrichment of upregulated differentially expressed lncRNAs coexpressed upregulated mRNAs. (C) Kyoto Encyclopedia of Genes and Genomes enrichment of upregulated differentially expressed lncRNAs coexpressed upregulated mRNAs. The top 20 enriched GO terms of the upregulated target mRNA were mostly involved in the immune response, especially those for T‐cell activation (Figure [84]2B). Three GO terms were related to T‐cell activation, and eight were related to lymphocyte and leukocyte activation. KEGG pathway enrichment analysis was performed to reveal the functions of the coexpressed regulated target mRNAs (Figure [85]2C). The four most enriched pathways were systemic lupus erythematosus, rheumatoid arthritis, Chagas disease, and the Toll‐like receptor signaling pathway. The tumor necrosis factor, interleukin 17, Th17, and T‐cell receptor signaling pathways were also involved. 3.4. The lower DNA methylation status in some gene regions in patients with inhibitors compared with patients without inhibitors by whole‐genome bisulfite sequencing Similar CG methylation levels were observed between samples of the inhibitor‐positive and inhibitor‐negative groups at the genome‐wide scale. However, the inhibitor‐positive samples exhibited a slightly lower methylation status in the CpG islands, (p < 0.001), 5′ untranslated region (p < 0.001), and exon regions (p < 0.001) (Figure [86]3A). A total of 3589 DMRs in CG were identified, including 1440 hypo‐DMRs and 2149 hyper‐DMRs within 1372 annotated differentially methylated genes. Among the GO enrichment terms, the top 20 immune response‐related GO terms are shown in Figure [87]3B. Nine of these terms were related to T‐cell or lymphocyte differentiation or activation. FIGURE 3. FIGURE 3 [88]Open in a new tab Methylation analysis of patients. (A) Comparison of methylation level between IP and IN at the genome‐wide scale. (B) Immune response‐related gene ontology term enrichment of genes with differentially methylated regions in CG. (C) Comparison of the gene body methylation level of upregulated and downregulated genes between IP and IN. (D) Comparison of the promoter methylation level of upregulated and downregulated genes between the IP and IN groups. CGI, CpG island; down, downregulated genes; IN, inhibitor negative; IP, inhibitor positive; Mc, methylation; up, upregulated genes; Utr3, 3′ untranslated region; Utr5, 5′ untranslated region We compared the methylation levels of the differentially methylated genes between the two groups and found no differences in the promoter and gene body regions (Figure [89]3C,D). 4. DISCUSSION This study compared the gene expression level and methylation level between hemophilia A patients with and without inhibitors. All patients carried the same F8 mutation with intron 22 inversion. GO terms and KEGG pathways enrichment analyses showed that differentially expressed genes were prominently enriched in the regulation of immune responses. The 20 most enriched KEGG pathways of upregulated differentially expressed mRNAs, such as the Toll‐like receptor and tumor necrosis factor signaling pathways, were considered to be related to inhibitor formation in hemophilia A patients.[90] ^8 , [91]^9 Six of the 20 most enriched GO terms of up‐ differentially expressed mRNAs were related to T‐cell activation. This was consistent with the theory that the production of antibodies against exogenous FVIII is a T‐cell‐dependent event.[92] ^10 After recognizing an FVIII epitope, T cells are activated and help B cells to proliferate and differentiate into plasma cells, thus leading to inhibitor development. In our patients, the time point of conducting RNA‐seq analysis was more than 2 years after the diagnosis of inhibitors, indicating a persistent T‐cell‐activated state after the inhibitor production. Treatment‐related factors such as surgery may challenge the immune system as a “danger signal,”[93] ^11 thus promoting the immune reaction to FVII for patients with those highly expressed genes. Evidence suggests that CD4^+ T cells are necessary for FVII inhibitor development. Th17, Th1, and Th2 cells have also been identified as important components of FVIII inhibitor development.[94] ^12 , [95]^13 In our study, Th17 differentiation, interleukin 17 signaling pathway, and Th1‐ and Th2‐cell differentiation pathways were also enriched in patients with inhibitors. T‐follicular helper cells are also related to inhibitor formation.[96] ^14 The expression level of the transcription factor B cell lymphoma‐6, expressed by T‐follicular helper cells, was also higher in patients with inhibitors in our study (Table [97]S2), indicating the activation of T‐follicular helper cells. The role of lncRNAs in inhibitor formation has not yet been studied. GO terms and KEGG enrichment analyses showed that the differentially expressed lncRNAs were also predominantly related to regulating multiple immune response associated genes, especially T‐cell activation. This result suggested a novel regulation role of lncRNA in the T‐cell response against exogenous FVIII. The nine differentially expressed lncRNAs identified in this study may serve as biomarkers and provide therapeutic target for the opportunities of inhibitor elimination. DNA methylation changes can regulate myeloid and lymphoid lineage development and inflammatory gene expression in response to environmental stressors.[98] ^15 Here, the CG methylation level was lower in inhibitor patients in several gene regions. GO analysis showed that differentially methylated genes were also related to immune response, especially T‐cell activation. We speculated that hypomethylation at the DNA level in inhibitor patients might cause T cells to be more easily activated when exposed to FVIII than noninhibitor patients. Many downregulated genes reportedly inhibit immune responses or inflammatory reactions. DNA methyltransferase 1 maintains global methylation after DNA replication and induces Treg cell differentiation.[99] ^16 Treg cells are associated with the induction of FVIII tolerance in patients with hemophilia.[100] ^17 SET domain bifurcated histone lysine methyltransferase 2 suppresses the expression of many genes with proinflammatory functions.[101] ^18 Plant homeodomain finger protein 8 silences genes involved in inflammatory responses.[102] ^19 Deficiency in Nei‐like DNA glycosylase 3 is associated with increased lymphocyte apoptosis, autoantibody production, and autoimmunity.[103] ^20 Adenosine monophosphate deaminase 2 is likely to be involved in the proinflammatory pathway by increasing adenosine triphosphate levels.[104] ^21 Stratifin plays a role in balancing the host inflammatory response to the virus.[105] ^22 Decreased expression of these genes may enhance the immune response or impair the immune tolerance to exogenous FVIII. Thus, these results indicate a decreased ability to suppress the inflammatory response in inhibitor patients. Our study has some limitations. We did not sort the immune cells and could not determine in which cell type the differentially expressed genes were enriched. Single‐cell sequencing may provide more information regarding these cell types. Additionally, our sample size was relatively small. We only selected patients with intron 22 inversions. It is unclear whether research on other types of F8 mutations would yield the same results. In summary, the results of our study reveal that there is an upregulation of genes involved with activation of the immune system in hemophilia A patients with inhibitors. This study is the first one that we are aware of to report the role of lncRNA and epigenetic modifications for inhibitor production in hemophilia A patients. These findings have the potential to reveal novel therapeutic targets for prevention and treatment of inhibitors. Our study also suggests that FVIII replacement therapy should be avoided when there is active inflammation or immune response in the patients. AUTHOR CONTRIBUTIONS W.L. and C.L. conceptualized the work, collected the data, wrote the original draft, and edited the draft. W.W. was involved in overall supervision, critical review, and editing of drafts. F.X., L.C., H.L., Y.C., R.W., and Y.F. collected the data. Y.M. contributed to the acquisition and interpretation of the data. L.Z. and R.Y. conceptualized and planned the work, and were involved in overall supervision, critical review, and editing of drafts. RELATIONSHIP DISCLOSURE L.Z. has received speaker/consultancy fees from Bayer, Novo Nordisk, Pfizer, Roche, and Takeda. R.Y. has received speaker/consultancy fees from Bayer, Novo Nordisk, Pfizer, Roche, and Takeda. The other authors declare no competing interests. Supporting information AppendixS1 [106]Click here for additional data file.^ (495.5KB, pdf) ACKNOWLEDGMENTS