RNA probes were designed against the entire HLA class II locus sequence (GRCh38 coordinates: chr6:32,426,802-34,167,129:1) retrieved from the UCSC (GENCODE version 38lift37) genome browser ( Figure 1A ) with a density of 3X probe tiling, i.e., each base of a target region was on average covered by three overlapping probes ( Supplementary Figure 1A ). The size of each probe is 120 nucleotides and tolerate mismatches and even smaller insertions/deletions. Probes were also designed against proximal promoters of 2,157 immune-related genes ( Supplementary Table 1 ) selected based on six publicly available transcriptomic data (38–42) from the monocytic/macrophage THP-1 cell line. We also utilized Hi-C data to include the interactions of the HLA with other genes using public database for THP-1 cells (43) ( Supplementary Table 2 ). Genes included in the final selection displayed transcriptional changes after inflammatory stimulation in at least three independent studies, complemented by additional genes from commercial myeloid panel genes (Nanostring) ( Supplementary Table 2 ). These non-HLA probes were designed based on promoter IDs obtained from the Eukaryotic Promoter Database (EPDnew released 2019, assembly GRCh38/hg38, https://epd.epfl.ch//index.php) with a 1-3X density. To improve the capture of regions with high GC%, we applied a so-called probe boosting ( Supplementary Figure 1 ), i.e., a dynamic process where each designed probe is replicated up to 16 times based on the GC% of a region. For the HLA sub-design masking strategy, we selected the Least Stringent option, (https://earray.chem.agilent.com/suredesign/help/Impact_of_parameters_on_probe_selection.htm). This approach allowed for the identification of repetitive elements (such as DNA transposons, LTR, non LTR-retrotransposons, LINE and SINE elements) by all three masking tools (RepeatMasker, WindowMasker, and the Duke Uniqueness 35 track) and the subsequent masking from the probe design. Of note, the pseudogenes of the HLA locus were not masked and the probe design enables their capture. The final panel design (Design ID: S3381502) comprised 120,990 probes covering 2,223 kilobase (kb), of which 65,384 probes cover 1,621 kb of the HLA class II locus, while 55,606 probes cover 602 kb of promoter sequences ( Figure 1A ; Table 1 ).

A total of 3.3 μg of genomic DNA was diluted in Low TE Buffer (10 mM Tris-HCl pH 8.0, 0.1 mM EDTA) in microTUBE-50 AFA Fiber Screw-Cap (Covaris) for fragmentation on Covaris M220. The DNA was sheared at 20°C for 160 seconds and 200 cycles per burst (Covaris M220) to achieve 250 bp-long fragmentation ( Figure 1B ). The quality of DNA was assessed using the Agilent 1000 DNA kit (Agilent Technologies) in a 2100 Bioanalyzer (Agilent Technologies).

The libraries were prepared using the SureSelect^XT Methyl-Seq kit (Agilent, part# G9651B) and beads from the Agencourt AMPure XP kit (Beckman Coulter) were used for purification of the DNA fragments, followed by end-repairing and A-tailing ( Figure 1B ). After ligation of the methylated adapters the samples were split into two reactions in a 2:1 volume ratio in order to accomodate paralelle Bisulfite (BS) and Enzymatic (EM) conversion. Two different library preparations followed: (i) BS-Seq & Variant Library and (ii) EM-Seq Library ( Figure 1C ). For the Bisulfite Conversion & Variant Library the hybridization of the samples was performed using the Custom Probe Design (Design ID: S3381502) for 16 hours at 65°C. After the hybridization, capture of the DNA-RNA formed complexes was performed using Dynabeads MyOne Streptavidin T1 beads (Thermo Fischer Scientific) under 30 min shaking at 1800 rpm on a mixer. In the last washing step, the 200 μl sample was split into two equal aliquots for variant genomic detection and BS methylation analysis methylation analysis.

For the construction of the BS Converted library, the DNA was eluted in 20 μl of 0.1 M NaOH and processed immediately using the EZ-DNA Methylation kit (Zymo Research). BS conversion lasted 150 minutes at 64°C and the BS converted DNA was eluted in a 20 μl final volume. Subsequently, the library was amplified using the SureSelect^XT Methyl reagent kit for 8 cycles. For the construction of the Variant Library, 30 μl of nuclease-free water was added to the beads, the DNA was amplified using the Herculase II Fusion DNA Polymerase (Agilent Technologies) for 8 PCR cycles and the libraries were cleaned using Agencourt AMPure XP beads.

For the construction of the EM Converted library, the DNA was oxidized using the NEB Next Enzymatic Methyl-seq Conversion Module (New England BioLabs). More specifically, the DNA was diluted in a 28 μL volume and incubated with the oxidation enzymes for 1h at 37°C in a thermocycler. Enzymes were inactivated using the Stop Reagent according to the manufacturer’s instructions and the converted DNA was purified using the Agencourt AMPure XP beads. The DNA was hybridized with the custom-made probes (Design ID: S3381502) for 16 hours at 65°C. The same hybridization and capture steps were followed as mentioned previously and the DNA was eluted in 20 μl of 0.1 M NaOH. Deamination of the cytosines was performed for 3 hours at 37°C and then the samples were cleaned with AMPure XP beads. The library was amplified for 8 cycles and indexed.

The quality and quantity of the libraries were assessed using the BioAnalyzer High Sensitivity DNA Kit (Agilent Technologies) before and after pooling in equal amounts ( Supplementary Figure 2 ) prior to sequencing with an Illumina NovaSeq sequencer and a 2x150bp protocol.

The THP-1 acute monocytic leukemic cell line (ATCC: The Global Bioresource Center) was expanded in complete RPMI medium composed of 10% Fetal Bovine Serum, supplemented with 1% Penicillin-Streptomycin (all from Merck KGaA) with the addition of L-glutamine, sodium pyruvate and 2-mercaptoethanol (1.7 mM, 0.87 mM and 0.02 mM, respectively). Cells were plated in a 6-well plate in complete RPMI medium and stimulated with 20 ng/ml of interferon-gamma (IFNγ) and/or 5 ng/ml lipopolysaccharide (LPS) for 24 hours. Following stimulation, cells were centrifuged at 350 g for 5 min, lysed in the RLT plus buffer (Qiagen) and kept at -80C until further use. Genomic DNA was extracted using the DNA/RNA Mini Kit (Qiagen) according to manufacturer’s instructions. The DNA was eluted in Elution Buffer (10 mM Tris-Cl, pH 8.5) and the DNA concentration was measured using the Qubit dsDNA Broad Range kit (Thermo Fischer Scientific). The purity of the extracted DNA was verified using the QIAexpert spectrophotometer (Qiagen), considering the OD260/280 absorbance ratio.

All analysis in R was conducted in version 4.0.3 (44). The sequencing data were deposited in the Gene Expression Omnibus (GEO) with accession number GSE238121.

The genes closest to CpGs that displayed methylation change > 5% between the investigated conditions with the same direction of the effect in BS-seq and EM-seq analysis, were selected for downstream analysis. Different networks with the different genes were created using STRING and community clustering with 10 different clusters was performed (54). The pathways were imported to Cytoscape and were graphically segregated depicting the differentially methylated positions in each network, applying Gene Ontology analysis through BiNGO (55, 56). GO overrepresentation analysis (https://www.webgestalt.org/) was performed using default parameters and the genes targeted with custom designed probes as a background. Representative GO terms with significant FDR and adj.p-value < 0.05 were graphically depicted in the pathway analysis.

We developed a probe-based DNA capture approach to investigate the genetic variation and DNA methylation in the HLA class II locus and proximal promoters of 2,157 immune genes ( Supplementary Table 1 ). The HLA class II genes and the extended class II region ( Figure 1A , hg38, chr6:32,426,802-34,167,129:1) were covered using the density of 3X probe tiling, which provided the largest in-silico coverage (96.66%) of the target locus in comparison to higher probe tiling, i.e.4X (96.62%) and 5X (96.56). The selection of the additional 2,157 immune genes relevant for myeloid and antigen presenting cell function was based on the myeloid panel genes from Nanostring as well as genes found to be affected by inflammatory stimulation of THP-1 myeloid cell line ( Supplementary Table 2 ). The genes that displayed differential expression following stimulation in several independent RNA-sequencing studies were prioritized (38–42). Each proximal promoter (average size between 100-200 bp) was covered according to regional specific parameters yielding the density of 1-3X probe tiling. The final design of the capture probes (Design ID: S3381502) comprised 120,990 probes covering 2,223 kb, of which 65,384 (54%) cover 1,621 kb of the HLA class II locus, while 55,606 (46%) probes cover 602 kb of promoter sequences ( Table 1 ).

We adapted the Agilent SureSelect Methyl-Seq platform (G7530-90002, Version E0, April 2018) ( Figure 1B ) (57, 58) for simultaneous DNA and methylome analysis producing a Variant Library (DNA-seq) and Bisulfite Conversion Library (BS-seq), respectively. Additionally, we tested the compatibility of the protocol for methylome analysis using enzymatic reaction by generating an Enzymatic Methylation Library (EM-seq) (59), which is becoming increasingly popular due to less degradation of the DNA compared to BS ( Figure 1C ).

We addressed the performance of the capture probe design using the human THP-1 immortalized monocytic cell line commonly used to study monocyte/macrophage differentiation and function. Following stimulation with LPS and IFNγ, THP-1 cells upregulate expression of HLA class II genes, co-stimulatory molecules and proinflammatory cytokines (60, 61), thus providing a suitable model for the study of APC cell functions. In addition, THP-1 cells are heterozygous for the HLA class II locus, carrying the HLA-DR1 and HLA-DR15 haplotypes (62), allowing us to test the effect of haplotype sequence variation on the efficiency of the DNA capture. The cells were stimulated with IFNγ/LPS or IFNγ only for 24 hours. Three libraries for each sample (six libraries in total, hereafter referred to as: DNA-IFNγ/LPS, BS-IFNγ/LPS and EM-IFNγ/LPS and DNA-IFNγ, BS-IFNγ, EM-IFNγ) were prepared, amplified, indexed and pooled for sequencing with an average of 4.1 Gb reads per sample and 1000X coverage.

DNA sequencing reads from IFNγ/LPS- and IFNγ-derived samples, which have identical genomes and thus represent technical replicates, were quality and adapter trimmed and subsequently mapped to the human reference genome hg38 ( Figure 1D ). Noticeably, approximately 42% (15,771/37,336 and 15,622/37,0700) of the reads were filtered in the deduplication step for both samples. After deduplication, 99.2% and 99.4% of the reads mapped with an average insert size of 275.7 and 287.3, respectively, and an average read quality for both samples of 33.1, indicating good reliability of the sequencing data.

Next, we assessed coverage across the 4,805 captured regions in DNA-IFNγ/LPS and DNA-IFNγ samples after applying a coverage filter of 10X ( Figure 2A ). Overall, 98.7% (4,744 of 4,805) and 98.5% (4,734 of 4,805) of the designed regions were successfully captured for each sample. We estimated a median coverage of 122X (range 11X - 4759X) across the genomic regions with little variation between the samples ( Table 2 ). A locus on chromosome 1 (hg38, chr1: 40,691,614-40,692,158) displayed surprisingly high coverage in both samples ( Figure 2A ). As expected, the extended HLA class II region (h38, chr6:32,426,802-34,167,129:1) had a larger median coverage of 202X (range 12X-789X) ( Table 2 ; Figure 2B ).

On average, 99.2% of all the sites were covered with at least 10X, 96.55% with > 30X and 83.2% with > 100X, with little variation among the replicates ( Figure 2C ; Supplementary Table 3 ). Noticeably, sites located upstream and downstream of the target regions were also captured ( Figure 2C ), with for example ~69% and 52% of the sites within ± 100 bp and 200 bp windows proximal to the targeted region, respectively, being covered with > 30X. Again, the coverage at the extended HLA class II region was higher, as more 99.7%, 98.7% and 93.3% of the sites) were covered with than 10X, 30X and 100X, respectively ( Figure 2D ; Supplementary Table 3 ).

Taken together, the data shows that the designed probes successfully capture the genomic regions of interest, and that the data are reproducible between technical replicates and provide sufficient coverage for downstream analysis such as variant calling.

We next called genotypes and inferred HLA haplotypes from the captured regions using bcftools (46) and HLA-HD (48), the latter method has previously been benchmarked for its high accuracy of HLA class II typing (63, 64). Genotypes could be called for 96.8% (2,154,496/2,221,666) and 96.6% (2,151,062/2,221,666) of the sites, respectively in each sample, at a minimum coverage of 30X. Variance, called as alternative alleles compared to the hg38 reference, was detected at 10,025 and 10,030 individual sites, respectively, the vast majority (98%, 9,822) of the variants overlapping between replicates.

As previously mentioned, the THP-1 cell line is heterozygous for the HLA class II locus, carrying both the HLA-DR1 and the HLA-DR15 haplotypes. Consistent with this, we called heterozygous genotypes within HLA-DRB1 as well as -DRB6, -DQA1, -DQB1 loci, whereas HLA-DQA2, -DQB2 and HLA-DRA were homozygous based on the bcftools-derived variants ( Figure 2E ). HLA-specific typing based on the HLA-HD algorithm and dictionary, which allows for up to 6-digit calling, confirmed heterozygosity at HLA-DRB1 (*01:01:01/*15:01:01), HLA-DQA1 (*01:01:01/*01:02:01), HLA-DQB1 (*05:01:01/*06:02:01) and HLA-DRB6 (*01:01/*01:02) haplotypes (alleles). Furthermore, HLA-DMB (*01:01:01/01:03:01) and HLA-DOA (*01:01:02/*01:01:01), HLA-DPA1 (*01:03:01/02:02:02), HLA-DPB1 (*04:02:01/02:01:02) were typed as heterozygous, whereas HLA-DMA*01:01:01, HLA-DOB*01:01:01, HLA-DRA*01:01:01, HLA-DRB5*01:01:01 and HLA-DRB9*01:01 were homozygous, as expected ( Table 3 ). HLA-DRB2, -3, -4, -7 and -8 pseudogenes were not typed. Of note, the HLA-HD-based typing conducted on quality-processed and adapter-trimmed fastq files yielded consistent results before and after deduplication.

In conclusion, the designed array efficiently captured, amongst others, the HLA-DR1 and HLA-DR15 haplotypes characterizing THP-1 cells and combined with the HLA-HD tool, allows for up to 6-digit resolution of HLA haplotypes.

Genomic DNA samples derived from cells stimulated with either IFNγ or IFNγ/LPS were simultaneously assessed for DNA methylation after BS or EM conversion ( Figure 1C ), both methods resulting in the conversion of unmethylated cytosines to thymines, with methylated cytosines left unchanged.

Quality and adapter filtered reads were mapped against an in silico converted human reference genome (hg38), which accounts for C>T conversion at non-methylated sites. The total number of reads differed slightly between conversion methods but not between conditions ( Supplementary Table 4 ). Accordingly, after deduplication, ~85% of the BS-derived and ~77% of the EM-derived reads mapped uniquely to hg38 reference ( Supplementary Table 4 ), suggesting that the SureSelect protocol performs slightly better with downstream analysis using BS compared to EM-based conversion. The difference was less pronounced in the extended HLA class II region (GRCh38, chr6:32,426,802-34,167,129:1), where 90.0% (825) and 87.7% (804) of the 917 designed regions were captured by each method, respectively. The median coverage across BS-derived libraries was on average 91.5X (range 11X - 3263X) and reached 51X (range 11X - 1851X) across EM-derived libraries ( Table 4 ). A higher coverage was generally estimated in the extended HLA class II region with an averaged median coverage for BS-seq of 151X (range 11X - 742X) and of 110X (range 11X - 446X) for EM-seq. Of note, all coverage estimations were conducted based on strand-merged CpGs and after applying a filter of 10X. The difference in sequencing depth and fraction of uniquely mapped reads between BS- and EM-derived libraries is further reflected in the number of designed regions that were captured for the DNA methylation libraries with an average of 89.8% (4,313/4,805) for BS-seq and 82.0% (3,941/4,805) for EM-seq, respectively ( Table 4 ). The median coverage was also visually assessed with Manhattan plots across chromosomes, which showed an expected peak at chromosome 6 and consistency across libraries ( Figure 3A ). The peak on chromosome 6 was specific to the targeted class II region (GRCh38, chr6:32426802-34167129:1) and drastically decreased in up- and downstream regions ( Figure 3B ).

We next asked whether DNA sequences from the BS- and EM-seq methylation data could be used for genotype calling using bcftools. The analysis yielded identical results compared to DNA-sequencing data, with HLA-DRB1,-DRB6, -DQA1 and -DQB1 loci being heterozygous, while homozygosity characterized HLA-DQA2, -DQB2 and HLA-DRA loci, based on the bcftools-derived variants ( Figure 3C ).

In conclusion, the designed regions are also enriched with BS and EM conversion, this enrichment is reproducible between technical replicates, and ensures sufficient coverage for estimation of DNA methylation levels.

As a first step to accurately estimate DNA methylation, we assessed the coverage at individual CpG sites within the captured regions. The higher the coverage, the more reliable the estimation. An average of 65,908 and 62,813 CpG sites were covered with > 10X within the designed regions of BS- and EM-seq libraries, respectively, with approximately a third of them, i.e. 33.8%, 22,264/65,908 for BS-seq and 35.3% 22,179/62,813 for EM-seq, mapping to the extended HLA class II region. The relatively large percentage of CpG sites located outside of the HLA class II region is likely due to the enrichment of promoters within our design, which are generally enriched in CpG sites compared to other genomic features (65). Importantly, the majority of the total CpG sites, 98.9% (64,772/65,501) for BS-seq and 97.8% (60,625/61,996) for EM-seq, overlapped between conditions at 10X coverage, which allows for assessing changes in DNA methylation across samples. As expected, the number of CpG sites dropped with an increasing coverage threshold, resulting to an average of 76% and 57% (> 50X) and 53% and 36% (> 100X) of initial CpGs found using BS-seq and EM-seq, respectively ( Figure 4A ; Supplementary Table 5 ). This was likely affecting non-HLA regions since, in contrast, the number of CpGs at the HLA class II locus remained above 93% for both sequencings at 50X and as high as 88% and 84% for BS-seq and EM-seq, respectively, at a 100X coverage depth ( Figure 4B ; Supplementary Table 5 ).

We then compared the conversion efficiency of the two methods. DNA methylation levels were estimated as the number of Cs divided by the total number of sites. Apart from neurons and stem cells, DNA methylation generally occurs in a CpG context. Hence, non-CpG sites are commonly used for estimating conversion rates. The analysis of conversion efficiency showed that EM-seq seemed to perform better as only 0.6% of total non-CpGs were methylated, compared to an average of 1.5% for BS-derived samples ( Supplementary Table 4 ).

In the next step, we compared the methylation levels detected by each method in the same samples, i.e. cells treated with either IFNγLPS or IFNγ at 10X coverage. We observed a linear correlation between methylation levels at all CpGs between EM-seq and BS-seq at all CpGs (Spearman r = 0.79, for IFNγ/LPS, r = 0.79 for IFNγ) as well as CpG sites located in the HLA class II locus (Spearman r = 0.93 for IFNγ/LPS, r = 0.94 for IFNγ) ( Figure 4C ).

We next functionally explored the differentially methylated positions (DMPs) identified in the samples exposed to IFNγ alone or in combination with LPS. We focused on the DMPs that displayed > 10% difference in one method and at least 5% in another with the same directionality in both methods for methylation analysis (BS- and EM-seq). We performed STRING pathway analysis followed by community clustering to subdivide the number of pathways and graphically illustrate the differential methylation observed in core genes of each network ( Figure 5A ). In total, 329 DMPs mapping either promoters or enhancers in close proximity of promoter/enhancers of 150 genes ( Supplementary Table 6 ) were selected for the analysis and segregated in the different pathways ( Figure 5A ). From these genes, some displayed hypomethylated and hypermethylated DMPs with or without hypomethylated/hypermethylated DMRs. In a next step we performed GO analysis of the differentially methylated genes and found significant (FDR < 0.05) enrichment of ‘biological processes’ related to MHC complex assembly, antigen processing and presentation via the MHC class II and antigen processing and presentation of exogenous peptide antigen pathways ( Figure 5B ). Moreover, we detected significant enrichment of the altered genes in the interferon-gamma mediated signaling pathway and response to interferon-gamma. Consistent with this, terms related to ‘molecular functions’ showed enrichment for the MHC class II receptor activity, peptide antigen binding and antigen binding ( Figure 5B ). GO categories linked to ‘cellular component’ highlighted enrichment in the MHC class II protein complex as well as pathways involved in ER and Golgi transportation and endosomal processing ( Figure 5B ).

These data demonstrate that with the current design we can simultaneously interrogate both the HLA regulation as well as the functional state of the myeloid immune lineage.