Genome-wide DNA methylation profiles analysis in primary warm autoimmune hemolytic anemia patients

ABSTRACT

Background

Autoimmune hemolytic anemia (AIHA) is caused by auto-antibodies, secreted by overactivated B cells, directed against self-red blood cells, resulting in hemolysis. It found that aberrant DNA methylation in B cells can induce the production of autoantibodies. Therefore, we attempted to explore if similar aberrant DNA methylation occur in AIHA patients.

Methods

A 49-year-old female wAIHA patient and a 47-year-old female healthy control (HC) were enrolled. Peripheral blood (PB) B cells DNA was extracted. After constructing genomic libraries, bisulfite genomic sequencing (BSP) and DNA methylation profiles were analyzed. BSP was verified using PB B cells from 10 patients with hemolysis, 10 patients with hemolytic remission, and 10 healthy controls (HCs) by Methylation-specific PCR.

Results

Total DNA methylation of whole-genome C bases (4.8%) and CG type bases (76.8%) in wAIHA patient were lower than those in the HC (5.3 and 82.5%, respectively) (p = 0.022 and p < 0.001). DNA methylation of C bases and CG type bases in whole-genome regulatory elements, such as coding sequence, up2Kb and down2Kb in the patient were also lower than those in the HC (p = 0.041, p = 0.038, and p = 0.029). 30,180 DNA-methylated regions (DMRs) on all 23 chromosomes were identified. DMR-related genes were mainly involved in the Rap1, phospholipase D, HIF-1, calcium, vascular endothelial growth factor (VEGF) and Ras signaling pathways.

Conclusion

The DNA methylation spectrum of B cells in AIHA patients is different from that of HC, and the proportion of hypo-methylation regions is higher than that of HC. DMR-related genes are mainly related to some signaling pathways.

KEYWORDS:

• Warm autoimmune hemolytic anemia
• B lymphocytes
• DNA methylation
• Methylation-specific PCR
• Bisulfite genomic sequence
• Autoimmunity disease
• DMR-related genes
• Primary

1. Introduction

AIHA is an autoimmune disease (AID) characterized by autoantibodies against RBC surface antigens, which are generated and secreted by over-reactive B cells. These autoantibodies, accompanied by abnormal activation of the complement system in some patients, lead to the destruction of RBCs, causing hemolytic anemia [1]. Based on the binding properties of autoantibodies, AIHA can be further divided into warm antibody AIHA (wAIHA), cold antibody AIHA, and mixed-type AIHA; wAIHA accounts for approximately 70% of AIHA cases [2]. Although the use of glucocorticoids is effective in most patients, only approximately 30% of achieve complete remission, with hemolysis relapse in most [3]. B cells play a pivotal role in the pathogenesis of autoimmune diseases and altered B-cell signaling and hyperactivated humoral responses are key drivers of human autoimmunity [4]. In AIHA patients absolute B cell numbers were significantly increased with upregulation of activation markers [5]. Duan N found that bruton’s tyrosine kinase (BTK) and phosphorylated BTK (p-BTK) increases in the B lymphocytes of AIHA patients [6]. Therefore, it is necessary to understand the mechanisms that lead to B-cell activation in patients with wAIHA to identify potential therapeutic targets.

The term ‘epigenetics’ can be outlined as the meiotically/mitotically heritable alterations in gene expression without changes to the sequence of bases in the DNA [7]. Epigenetic dysregulation directly influences the development of autoimmunity by regulating immune cell functions [8]. Tanaka S provided clear evidence that an epigenetic regulator has a critical role in peripheral B cell tolerance and prevention of autoimmune disease [9]. Many studies have highlighted the essential roles of epigenetic regulation in different B-cell subsets and this epigenetic regulation contributes to the pathogenesis of autoimmune diseases, including systemic lupus erythematosus (SLE), primary Sjögren’s syndrome (pSS), and rheumatoid arthritis (RA) [10–12]. DNA methylation is one of the most studied epigenetic marks and methylation of the 5’ position of cytosines (5mC), in the context of CpG dinucleotides, is the main form of DNA methylation in mammals and is usually associated with gene suppression [13]. Abnormal DNA methylation in B cells can induce the production of the autoantibodies and mediate the occurrence of AID [14]. Recent findings underlined the critical role of DNA methylation in numerous autoimmune diseases by altering gene expression profiles [15,16]. We hypothesized that DNA methylation in peripheral blood B cells, which participate in the pathogenesis of AIHA, is decreased. We constructed and analyzed a wAIHA DNA methylation profile to further explore the mechanisms that lead to B cell overreaction and autoantibody secretion, in a patient with wAIHA who suffered a hemolytic episode, with an age- and sex-matched healthy control (HC). The peripheral B lymphocytes of 20 patients with wAIHA were used for verification.

2. Materials and methods

2.1. Study subjects

A 49-year-old female with wAIHA hospitalized in the Hematology Department in December 2020 and a 47-year-old female HC were enrolled for bisulfite genomic sequencing (BSP). 20 patients, including 10 with hemolysis and 10 patients in remission, and 10 HCs were enrolled for methylation-specific PCR (MSP) to verify the sequencing results of patients with hemolytic episodes and healthy control. All AIHA patients were diagnosed in accordance with Chinese guideline criteria [17]. The clinical characteristics of the 20 patients are shown in Table 1.

Table 1. Clinical characteristics of wAIHA patients.

Characteristics	Patients (N = 20)
Median age, years	63 (18∼80)
Gender
Female, N (%)	16 (80%)
Male, N (%)	4 (20%)
Diagnosis	wAIHA
RBC (×10^12/L), median (range)	2.015 (0.99∼5.82)
Hb (g/L), median (range)	79.5 (51∼166)
Ret (×10^12/L), median (range)	232.75 (34.9∼563.9)
Ret (%), median (range)	10.315 (1.77∼40.25)
WBC (×10^9/L), median (range)	6.075 (2.94∼14.36)
PLT (×10^9/L), median (range)	176 (86∼322)
TBIL (umoL/L), median (range)	26.15 (8.3∼86.4)
IBIL (umoL/L), median (range)	17.3 (6.1∼68.1)
LDH (U/L), median (range)	302 (149∼1162)
Hp (g/L), median (range)	0.1875 (0.06∼3)
FHb (mg/L), median (range)	49 (25∼239)
CD5 ^+CD19^+/CD19^+ (%), median (range)	15.46 (2.86∼44.03)
Therapies
Steroids, N (%)	20 (100%)
IVIGs, N (%)	3 (15%)
Rituximab, N (%)	2 (10%)
Other Immunosuppressors, N (%)	16 (80%)
Complications
Thrombosis, N (%)	3 (15%)
Infection, N (%)	7 (35%)

wAIHA: warm autoimmune hemolytic anemia; RBC: red blood cell; Hb: hemoglobin; Ret (%): reticulocytes (%); WBC: white blood cell; PLT: platelet; TBIL: total bilirubin; IBIL: indirect bilirubin; LDH: lactic dehydrogenase; Hp: haptoglobin; FHb: free hemoglobin; IVIG: intravenous immunoglobin.

2.2. Separation and purification of peripheral blood B lymphocytes

Peripheral venous blood (10 mL) was obtained from each subject, EDTA was added, and peripheral blood mononuclear cells were isolated using a Lymphocyte Separation Medium (Beijing Solarbio Science & Technology Co., Ltd., China). Anti-human CD19 magnetic beads (20 µl) and 80 µl separation buffer (Miltenyi Biotec, Germany) were added to 10⁷ mononuclear cells. The suspension was fully mixed in the dark and incubated at 4°C for 15 min. Separation buffer (1 mL) was added, followed by pelleting the cells at 1500 rpm for 10 min. After discarding the supernatant, 500 µl separation buffer was added to resuspend the PBMCs. A magnetic cell sorting device was installed, and CD19⁺ B cells were separated according to the manufacturer’s instructions. Sorting purity was assayed by flow cytometry and was above 90% (Figure 1).

Figure 1. The purity of CD19⁺ B lymphocytes separated by magnet-activated cell sorting (MACS) and the purity of CD19 positive B cells is greater than 95%.

2.3. DNA extraction, genomic library construction, and BSP

B-lymphocyte DNA was extracted according to the instructions of the DNA Isolation Mini Kit (Vazyme Biotech, China). DNA concentration and purity were measured by spectrophotometry. Genomic DNA samples were sheared into 100∼300 bp fragments using ultrasound [18]. Their ends were repaired and connected with a methylation sequencing adapter. Fragments were modified with a DNA methylation modification kit (EZ DNA Methylation-Gold Kit; Zymo Research, CA, USA). After 2% agarose gel electrophoresis, fragments were selected based on their lengths. Appropriate DNA fragments were purified using a QIAquick Gel Extraction Kit (Qiagen, Germany), and PCR was used to amplify them.

2.4. BSP data analysis

Clean data were obtained after filtering the sequencing data, including decontamination and low-quality data removal. The clean data were mapped to the human reference genome using BSMAP software. Mapping and bisulfite conversion rates of each sample were calculated. The methylation level of site C in the whole genome was calculated using the formula: $\mathrm{R}{\mathrm{m}}_{\mathrm{a}\mathrm{v}\mathrm{e}\mathrm{r}\mathrm{a}\mathrm{g}\mathrm{e}}=[\mathrm{N}{\mathrm{m}}_{\mathrm{a}\mathrm{l}\mathrm{l}}/(\mathrm{N}{\mathrm{m}}_{\mathrm{a}\mathrm{l}\mathrm{l}}+\mathrm{N}\mathrm{n}{\mathrm{m}}_{\mathrm{a}\mathrm{l}\mathrm{l}}\times 100\left)\right]$where Nm represents the number of reads indicating methylated sites, while Nnm is the number of reads indicating non-methylated sites.

DMRs in B cells of the patient with AIHA were identified based on methylation levels of C sites in the whole genome compared with the HC. Each DMR contained at least five CG sites and the average methylation levels in this region were compared between the two samples. If the fold change of the average DNA methylation level in this region was ≥2 (p ≤ 0.05) by Fisher's exact test compared with the HC, then this region was to be a DMR. Gene ontology (GO) functional enrichment and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses were performed on DMR-related genes [19].

2.5. Validation of BSP data by MSP

Five DMRs identified by BSP were selected randomly, including two genes in which the DNA methylation of the promoter was significantly decreased (TCF4 and IL-18) and three genes with significantly elevated promoter DNA methylation (PTPN7, NOD2, and WT1). Based on the sequencing data, MSP primers were designed using online primer design tools (http://www.urogene.org/methprimer/index.html) (Table 2). An additional 10 patients with hemolytic wAIHA, 10 with wAIHA in remission, and 10 HCs were enrolled. B-lymphocytes from each group were isolated and purified as described above. DNA was extracted and modified using an EpiArt^TM DNA Methylation Bisulfite Kit (Vazyme Biotech). After PCR amplification and agarose gel electrophoresis of the bisulfite transformation DNA products, MSP data were obtained using a gel imaging analysis system (Guangzhou Biolight Biotechnology Co., Ltd., China).

Table 2. The primer sequences for methylation specific PCR.

	log2 (fold change of the DNA methylation level)	Methylation primer sequences	Unmethylation primer sequences
NOD2	−2.017	F-TGGGTTTTAGGGAGGGTTTAC R-AACTTCCCCTCTTCAAACGA	F-GGTTTTAGGGAGGGTTTATGT R-AACTTCCCCTCTTCAAACAAA
WT1	−2.263	F-TTGGGTAAGTATATATGAAGGGGC R-TCTCTACCTACAAAATATACGACGTA	F-TTTGGGTAAGTATATATGAAGGGGT R-TCTCTACCTACAAAATATACAACATA
PTPN7	−2.888	F-TATTTTTTACGTAGTTTACGTTCGG R-ACAAAAAAATCATTAAAAACTCACG	F-ATTTATTTTTTATGTAGTTTATGTTTGG R-AAAAAAATCATTAAAAACTCACACA
IL-18	3.322	F-GTAGATTGTAGTAGGTGGTAGTCGT R-AAAAAAAACCTATTTAAAAACGAA	F-GTAGATTGTAGTAGGTGGTAGTTGT R-AAAAAAAACCTATTTAAAAACAAA
TCF4	3.369	F-GGAATTAGAGAGGTGAGTAGGAATC R-CTAAATATTAAATATTTTCGCCGAA	F-AATTAGAGAGGTGAGTAGGAATTGG R-ACTAAATATTAAATATTTTCACCAAA

2.6. Statistical analysis

All analyses were performed using SPSS 22.0. Data are expressed as means ± standard deviation. Comparisons between two groups were analyzed using Student’s t-test. One-way analysis of variance was used for multiple independent groups. Fisher's exact test was used to analyze the data. Statistical significance was set at p < 0.05. GraphPad Prism 6 software was used to visually depict the statistical analyses.

3. Results

3.1. Sequencing library construction

Sequencing data were more than 90% clean, which met the requirements for further analysis (Table 3). Figure 2A shows the base content distribution of the clean reads. The left part represents strand 1, which was directly obtained during sequencing. Because of bisulfite modification, the unmethylated cytosine (C) was converted to uracil (U) and then to thymine (T). The T content was relatively high, whereas C content was relatively low in strand 1. The right part shows the base content distribution of complementary strand 2. The adenine (A) content was relatively high, and the guanine (G) content was relatively low, which corresponded to strand 1. Figure 2B shows the sequencing quality distribution of bases in the clean reads. During sequencing, the quality of most of the bases was greater than 35 (Q > 35), while the proportion of low-quality (Q < 20) bases was extremely low. Clean read base content distribution and base sequencing quality distribution indicated the reliability of the BSP data. The base content distribution and base mass distribution maps initially verified the reliability of the sequencing results. In summary, the data used for database construction meets the standards and can be used for methylation analysis.

Figure 2. (A) After getting the original data, first filter the data, including decontamination and low-quality data, to get ‘clean read’. The base content distribution of clean reads. The horizontal axis indicates the base position on the reads and vertical axis indicates the proportion of the bases. The left part represents strand 1 which was obtained directly during the sequencing process. The content of T bases was relatively high and the content of C bases was relatively low. The right part represents the base contents distribution of complementary strand 2, and the content of adenine (A) was relatively high, the content of guanine (G) was relatively low which was corresponding to strand 1. (B) The base quality distribution of clean reads. The horizontal axis indicates the base position on the reads and vertical axis indicates the base sequencing quality. Each point in the figure represents the sequencing quality of the corresponding base in reads. The figure indicates that the distribution of base content and quality distribution meet the standards for further analysis.

Table 3. The basic processing and quality control of the sequencing data.

Sample ID	Fragment length	Sequencing strategy	Clean Reads number	Clean data size (bp)	Clean rate (%)
AIHA patient	100∼500 bp	PE100	757,483,644	113,622,546,600	91.75
Healthy Control	100∼500 bp	PE100	805,140,360	120,771,054,000	90.98

3.2. Mapping, depth, and coverage of C sites

After clean read acquisition, BSMAP was used to map the clean reads to the human genome. Mapping rates were greater than 95%. The bisulfite conversion rate was approximately 99%. The sequencing depth of C-sites for each sample was more than 18X and coverage were both greater than 85% (Table 4).

Table 4. Mapping results of the sequencing data and QC control.

Sample ID	Clean reads Q20 Rate (%)	Clean data size (bp)	Mapping rate (%)	Bisulfate conversion rate (%)	Duplication rate (%)	Average depth (X)	Coverage (%)
AIHA patient	93.05	757,483,644	95.82	99.06	5.77	18.69	85.753
Healthy Control	93.07	805,140,360	95.38	98.99	4.01	20.18	85.809

Sequencing depth was approximately evenly distributed (Figure 3). The depth cumulative distribution curves for different types of C bases in the patient with AIHA and the HC roughly coincided. The proportion of C bases with sequencing depths of no less than 20X accounted for approximately 80% of all C bases in the whole genome. The coverages of C sites within the range of all chromosomes and regulatory elements in whole genomes of the patient with AIHA and the HC were approximately 80%. These data met the criteria for further analysis (Figure 4).

Figure 3. The sequencing depth distribution of (A) autoimmune hemolytic anemia (AIHA) patient and (B) healthy control. The horizontal axis represents the sequencing depth and the vertical axis is the percentage of the whole genome bases corresponding to the sequencing depth. The sequencing depth corresponding to the highest point of the curve is close to or consistent with the average coverage depth of the entire genome, and the distribution map shows that the sequencing data results are roughly evenly distributed.

Figure 4. The cumulative sequencing depth distribution of C bases in (A) autoimmune hemolytic anemia (AIHA) patient and (B) Health Control. The horizontal axis represents the sequencing depth and the vertical axis represent the percentage of covered C bases which sequencing depth is no less than the corresponding sequencing depth. The proportion of C bases in the entire genome is approximately 80%, and the cumulative distribution curves of sequencing depth for different types of C bases roughly overlap.

The above suggests that the coverage of C sites within the range of chromosomes and various types of regulatory elements in the entire genome meets the analysis standards.

3.3. Whole-genome DNA methylation in AIHA

Methylation level of the C base was calculated as the number of sequences covering this site, which supported methylated C sites divided by the total number of sequences effectively covering this site during sequencing. The total DNA methylation levels of whole-genome C bases (4.848%) and CG type bases (76.793%) in the patient were lower than those in the HC (5.335% and 82.504%) (p = 0.022 and p < 0.001). DNA methylation levels of C bases and CG type bases in regions of the whole genome regulatory element, such as coding sequence, up2Kb and down2Kb in the patient with AIHA were lower than those in the HC (p = 0.041, p = 0.038, and p = 0.029) (Tables 5 and 6). This result shows the degree of methylation of the c base contained on each chromosome of AIHA patients and healthy controls.

Table 5. The average DNA methylation level of all chromosomes and regulatory elements in the whole genome of wAIHA patient.

Chr	C (%)	CG (%)	CHG (%)	CHH (%)
NC_000001.11	5.213	77.408	1.187	1.23
NC_000002.12	4.535	77.998	1.165	1.204
NC_000003.12	4.341	77.918	1.158	1.204
NC_000004.12	4.373	76.215	1.195	1.219
NC_000005.10	4.451	77.106	1.165	1.221
NC_000006.12	4.467	77.738	1.163	1.203
NC_000007.14	4.842	78.126	1.175	1.209
NC_000008.11	4.432	77.548	1.154	1.201
NC_000009.12	4.798	77.421	1.174	1.207
NC_000010.11	5.128	77.585	1.191	1.221
NC_000011.10	4.542	76.223	1.165	1.204
NC_000012.12	4.713	77.722	1.166	1.203
NC_000013.11	4.648	77.184	1.162	1.219
NC_000014.9	4.608	77.392	1.164	1.199
NC_000015.10	4.884	77.555	1.177	1.205
NC_000016.10	7.37	77.603	1.212	1.312
NC_000017.11	5.947	76.34	1.216	1.215
NC_000018.10	4.71	78.006	1.178	1.211
NC_000019.10	6.722	75.352	1.239	1.205
NC_000020.11	5.22	75.929	1.19	1.22
NC_000021.9	5.145	76.297	1.187	1.221
NC_000022.11	6.277	77.285	1.232	1.214
NC_000023.11	4.028	72.413	1.138	1.189
Whole Genome	4.848	76.793	1.178	1.216
Different Genome Regions	C (%)	CG (%)	CHG (%)	CHH (%)
CDS	8.656	64.778	1.258	1.216
Down2k	5.514	75.396	1.2	1.199
Up2k	5.32	48.834	1.191	1.186
mRNA	5.026	78.349	1.175	1.2
repeat	5.336	82.45	1.208	1.22
CpG Island	13.671	44.803	1.577	1.261

Table 6. The average DNA methylation level of all chromosomes and regulatory elements in the whole genome of HCs.

Chr	C (%)	CG (%)	CHG (%)	CHH (%)
NC_000001.11	5.739	82.967	1.275	1.322
NC_000002.12	5.001	83.766	1.25	1.294
NC_000003.12	4.785	83.86	1.244	1.295
NC_000004.12	4.841	82.171	1.286	1.314
NC_000005.10	4.901	83.066	1.252	1.315
NC_000006.12	4.913	83.419	1.245	1.293
NC_000007.14	5.333	83.785	1.262	1.298
NC_000008.11	4.88	83.574	1.241	1.287
NC_000009.12	5.278	82.943	1.262	1.292
NC_000010.11	5.594	83.313	1.279	1.309
NC_000011.10	5.009	82.023	1.25	1.291
NC_000012.12	5.171	83.087	1.249	1.292
NC_000013.11	5.092	83.066	1.248	1.313
NC_000014.9	5.063	82.715	1.25	1.288
NC_000015.10	5.352	82.853	1.256	1.284
NC_000016.10	8.162	83.352	1.312	1.427
NC_000017.11	6.462	81.363	1.305	1.296
NC_000018.10	5.191	83.972	1.262	1.298
NC_000019.10	7.298	79.698	1.333	1.29
NC_000020.11	5.68	81.588	1.271	1.297
NC_000021.9	5.689	82.615	1.281	1.31
NC_000022.11	6.821	82.211	1.317	1.287
NC_000023.11	4.497	79.446	1.221	1.279
Whole Genome	5.335	82.504	1.265	1.306
Different Genome Regions	C (%)	CG (%)	CHG (%)	CHH (%)
CDS	9.351	66.688	1.36	1.287
Down2k	6.036	79.7	1.289	1.279
Up2k	5.8	50.241	1.276	1.267
mRNA	5.519	83.215	1.262	1.287
repeat	5.823	88.262	1.296	1.309
CpG Island	13.613	44.134	1.694	1.337

3.4. GO functional enrichment and KEGG pathway enrichment analyses

Table 7 shows the number of CG-base DMRs detected and the total length of DMRs on each chromosome in the patient with AIHA compared with the HC. DMR-related genes and genes corresponding to DMR-related promoters were identified. Compared with the HC, the top 10 DMR-related genes and the top 10 genes corresponding to DMR-related promoters, in which DNA methylation levels were significantly increased or decreased in B cells of the patient with AIHA are shown (Figures 5 and 6; Tables 8 and 9). We have identified genes and promoter regions with significant differences in methylation degree between AIHA patients and healthy controls.

Figure 5. Compared with healthy control, top 10 DNA-methylated regions (DMRs)-related genes which the DNA methylation level significantly increased or decreased in B cells of AIHA patient.

Figure 6. Compared with healthy control, the top 10 genes corresponding to DNA-methylated regions (DMRs)-related promotor regions which the DNA methylation level significantly increased or decreased in B cells of AIHA patient.

Table 7. The differential methylation regions between AIHA patients and HCs.

Chr	DMRs number	DMRs length
NC_000001.11	2724	889,332
NC_000002.12	2125	706,958
NC_000003.12	1562	527,986
NC_000004.12	973	311,225
NC_000005.10	1375	441,706
NC_000006.12	1366	471,583
NC_000007.14	1543	479,588
NC_000008.11	1093	354,673
NC_000009.12	1433	452,855
NC_000010.11	1517	495,155
NC_000011.10	1493	488,337
NC_000012.12	1381	450,669
NC_000013.11	563	187,522
NC_000014.9	883	297,286
NC_000015.10	915	315,794
NC_000016.10	1493	476,881
NC_000017.11	1770	573,758
NC_000018.10	604	203,224
NC_000019.10	1749	530,270
NC_000020.11	1133	366,086
NC_000021.9	513	166,273
NC_000022.11	1076	343,951
NC_000023.11	896	257,372

Table 8. Top 10 significantly altered DMRs related genes in B cells from AIHA patient.

Gene ID	Gene name	log2Ratio (AIHA/HC)	Type of DNA methylation level differences
114569	MAL2	−20	Reduced
114821	ZBED9	−20
3054	HCFC1	−20
54880	BCOR	−20
9071	CLDN10	−20
9249	DHRS3	−20
221	ALDH3B1	−5
26166	RGS22	−4.755
1824	DSC2	−4.459
728780	ANKDD1B	−4.459
1822	ATN1	3.524	Elevated
118	ADD1	3.544
6886	TAL1	3.585
91050	CCDC149	3.585
753	LDLRAD4	3.954
107987075	LOC107987075	4.17
55786	ZNF415	4.392
5937	RBMS1	4.459
223117	SEMA3D	4.644
58480	RHOU	4.644

Table 9. Top 10 significantly altered genes corresponding to DMRs related promotor regions in B cells of AIHA patient.

Gene ID	Gene name	log2Ratio (AIHA/HC)	Type of DNA methylation level differences
3054	HCFC1	−20	Reduced
8269	TMEM187	−20
9249	DHRS3	−20
221	ALDH3B1	−5
26166	RGS22	−4.755
1824	DSC2	−4.459
3831	KLC1	−3.907
3981	LIG4	−3.907
6239	RREB1	−3.624
25851	TECPR1	−3.459
55669	MFN1	3.17	Elevated
3606	IL18	3.322
6925	TCF4	3.369
219844	HYLS1	3.392
640	BLK	3.402
91050	CCDC149	3.585
55757	UGGT2	3.7
124842	TMEM132E	4.392
400591	C17orf102	4.392
574029	DUSP5P1	4.644

GO functional enrichment analysis was performed for DMR-related genes. The top 10 genes with the most significant enrichment were selected for each of the three types of gene ontology (Figure 7A). KEGG pathway enrichment analysis was performed for DMR-related genes; the top 20 pathways with the most significant enrichment are shown in Figure 7B. Through analysis, it was found that DMR-related genes are mainly involved in Rap1, phospholipase D, HIF-1, calcium, VEGF, and Ras signaling pathways, regulation of transient receptor potential ion channels in the inflammatory response, and the process of neurotransmitter transmission.

Figure 7. (A) Gene ontology (GO) functional enrichment analysis of DMRs-related genes. (B) Kyoto encyclopedia of genes and genomes (KEGG) pathway enrichment analysis of DNA-methylated regions (DMRs)-related genes. (C) GO functional enrichment analysis of genes corresponding to DNA-methylated regions (DMRs)-related promoter regions. (D) KEGG pathway enrichment analysis of genes corresponding to DNA-methylated regions (DMRs)-related promoter regions.

GO functional enrichment analysis was performed for DMR-related promoters. The top 10 genes with the most significant enrichment were selected for each of the three types of gene ontology (Figure 7C). KEGG pathway enrichment analysis was also performed, and the top 20 pathway items with the most significant enrichment are shown in Figure 7D. It found that genes corresponding to DMR-related to promoter regions are mainly involved in Rap1, cytokines, Ras, HIF-1, NF-κB, and phospholipase D signaling pathways, biosynthesis of fatty acids, metabolism of fructose, mannose, fatty acid, phosphonate, phosphinate, and glycerol-phospholipid, primary immune deficiency, and the interaction of cytokines and their receptors.

3.5. Verification of BSP using MSP

An additional 10 patients with hemolytic wAIHA, 10 with remission wAIHA, and 10 HCs were enrolled in the study. B lymphocytes from each were isolated and purified. Their DNA was extracted and modified. After PCR amplification and agarose gel electrophoresis of the bisulfite transformation DNA products, MSP data were quantified using a gel imaging system. DNA products that could be amplified by methylation primers were to be methylated in the promoter, and those that could be amplified by non-methylation primers were to be unmethylated in the promoter. DNA products that could be amplified by both methylation and non-methylation primers were to be partially methylated in the promoter.

The methylation status of the NOD2 promoter region in B cells was assessed. In one case, no target product was amplified by non-methylation or methylation-primers. Among the remaining 29 cases, there were 10 cases in the hemolytic group, 10 cases in the remission group, and 9 cases in the HC group. In the hemolytic group, methylation and partial methylation were present in six and four cases, respectively, while non-methylation was absent in the NOD2 promoter in peripheral blood B cells. In the remission group, methylation and partial methylation were present in six and four cases, respectively, while non-methylation was absent. In the HC group, methylation and partial methylation were present in eight and one case, respectively, while non-methylation absent. The percentage of cases in which the NOD2 promoter was methylated in the B cells of patients with AIHA (60%) was less than that in HCs (88.9%). However, methylation of the NOD2 promoter was not significantly different among the three groups according to Fisher’s exact test (p = 0.311) (Figure 8A and Table 10).

Figure 8. (A) Methylation status of NOD2 gene promoter region in peripheral blood B cells. The proportion of methylation in the NOD2 gene promoter region of peripheral blood B cells in AIHA patients (60%) is lower than that in healthy controls (88.9%), but there was no statistical difference. (B) Methylation status of TCF4 gene promoter region in peripheral blood B cells. The proportion of unmethylated TCF4 gene promoter region in peripheral blood B cells of AIHA patients (21.1%) is lower than that of healthy controls (40%), which also lacks statistical significance.

Table 10. Methylation status of NOD2 gene promoter region in peripheral blood B cells of AIHA patients and HCs.

	Cases of methylation (%)	Cases of partial methylation (%)	Cases of unmethylation (%)	p
Hemolytic Group	6 (60%)	4 (40%)	0 (0%)
Remission Group	6 (60%)	4 (40%)	0 (0%)
Healthy Controls	8 (88.9%)	1 (11.1%)	0 (0%)	0.311

The methylation status of TCF4 promoter in B cells was assessed. In one case, no target product was amplified by non-methylation or methylation primers. Among the remaining 29 cases, there were nine cases in the hemolytic group, 10 cases in the remission group, and 10 cases in the HCs. In the hemolytic group, methylation was absent, partial methylation was present in seven cases, and non-methylation was present in two cases for the TCF4 promoter. In the remission group, methylation was absent, while partial methylation and non-methylation were present in eight and two cases, respectively. In HCs, methylation was absent, partial methylation in six cases, and non-methylation in four cases. The percentage of cases in which the TCF4 promoter was non-methylated in peripheral blood B cells of patients with AIHA (21.1%) was lower than that in HCs (40%). However, methylation of the TCF4 promoter was not significantly different among the three groups according to Fisher’s exact test (p = 0.563) (Figure 8B and Table 11).

Table 11. Methylation status of TCF4 gene promoter region in peripheral blood B cells of AIHA patients and HCs.

	Cases of methylation (%)	Cases of partial methylation (%)	Cases of unmethylation (%)	p
Hemolytic Group	0 (0%)	7 (77.8%)	2 (22.2%)
Remission Group	0 (0%)	8 (80%)	2 (20%)
Healthy Controls	0 (0%)	6 (60%)	4 (40%)	0.563

The degree of NOD2 gene promoter methylation in peripheral blood B cells of AIHA patients was lower than that of healthy controls, while the degree of TCF4 gene promoter methylation was higher than that of healthy controls, consistent with BSP sequencing results, but lacking statistical differences.

The promoters of the WT1 and PTPN7 genes were methylated in all DNA samples analyzed by MSP. In our methylation analysis of the IL-18 promoter, we found no amplification by either non-methylation or methylation primers designed using online tools in any DNA sample. It showed that the WT1 and PTPN7 gene promoters of AIHA patients and healthy controls were all methylated, and there was no difference between the two groups, and there was no difference in the methylation of the IL-18 promoter region between the two groups.

4. Discussion

BSP is the gold standard for DNA methylation studies, as it provides accurate data for single cytosine bases (C-bases) in the whole genome [20]. After our construction of the B-cell BSP methylation profile in AIHA, more DMRs and DMR-related downstream genes were identified. The specific mechanisms by which abnormal B cell DNA methylation participates in the pathogenesis of AIHA needs to be further explored.

In this study, we found that total DNA methylation levels of both C and CG type bases in the patient with AIHA were lower than those in the HC. In the patient with AIHA, DNA methylation levels of C bases and CG type bases in regions of the whole genome regulatory element, such as coding sequence, up2Kb, down2Kb, and mRNA regions, were lower than those in the HC. These data preliminarily demonstrate that DNA methylation profiles differ between patients with AIHA and HC. A total of 30,180 DMRs were identified on 23 chromosomes of the patient with AIHA. There are two types of corresponding relationships between DMRs and their downstream genes. First, any regions identified as DMRs can be screened into the Gene List as DMR-related genes. Second, only promoters identified as DMRs can be screened into the Gene List as genes corresponding to DMR-related promoters. Studies have found that the pathways enriched by the KEGG public database for downstream genes screened by two different corresponding relationships overlap. Phospholipase D (PLD), HIF-1, Rap1, and Ras signaling pathways were enriched in two different analyses. Many studies have found that abnormal activation of these signaling pathways is correlated with the occurrence of multiple AIDs.

PLD on the surface of mammalian cells catalyzes the hydrolysis of phosphatidylcholine to produce choline and phosphatidic acid, which function as secondary messengers to transmit signals and promote the transport of vesicles [21]. Among the PLD protein family, PLD1 and PLD2 are two major subtypes that have been extensively studied. Lymphocyte adhesion and transport through the endothelial barrier are key factors that mediate the development of central nervous system immune diseases such as multiple sclerosis. PLD1 expression is elevated in the lymphocytes of patients with multiple sclerosis, where it plays a central role in immune cell migration [22]. PLD-2 signaling promotes chemotaxis, migration, and cytoskeletal remodeling of immune cells, which participate in the pathogenesis of multiple AIDs [23]. It has been reported that PLD2 mutations may be involved in familial systemic lupus erythematosus by inhibiting PLD2-RAS signaling and downregulating Ras activity [24].

Recent studies have shown that HIF-1 signaling regulates the development and differentiation of multiple immune cells. HIF-1 induces the transformation of macrophages into M1 macrophages that participate in the anti-infection immune response by producing IL-1β, inducible nitric oxide synthase, and other proinflammatory cytokines [25]. It also induces the expression of the co-stimulatory factors CD80 and CD86 on dendritic cells to promote T-cell proliferation [26]. It also inhibits Treg development by binding Foxp3 and promoting its degradation. Studies have also indicated that HIF-1 directly activates the transcription factor RORγt and glycolysis pathway to promote the development of Th17 cells that are involved in innate immunity and inflammation, resulting in the overexpression of IL-17, which mediates the early inflammatory response and histopathological damage [27]. Additionally, HIF-1 plays an important role in the proliferation and differentiation of progenitor B cells, regulating the development of B1 and B2 cell subsets. It has been shown that the expression level of CD45 (B220) on the surface of B1 cells in the peritoneal cavity of HIF1A-deficient murine Hif1a^-/-Rag2^-/- model was significantly increased, and B2 cells matured abnormally with increased autoreactivity [28].

Abnormal activation of the HIF-1 signaling pathway correlates with the occurrence of many AIDs, including rheumatoid arthritis. The expression of HIF-1 is significantly increased, and its level is significantly associated with the number of new blood vessels in the synovial tissue of patients with rheumatoid arthritis [28]. The expression of HIF-1 promotes the proliferation of Th1 and Th17 cells, which is mediated by synovial fibroblasts, resulting in the production of IFN-γ and IL-17 [29]. In addition, the inflammatory response in synovial tissues is promoted by activation of the HIF-1/IL-33 signaling pathway [30]. HIF-1 also participates in the proliferation and migration of fibroblast-like synovial cells mediated by matrix metalloproteinases [31]. HIF-1 promotes glycolysis in osteoclasts and synovial cells to produce a large amount of ATP, leading to pathological bone damage in patients with rheumatoid arthritis [32]. It has been shown that the expression of HIF-1 is elevated in the skin and serum of patients with psoriasis [33]. It promotes angiogenesis by promoting the expression of VEGF and inducing keratinocyte proliferation through glucose transporter member-1 [34]. HIF-1 directly activates the transcription factor RORγt to promote Th17 cell differentiation and IL-17 production [35]. Th17 cells promote the secretion of inflammatory cytokines such as IL-6 and IL-8 by keratinocytes, which play an important role in the maintenance of chronic inflammation in patients with psoriasis. In conclusion, HIF-1 is a potential therapeutic target for numerous AIDs.

The Rap1 signaling pathway also regulates the development and differentiation of multiple immune cells and responses [36]. Overactivation of Rap1 signaling leads to inhibition of double-positive T cell proliferation in the thymus [37]. Meanwhile, activation of Rap1 signaling induces the expression of lymphocyte function-associated antigen-1 and mediates its binding with its ligand intercellular adhesion molecule-1 to promote the formation of immune synapses between T cells and antigen-presenting cells. However, overactivation of Rap1 signaling results in T-cell anergy through inhibition of RAS-ERK signaling [38]. The Rap1 signaling pathway also regulates the differentiation and development of B cells and the development of immune tolerance [39]. After developing into immature B cells in the bone marrow, some B cells recognize self-antigens and have a tendency for autoimmunity. These immature B cells establish central immune tolerance through clonal deletion, receptor editing, and anergy ^[40]. During receptor editing, the light chain Vκ gene is expressed and rearranged, resulting in the synthesis of a new light chain to replace the self-reactive light chain and allowing further development of the B cells. Overactivation of Rap1 signaling promotes the expression of the transcription factor OcaB through downstream p38MAPK signaling, which leads to the gradual transition from light chain Vκ gene expression and rearrangement to Vλ gene expression and rearrangement, resulting in the allelic fusion of Vκ and Vλ genes [41]. B cells arising from incomplete receptor editing have a tendency to self-react and migrate to the peritoneal cavity to become B1a cells. A large number of autoantibodies are produced and secreted from B1a cells after receiving peripheral self-antigen stimulation through class switching and affinity maturation, which mediates the occurrence of AIDs [42]. Minato et al. found that leukemia occurred in some Spa1^-/- mice, characterized by overactivation of the Rap1 signaling pathway, in which the phenotype of malignant cells was identical to that of B1a cells (CD5 ⁺Mac-1 ⁺B220⁺) [43]. The production of anti-RBC autoantibodies and the onset of hemolytic anemia were also observed in a few Spa1^-/- mice, with clinical features similar to those of chronic lymphocytic leukemia in humans.

The Rap1 signaling pathway was significantly enriched in the KEGG pathway enrichment analysis for DMR-related genes. Its activation may promote the expression of transcription factor OcaB via p38MAPK signaling, induce the allelic fusion of Vκ and Vλ genes during receptor editing, and promote the proliferation of self-reactive B1a cells, which mediate the occurrence of AIHA.

In this study, we used MSP with many patients to verify our BSP data from a single patient. The methylation of the NOD2 promoter in peripheral blood B cells of patients with AIHA was lower than that in HCs. The methylation level of the TCF4 promoter region was higher than that of HCs. These data were consistent with those obtained by BSP. However, Fisher's exact test indicated that there were no significant differences in methylation of the NOD2 or TCF4 promoters among the groups, which may be related to small sample size. BSP can only roughly describe and analyze the DNA methylation profile of B cells in patients with AIHA, so our results need to be further verified by MSP with large sample sizes. However, GO functional enrichment and KEGG pathway enrichment analyses for DMR-related genes provides a direction for further research on the mechanism of B cell hyperimmunity caused by abnormal DNA methylation in AIHA.

5. Conclusions

The total DNA methylation level of whole genome C bases and CG type bases in peripheral blood B cells of patients with wAIHA was reduced, and the DNA methylation profile of B lymphocytes of the patient with AIHA was different from that of the HC, but the studies with larger samples will clarify if these differences are statistically significant. Functional enrichment and pathway enrichment analyses of DMR-related genes suggest that activation of the Rap1 signaling pathway may lead to B lymphocyte over-immunity and autoantibody production in patients with AIHA.

Ethics approval and consent to participate

All study participants provided informed consent. The study design was approved by the Ethics Committee of Tianjin Medical University General Hospital and was performed in accordance with the Declaration of Helsinki.

Acknowledgements

The authors thank Editage (www.editage.cn) for English language editing.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

This project was funded by the National Natural Science Foundation of China (81170472, 81970116, and 81770118).