Ninety-five Chinese Han ethnicity adult outpatients with active CSU diagnosed according to the criteria of the EAACI/GA²LEN/EDF/WAO guideline were enrolled from June 2014 to April 2018 (8). Recurrent CSU was defined as occurrence of CSU after absence of symptoms for at least 1 year without administering any medications. Patients with chronic inducible urticaria, urticarial vasculitis, or CSU patients who had been treated by systemic corticosteroids and/or immunosuppressants were excluded. CSU patients who simultaneously suffered from active autoimmune diseases or allergic diseases were excluded. Demographic data, medical histories, and clinical data were obtained using a questionnaire on presentation. Physical examinations and urticaria activity score of 1 day (UASday) were performed and evaluated by a dermatologist and a postgraduate. Ninety-five ethnicity-, age-, and sex-matched healthy controls with serum total IgE (tIgE) <100 IU/ml, negative specific IgEs (sIgEs), negative anti-thyroid peroxidase (TPO) IgG autoantibodies (AAbs), and negative anti-thyroglobulin (TG) IgG AAbs were recruited. This study was approved by the institutional ethics committee of The First Hospital of China Medical University. Each participant signed a written informed consent.

We used a therapeutic regimen as previously described (9). Second-generation H₁-antihistamines (sgAHs) in licensed doses were used as the first-step treatment. An increase in the dose of the sgAHs to a maximum of fourfold dose or combinations of sgAHs up to a fourfold equivalent dose comprised the second-step treatments. H₂-antihistamines (ranitidine or famotidine) and/or leukotriene receptor antagonist (montelukast) was added as the third-step treatments. Before complete control was achieved, or during the dose-increment stage, patients were assessed weekly and UASday was recorded at each visit. CSU cases failing to be completely controlled by fourfold doses or fourfold equivalent doses of sgAHs for at least 2 weeks were designated as sgAHs-refractory cases.

Serum levels of tIgE (Euroimmun, Lübeck, Germany) were measured using ELISA kits according to the manufacturer’s protocols (normal range: 0–100 IU/ml). Serum sIgEs were detected using immunoblotting (Euroimmun) as previously described (10). The detected sIgEs included 8 inhalant sIgEs and 10 food sIgEs. A test result of < 0.70 kU/L was considered negative. Serum levels of anti-TPO (normal range: 0.00–5.61 IU/ml) and anti-TG (normal range: 0.00–4.11 IU/ml) IgG AAbs were measured according to the manufacturer’s protocols (Abbott Park, Middletown, USA).

Genomic DNA was isolated from whole blood using QIAamp DNA MiniKit (Qiagen, Hilden, Germany). The purity and concentration of DNA was estimated using Nanodrop 2000 (Thermo)/Quibt 3.0. Then, 500 ng DNA of each sample was used to bisulfite converted using EZ DNA Methylation Kits (Zymo Research, USA), and the converted products were put into Illumina Infinium Human Methylation 850K BeadChip (Illumina Inc, CA, USA) in accordance with the manufacturer’s guidelines and protocol.

The array data of Illumina methylation chip (.IDAT files) were analyzed using ChAMP package in the R software for deriving the methylation level. Firstly, CpG probe filtering of raw data was performed by removing (1) CpG probes with a detection p ≥ 0.01 in 1% of samples; (2) CpG probes with bead count < 3 in 5% of samples; (3) non-CpG probes; (4) CpG probes with single-nucleotide polymorphisms (SNPs) relation; and (5) CpG probes that aligned to multiple locations. Secondly, beta-mixture quantile normalization adjustments for correcting type I and type II probe design bias were used to standardize the methylation data (11). Thirdly, we used singular value decomposition (SVD) analysis to analyze the batch effect caused by BeadChip slide and array, then applied Combat to correct this batch effect (12). SVD analysis was also used to analyze the effect of confounders including age, sex and smoking status on the DNA methylation of the samples.

As genomic DNA was extracted from clustering leukocytes, we estimated leukocyte compositions including CD4⁺ T lymphocytes, CD8⁺ T lymphocytes, natural killer (NK) cells, B lymphocytes, monocytes, and granulocytes, by the RefbaseEWAS method for inferring changes in the distribution of leukocytes between the CSU patients and the healthy controls using DNA methylation signatures, in combination with a previously obtained external validation set consisting of signatures from purified leukocyte samples (13).

Differential methylated CpG positions were calculated by differentially methylated position (DMP) analysis of ChAMP, and the corrected p values were computed using the Benjamini–Hochberg method (14). The methylation level of each CpG probe was denoted as a β value. The Δβ of each CpG site represents the difference in average β values of the CSU patients and the healthy controls. A CpG site with |Δβ| ≥ 0.06 and p < 0.01 was considered as a DMP or a differentially methylated CpG site. A CpG site was considered hypermethylated if Δβ ≥ 0.06 or hypomethylated when Δβ ≤ −0.06.

Since there are many different cell types in the whole blood, the difference caused by the change of a single cell type is not significant in the whole blood. Therefore, the cutoff value for β of the whole blood sample is generally low. In previous literature, 0.1 was commonly used as a cutoff value for β in peripheral blood samples. In this study, we used 0.06 as the cutoff value for β (simultaneously p < 0.01) for bioinformatic analysis and 0.1 as the cutoff value for β for single CpG site analysis. Moreover, 0.06 was used as the cutoff value for β to avoid the chip-based error of 0.05 in this study.

To explore the potential biological functions of target genes of DMPs or differentially methylated CpG sites, we performed gene ontology enrichment analysis including molecular function, cellular component, and biological process using the Annotation, Visualization, and Integrated Discovery (DAVID) database as well as pathway enrichment analysis using the Kyoto Encyclopedia of Genes and Genomes (KEGG) database. A p value < 0.05 was used as the threshold to select the significantly enriched gene ontology terms and pathways.

Data were processed using GraphPad Prism 8.0 and SPSS 23.0. Normally distributed continuous variables were compared by using Student’s t-test while Mann–Whitney U-test was used for comparison of non-normally distributed continuous variables. The associations between categorical variables were analyzed using Pearson’s chi-squared test. p < 0.05 was considered statistically significant. We use p-value in the test of the entire article rather than corrected p-value.

Demographic, clinical, and laboratory data of the 95 CSU patients and the 95 healthy controls are summarized in Supplementary Table 2 . The 95 CSU patients achieved complete control using different treatment regimens. Among the 95 patients, 36 were refractory. The therapeutic regimens and characteristics of the 36 refractory cases and 59 non-refractory cases are listed in Supplementary Tables 3–5 .

A total of 832,347 CpG sites were analyzed following quality filtering and data normalization. The raw data will be available at https://www.biosino.org/node with accession number OEP002482 upon approval by Ministry of Science and Technology, China. Density plot showed relatively similar methylation levels in the samples from the CSU patients and the healthy controls ( Figure 1A ). Principal component analysis (PCA) based on all 832,347 CpG sites did not reveal any discernable separation of the CSU patients from the healthy controls ( Figure 1B ). Leukocyte composition analysis showed that the percentages of CD8⁺ T lymphocytes, NK cells, and monocytes were significantly lower in the CSU patients than those in the healthy controls, whereas the percentage of granulocytes was significantly higher in the CSU patients than that in the healthy controls ( Table 1 ).

The methylation levels of the CpG sites in the CSU patients and the healthy controls were strongly correlated (R ² = 0.999, p < 0.0001, Figure 2A ). Compared with the 95 healthy controls, 439 DMPs (p < 0.01 and |Δβ| ≥ 0.06) were identified in the 95 CSU patients. Among the 439 DMPs, 380 (86.56%) were hypomethylated (p < 0.01 and Δβ ≤ −0.06) and 59 (13.44%) were hypermethylated (p < 0.01 and Δβ ≥ 0.06) ( Supplementary Table 6 ). Volcano plot of the differential DNA methylation analysis showed obviously that the number of hypomethylated DMPs were far more than that of hypermethylated DMPs ( Figure 2B ). Although there was some overlap between the CSU patients and the healthy controls, the PCA based on the 439 DMPs showed two different clusters between the CSU patients and the healthy controls ( Figure 1C ).

The average global DNA methylation levels of the 832,347 CpG sites between the 95 CSU patients and the 95 healthy controls showed no significant difference ( Figure 3A ). The average global DNA methylation levels of the 439 DMPs in the 95 CSU patients were significantly lower than those in the 95 healthy controls (p < 0.0001) ( Figure 3B ). The average methylation levels of the 59 hypermethylated DMPs in the 95 CSU patients were significantly higher than those in the 95 healthy controls (p = 0.02) and the average methylation levels of the 380 hypomethylated DMPs in the 95 CSU patients were significantly lower than those in the 95 healthy controls (p < 0.0001) ( Figures 3C, D ).

We used SVD analysis to analyze the effect of confounders including age, sex, and smoking status on the DNA methylation of the samples. From the statistical p values corresponding to the color blocks ( Supplementary Figure 1 ), age, sex and smoking status (28 CSU patients and 4 healthy controls were active smokers) showed significant correlation with DNA methylation of the samples. Then, we searched in PubMed for previously reported age-, sex-, and smoking-related CpGs and used the online tool of Venny 2.1 assisted by manual checking to identify overlapped DMPs. No overlapped DMP was found between the 439 DMPs discovered in our study and previously reported age-related CpGs (15) or smoking-related CpGs (16). In the 439 DMPs, 5 CpGs (cg04732279, cg12873598, cg16158408, cg11887420, and cg18105467) were located on chromosome X. However, none of the 5 CpGs was overlapped with previously reported sex-related CpGs (17). Therefore, the effect of confounders including age, sex, and smoking status on the DNA methylation of the 439 DMPs discovered in our study was excluded.

The 439 DMPs in the 95 CSU patients were distributed widely on chromosome 1 to 22 and chromosome X ( Figure 4A ). Among the 23 chromosomes, chromosome 6 contained the largest number of DMPs (n = 51), among which 15 DMPs-annotated genes TRIM38, MAP3K5, ARG1, C6orf106, SRPK1, MAPK14, COL11A2, HLA-DRB1, HLA-DPB2, HLA-DQB1, HLA-C, C6orf10, BAT3, HLA-DPB1, and COL21A1 were definitely involved in autoimmunity in the literature. Moreover, the comparisons between the distribution intensity of the 439 DMPs on each chromosome and the corresponding distribution intensity of the 832,347 CpG sites on each chromosome showed that only the ratio of chromosome 6 had a significant elevation ( Table 2 ). With regard to corresponding gene regions ( Figure 4B ), 79 (16.70%) DMPs were located in proximal promoter regions including TSS1550 (200–1,500 bp upstream of the transcription starting site), TSS200 (200 bp upstream of the transcription starting site), 5’UTR, and 1st exon. We observed significant differences in the distribution percentages of the 439 DMPs on the TSS1500, the TSS200, the 1st Exon, the gene body, and the 3’ UTR as compared with the distribution percentages of the 832,347 CpG sites on each gene region, respectively. While the distribution percentages of the 439 DMPs and the 832,347 CpG sites on the 5’UTR and the intergenic region were not significantly different ( Supplementary Table 7 ). Concerning CpG island regions in relation to the 439 DMPs ( Figure 4C ), 350 (79.73%) DMPs were located in the Open Sea (DNA outside the CpG island regions), 15 were located in the CpG island, 43 were located on the CpG shores (0–2 kb upstream and downstream of the CpG island, also called North and South Shore), and 31 were located on the Shelves (2–4 kb upstream and downstream of the CpG island, also called North and South Shelf). Significant differences were observed between the distribution percentages of the 439 DMPs and the 832,347 CpG sites on the island, S shore, and the open sea ( Supplementary Table 8 ).

According to the Illumina annotated CpGs and transcript pairs, all CpGs in the Open Sea were excluded from the 439 DMPs. Finally, 313 DMPs including 273 hypomethylated DMPs and 40 hypermethylated DMPs were annotated to genes. As some CpG sites were annotated to more than one gene, 304 DMGs were identified after the duplicate genes were removed. In the 18 DMGs that mapped from more than one CpG site, 15 DMGs were hypomethylated, 2 DMGs were hypermethylated, and only 1 DMG showed mixed methylation ( Supplementary Table 9 ). Of the 18 DMGs, when Δβ ≥ 0.1, HLA-DPB2, HLA-DRB1, PPP2R5C, and LTF were identified.

Across the whole genome, the CpG site cg07052231 located on the PEX5 gene had the highest association with CSU (p = 5.16E−13, Δβ = −0.08). Compared with the 95 healthy controls, 41 (22 hypomethylated, 19 hypermethylated) DMPs annotated to 28 DMGs (p < 0.01 and |Δβ| ≥ 0.1) were identified in the 95 CSU patients ( Table 3 ). The heatmap of the 41 DMPs between the 95 CSU patients and the 95 healthy controls (p < 0.01 and |Δβ| ≥ 0.1) is shown in Figure 5 . The 41 DMP-level differential methylation analysis results (DMP specific violin plots) are shown in Supplementary Figure 2 .

We divided the 95 CSU patients into eight pairs of subgroups according to eight aspects: elevated/normal tIgE, positive/negative anti-TPO IgG, positive/negative anti-TG IgG, with/without angioedema, refractory/non-refractory, UASday > 4/UASday ≤ 4, disease duration > 6 months/disease duration ≤ 6 months, and recurrent/non-recurrent. Except for the refractory/non-refractory subgroup and the disease duration > 6 months/disease duration ≤ 6 months subgroup, phenotype-specific DMPs were generated in the 439 DMPs in the other six subgroups ( Table 4 ).

Biological process analysis of the 304 DMGs showed that these genes were mainly enriched in regulation of cell signaling, cellular and intracellular compounds’ catabolic and metabolic processes, and cellular response to and regulation of endogenous and exogenous stimulus. Cellular component analysis of the 304 DMGs showed that these genes were significantly enriched in the intracellular part including cytoplasm and organelles. Molecular function analysis of the 304 DMGs showed that these genes were significantly enriched in catalytic activity, protein binding, pyrophosphatase, and hydrolase activity. Gene ontology enrichment analyses in above three aspects were displayed with the top 10 significant p values, respectively ( Supplementary Table 10 ).

The top 30 prominently enriched KEGG pathways of the 304 DMGs are shown in Table 5 . The DMGs were enriched in the pathways related to autoimmune diseases including IBD, ATD, and TID, as well as important immune-related pathways including antigen processing and presentation, Th17 cell differentiation, Th1 and Th2 cell differentiation, PI3K-Akt signaling pathway, TNF signaling pathway, and AMPK signaling pathway, which play important roles in autoimmune diseases. Among these pathways, the most significant one was the sphingolipid signaling pathway (p = 9.08E−05). The pathway with the most number of DMGs was the PI3K-Akt signaling pathway (13 DMGs, p = 3.29E−03).