Risk factors for inhibitors in hemophilia A based on RNA-seq and DNA methylation.

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.

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.

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. 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, affecting treatment effectiveness, increasing morbidity, and raising the management cost. 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. 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. , 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, 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.

METHODS

Ethics statement

Ethical approval was granted by the Ethics Committee of Blood Diseases Hospital, Chinese Academy of Medical Sciences (reference no. IIT2017005‐EC‐1).

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.

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.

Low‐titer inhibitor was defined as FVIII inhibitor titer <5 BU/ml.

High‐titer inhibitor was defined as inhibitor ≥5 BU/ml.

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 S1).

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.

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.

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.

RESULTS

General information

This study included seven hemophilia A patients with inhibitors and 10 without. Table 1 presents the general information of the enrolled patients.

TABLE 1

The general information of enrolled patients

Patient number a	FVIII:C c (IU/ml)	Inhibitor titer (BU) e	ED (d)	Age (y) f	Age (y) g
IP1	<0.01	0.66	<75 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 d	56	62
1P6	<0.01	200.00	53	1	3
1P7 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

Abbreviations: ED, exposure day of factor VII at inhibitor development; IN, inhibitor negative; IP, inhibitor positive.

All patients are ethnic Chinese.

This patient number is listed as IP12 in Figure 1B.

The factor VIII:C level before inhibitor formation, for patients with inhibitor.

Patients were unable to recall the specific exposure days.

Inhibitor titer at the time of blood sample taken.

Age at the time of inhibitor formation.

Age at the time of blood sample taken.

FIGURE 1

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

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.

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 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 S2–S4).

GO and KEGG enrichment analyses were performed. The 20 most overrepresented GO terms for the upregulated differentially expressed mRNAs are shown in Figure 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 1D). The upregulated and downregulated genes involved in the immune response are shown in Tables 2 and 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 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 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 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

Immune polymorphism of IL‐1ß gene has been reported to be associated with inhibitor development.

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

The results of the KEGG enrichment analysis are shown in Figure 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 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 S1). The results were consistent with those of the RNA‐seq findings.

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 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 2A).

FIGURE 2

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 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 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.

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 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 3B. Nine of these terms were related to T‐cell or lymphocyte differentiation or activation.

FIGURE 3

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 3C,D).

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. , 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. 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,” 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. , 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. 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 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. 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. Treg cells are associated with the induction of FVIII tolerance in patients with hemophilia. SET domain bifurcated histone lysine methyltransferase 2 suppresses the expression of many genes with proinflammatory functions. Plant homeodomain finger protein 8 silences genes involved in inflammatory responses. Deficiency in Nei‐like DNA glycosylase 3 is associated with increased lymphocyte apoptosis, autoantibody production, and autoimmunity. Adenosine monophosphate deaminase 2 is likely to be involved in the proinflammatory pathway by increasing adenosine triphosphate levels. Stratifin plays a role in balancing the host inflammatory response to the virus. 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.

AppendixS1

Click here for additional data file.