ARTICLE https://doi.org/10.1038/s41467-022-29450-x OPEN Single-cell Atlas of common variable immunodeficiency shows germinal center- associated epigenetic dysregulation in B-cell responses Javier Rodríguez-Ubreva 1,2,20 ✉, Anna Arutyunyan3,4,20, Marc Jan Bonder5,6,7, Lucía Del Pino-Molina 8, Stephen J. Clark 9, Carlos de la Calle-Fabregat1,2, Luz Garcia-Alonso3, Louis-François Handfield3, 1234567890():,; Laura Ciudad1,2, Eduardo Andrés-León 10, Felix Krueger 11, Francesc Català-Moll1,2, Virginia C. Rodríguez-Cortez2, Krzysztof Polanski3, Lira Mamanova3, Stijn van Dongen3, Vladimir Yu. Kiselev3, María T. Martínez-Saavedra 12, Holger Heyn 13,14, Javier Martín10, Klaus Warnatz 15,16, Eduardo López-Granados8, Carlos Rodríguez-Gallego 12,17, Oliver Stegle 5,6,7,18, Gavin Kelsey 4,9, Roser Vento-Tormo 3,4,21 ✉ & Esteban Ballestar 1,2,19,21 ✉ Common variable immunodeficiency (CVID), the most prevalent symptomatic primary immunodeficiency, displays impaired terminal B-cell differentiation and defective antibody responses. Incomplete genetic penetrance and ample phenotypic expressivity in CVID sug- gest the participation of additional pathogenic mechanisms. Monozygotic (MZ) twins dis- cordant for CVID are uniquely valuable for studying the contribution of epigenetics to the disease. Here, we generate a single-cell epigenomics and transcriptomics census of naïve-to- memory B cell differentiation in a CVID-discordant MZ twin pair. Our analysis identifies DNA methylation, chromatin accessibility and transcriptional defects in memory B-cells mirroring defective cell-cell communication upon activation. These findings are validated in a cohort of CVID patients and healthy donors. Our findings provide a comprehensive multi-omics map of alterations in naïve-to-memory B-cell transition in CVID and indicate links between the epigenome and immune cell cross-talk. Our resource, publicly available at the Human Cell Atlas, gives insight into future diagnosis and treatments of CVID patients. A full list of author affiliations appears at the end of the paper. NATURE COMMUNICATIONS | (2022)13:1779 | https://doi.org/10.1038/s41467-022-29450-x | www.nature.com/naturecommunications 1 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-022-29450-x C ommon variable immunodeficiency (CVID), the most methylomes15 of around 200 single cells (Supplementary Data 2) frequent symptomatic primary immunodeficiency1, is corresponding to naïve, unswitched memory (US-mem) and swit- represented by a heterogeneous group of patients with low ched memory (S-mem) B cells from the control and CVID twins serum immunoglobulin concentrations, defective specific anti- (Fig. 1a and Supplementary Figs. 1 and 2). We then imputed body production and increased susceptibility to bacterial infec- missing methylation values using DeepCpG16 per cell type, and tions of the respiratory and gastrointestinal tracts2. A high generated merged data from the same subpopulations from each proportion of CVID patients have a low frequency of memory B donor (pseudo-bulks) to outline their methylomes. cells3. In fact, CVID is mainly considered to be a humoral In parallel, we performed whole-genome sequencing (WGS) of immunodeficiency, although several immune cell types might be the MZ twins to study a potential genetic origin of the CVID affected4–7. discordance observed in these individuals. WGS analysis revealed CVID is usually detected sporadically in patients with no the existence of 1,400 somatic variants in the CVID twin family history of immunodeficiency. However, for about 10% of (Supplementary Data 3), which is in the same range as the subjects, other first-degree relatives may either be hypo- variants observed in twin cohorts17,18. According to ClinVar19 gammaglobulinemic or present selective IgA deficiency8. The there are only two differential SNPs and two differential indels existence of a genetic component in CVID has been recognized between the twins in genes potentially related to a disease for decades and next-generation sequencing has revealed a (VCV000332641.2, VCV000670396.1, VCV000801645.1 and number of pathogenic genes9. However, CVID pathogenesis VCV000402355.1), although none of these genetic variants have remains largely unexplained, given that only 20% of CVID cases previously been linked to CVID. We also analyzed the effect of can be accounted for by monogenic gene defects, and genetic the identified genetic variants using SNPEff20 and observed that mutations remain elusive in the majority of CVID cases. The only one of them induced a relevant effect on the gene TAS2R31. absence of mutations for a high proportion of CVID patients The link of this gene to CVID or to any other disease has not together with the incomplete disease penetrance for those har- been described (Supplementary Data 3). Of note, we identified boring mutations suggest the operation of additional, as yet genetic variants shared by both CVID-discordant twins affecting undefined, mechanisms (e.g., polygenic, epigenetic, and envir- genes previously associated with CVID, such as CR2, NFKB2, onmental contributors) that help determine the CVID clinical CD19 and TNFRSF13B. This observation reinforces the hypoth- phenotype. esis that additional alterations, such as epigenetic changes, might Recently, it was reported the occurrence of DNA methylation be contributing factors to the differential phenotype of these alterations associated with CVID, using a pair of monozygotic siblings to the development of CVID. (MZ) twins discordant for the disease and a cohort of patients In relation to the single-cell methylation analysis, in the control and controls, providing a concrete proof-of-principle of an epi- twin, DNA methylation levels of both US-mem and S-mem B cells genetic dimension to CVID10. That work provided the first evi- decreased relative to naïve B cells (Fig. 1b, top, and Supplementary dence of epigenetic dysregulation in primary antibody Fig. 3a), consistent with the progressive demethylation occurring deficiencies (PADs) and has spearheaded interest in epigenetic during B-cell differentiation12,13. However, we found that S-mem B aberrancies in the broader field of primary immunodeficiency cells in CVID displayed impaired DNA demethylation throughout research11. However, that study only interrogated a limited set of the entire genome (Fig. 1b, bottom, and Supplementary Fig. 3a). The CpG sites and did not explore links with the functional impact of methylomes of naïve and US-mem B cells from the CVID sibling such epigenetic defects. DNA methylation has an established were virtually identical to those of their control counterparts, function in B-cell differentiation and biology12,13. It is likely that whereas the methylomes of S-mem B cells displayed profound DNA methylation alterations are accompanied by other epige- differences between the two twins, as ascertained by Pearson netic changes and by transcriptional alterations in the B-cell correlation and principal component analysis (PCA) (Fig. 1c and compartment of CVID patients. In this regard, single-cell tran- Supplementary Fig. 3b, c). scriptomics and epigenomics emerge as appropriate approaches To identify the functional loci potentially involved in transcrip- for studying the diversity, plasticity and adaptability of immune tional regulation, we defined differentially methylated regions cell subsets14. (DMRs) using Metilene21 (Supplementary Data 4 and Supplemen- In this work, we generate and integrate single-cell multi-omics tary Fig. 1). In the transition from naïve-to-US-mem B cells, we datasets (DNA methylation, chromatin accessibility, and tran- detected 1537 DMRs that underwent hypomethylation in both scriptomics) corresponding to the aforementioned pair of CVID- control and CVID twins. However, of the 628 DMRs that became discordant MZ twins and a cohort of patients and healthy donors hypermethylated in control, 192 were not hypermethylated in to understand the heterogeneity and impact of epigenetic dysre- CVID (Fig. 1d). The most dramatic alterations were found in the gulation in CVID. Our data indicate prevalent and widespread transition from naïve-to-S-mem B cells. Specifically, 3745 of the alterations affecting DNA methylation and chromatin accessi- 9204 hypomethylated DMRs (40.7%) of the naïve-to-S-mem B cells bility of memory B cells, and provide evidence that these changes comparison of the healthy twin, did not undergo demethylation in impact the transcriptome of these cells following activation. These his CVID sibling. Of the 1555 hypermethylated DMRs in the naïve- epigenetic and transcriptional alterations show defective immune to-S-mem B-cell differentiation in the healthy control, 404 (26%) cell–cell communication, which may explain the impaired were not hypermethylated in CVID (Fig. 1d, e). No DMRs were immune responses observed in these patients. found in the comparison of naïve B cells from CVID and healthy control. FACS analysis showed that the proportions of IgG+ and IgA+ B Results cells within the memory compartment of the twins are comparable DNA methylation alterations in CVID in memory B-cell dif- (around 50% of each in both individuals), eliminating the possibility ferentiation. To define the DNA methylation profiles in CVID, we that altered methylation observed in the CVID twin was a performed single-cell whole-genome bisulfite sequencing of B cells consequence of a different memory B-cell subtype abundance. isolated from a pair of CVID-discordant MZ twins (Fig. 1a) without Our results indicate that CVID patients display methylation identified pathogenic mutations. Clinical information, vaccination defects that do not preexist in naïve B cells but are established response and B-cell phenotype of these individuals are displayed in during memory B-cell differentiation, and affect only specific Supplementary Data 1. Specifically, we generated the whole DNA genomic regions. Since most of the alterations detected in CVID 2 NATURE COMMUNICATIONS | (2022)13:1779 | https://doi.org/10.1038/s41467-022-29450-x | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-022-29450-x ARTICLE a CH3 DNA methylation (scBS-seq) CH3 IgD CD19 Naïve BCs Chromatin accessibility DMRs identification (scATAC-seq) IgD CD19 CD27 US-mem BCs Transcriptomics AAAAA Healthy CVID (scRNA-seq SS2) CVID discordant-twins Single cell omics CD19 CD27 S-mem BCs Healthy CVID Peripheral blood B cells sorting Additional cohort Amplicon-seq of DMRs AA AA AA A AAAAA PBMCs stimulation Transcriptomics Cell-cell communications CITE-seq (scRNA-seq 10x) b c d Pearson correlation coefficient 22 21 20 1 10000 19 18 0.7 1.0 17 Similar hypo 2 16 8000 Altered hypo Number of DMRs S-mem BCs 0.81 15 Similar hyper 3 CVID.US-mem BCs 0.74 0.9 6000 Altered hyper 14 Control 4 US-mem BCs 13 twin 0.97 0.73 0.9 4000 12 5 CVID.Naïve BCs 0.94 0.94 0.72 0.89 2000 11 Naïve BCs 0.98 0.95 0.94 0.72 0.9 6 0 10 Naïve BCs 7 DNA methylation Naïve Naïve 9 8 s s s s s vs. vs. BC BC BC BC BC US-mem BCs US-mem S-mem ve em em em em <0.65 >0.95 aï m m m -m S-mem BCs 22 .N 21 S- S- S- .S 20 D 1 U .U 19 D VI D VI C 18 VI C 17 C 2 16 15 14 3 e DMR chr12: 31,100,759 - 31,101,148 DMR chr19: 605,091 - 605,988 CVID DDX11 HCN2 DNA Methylation 4 twin DNA Methylation 13 1.0 1.0 0.8 0.8 12 5 0.6 0.6 Naïve BCs 0.4 0.4 CVID.Naïve BCs 11 S-mem BCs 0.2 0.2 CVID.S-mem BCs 6 0.0 10 0.0 7 9 0.5 Kb 0.5 Kb 8 f g h i 14 H3K4me1 Naïve BCs 6 GC BCs Control twin 2 Odds Ratio PMDs-DMRs S-mem BCs Normalized gene expression 12 CVID twin 4 Plasma cells (25.4%) Total number of cell divisions 2 10 1 0 H3K27ac 8 6 Odds Ratio 0 4 6 2 4 −1 0 H3K27me3 2 Non-PMDs-DMRs 6 Odds Ratio (74.6%) −2 4 0 s s s lls 2 BC BC BC ce s s s BC BC BC ve C em a m G 0 aï ve em em m as -5Kb DMR +5Kb N S- aï Pl m m N S- S- Distance to DMR U affected DNA methylation loss and occurred in S-mem B cells, we the impaired DNA demethylation observed in CVID S-mem B focused on those regions and referred to them as CVID.no- cells results from alterations in passive cell-division-associated demeth DMRs hereafter for brevity (Supplementary Data 5). demethylation, we analyzed their replicative history using the A high proportion of the DNA demethylation taking place KREC assay24 (Formula 1). Consistent with our previous during cell differentiation occurs in genomic regions named findings25, S-mem B cells of the CVID twin displayed a lower partially methylated domains (PMDs), which are characterized by cell-division rate (Fig. 1f), which would be compatible with late replication and passive demethylation22,23. To test whether reduced levels of passive PMD-associated demethylation. NATURE COMMUNICATIONS | (2022)13:1779 | https://doi.org/10.1038/s41467-022-29450-x | www.nature.com/naturecommunications 3 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-022-29450-x Fig. 1 CVID methylation defects take place in the memory compartment at genomic regions involved in B-cell function. a Scheme depicting sorter strategy for B-cell isolation from the CVID-discordant twins and the cohort of CVID patients and healthy controls, as well as the multi-omics single-cell approaches used. b Circular representation of DNA methylation levels in naïve B cells, US-mem B cells, and S-mem B cells from control and CVID twins. Histogram tracks represent the average methylation levels over 5 Mb windows throughout the genome. Autosomal chromosomes from 1 to 22 are represented. c Heatmap showing Pearson correlation coefficients from the comparison between control and CVID within the various B-cell compartments. d Bar plot representing the number of significant DMRs identified (q < 0.05 and meth.diff > 20%) in the transition from naïve to US-mem B cells or from naïve to S-mem B cells. e Selected examples of smoothed DNA methylation data in altered DMRs in the CVID twin. f Bar plot indicating the total number of cell divisions in naïve, US-mem and S-mem B cells in control and CVID twins. g Overlap of CVID.no-demeth DMRs with partially methylated domains (PMDs, light blue; non-PMDs, dark blue). h Violin plots representing normalized expression levels of the genes nearest to CVID.no-demeth DMRs in naïve, GC and S-mem B cells, as well as in PCs. Gene expression data (n = 1 biological independent sample for each B-cell population) were obtained from publicly-available RNA-seq datasets (see Methods). Whiskers correspond with the minimum and maximum values of the dataset (excluding any outliers). The box is drawn from Q1 to Q3 with a horizontal line to indicate the median. i Plots showing odds ratios in CVID.no-demeth DMRs ± 5 Kb flanking regions in naïve, GC and S-mem B cells, as well as PCs for different enhancer-associated histone modifications (H3K4me1, H3K27ac, and H3K27me3). Active enhancer regions were defined for the concurrence of H3K4me1 and H3K27ac, repressed enhancers by H3K4me1 and H3K27me3, whereas primed enhancers by H3K4me1 exclusively. Histone modifications data were obtained from publicly-available ChIP-seq datasets (see Methods). Source data for panels in this figure are provided in the Source Data Fig. 1 file. However, the majority of the CVID.no-demeth DMRs overlapped involved in B-cell differentiation and function, including members with non-PMDs and with early replication genomic regions of the bZIP family (e.g., BATF, JUNB, Fosl2, and Fra2), NF-kB (p65 (Fig. 1g and Supplementary Fig. 3d), indicating a prevalence of and p50), CTCF, IRFs and PU.121 (Supplementary Fig. 4a). Using alterations in active demethylation mechanisms during S-mem B- publicly-available ChIP-seq data from human primary B cells and cell differentiation in CVID. We found no significant differences immortalized B cells (accession numbers indicated in Methods), we in the expression of methylcytosine dioxygenase Ten-eleven observed a greater degree of binding for several of these TFs, such as translocation (TET) enzymes TET1, TET2, and TET3, involved in CTCF, BATF and JUNB in CVID.no-demeth DMRs, which implies active demethylation, in the different B-cell compartments that the regions are indeed bound to these TFs in vivo in B cells (Supplementary Fig. 3e), suggesting a defect in the enzyme (Fig. 2d and Supplementary Fig. 4b). Motif-enrichment analysis of activity and/or its recruitment to genomic sites rather than an CVID.no-demeth DMRs at the single-cell level showed similar motif alteration in the levels of these proteins. enrichment for TFs, such as CTCF, NF-kB-p65, IRF8, and members To evaluate the functional relevance of the CVID.no-demeth of the bZIP family, in almost all of the individual cells, regardless of DMRs, we analyzed publicly-available gene expression data from the number of altered DMRs they presented (Fig. 2e). various B-cell populations (accession numbers indicated in Our results indicate that, despite the heterogeneous demethy- Methods). We observed that a large number of the nearest genes lation defects affecting CVID S-mem B cells, the entire to these CVID.no-demeth DMRs were specifically upregulated in compartment displayed DNA methylation alterations in common germinal center (GC) B cells compared with naïve B cells, S-mem TF binding sites that might compromise the proper binding of B cells or plasma cells (PCs) from healthy individuals (Fig. 1h) those TFs. and corresponded to functional categories related to B-cell biology (Supplementary Fig. 3f). Additionally, the vast majority Chromatin accessibility changes in regulatory regions in CVID. of the CVID.no-demeth DMRs were located at intronic and In addition, we performed single-cell assay for transposase-accessible intergenic regions (Supplementary Fig. 3g), suggesting a function chromatin with sequencing (scATAC-seq) of the B-cell compart- in modulating regulatory elements. In this regard, using publicly- ment to determine chromatin accessibility profiles and explore their available ChIP-seq data of histone marks (accession numbers relationship with the previously identified DNA methylation indicated in Methods), we observed that CVID.no-demeth DMRs alterations in CVID. Uniform Manifold Approximation and Pro- were more highly enriched in active enhancer histone marks jection (UMAP) visualization showed that the three B-cell subsets (H3K4me1 and H3K27ac) in GC B cells, S-mem B cells and PCs analyzed (naïve, US-mem and S-mem B cells) were clustered on the than in naïve B cells (Fig. 1i). All these results suggest that CVID basis of their chromatin accessibility status (Fig. 3a). aberrant hypermethylation in S-mem B cells occurs at genomic We then used cisTopic26 to identify chromatin accessibility regions that control B-cell function, comprising genes and modules (groups of co-occurring accessible chromatin regions enhancer regions that become activated during GC reaction. referred to hereafter as topics) in scATAC-seq data for the three B-cell subpopulations (Supplementary Fig. 1). A topic may be Heterogeneous DNA methylation alterations define key TFs in linked to a particular biological feature in the cell. As in the CVID. We assessed the heterogeneity of DNA methylation pro- previous DNA methylation analysis, most CVID-associated files in the S-mem B-cell compartment from the CVID- differences in chromatin accessibility occurred in the S-mem discordant twins. PCA based on hypomethylated DMRs in the B-cell compartment (Fig. 3b). Among the 18 regulatory topics naïve-to-S-mem B-cell transition revealed greater methylation identified using cisTopic, we selected six of them with the greatest heterogeneity in S-mem than in naïve B cells (Fig. 2a). Some of differences (above the 95th and below the 5th percentiles). Topics the cells within the CVID S-mem B-cell compartment displayed 2, 7, and 8 corresponded to less accessible regions in CVID similar DNA methylation profiles to those of control cells, S-mem B cells (CVID-depleted topics), whereas topics 4, 5, and 6 whereas others clustered separately and showed aberrant hyper- defined more accessible regions in CVID S-mem B cells (CVID- methylation (Fig. 2a, b). In-depth analysis of the data (see enriched topics) (Fig. 3b). TF binding motif analysis of these Methods) showed that the entire set of S-mem B cells from the topics indicated the significant enrichment of several members of CVID twin displayed methylation alterations, although the fre- the bZIP, bHLH and ZF TF families in CVID-depleted topics, and quency of these alterations was highly variable (Fig. 2c). enrichment of CTCF in CVID-enriched topics (Fig. 3c), coin- Transcription factor (TF) motif-enrichment analysis of CVID.no- cident with several of the TF binding motifs enriched in the demeth DMRs revealed significant enrichment of several TFs CVID.no-demeth DMRs. 4 NATURE COMMUNICATIONS | (2022)13:1779 | https://doi.org/10.1038/s41467-022-29450-x | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-022-29450-x ARTICLE a b naïve control Principal Component 2 (1.9%) naïve CVID S-mem control S-mem CVID S-mem control twin % DNA methylation Principal Component 1 (82.3%) S-mem CVID twin 0 100 c DMRs number e 2000 1 #1 >1000 2 #1 500-1000 4 #1 IRF1 (IRF) 0 #1 <500 Number of altered DMRs Ets1−distal (ETS) 2 #2 NFkB−p50 (RHD) 1500 IRF2 (IRF) IRF8 (IRF) NFkB−p65 (RHD) 7 #2 #5 #8 #1 ETS:E−box (ETS,bHLH) 8 1000 #1 0 SpiB (ETS) #4 #2 GATA3 (Zf) 9 #6 #1 TFE3 (bHLH) #3 IRF:BATF (IRF:bZIP) 6 #7 500 #1 #1 Fosl2 (bZIP) #9 5 3 Jun−AP1 (bZIP) #1 #1 1 3 #2 #2 CTCF (Zf) Fra2 (bZIP) 0 JunB (bZIP) Atf3 (bZIP) single S.mem BCs from CVID twin Fra1 (bZIP) d DMRs bg CTCF binding RORa (NR) Tcfcp2l1 (CP2) Nrf2 (bZIP) 3 CTCF PAX6 (Paired,Homeobox) bg PAX5 (Paired,Homeobox) CTCF ChIP-seq signal CVID.no-demeth DMRs AP−1 (bZIP) Log2(FC) BATF (bZIP) 2 4 Oct11 (POU,Homeobox) ZNF768 (Zf) Brn1 (POU,Homeobox) 3 CTCF ChIP-seq signal 3.5 Oct6 (POU,Homeobox) 2.5 1 ZNF528 (Zf) 2 ZFP3 (Zf) 1.5 CTCF (Zf) 1 0.5 ZNF382 (Zf) 0 #17 #3 #16 #15 #23 #21 #9 #18 #20 #19 #5 #11 #12 #14 #6 #4 #10 #22 #1 #2 #8 #13 #7 0 -2.5 0 2.5 -2.5 0 2.5 -2.5 0 2.5 Distance to DMR (Kb) Distance to DMR (Kb) Fig. 2 Heterogeneous DNA methylation alterations in CVID S-mem B cells. a Dot plot showing principal component analysis based on the DMRs at the single-cell level of naïve and S-mem B cells in the control and CVID twins. The elliptical shadows correspond to the confidence interval (90% CI) for each subpopulation. b Heatmap representing the DNA methylation percentage of the 300 most variable CVID.no-demeth DMRs at the single-cell level. Dark blue samples correspond to S-mem B cells from the control twin; light blue samples correspond to S-mem B cells from the CVID twin. c Plot showing the number of altered hypomethylated DMRs for each individual S-mem B cell for the entire set of cells analyzed at the single-cell level in the CVID twin. d Heatmap and histogram depicting CTCF binding signal in human primary B cells at CVID.no-demeth DMRs ± 2.5 Kb flanking regions. CTCF binding data was obtained from publicly-available ChIP-seq dataset (see Methods). e Heatmap representing the log2(FC) enrichment for several TF binding motifs at the single-cell level over the altered hypomethylated DMRs in S-mem B cells from the CVID twin. Source data for panels in this figure are provided in the Source Data Fig. 2 file. Subsequently, we defined active (H3K4me1 and H3K27ac), nine identified clusters, cluster 3 (in red) was significantly primed (H3K4me1) and repressed (H3K4me1 and H3K27me3) enriched in CVID B cells compared with control B cells (Fisher’s enhancers occurring in naïve, S-mem and GC B cells as well as in test p = 0.0001). This cluster also displayed the highest number PCs, and investigated how the identified topics were associated (384) of significant differentially accessible regions (DARs, with those regulatory elements. We observed a stronger calculated with limma, p < 0.05). Among those regions, 356 association of CVID-depleted topics with active enhancers in DARs were less accessible in the CVID sample (CVID-depleted), GC B cells and PCs (Fig. 3d). Our results indicate the existence of and 28 were more accessible (CVID-enriched). To further give CVID chromatin accessibility defects in regulatory elements that those DARs meaning in terms of the previous features that we undergo activation upon antigen encountering by naïve or had inspected (Fig. 3d, e), we examined the overlap of TFs memory B cells. Consistent with this observation, we detected binding or enhancer regions (see Methods, using the same greater motif enrichment of several TFs involved in enhancer publicly-available files as for Fig. 3d, e) with those DARs regulation, such as BATF, IRF4, and CTCF27–29 in CVID- (Supplementary Fig. 4d, e). The analysis of this specific cluster depleted topics (Fig. 3e). 3 showed that CVID-depleted DARs were enriched in active In order to validate the findings obtained with cisTopic, we also enhancers in GC or post-GC states (memory B cells and plasma performed an unsupervised clustering of the cells based on cells), which is in line with our previous findings using the scATAC-seq data (Supplementary Fig. 4c). Of note, among the cisTopic tool (Supplementary Fig. 4d). It also indicates that all the NATURE COMMUNICATIONS | (2022)13:1779 | https://doi.org/10.1038/s41467-022-29450-x | www.nature.com/naturecommunications 5 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-022-29450-x a BC b BC Cell subset Donor BC em em ïve -m m Naïve BCs Na US C Control twin S- US-mem BCs CVID twin C S-mem BCs 1 2 3 UMAP2 UMAP2 4 UMAP1 UMAP1 5 6 Topics c CVID depleted CVID enriched 7 -log10(p-value) 8 >10 9 5 Topics: 4 5 6 2 10 -log10(p-value) 0.1 100 11 50 20 12 Diff. Normalized AUC score 10 0.02 13 (CVID-control) 14 log2(FC) in topic vs background 2.0 15 0 1.5 16 1.0 17 -0.02 0.5 18 0.0 f control twin CVID twin control twin CVID twin Naïve BCs Naïve BCs S.mem BCs S.mem BCs Topics: 2 7 8 60 d e 40 CVID CVID 20 CVID CVID depleted enriched depleted enriched 0 Naïve BCs ATF2_GM12878 enhancers S-mem BCs ATF7_GM12878 Active GC BCs BATF_GM12878 PCs CTCF_Bcells CVID.no-demeth DMRs Naïve BCs CTCF_GM12878 enhancers E2F4_GM12878 Primed S-mem BCs GC BCs EBF1_GM12878 PCs ETS1_GM12878 70 IRF4_GM12878 60 Naïve BCs Repressed enhancers NFATC1_GM12878 50 S-mem BCs 40 GC BCs PAX5_GM12878 30 PCs SP1_GM12878 20 Topic 2 Topic 8 Topic 7 Topic 4 Topic 5 Topic 6 Topic 2 Topic 8 Topic 7 Topic 4 Topic 5 Topic 6 10 0 Normalised AUC Score DMR DMR DMR DMR Normalised AUC Score -3.0 -1.5 0.0 1.5 3.0 -3.0 -1.5 0.0 1.5 3.0 center center center center Fig. 3 Chromatin accessibility defects in CVID. a UMAP visualization of unsupervised analysis of chromatin accessibility in naïve, US-mem and S-mem B cells from the control and CVID twin (1179 cells) indicating the B-cell subset and the origin of the cell. b Heatmap representing the difference in normalized AUC scores for the different regulatory topics defined by CisTopic between CVID and control twins in the B-cell subsets analyzed. Red color corresponds to a stronger association of the topic with the CVID twin; blue corresponds to a stronger association of the topic with the control twin. Circle size indicates the significance of the association (Mann–Whitney test). c Dot plot representing the significantly enriched TF binding motifs in selected topics that are most closely associated with chromatin accessibility in control twin (CVID-depleted: topics 2, 7, and 8) or most associated with chromatin accessibility in the CVID twin (CVID-enriched: topics 4, 5, and 6). d Heatmap representing the association of the CVID-depleted and CVID-enriched topics with different enhancer stages during B-cell activation (active, primed and repressed enhancers). Active enhancer regions were defined for the concurrence of H3K4me1 and H3K27ac, repressed enhancers by H3K4me1 and H3K27me3, whereas primed enhancers by H3K4me1 exclusively. e Heatmap showing the association of CVID-depleted and CVID-enriched topics with several TFs binding signal in human immortalized B cells and CTCF in human primary B cells. f Heatmaps and histograms representing chromatin accessibility signals in naïve and S-mem B cells at CVID.no-demeth DMRs ± 2.5 Kb flanking regions in control and CVID twins. Source data for panels in this figure are provided in the Source Data Fig. 3 file. 6 NATURE COMMUNICATIONS | (2022)13:1779 | https://doi.org/10.1038/s41467-022-29450-x | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-022-29450-x ARTICLE TFs peaks examined were more enriched in CVID-depleted memory B cells included B-cell-related upregulated genes such as DARs which is consistent with our previous observations for FCGR2B, CD72, PTPRC, CD79A, CD22, and CCR6, as well as topics 2 and 8 in cisTopic analysis (Supplementary Fig. 4e). These downregulated genes such as CCL22, CD40, ICAM1, CCR7, results may suggest that chromatin accessibility dysregulation in CCL17, CD80, and CD86 (Fig. 4e). the CVID twin is taking place in a specific subset of memory Separate differential gene expression analysis of the two B cells. activation conditions (CD40L + IL-21 or anti-CD3/CD28), iden- We also noted an increase in chromatin accessibility in CVID.no- tified that relevant genes such as CCL22, IL4R, IRF4, ATF5 and demeth DMRs in the naïve-to-S-mem B-cell transition. However, IFNGR2 were specifically dysregulated upon CD40L + IL-21 those regions became accessible despite not being properly activation, whereas CD74, CD79A, CD79B, PTPN6, CD72, demethylated in the CVID twin (Fig. 3f). This analysis suggests CD22, or PTPN1 were dysregulated specifically upon anti-CD3/ that increased chromatin accessibility might be a prior process that CD28 stimulation. Interestingly, other genes such as CD19, is necessary but not sufficient for DNA demethylation and that, in LGALS9, TGFB1, PTPRCAP, and LY6E were similarly dysregu- the CVID twin, defects in the proper recruitment/activation of the lated regardless of the activation condition (Supplementary demethylating machinery is impaired. Data 7). However, we noticed that separate analysis causes a Overall, our results indicate that chromatin accessibility drastic decrease in the statistical power of the analysis with a alterations in CVID S-mem B cells occur in regulatory elements reduction of 50% in the number of detected differentially that become active during the B-cell response. Additionally, the expressed genes. In this regard, a large group of the previously impaired DNA demethylation observed in CVID proves not to be identified DEGs such as PTPRC, NFKB2, TNF, CD80, CD86, a consequence of changes in chromatin accessibility. CD40, CCR6, and CD70, among others, were not detected as differentially expressed now (Supplementary Data 7). The dysregulation of relevant genes for B-cell activation such as CVID epigenetic defects perturb transcription during B-cell CD79A, CD79B, CD22, CD72, and PTPN6, all of them genes activation. DNA methylation and chromatin accessibility defects involved in the modulation of the B-cell receptor pathway, in CVID S-mem B cells are likely to affect gene expression. We specifically upon anti-CD3/CD28 treatment (indirect B-cell therefore performed plate-based Smart-seq230 single-cell RNA- activation via T-cell stimulation) and not upon CD40L + IL-21 sequencing (RNA-seq) in naïve, US-mem and S-mem B cells treatment (direct B-cell stimulation), would suggest the existence from the aforementioned CVID-discordant twins. Unsupervised of extrinsic B-cell defects that might contribute to the impaired dimensionality reduction, clustering and differential expression B-cell responses observed in CVID patients. analysis in each of the subsets revealed virtually identical tran- As shown above, most of the DNA methylation alterations in scriptional profiles in both twins (Fig. 4a, Supplementary Fig. 1 the CVID twin occur in enhancers active during the B-cell and Supplementary Data 6). response. However, we identified a group of altered genes in The transcriptional signatures of naïve and memory B cells are CVID activated memory B cells that were also in the proximity of very similar despite the profound epigenetic reprogramming of aberrantly hypermethylated regions (Fig. 4f). 907 (18%) of the naïve B cells following antigen encounter31. We therefore upregulated genes and 53 (30%) of the downregulated genes in speculated that epigenetic alterations in CVID might influence S-mem BCs of the CVID twin correspond to genes near CVID.no- the transcriptome of activated B cells, instead of the transcrip- demeth DMRs. Interestingly, some of these hypermethylated and tome of steady-state B cells (naïve and memory B cells without downregulated genes were specifically upregulated in memory B any stimulus). This would also be consistent with our analysis cells upon cell activation in comparison with naïve B cells. For showing that DNA methylation and chromatin accessibility instance, selected members of the CD40 signaling pathway (e.g., alterations in CVID are both associated with activated B-cell CD40, TRAF1, and NFKB2), and other genes related to B-cell features (GC upregulated genes and GC active enhancers). response (e.g., CD70, ICAM1 (CD54), CCL17) were specifically Therefore, we studied the transcriptome of activated B cells from hypermethylated and downregulated in activated memory B cells the CVID-discordant twins. To this end, we conducted in vitro in the CVID twin (Fig. 4f). activation using peripheral blood mononuclear cells (PBMCs). To integrate the multiple single-cell-omic data generated in this Specifically, we stimulated PBMCs using a combination of CD40L study (i.e., scBS-seq, scATAC-seq and scRNA-seq data), we firstly and IL-21 (CD40L/21) (to directly activate B cells) or, alternatively, analyzed the expression of relevant TFs. We then related the using α-CD3/CD28 (CD3/28) (to define the influence of the T-cell concurrence of the dysregulated TF with the presence of TF compartment on B cells). We then used droplet-based techniques32 motifs in the regions with DNA methylation and chromatin to create a transcriptomic census of different immune cells of the accessibility alterations. We identified that the expression levels of CVID-discordant twin pair and performed gene expression analysis members of the IRF, ETS and bZIP, TF families such as IRF1, (Fig. 4 and Supplementary Fig. 1). SPIB, JUN and FOS were significantly altered in the activated Under these experimental conditions, cell deconvolution using memory B cells of the CVID twin (Supplementary Data 7). logistic regression and cell markers expression analysis allowed us Interestingly, several members of those TF families were also to identify naïve and memory B-cell subsets (Fig. 4c, d and significantly enriched in the genomic regions displaying epige- Supplementary Fig. 5a). Despite the high similarities between the netic alterations (Supplementary Fig. 5d), both in DNA control and the CVID individual in the expression of B-cell methylation (previously shown in Fig. 2e) and in chromatin activation marker genes (Supplementary Fig. 5b), differential accessibility (previously shown in Fig. 3e). expression analysis identified 870 upregulated and 11 down- Our findings demonstrate that the transcriptional program of regulated genes (FDR < 0.05) in the activated naïve B-cell cluster activated B cells is compromised in the CVID twin and that of the CVID twin (Supplementary Fig. 5c and Supplementary CVID-associated epigenetic alterations might contribute to Data 7). However, in line with the previously identified epigenetic transcriptional perturbations that occur during B-cell response. defects occurring specifically in the memory B-cell compartment, we detected 6-fold more dysregulated genes in activated memory B cells (5116 upregulated and 179 downregulated genes) when the Alterations in immune cell communication in B-cell responses CVID twin was compared with his healthy control sibling (Fig. 4e in CVID. Our single-cell transcriptomics census allowed us to and Supplementary Data 7). Dysregulated genes in activated generate the transcriptional profiles of immune compartments NATURE COMMUNICATIONS | (2022)13:1779 | https://doi.org/10.1038/s41467-022-29450-x | www.nature.com/naturecommunications 7 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-022-29450-x a b Cell type Disease condition Naïve BCs Healthy US-mem BCs CVID S-mem BCs Doublets UMAP2 UMAP2 UMAP1 UMAP1 c d TXN B cell cluster Activation condition Disease condition CD37 CXCR4 CD52 CD82 TCL1A TNFRSF13B CD27 UMAP2 UMAP2 UMAP2 mem BCs Un-stimulated CVID Act_mem BCs CD40L Healthy IGHV4-31 Naïve BCs aCD3 IGHV4-28 UMAP1 Act_Naïve BCs UMAP1 UMAP1 ve s t_ m Cs O mB s o. s s Oligoclonal BCs aï BC e C lig C BC Ac me B m B t_ ïve Ac Na N Naïve BCs Memory BCs e f (activated) (activated) Healthy CVID Healthy CVID twin twin twin twin TRAF1 40 ICAM1 CD19 IRF4 CD40 ATF5 CCL17 -Log(FDR) PTP1N CD22 20 NFKB2 CD79B CD70 CCL22 CCR7 FCGR2B PTPN6 CD72 CD79A ICAM1 TNF CCR6 PTPRC IFNGR2 CD80 CD81 CD74 CD40 CD86 CCL17 IL4R 0 IL2RA z-scores -0.5 0 0.5 −1 0 1 Log(FC CVID/healthy) Fig. 4 CVID transcriptional defects in activated B cells. a Unsupervised UMAP visualizations of SS2 single-cell transcriptome data based on highly variable genes and indicating B-cell subset and cell origin (1324 cells). b Scheme depicting immune cell activation with CD40L + IL-21 (B-cell activation) or α-CD3/CD28 (T-cell activation), the establishment of cell–cell communications after cell stimulation and the use of 10× Chromium technology for single- cell transcriptomic analysis. Graph created with BioRender.com. c Unsupervised UMAP visualizations (5119 cells) showing the predicted classification of B cells into steady state or activated naïve and memory B cells using logistic regression, the stimulation conditions and disease distributions. d Violin plots representing the expression of selected marker genes used to annotate the B-cell subpopulations in the CVID and healthy twins. e Volcano plot representing differentially expressed genes in activated memory B cells comparing control and CVID twins. Some relevant upregulated and downregulated genes are labeled. f Heatmap representing the expression levels of selected dysregulated genes in CVID activated memory B cells corresponding to the nearest gene to a CVID.no-demeth DMR. Activated naïve B cells and activated memory B cells from CVID-discordant twins are represented. Source data for panels in this figure are provided in the Source Data Fig. 4 file. additional to B cells from the CVID-discordant twin pair. Clus- To systematically analyze the effect of cell–cell communication tering and cell type-specific marker gene expression allowed us to on B-cell responses in the healthy and the CVID twins, we used identify several populations of immune cells: B cells (including CellPhoneDB33 (Supplementary Fig. 1), our recent database of naïve and memory B cells), CD4+ and CD8+ T cells (naïve and ligands, receptors and their interaction. This approach, available memory), CD8+ MAIT cells, Treg cells, γδ-T cells, two subsets of on GitHub (https://github.com/Teichlab/cellphonedb), allowed us NK cells, myeloid cells and their corresponding activated subsets to generate a potential cell–cell communication network upon (Fig. 5a, b and Supplementary Fig. 6a). Differential gene expres- immune cell challenge in the CVID-discordant twins. sion analysis between the two CVID-discordant siblings within CellPhoneDB analysis revealed defects in several ligand- each compartment revealed transcriptional dysregulation in receptor (L/R) pairs relevant to B-cell regulation (Fig. 5c and immune subsets other than B cells (Supplementary Data 8). Supplementary Data 8). Remarkably, activated memory B cells 8 NATURE COMMUNICATIONS | (2022)13:1779 | https://doi.org/10.1038/s41467-022-29450-x | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-022-29450-x ARTICLE a act_NK16 NK cells b NK56 CD3D NK16 CD8A act_NK56 CCL5 B cells IL7R T8mem GZMB Myeloid cells Bnaïve Tγδ KLRB1 GZMK Bmem Oligo BCs ANXA1 T8MAI IL2RA act2_T8 TRDV2 act_T8MAI act_Bnaïve PRF1 FCGR3A GZMH T4mem XCL1 act_T4 CCND2 T4reg CD19 act_Bmem HLA-DRA T4naïve IGHV4-28 T cells CD14 UMAP 2 T8naïve act1_T8 CD36 T4mem T4naive T8naive NK56 NK16 act_NK16 T8MAI T8mem Tγδ act_NK56 act_T8MAI act2_T8 Tregs act_T4 act1_T8 act_Bmem Myeloid act_Bnaive Bnaive Bmem Oligo BCs UMAP 1 c DE gene : partner gene Naïve B cells (steady state) Naïve B cells (activated) Memory B cells (steady state) Memory B cells (activated) IL7:IL2RG IL7 receptor IL7:IL7R LogFC (CVID/healthy) ITGAL:ICAM3 1.0 Integrin aLb2 ITGB2:ICAM3 ITGAL:ICAM2 0.5 Integrin aLb2 ITGB2:ICAM2 0.0 ITGAL:ICAM1 Integrin aLb2 ITGB2:ICAM1 -0.5 CD22:PTPRC TNFRSF14:LTA -1.0 CLEC2D:KLRB1 ICAM3:ITGAL Integrin aLb2 ICAM3:ITGB2 CD72:SEMA4D expressing the partner CD74:APP Percentage of cells TNFRSF13B:TNFSF13B 25% LGALS9:HAVCR2 50% CCR6:CCL20 75% CCL17:CCR4 CXCL9:DPP4 100% CXCL9:CXCR3 CXCL10:DPP4 CCL22:DPP4 CCL22:CCR4 CD40:CD40LG ICAM1:ITGAL Integrin aLb2 ICAM1:ITGB2 TNF:TNFRSF1A TNF:TNFRSF1B LTA:TNFRSF14 LTA:TNFRSF1B LTA:TNFRSF1A HLA-E:KLRC1 FAS:FASLG CD38:PECAM1 CD86:CTLA4 CD86:CD28 ICAM1:SPN NK16 act_NK16 NK56 act_NK56 T4naive T4mem act_T4 Treg γδTcells T8MAI T8naive act1_T8 act2_T8 NK16 act_NK16 NK56 act_NK56 T4naive T4mem act_T4 Treg γδTcells T8MAI T8naive act1_T8 act2_T8 NK16 act_NK16 NK56 act_NK56 T4naive T4mem act_T4 Treg γδTcells T8MAI T8naive act1_T8 act2_T8 NK16 act_NK16 NK56 act_NK56 T4naive T4mem act_T4 Treg γδTcells T8MAI T8naive act1_T8 act2_T8 NK cells T cells NK cells T cells NK cells T cells NK cells T cells displayed the highest number of dysregulated L/R pairs (Fig. 5c might constitute a repressive autocrine/paracrine mechanism and Supplementary Fig. 6b). For instance, several genes that code affecting B-cell activation in the CVID twin (Supplementary for negative regulators of the BCR function, such as CD22, CD45 Fig. 6b). Conversely, we also observed that several genes, like (PTPRC), Galectin-9 (LGALS9), and CD72, were significantly CD86 and the chemokine genes CXCL10, CCL22, and CCL17, upregulated in CVID B cells, which might explain the defective were significantly downregulated in CVID B cells. These B-cell response in CVID twins. The interacting partners of these cytokines presented their receptor partners in both NK and upregulated genes in B cells were expressed in several immune T-cell subsets (Fig. 5c). We also found expression changes in cells that are present both in the NK and T-cell compartments other cell compartments that interact with B cells (Supplementary (Fig. 5c). Furthermore, we observed that the negative regulators Fig. 6c and Supplementary Data 8). For instance, in the CD4+ CD22 and CD45 were also upregulated in CVID B cells, which T-cell compartment, we observed significant downregulation of NATURE COMMUNICATIONS | (2022)13:1779 | https://doi.org/10.1038/s41467-022-29450-x | www.nature.com/naturecommunications 9 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-022-29450-x Fig. 5 Cell–cell communication alterations in the CVID twin. a UMAP visualization showing different immune cell populations identified from Louvain clustering and cell-specific marker gene expression (16,901 cells). The B-cell compartment is outlined with a blue dotted line and includes naïve BCs (Bnaïve), memory BCs (Bmem) activated naïve BCs (act_Bnaïve), activated memory BCs (act_Bmem), and oligoclonal BCs (oligo BCs). The T-cell compartment is outlined with a green dotted line and includes CD8+ naïve T cells (T8naïve), CD8+ memory T cells (T8mem), CD8+ activated T cells (act1_T8 and act2_T8), CD8+ MAIT cells (T8MAI), CD8+ activated MAIT cells (act_T8MAI), CD4+ naïve T cells (T4naïve), CD4+ activated T cells (act_T4), CD4+ regulatory T cells (Treg) and γδ-T cells (Tγδ). The NK cell compartment is outlined with a purple dotted line and includes NK CD16+ cells (NK16), activated NK16 cells (act_NK16), NK CD56+ cells (NK56), and activated NK56 cells (act_NK56). Myeloid cells are also represented (b) Violin plots representing the expression of selected marker genes in the cell populations identified. c Overview of selected dysregulated L/R interactions between B cells and the other immune cell compartments in the CVID twin in naïve and memory B cells in steady state and following activation. The scale indicates the log2(FC) gene expression of the defined B-cell subsets in the CVID vs. control comparison. Only differentially expressed genes with an FDR < 0.05 were considered in the analysis. The percentage of other immune cells expressing the partner molecule is indicated by the circle size. Molecules of the L/R pairs expressed in B cells are shown in blue; molecules of the L/R pairs expressed in other non-B immune cells are shown in red. Assays were carried out at the mRNA level, but were extrapolated to protein interactions. Source data for panels in this figure are provided in the Source Data Fig. 5 file. CD40L, the ligand of CD40, which is expressed in B cells and is response such as CCL22, CD70, ICAM1, LGALS9, CD72 or required for B-cell activation. This was also observed at the PTPRC. Inspection of the CVID Ib subcohort led us to identify protein level by FACS (Supplementary Fig. 6d), although due to additional dysregulated genes shared with the CVID twin. the inherent limitations of these particular rare samples, this is a Examples of these genes include CCR6, PTPN1, CD81 and single experiment that should not be taken as definitive proof. CD79A (upregulated) and CCR7, CD80, CD86, and CCL17 All these results indicate that, in addition to B-cell-intrinsic (downregulated) (Fig. 6f and Supplementary Data 11). Some of alterations, defects in other immune cell compartments might these dysregulated genes were validated at the protein level using compromise the correct B-cell response in the context of this CITE-seq analysis (Fig. 6g and Supplementary Data 12, 13, and primary immunodeficiency. 14). For example, we detected protein downregulation of CD54 (ICAM1), CD86, CD80, CD70, CD71 and CD307d (FCRL4), Epigenetic, transcriptional, and communication changes in a together with protein upregulation of IgB (CD79B), CD196 CVID cohort. We performed DNA methylation and tran- (CCR6), CD45 (PTPRC), CD32 (FCGR2), CD19, and CD35 (CR1) scriptomic analysis in a cohort consisting of 10 CVID patients in the activated memory B cells of the CVID patients compared and an equivalent number of healthy controls and compared the with the healthy donors. In addition, the analysis of a larger results obtained in our twin-study. The CVID patients—without number of cells in the CVID cohort allowed us to annotate any known monogenic defect—represented the three CVID additional immune cell populations such as progenitors, DCs and subtypes (Ia, Ib, and II) according to the Freiburg classification34. monocytes among others, besides the ones identified in the twin Clinical information and B-cell phenotype of these individuals are analysis (Fig. 6h and Supplementary Fig. 7d). We also observed provided in Supplementary Data 1. shared dysregulated L/R pairs in the cohort analysis in We first performed amplicon sequencing to interrogate 1058 comparison to the twins dataset not only in the B-cell CpG sites corresponding to 162 selected CVID.no-demeth DMRs compartment (eg. CD72:SEMA4D, CD74:APP, CCL22:CCR4, (Supplementary Fig. 1 and Supplementary Data 9). We observed LGALS9:HAVCR2, and CD40:CD40LG among others) but also that the main methylation alterations observed in the CVID- in compartments beyond B cells such as activated CD4+ T cells discordant twin are also present in the validation cohort (Fig. 6a). (Fig. 6i, left panel, Supplementary Figs. 7e and 8a, b) As Selected DMRs (Fig. 6b), confirmed the occurrence of hypermethy- previously observed for DEGs, we detected additional dysregu- lation in S-mem B cells. Methylation levels in S-mem B cells were lated L/R pairs when focusing in CVID Ib subtype patients (eg. significantly higher for all 3 CVID subtypes represented in our CCR6:CCL20, CXCL9:CXCR3, CXCL10:DPP4, and CD86:CTL4 cohort (Fig. 6c). Interestingly, the effects were most accentuated for among others) (Fig. 6i, right panel and Supplementary Data 15). those CVID patients of the Ib subtype (the subtype to which the CVID twin belongs). Furthermore, in agreement with our previous Discussion works10,25, we observed aberrant hypermethylation in genes such as In this study, the use of multi-omics single-cell atlas technologies BCL2L1, TCF3, BCL6, BCL10, and AICDA among others allowed us to identify functional alterations occurring in specific (Supplementary Fig. 7a). As in the case of the twins (Fig. 1f), scarce immune cell populations that would have gone undetected CVID patients showed lower proliferation rates in both naïve and if conventional bulk-based approaches had been used. We show memory B cells after CD40L + IL-21 stimulation (Supplementary that CVID B cells display widespread heterogeneous DNA Fig. 7b). The existence of these proliferative alterations in the naïve methylation and chromatin accessibility changes that are B-cell compartment raises the possibility of additional primed restricted to the memory compartment. Similarly to other epi- defects established in previous steps of B-cell development. genomic studies for congenital syndromes35,36, the DNA We also obtained the single-cell transcriptomes of steady-state methylation patterns identified in this study might lead to the and activated B cells in the cohort of CVID patients and healthy identification of a diagnostic “episignature” for CVID. Epigenetic controls. We obtained the transcriptome of ~100 K PBMCs signatures in the S-mem B-cell compartment can be used as a stimulated with (i) either CD40L and IL-21 (CD40L/21), or (ii) α- historical record of epigenetic processes occurring during the GC CD3/CD28 (CD3/28) (Supplementary Fig. 1 and Supplementary reaction. In this regard, we have demonstrated that epigenetic Fig. 7c). We observed shared transcriptomic signatures between alterations in CVID memory B cells affect genes and regulatory the CVID-discordant twin pair and the CVID cohort. Within the elements that undergo activation in GC B cells and subsequent B-cell compartment, activated memory B cells displayed the PC differentiation. highest number of dysregulated genes, as previously observed for DNA demethylation changes in CVID S-mem B cells are the CVID twin (Fig. 6d, e and Supplementary Data 10). Several of mainly associated with altered active demethylation mechanisms, these differentially expressed genes in both the CVID cohort and as shown by the predominance of alterations outside PMDs. We the CVID twin, include relevant genes for B-cell immune did not detect any significant changes in the expression levels of 10 NATURE COMMUNICATIONS | (2022)13:1779 | https://doi.org/10.1038/s41467-022-29450-x | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-022-29450-x ARTICLE a b Chr1: 1167415 - 1167510 Chr20: 37438649 - 37438726 c DNA methylation DNA methylation **** **** 100 100 80 80 60 60 **** **** 40 40 **** 20 20 **** 2 2 Normalized DNA methylation Normalized DNA methylation 0 0 Chr21: 29112176 - 29122480 Chr22: 42769446 - 42769520 DNA methylation DNA methylation 100 100 80 80 Healthy 0 0 Healthy 60 60 CVID II 40 40 CVID 20 20 CVID Ia 0 0 CVID Ib Chr4: 2320016 - 2320264 Chr5: 1294267 - 1294373 −2 DNA methylation DNA methylation −2 100 100 80 80 60 60 40 40 20 20 −4 0 0 Naïve Healthy S-mem Healthy S-mem BCs ve em em aï Naïve CVID m m S-mem CVID N S- S- U d e f g Ib D D VI VI CCL22 LGALS9 C C HLA−DRA 12 7 PTPN1 mRNA levels CD74 N.S CD54 (ICAM1) Downregulated Upregulated 6 IL4R 9 5 HLA−DPA1 N.S in CVID in CVID CD52 CD86 (CD86) 6 4 CCR6 3 PTPRC CD80 (CD80) 15 3 2 LGALS9 CD79B CD72 CD71 (CD71) CD72 CD79A CD70 FCRLA CD19 CD70 (CD70) -log10(adj.P.Val) 5 mRNA levels 2.0 CD81 FCGR2B CD307d (FCRL4) 10 1.5 4 CD37 logFC CR1 3 PTPRCAP IgB (CD79B) CD79B 1.0 PTPN6 0.25 2 IFNGR1 0.5 HLA−DQB2 0.00 CD196 (CCR6) CD70 IL2RA CCR7 -0.25 CD45 (PTPRC) CCL22 CL L CD72 CCL17 5 FCRL4 ICAM1 PTPRC CCL22 5 16 LTA -logFDR CD32 (FCGR2) TRAF4 AF4 F F4 PTPRC mRNA levels FCGR2 FCGR2B CGR2B TNF IFIT IFITM IT T 3 4 14 ICAM1 2 LTA T CD19 PTPRC P AP CXCL10 CD19 (CD19) ICAM1 AM IFNGR 12 CXCL9 20 CD40 0 FNGR1 FNGR 1 CD37 3 CD80 CD35 (CR1) 10 CD70 0 2 CD86 40 IFITM3 8 FCRL4 -0.2 -0.1 0.0 0.1 0.2 TRAF4 y D y D lth lth VI VI ea ea C C logFC z-score Healthy H H CVID Ib -2 0 +2 h i Memory B cells CVID (activated) Memory B cells CVID Ib (activated) Bnaïve CD72:SEMA4D CD74:APP : act_Bmem act_Bnaïve * LGALS9:HAVCR2 * CCR6:CCL20 * expressing the partner PTPRC:CD22 Percentage of cells Bmem CXCL9:DPP4 cDC1 act_T8 PCs * CD70:CD27 CXCL9:CXCR3 25% 50% Treg CXCL10:CXCR3 pDC CXCL10:DPP4 75% prog mono16 CCL22:DPP4 100% act_T4 mΦ CCL22:CCR4 cDC1 TGFB1:TGFBR1 TGFB1:TGFBR2 iNKT act_NK56 CD40:CD40LG LogFC (CVID/healthy) NK56 mono14 * ICAM1:ITGAM 0.2 Tγδ * ICAM1:ITGB2 0.1 act_NK16 * ICAM1:ITGAX * ICAM1:ITGB2 0.0 T4mem mono_int * ICAM1:ITGAL * ICAM1:ITGB2 -0.1 T4naïve NK16 * ICAM1:SPN TNF:TNFRSF1A -0.2 UMAP 2 TNF:TNFRSF1B T8MRA LTA:TNFRSF14 T8MAI T8naïve T8EM LTA:TNFRSF1B T8CM LTA TNFRSF1A LTA: UMAP 1 *CD86:CTLA4 CD86:CD28 * Bnaïve act_Bnaïve Bmem act_Bmem PCs T4naïve T4mem act_T4 Treg Tγδ T8naíve act_T8 T8MAI T8CM T8EM T8MRA iNKT NK16 act_NK16 NK56 mΦ mono14 mono16 mono_int cDC1 cDC2 prog Bnaïve act_Bnaïve Bmem act_Bmem PCs T4naïve T4mem act_T4 Treg Tγδ T8naíve act_T8 T8MAI T8CM T8EM T8MRA iNKT NK16 act_NK16 NK56 mΦ mono14 mono16 mono_int cDC1 cDC2 prog B cells T cells NK myeloid B cells T cells NK myeloid cells cells cells cells ten-eleven translocation (TET) demethylating enzymes. However, differentiation into PCs41. In addition, IRF8 is crucial to the given that several TFs interact with TETs37–40, it is likely that development of GC B cells and to the maintenance of the B-cell dysregulation of certain TFs might alter the recruitment of program42,43, and CTCF is required to sustain the GC tran- demethylating activities to specific genomic loci. scriptional program to avoid premature PC differentiation44. Regardless of the heterogeneous degree of CVID-associated Interestingly, several studies have shown how DNA methylation methylation defects observed in the S-mem B-cell compartment, may affect the binding of some TFs, including CTCF, AP-1 and we detected a common enrichment in the binding motif of TFs BATF45–47. Therefore, the observed methylation defects may associated with B-cell function in those altered regions. For compromise the functionality of the entire memory B-cell com- instance, Fra1 is involved in the regulation of follicular B-cell partment in CVID. NATURE COMMUNICATIONS | (2022)13:1779 | https://doi.org/10.1038/s41467-022-29450-x | www.nature.com/naturecommunications 11 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-022-29450-x Fig. 6 DNA methylation, transcriptomic and cell–cell communication defects in a CVID cohort. a Box plot depicting DNA methylation levels of selected DMRs in different B-cell subsets from CVID patients (n = 11) and healthy individuals (n = 10). Whiskers correspond with the minimum and maximum values of the dataset (excluding outliers). The box is drawn from Q1 to Q3 with a horizontal line to indicate the median. Two-sided Wilcoxon test with no multiple test correction, where **** represents p value ≤ 0.0001. b Plot showing smoothed DNA methylation data in altered DMRs. c Box plots depicting DNA methylation levels of the inspected CpG sites in S-mem B cells from CVID patients (n = 11) and healthy controls (n = 10). Two-sided Wilcoxon test with no multiple test correction, where **** represents p value < 0.0001. d Volcano plot representing DEGs in activated memory B cells comparing healthy donors (n = 8) and CVID patients (n = 10). Downregulated (green) and upregulated (orange) genes are plotted. Dysregulated genes in the CVID- discordant twin are highlighted with a darker color. Relevant genes are labeled. e Box plots representing selected dysregulated genes in the CVID cohort (n = 10) compared to healthy donors (n = 8). f Heatmap representing DEGs from the comparison of healthy controls and CVID Ib patients. Genes also dysregulated in the discordant twins are labeled in red. g Dot plot representing selected differentially expressed proteins obtained from CITE-seq analysis (gene name in brackets). Scale indicates the log2(FC) protein expression of the activated memory B-cell subset in the CVID or CVID Ib vs. control comparison, logFDR is indicated by the circle size. N.S means not significant. h UMAP visualization showing different immune cell populations identified (95,064 cells). i Overview of selected dysregulated L/R interactions. Scale indicates the log2(FC) gene expression of memory activated B cells in the CVID or CVID Ib vs. control comparison. The percentage of other immune cells expressing the partner molecule is indicated by the circle size. Molecules of the L/R pairs expressed in activated memory B cells are shown in blue; molecules of the L/R pairs expressed in the immune cell partner appear in red. Validated gene expression at protein level is indicated with an asterisk. Source data for panels in this figure are provided in the Source Data Fig. 6 file. CVID-associated epigenetic alterations are not reflected in the alterations taking place in other immune subsets distinct to B transcriptomes of the few steady-state S-mem B cells that are cells, together with the specific dysregulation in the B-cell com- present in peripheral blood and which are virtually identical to partment upon indirect B-cell activation through T-cell stimula- those in the healthy controls. Instead, B cells display transcrip- tion, would constitute B-cell-extrinsic defects that also might tional differences upon activation. Since several donors included contribute to the failure of B-cell activation in CVID patients. in this study are healthy individuals, any invasive approach One of the major constraints of our study consists in the limited including a lymph node biopsy would not be ethically justifiable size of the CVID patients’ cohort. In this sense, we mainly focused on for research purposes. Given these constraints, we used in vitro- CVID Ib patients, being the other two subgroups (CVID Ia and activated B cells in order to partially emulate the GC reaction. CVID II) barely represented. Moreover, our analyses are based on The GC reaction is a tightly regulated process in which the bal- immune cells isolated from peripheral blood, therefore CVID- ance between survival and apoptotic signals is critical to promote associated immune features taking place in secondary lymphoid antigen affinity while simultaneously avoiding self-antigen organs or in the bone marrow still need to be determined. In addition recognition. In this context, we identified significant down- to the cohort of CVID patients and healthy controls used, one of the regulation of several chemokine genes in CVID B cells such as strong points of this study is the use of monozygotic twins, which are CCL22 and CCL17, which have been described as fine-tuning genetically identical. It is likely that the exposure to different envir- antibody affinity maturation and positive selection in GC pro- onments at prenatal or post-delivery time, or even their differential moting B-cell interactions with follicular helper T cells48. Other early infection history, may have induced epigenetic alterations that relevant downregulated genes in the B cells of CVID patients are ultimately lead to a difference in gene expression and twin dis- those encoding CD70, which has been involved in B-cell cordance in CVID. The reduction of both memory B cell and PC activation49,50, and ICAM1 (CD54), whose interaction with the counts in CVID patients, as well as the severe hypogammaglobuli- integrin LFA-1 lowers the threshold for B-cell activation by nemia in one of the twins could be a consequence of the defects in promoting B-cell adhesion and immunological synapse the proper activation of naïve B cells, alterations during memory and formation51. Additionally, our data indicate that genes encoding PC generation, or even both. However, our results indicate that naïve Gal-9, CD45, and CD22 are upregulated in the B-cell compart- B-cell populations are almost identical in both twins in terms of their ment of the CVID twin. Gal-9 hampers BCR activation through DNA methylation, chromatin accessibility and transcription profile, its interaction with the N-glycan repertoire of CD45 molecules, indicating that CVID epigenetic and gene expression defects are which ultimately inhibits BCR signaling via CD2252. Another established later on during naïve B-cell activation and memory B-cell study showed that Gal-9 facilitates the interaction of IgM-BCR generation at secondary lymphoid organs, where the establishment of with the inhibitory molecules CD45 and CD2253. Additionally, appropriate cell–cell communication is crucial for mounting efficient we observed significant upregulation of the CD72 gene in CVID B immune responses. cells, which is another molecule that induces cell-cycle arrest and apoptosis in mature B cells and that is involved in the inhibition Methods of BCR downstream signaling pathways54,55. The aberrant Patients and ethics approval. Human blood samples used in this study were upregulation of these inhibitory molecules might interfere with obtained from a pair of monozygotic twins discordant for CVID (the same pair the appropriate activation of the BCR, and might be one of the participating in a previous study10) and an additional cohort of ten CVID patients and ten healthy donors (Supplementary Data 1). CVID patients were diagnosed causes of the deficient immune responses and the defective B-cell according to European Society for Immunodeficiencies (ESID) criteria59. They counts that characterize CVID patients. were collected at the University Hospital Dr Negrín, Gran Canaria, Spain and at the It is of particular note that we also detected transcriptional Hospital La Paz, Madrid, Spain. All donors received oral and written information defects beyond the B-cell compartment. In this regard, we found about the possibility that their blood would be used for research purposes, and any questions that arose were then answered. Before giving their first blood sample the several dysregulated genes in the activated CD4+ T cells of the donors signed a consent form approved by the Ethics Committee at their corre- CVID patients. For instance, we found downregulation of genes sponding hospital (Hospital La Paz PI-2833), which adhered to the principles set such as SEMA4D (CD100), which is an interacting partner of out in the WMA Declaration of Helsinki. The protocol used to isolate B cells from CD72, as well as upregulation of the B-cell inhibitory molecule these donors was approved by the Ethics Committee of the Bellvitge University Hospital (CEIC) on 9 March 2017 (PR053/17). PTPRC (CD45). Additionally, we found a significant down- regulation of CD40LG (CD40L) in this T-cell compartment as Sample collection and immune cell activation. PBMCs were obtained from previously described by others56,57, although the dysregulation of peripheral blood by Ficoll gradient using Lymphoprep (Stemcell Technologies, Cat. this molecule is still a contentious topic in the field58. All these No. # 07801). For the isolation of CD19+CD27negIgD+ naïve, CD19+CD27+IgD+ 12 NATURE COMMUNICATIONS | (2022)13:1779 | https://doi.org/10.1038/s41467-022-29450-x | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-022-29450-x ARTICLE unswitched memory (US-mem) and CD19+CD27+IgDneg switched memory (S- KREC assay. The replication history of naïve, US-mem and S-mem B cells from mem) B cells, PBMCs were stained with 2 μL/million cells of anti-CD19-FITC (BD CVID MZ discordant twins, was inferred from isolated DNAs by the κ − deleting Biosciences, clone: 4G7, Cat. No. 15856028), 2 μL/million cells of anti-CD27-APC recombination excision circle (KREC) assay, as previously described24. In brief, the (Miltenyi Biotec, clone: M-T271, Cat. No. 130-113-630) and 2 μL/100 million cells replication history is estimated by the ratio between genomic coding joints and of anti-IgD-PE (Southern Biotech, Birmingham, AL, USA, clone: NA, Ref: 2032-09) signal joints on KRECs of the Igκ intron RSS-K deleting rearrangement, which are in MACS buffer (PBS + 2% FBS + 2 mM EDTA). quantified by Real-Time PCR, applying the formula: ΔCT sample = Ct (signal For immune cell activation, PBMCs were resuspended in complete medium joint) − Ct (coding joint) (formula 1). (RPMI Medium 1640 + GlutaMAX (Gibco, Life Technologies, Cat. No. 61870036) containing 20% fetal bovine serum (Gibco, Life Technologies, Cat.No. 10270-106), 100 units/mL penicillin, and 100 μg/mL streptomycin). PBMCs were then activated Single-cell library preparations. For single-cell bisulfite sequencing (scBS-seq) with a B-cell stimulus consisting of a combination of 0.1 µg/mL of CD40L assays, naïve, US-mem and S-mem B cells were collected by sorting in 12 μl of lysis (MEGACD40L Protein, ENZO, Cat. No. ALX-522-110-C010) and 50 ng/mL of IL- buffer (10 mM Tris-Cl pH 7.4, 0.6% SDS, 0.5 μl proteinase K) using a Beckman 21 (Tebu-bio, Cat. No. 200-21), or with a T-cell stimulus consisting of α-CD3/ Coulter MoFlo Astrios EQ Cell Sorter instrument in single-cell 1 drop mode. CD28 Dynabeads 20 μL of beads/million PBMCs (Invitrogen, Cat. No. 11131D) for ToPro-3 and Hoechst 33342 staining were used to select for live cells with low 48 h. As a control, some PBMCs were also cultured with no stimulus. After DNA content (i.e., in G0/G1). Negative controls were sorted using BD Accudrop stimulation, PBMCs were harvested and a fraction of the cells were stained with Beads on lysis buffer, and were prepared and processed concomitantly with all 2 μL/million of cells of anti-CD19-FITC for B-cell sorting in Beckman Coulter single-cell samples. Libraries for scBS-seq were generated as previously described15. MoFlo Astrios EQ Cell Sorter instrument. Finally, PBMCs were mixed with sorted Briefly, bisulfite conversion was performed on cell lysates using the Imprint DNA CD19+ B cells in a 2:1 ratio in order to obtain PBMCs samples enriched in B cells. Modification Kit (Sigma, Cat. No. MOD50-1KT) with the following modifications: This B-cell enrichment step was performed in order to avoid T cells and other all volumes were halved, and chemical denaturation was followed by incubation at immune cell types mask the B-cell compartment due to the drastically reduced 65 °C for 90 min, 95 °C for 3 min and 65 °C for 20 min. DNA was purified and levels of B cells observed in CVID patients. eluted in 10 mM Tris-Cl (pH 8.5) and combined with 0.4 mM dNTPs, 0.4 μM For the analysis of CD40L expression in CD4+ T cells, PBMCs were isolated oligo1 ([Btn]CTACACGACGCTCTTCCGATCTNNNNNNNNN) and 1× Blue from the CVID-discordant twins and cultured at a concentration of 200,000 cells/ Buffer (Sigma) (24 μl final) before incubation at 65 °C for 3 min followed by 4 °C 200 μL of complete medium in a 96-well plate. Cells were stimulated with α-CD3/ pause. 50U of Klenow exo- (Sigma) were added and the samples incubated at 4 °C CD28 Dynabeads 20 μL or 4 μL of beads/million of PBMCs for 48 h. As a control, for 5 min, +1 °C/15 s to 37 °C, 37 °C for 30 min. Samples were incubated at 95 °C some PBMCs were also cultured with no stimulus. After stimulation, PBMCs were for 1 min and transferred immediately to ice before addition of fresh oligo1 harvested and stained with 2 μL/million cells of anti-CD4-APC (BD Biosciences, (10pmol), Klenow exo- (25U), and dNTPs (1nmol) in 2.5 μl total. The samples clone: RPA-T4. Cat. No. 555349) and 2 μL/million cells of anti-CD40L-PE (BD were incubated at 4 °C for 5 min, +1 °C/15 s to 37 °C, 37 °C for 30 min. This Biosciences, clone: TRAP1, Cat. No. 555700) in MACS buffer. Cells were acquired random priming and extension was repeated a further 3 times (5 rounds in total). using Gallios Flow Cytometer (Beckman Coulter, CA, USA) and analyzed using the Samples were then incubated with 40U exonuclease I (NEB) for 1 h at 37 °C before Flowjo software. DNA was purified using 0.8× Agencourt Ampure XP beads (Beckman Coulter, Cat. No. A63881) according to the manufacturer’s guidelines. Samples were eluted in 10 mM Tris-Cl (pH 8.5) and incubated with washed M-280 Streptavidin Dyna- Whole-genome sequencing analysis. WGS of genomic DNA obtained from the beads (ThemoFisher Scientific, Cat. No. 11205D) for 20 min with rotation at room pair of monozygotic twins discordant for CVID. The quality of the isolated temperature. Beads were washed twice with 0.1 N NaOH, and twice with 10 mM genomic DNA was checked on 1% agarose gel. Besides, DNA concentration was Tris-Cl (pH 8.5) and resuspended in 48 μl of 0.4 mM dNTPs, 0.4 μM oligo2 measured using Qubit DNA Assay Kit (Cat. No. 10146592) in Qubit 2.0 Fluo- (TGCTGAACCGCTCTTCCGATCTNNNNNNNNN) and 1× Blue Buffer. Sam- rometer (Life Technologies, CA, USA). The genomic DNA was randomly frag- ples were incubated at 95 °C for 45 s and transferred immediately to ice before mented by sonication to the size of 350 bp, then DNA fragments were end addition of 100U Klenow exo- (Sigma, Cat. No. KEM0026) and incubation at 4 °C polished, A-tailed, and ligated with the full-length adapters of Illumina sequencing, for 5 min, +1 °C/15 s to 37 °C, 37 °C for 90 min. Beads were washed with 10 mM and followed by further PCR amplification. The PCR products as the final con- Tris-Cl (pH 8.5) and resuspended in 50 μl of 0.4 mM dNTPs, 0.4 μM PE1.0 forward struction of the libraries were purified with the AMPure XP system. Then libraries primer (AATGATACGGCGACCACCGAGATCTACACTCTTTC-CCTA- were checked for size distribution using Agilent 2100 Bioanalyzer (Agilent Tech- CACGACGCTCTTCCGATCT), 0.4 μM indexed iPCRTag reverse primer20, 1U nologies, CA, USA), and quantified by real-time PCR (to meet the criteria of 3 nM). KAPA HiFi HotStart DNA Polymerase (KAPA Biosystems, Cat. No. KAPA The 150 × 2 paired-end sequencing was conducted on a Novaseq sequencer KK2502) in 1× HiFi Fidelity Buffer. Libraries were then amplified by PCR as (Illumina, CA, USA) producing 612,180,910 raw paired reads on average. follows: 95 °C 2 min, 12–13 repeats of (94 °C 80 s, 65 °C 30 s, 72 °C 30 s), 72 °C Quality check (QC) and read alignment were performed using miARma-seq 3 min, 4 °C hold. Amplified libraries were purified using 0.8× Agencourt Ampure pipeline60. This software uses FastQC (version 0.11.5)61 to detect adapter XP beads, according to the manufacturer’s guidelines, and were assessed for quality contamination and low quality reads. Reads with a quality lower than 30 were and quantity using High-Sensitivity DNA chips on the Agilent Bioanalyser, and the stringently eliminated to reduce false positives from the sequencing output. After KAPA Library Quantification Kit for Illumina (KAPA Biosystems, Cat. No. QC, the sequence data were subjected to BWA (version 0.7.13)62 for alignment to KK4824). Pools of 12–14 single-cell libraries were prepared for 100 bp paired-end the human reference genome (GRCh38p12). The alignment data (BAM files) were sequencing on a HiSeq2500 in rapid-run mode (2 lanes/run). processed by GATK (version 4.1.7.0)63 according to their best practices workflows. The plate-based single-cell ATAC-seq assay was performed as previously In that sense markduplicates, base recalibration, covariate analysis and variant described68. Briefly, naïve, US-mem and S-mem B cells were collected by FACS in filtration (for SNPs: Variant Confidence/Quality by Depth (QD) < 2.0, Phred- MACS buffer, and cells were centrifuged at 500 × g at 4 °C for 5 min. Cell pellets scaled p value using Fisher’s exact test to detect strand bias (FS) > 60.0, RMS were resuspended in 50 μl tagmentation mix (33 mM Tris-acetate, pH 7.8, 66 mM Mapping Quality (RMQ) < 40.0, Symmetric Odds Ratio of 2 × 2 contingency table potassium acetate, 10 mM magnesium acetate, 16% dimethylformamide, 0.01% to detect strand bias (SOR) > 4.0, Z-score From Wilcoxon rank sum test of Alt vs. digitonin and 5 μl of Tn5 from the Nextera kit from Illumina, Cat. No. FC-121- Ref read mapping qualities (MQRankSum) < −12.5 and Z-score from Wilcoxon 1030). The tagmentation reaction was done on a thermomixer at 800 rpm, 37 °C rank sum test of Alt vs. Ref read position bia (ReadPosRankSum) < −8.0. For for 30 min. The reaction was then stopped by adding equal volume (50 μl) of indels: QD < 2.0, FS > 200.0, SOR > 10.0) was performed to optimize the final depth tagmentation stop buffer (10 mM Tris-HCl, pH 8.0, 20 mM EDTA, pH 8.0) and left and base-call quality for variant detection. on ice for 10 min. A volume of 200 μl MACS buffer with 0.5% BSA was added and The set of common variants between both individuals was obtained using the the nuclei suspension was transferred to a FACS tube. DAPI was added at a final intersectBed tool from the bedtools suite64. Furthermore, putative somatic concentration of 1 μg/μl to stain the nuclei. DAPI-positive single nuclei were sorted mutations were detected with the MuTect2 software by GATK (version 4.1.7.0) by into 96-well full-skirted Eppendorf plates chilled to 4 °C, prepared with 2 µl 2X lysis comparing BAM files from the twin and the co-twin. A base recalibrator process buffer (100 mM Tris.HCl, pH 8.0, 100 mM NaCl, 40 µg/ml Proteinase K (Ambion, was performed providing known sites of variation such as the allele frequency by AM2546, 20 mg/ml stock) and 0.4% SDS) and 2 µl of 10 µM S5xx/N7xx Nextera the Genome Aggregation Database (gnomAD)65, high confidence SNPs and indels Index Primer Mix (5 µM each). Plates were immediately sealed, spun down at from the 1000 Genomes project66 and the NCBI database of genetic variation 300 × g at 4 °C for 1 min and stored at −80 °C. (dbSNP)67, to label sequencing errors and distinguish poor base quality from For single-cell RNA-seq using the plate-based Smart-seq2 (SS2) protocol, naïve, variants. These fine-tuned parameters allowed to optimize the final depth and base- S-mem and S-mem B single cells were collected in 96-well full-skirted Eppendorf call quality for somatic mutation Identification. Then, the FilterMutectCalls tool plates chilled to 4 °C, prepared with lysis buffer consisting in 10 μl of TCL buffer from the GATK toolbox was applied to use allele fractions to distinguish somatic (Qiagen, Cat. No. 1070498) supplemented with 1% β-mercaptoethanol. Plates were variants from sequencing errors. sealed, vortexed, spun down at 300 × g at 4 °C for 1 min, immediately placed on dry The relationship of identified variants with phenotypes was performed with the ice and transferred for storage at −80 °C. The Smart-seq2 protocol was performed intersectBed utility using the data provided by the ClinVar database19. The impact basically as previously described69. Libraries were sequenced, aiming for an average of the sequence modification by these variants at the functional level was annotated depth of 1 million reads per cell, on an Illumina HiSeq 2000 with version 4 using SNPEff20. chemistry (paired-end, 75-bp reads). WGS results are not publicly-available due to ethical regulations, but are For the droplet scRNA-seq method, stimulated and unstimulated PBMCs available from the corresponding authors upon request. enriched with CD19+ B cells were counted using a Neubauer hemocytometer, and NATURE COMMUNICATIONS | (2022)13:1779 | https://doi.org/10.1038/s41467-022-29450-x | www.nature.com/naturecommunications 13 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-022-29450-x 5,000 cells per sample were loaded in the 10x-Genomics Chromium. The 10x- v4.10.3. Specifically, we used the findMotifsGenome.pl algorithm (with parameters Genomics 5’ libraries were prepared following the manufacturer’s instructions. -size given -cpg) to search for significant enrichment compared with background, Libraries were sequenced, aiming at a minimum coverage of 50,000 raw reads per adjusted to ensure an equal number of sequences with equal lengths and with equal cell, on an Illumina HiSeq 4000 (paired-end; read 1: 26 cycles; i7 index: 8 cycles, i5 CpG and GC content. index: 0 cycles; read 2: 98 cycles). For samples of the expanded cohort of CVID DMR annotation for genetic context location was performed using the patients and healthy controls, we processed PBMC samples for scRNA-seq with annotatePeaks.pl algorithm in the HOMER software v4.10.3. paired measurements of 192 surface proteins (CITE-seq, using TotalSeq reagents Replication timing data in the GM12878 and GM12801 lymphoblastoid cell from Biolegend). The list of the antibodies included in the TotalSeq custom panel lines were obtained from the UW Repli-seq track of the UCSC Genome Browser. are detailed in Supplementary Data 12. Briefly, 1 million cells pooled from several Replication timing values were binned in deciles to perform the overlap with the individuals were stained with TotalSeq-C custom human panel together with CVID.no-demeth DMRs. Human TruStain FcX (Biolegend, Cat. No. 422301) on ice for 30 min. Then cells GREAT software74 was used to enrich downstream pathways and gene were washed three washes with cold PBS + 4% FBS. After completing the last ontologies. We used the single nearest gene option for the association between washing steps, cells were resuspended in PBS + 0.04% BSA and filtered with a 40 genomic regions with genes. um Bel-art Flowmi strainer into a new 1.5 ml lo-bind eppendorf tube. Cells were To define the genomic coordinates corresponding to active, primed or repressed counted and concentration adjusted to load 50,000 cells on the 10X-Genomics enhancers, we used the following command via bedtools (version v2.28.0): Chromium Station. The 10x-Genomics 5’ libraries together with CITE-seq libraries “bedtools intersect file_A.bed file_B.bed”, in order to identify the regions with were prepared following the manufacturer’s instructions. specific histone marks: For active enhancers file_A.bed corresponded to H3K4me1 peak file of the cell type of interest, and file_B.bed to H3K27ac. For the repressed enhancers file_A.bed corresponded to H3K4me1 peak file of the cell type of Single-cell whole-genome bisulfite sequencing. For DNA methylation data, interest, and file_B.bed to H3K27me3. For defining primed enhancers, H3K4me1 single-cell bisulfite sequencing (scBS-seq) data was processed as previously bed file of the cell type of interest was used. All of the above was done separately for described70. Reads were trimmed with TrimGalore! (https://www.bioinformatics naïve, S-mem, GC B cells, as well as for plasma cells. .babraham.ac.uk/projects/trim_galore/), Cutadapt71 and FastQC (https://www. Inference of TF activities from expression values were calculated using bioinformatics.babraham.ac.uk/projects/fastqc/), using default settings for DNA DoRothEA75. The B_viperRegulon.rdata dataset was loaded in R and its Integrated methylation data and additionally removing the first 6 bp. Subsequently, Bismark72 Development Environment (RStudio) and used for the analysis. (v0.16.3) was used to map the bisulfite data to the human reference genome (GRCh38), in single-end non-directional mode, followed by de-duplication and DNA methylation calling using default settings. We removed cells with a library size of Alignment, quantification, and quality control of scRNA-seq data. Smart-seq2 fewer than 3 M reads, resulting in 196 cells with DNA methylation information. (SS2) sequencing data were aligned with STAR76 (version 2.5.1b) using the Specifically, there were 55 S-mem (23 CVID), 69 US-mem (31 CVID) and 72 naïve GRCh38 (Ensembl release 84) STAR index and annotation. Gene-specific read cells (37 CVID) (Supplementary Data 2). counts were calculated using HTSeq-count (version 0.10.0). Cells with fewer than 200 detected genes and more than 20% mitochondrial gene expression content Imputation of DNA methylation signals. To mitigate the incomplete coverage of were removed. Genes expressed in fewer than 3 cells were removed. scBS-seq profiles, [healthy twin → naïve BCs = 27.0% ± 3.5%, US-mem BCs = Droplet-based RNA-sequencing data for B-cell receptor and T-cell receptor 23.8% ± 5.9%, S-mem BCs = 21.8% ± 4.0%//CVID twin → naïve BCs = 23.0% ± (what we later on refer to as “BCR data” and “TCR data”) were aligned with Cell- 2.9%, US-mem BCs = 25.7% ± 8.5%, S-mem BCs = 24.2% ± 7.9%] (see Supple- Ranger VDJ SingleCell Software Suite (version 3.0.2, 10x-Genomics Inc.) using the mentary Data 2 for coverage percentages per single-cell), we applied DeepCpG16 to GRCh38 official Cell-Ranger V(D)J reference, version 2.0.0. Downstream analysis impute unobserved methylation states of individual CpG sites. DNA methylation was performed using scirpy (https://icbi-lab.github.io/scirpy/). profiles of the three cell types were imputed separately; the two donors were jointly Droplet-based RNA-sequencing data was aligned and quantified using the Cell- imputed to mitigate potential imputation artifacts between the samples. The cell Ranger SingleCell Software Suite (version 3.0.2, 10× Genomics Inc.) using the type-specific models were developed using CpG and genomic information GRCh38 human reference genome (official Cell-Ranger reference, version 1.2.0). according to DeepCpG’s setup of a joint model (see DeepCpG16 for details and Cells with fewer 200 detected genes or more than 20% mitochondrial gene default values). The accuracy of the models per cell, as measured using AUC, expression content were removed. Genes expressed in fewer than 3 cells were ranged between 0.78 and 0.93, see Supplementary Data 2 for imputation accuracy removed. For the SS2 and droplet-based datasets, the Scanpy (version 1.4.4) Python per sample. Of note, the DNA methylation patterns after imputation recapitulates package77 was used to load the cell-gene count matrix and perform downstream better the methylation profiles obtained from WGBS than those before imputation analysis. (when compared to data from the Blueprint consortium; EGAD00001002361 and To remove cell-cycle-associated variation, cell cycle-associated genes were EGAD00001002354), indicating that imputation is an appropriate approach to removed from SS2 and droplet-based datasets. To do so, we normalized mitigate the technical noise associated with scBS. (scanpy.pp.normalize_per_cell method, scaling factor 10000), log-transformed (scanpy.pp.log1p) and subsetted the highly variable gene (scanpy.pp.filter_gene_dispersion) count matrices of both datasets. Next, datasets DMR identification. After imputation we merged the data per cell type - donor were transposed to operate in gene space with cells as features (gene-centered combination to identify differentially methylated regions (DMRs), using Metilene21 analysis). Data-feature scaling (scanpy.pp.scale), PCA analysis (scanpy.pp.pca, from (v0.2.5) in de novo mode. We selected only DMRs with values of q < 0.05 and a variable genes), neighborhood graph building (scanpy.pp.neighbors) and Louvain methylation difference of >20% for subsequent analysis. graph-based clustering (scanpy.tl.louvain, clustering resolution manually tuned) To identify altered DMRs in the CVID twin, we analyzed the previously selected were performed. Previously known cell-cycle genes (CDK1, MKI67, CCNB2, DRMs using the limma package73, comparing the methylation values of those PCNA) were used, and all genes that cluster together with any of these were DMRs in single cells (FDR < 0.05 and methylation difference >10%). removed. To define altered DMRs at the single-cell level in CVID, we first calculated the methylation distributions of each DMR in healthy control S-mem B cells. We then considered as altered DMRs those DMRs in CVID S-mem B cells whose mean Deconvolution of donors in pooled CITE-seq samples. Pooled CITE-seq samples methylation was more than 1.5*IQR above the third quartile or below the first containing cells from multiple individuals were demultiplexed using souporcell78. quartile of the control DMR distributions. Briefly, the algorithm identifies genotypic differences between single cells by variant calling aligned reads using STAR76 and generating a VCF (Variant Call Format) file using Freebayes (Marth 2012). Souporcell was run with the following command for Amplicon sequencing for DNA methylation analysis. B-cell subsets (naïve, US- all samples (argument $1 corresponds to sample ID, and N is the number of mem and S-mem B cells) from 10 healthy donors and 11 CVID patients were multiplexed individuals in a sample): /software/singularity-v3.5.1/bin/singularity sorted as described in the Sample collection section. Genomic DNA was then exec./souporcell.sif./souporcell_pipeline.py -i./cellranger302_count_$1_GRCh38- isolated from each population using QIAamp DNA micro kit (QIAGEN, Cat. No: 1_2_0/possorted_genome_bam.bam -b./cellranger302_count_$1_GRCh38-1_2_0/fil- 56304) and was bisulfite converted using MethylCode Bisulfite Conversion Kit tered_feature_bc_matrix/barcodes.tsv -f./refdata-cellranger-GRCh38-1.2.0/fasta/ (Thermofisher, Cat. No. MECOV50) following manufacturer instructions. From 1 genome.fa -t 8 -o souporcell_result_$1 -k N --skip_remap True --common_variants./ to 10 ng of bisulfite converted DNA was used for amplicon generation using a filtered_2p_1kgenomes_GRCh38.vcf, where the last VCF file with common variants custom panel, as well as for library preparation. The list of analyzed CpG sites was downloaded as per instructions in https://github.com/wheaton5/souporcell (1058 CpG sites corresponding to 162 selected CVID.no-demeth DMRs) is available with the following command: wget --load-cookies /tmp/cookies.txt “https:// in Supplementary Data 9. Libraries were sequenced using IonTorrent technology docs.google.com/uc?export=download&confirm=$ (wget-quiet-save-cookies/tmp/ (Thermofisher) according to manufacturer instructions. For methylation measure, cookies.txt-keep-session-cookies-no-check-certificate ‘ https ://docs.google.com/uc? only those CpG sites with a coverage > 50 (methylated + un-methylated reads) export=download&id=13aebUpEKrtjliyT9rYzRijtkNJVUk5F_’-O-|sed-rn’s/. were considered. *confirm=([0-9A-Za-z_]+).*/\1\n/p’)&id=13aebUpEKrtjliyT9rYzRijtk NJVUk5F_”-Ocommon_variants_grch38.vcf && rm-rf/tmp/cookies.txt. Analysis of the different-omics datasets in relation with several biological Furthermore, to find out the exact donor identity of each donor’s barcode features. Transcription factor motifs were enriched for each set of identified cluster in souporcell results, the cells in parallel have undergone genotyping using CVID.hypo.DMRs in pseudo-bulk or in each single-cell using HOMER software Illumina Infinium Global Screening Array. In order to unambiguously identify 14 NATURE COMMUNICATIONS | (2022)13:1779 | https://doi.org/10.1038/s41467-022-29450-x | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-022-29450-x ARTICLE every individual in the pooled samples, each donor’s variants were separated from Alignment, quantification, and downstream analysis of single-cell ATAC-seq the pooled VCF and each single-donor VCF was matched to the genotype data data. Plate-based single-cell ATAC sequencing data were aligned and quantified as using PLINK79. This software matches each souporcell sample with the genotype previously described68 and further analyzed with the standard scanpy pipeline described data giving a concordance ratio (based on the similarity of the variants) that allows above (excluding cell-cycle genes removal and doublet detection). Further downstream us to distinctly identify each sample with each donor ID. analysis was performed using cisTopic24 approach following the basic tutorial (https:// rawcdn.githack.com/aertslab/cisTopic/f628c6f60918511ba0fa4a85366ebf52db5940f7/ vignettes/CompleteAnalysis.html). Doublet removal from scRNA-seq and CITE-seq data. To exclude doublets from Transcription factor motif-enrichment analysis was performed for each single-cell RNA and CITE sequencing data, a previously described two-step identified topic using HOMER software v4.10.3. We used the findMotifsGenome.pl approach80 was used. Firstly, we applied the Scrublet81 algorithm for each sample algorithm (with parameter -size given) to search for significant enrichment relative to calculate the scrublet-predicted doublet score per cell. Then, to ensure that small to the background. clusters with high doublet density were not grouped with large numbers of singlets, Differentially accessible region (DAR) calculation was performed in a scRNA- we overclustered the manifold (running a basic scanpy pipeline up to clustering, seq like manner using limma73 approach as described above. We used.bed files additionally clustering each cluster separately), computed per-cluster Scrublet corresponding to significant DARs and.bed files corresponding to TF and enhancer scores as the median of the observed values, and computed normal distribution Chip-seq data (same publicly-available data as used for Fig. 3e, d) to first calculate parameters, centered on the median of the scores, with a MAD-derived standard overlaps between them and then estimate the enrichment of target (TF or deviation (the score distribution is zero-truncated, so we used only above-median enhancer) peaks in DARs as the ratio of total overlap length to length of values to compute the MAD). We then multiplied the MAD value by 1.4826 to target peaks. obtain a literature-derived normal distribution standard deviation estimate and FDR-corrected the p values using the Benjamini–Hochberg method62. Secondly, Statistical analysis. Statistical analyses were performed in R (v4.0.3). Data dis- once the dataset had been processed and clustered with the Louvain algorithm as tributions were tested for normality. Normally distributed data were tested using described above, step 1 was repeated for each of the clusters identified, thus cap- two-tailed unpaired Student’s t tests; non-normal data were analyzed with the turing doublets within each cluster that were not detected with the initial doublet appropriate non-parametric statistical test. calling. Cells called as doublets in this last step were flagged as “Doublets”. Fur- DNA methylation correlations between different B-cell subsets were expressed thermore, we identified physical doublets of T and B cells using BCR data as cells in as Pearson correlation coefficients. The statistical significance of potential DMRs T-cell clusters that had BCR data and, conversely, using TCR data as cells in B-cell was assessed by a two-dimensional version of the Kolmogorov–Smirnov test to clusters that had TCR data. These cells were flagged as “B:T_doublets_by_TCR” or calculate P values during segmentation and to use it as a termination criterion “T:B_doublets_by_BCR”. In the zoomed in reanalysis of B-cell populations of the during recursive segmentations. The output of metilene provides the adjusted and twins we also labeled a cluster of cells expressing CD3E as “T_contaminants”. All unadjusted P value for the 2D-KS test and additionally the results of an the three types of detected doublet populations were removed from subsequent independent Mann–Whitney U test. downstream analysis. In addition, for pooled samples, any cells that were called as Transcription factor motif enrichments were calculated using the cumulative inter-individual doublets by souporcell were also removed from subsequent hypergeometric distributions (input vs. background sequences). Significance of downstream analysis. associations between topics from cisTopic and genomic features (enhancers, TF motifs or TF peaks) were calculated with a Mann–Whitney test. Significance of B cells enrichments in scATAC-seq clusters were calculated using a Fisher’s test. Denoising protein counts in CITE-seq data. To remove ambient signal from Differentially expressed genes were calculated using a two-sided limma protein counts of the CITE-seq data, we used SoupX[https://academic.oup.com/ moderated t-test with Benjamini-Horchberg correction for multiple testing. gigascience/article/9/12/giaa151/6049831] in the following way: first, we got the expression profile of the ambient signal from empty droplets (estimateSoup with soupRange = c (4, 100)), then we automatically calculated contamination fraction Reporting summary. Further information on research design is available in the Nature (autoEstCont with soupQuantile = 0.1 and tfidfMin = 0.05) and removed back- Research Reporting Summary linked to this article. ground contamination from the count matrix (adjustCounts). The corrected pro- tein counts were then merged with the gene expression part of CITE-seq data and further treated as raw counts and analysed. Protein features in the CITE-seq data Data availability The scBS-seq data generated in this study is available under ENA study accession were considered equal to gene features in subsequent downstream analysis. PRJEB50820. Smart-seq2 (SS2) sequencing data, scATAC-seq data and 10x-Genomics Chromium droplet-based scRNA-seq data can be accessed under EGA dataset accession EGAD00001008575. Interactive visualizations of single-cell transcriptomic datasets from Clustering and annotation of scRNA-seq data and CITE-seq data. Downstream this study can be accessed via https://www.immunodeficiencycellatlas.org/. WGS results analysis included data normalization (scanpy.pp.normalize_per_cell method, are not publicly available due to ethical regulations, but are available from the scaling factor 10000), log-transformation (scanpy.pp.log1p), variable gene detec- tion (scanpy.pp.filter_gene_dispersion), data-feature scaling (scanpy.pp.scale), PCA corresponding authors upon request. RNA-seq and ChIP-seq publicly-available datasets analysis (scanpy.pp.pca, from variable genes), batch-balanced neighborhood graph for naïve-, GC, S-mem B cells, as well as for PCs from healthy donors were obtained from building (scanpy.pp.bbknn) by experimental sample batch key ‘sample’ and Lou- the Blueprint Consortium (http://dcc.blueprint-epigenome.eu/#/datasets), with the vain graph-based clustering (scanpy.tl.louvain, clustering resolution manually following accession numbers: Naïve BCs RNA-seq (EGAD00001002315), GC BCs RNA- tuned). Cluster cell identity was assigned by manual annotation using known seq (EGAD00001002452), S-mem BCs RNA-seq (EGAD00001002476), Plasma cells marker genes (Fig. 5b). RNA-seq (EGAD00001002323), Naïve BCs ChIP-seq (EGAD00001002466), GC BCs While annotating B-cell populations, we perform more resolved clustering of ChIP-seq (EGAD00001002442), S-mem BCs ChIP-seq (EGAD00001002430), Plasma the coarse “B cells” cluster (in the main manifold with all cells) to tease apart cells ChIP-seq (EGAD00001002281). We used additional ChIP-seq publicly-available subtler B-cell states (represented in Fig. 4c left panel) via (1) examining genes datasets for several transcription factors with the following accession numbers: CTCF defining those states (as shown in Fig. 4d) and (2) examining logistic regression from human B cells (GSM1003474), as well as JUNB (GSE96455), BATF (GSM803538), predictions from (a) unstimulated B-cell scRNA-seq of twins (used in the first part CTCF (ENCSR184YZV) and PAX8 (GSE127505) from the GM12878 human of study and visualized in Fig. 4a) and (b) dataset from King et. al., Science lymphoblastoid cell line. Genome Reference Consortium Human Build 38 patch release Immunology, 2021 to better tell apart naïve and memory phenotype. 12 (GRCh38p12). Source data are provided with this paper. Differential gene expression analysis in scRNA-seq. To calculate differentially Code availability expressed genes in each comparison cell state between CVID and control twin we All the code used in the analyses can be found at https://github.com/ventolab/CVID83. adopted the limma73 approach (limma version 3.46.0, edgeR version 3.32.1). Received: 29 July 2021; Accepted: 8 March 2022; Cell–cell communication analysis. Cell–cell interactions were inferred using CellPhoneDB (www.CellPhoneDB.org), as previously described33,82. We modified the statistical method to calculate significant L/R interactions enriched in disease. To do this, instead of random shuffling of cells to select specific interactions, we used genes differentially expressed between the twins in the different identified immune cell subsets (FDR < 0.05) and only included an L/R pair in the analysis References when both molecules were present in at least 10% of the interacting cell clusters 1. Bousfiha, A. et al. Human inborn errors of immunity: 2019 update of the IUIS and at least one of the molecules was DE between CVID and healthy controls. phenotypical classification. J. Clin. Immunol. 40, 66–81 (2020). Additionally, since CVID is considered a B-cell-centered disease, we focused our 2. Bonilla, F. A. et al. International Consensus Document (ICON): common analysis on the potential interactions between B cells and the other immune cell variable immunodeficiency disorders. J. Allergy Clin. Immunol. Pract. 4, 38–59 compartments. (2016). NATURE COMMUNICATIONS | (2022)13:1779 | https://doi.org/10.1038/s41467-022-29450-x | www.nature.com/naturecommunications 15 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-022-29450-x 3. Blanco, E. et al. Defects in memory B-cell and plasma cell subsets expressing 32. Zheng, G. X. Y. et al. Massively parallel digital transcriptional profiling of different immunoglobulin-subclasses in patients with CVID and single cells. Nat. Commun. 8, 14049 (2017). immunoglobulin subclass deficiencies. J. Allergy Clin. Immunol. 144, 809–824 33. Efremova, M., Vento-Tormo, M., Teichmann, S. A. & Vento-Tormo, R. (2019). CellPhoneDB: inferring cell-cell communication from combined expression of 4. Ramesh, M., Hamm, D., Simchoni, N. & Cunningham-Rundles, C. Clonal and multi-subunit ligand-receptor complexes. Nat. Protoc. 15, 1484–1506 (2020). constricted T cell repertoire in Common Variable Immune Deficiency. Clin. 34. Warnatz, K. et al. Severe deficiency of switched memory B cells (CD27+IgM- Immunol. 178, 1–9 (2017). IgD-) in subgroups of patients with common variable immunodeficiency: a 5. Gao, Y. et al. Common variable immunodeficiency is associated with a new approach to classify a heterogeneous disease. Blood 99, 1544–1551 (2002). functional deficiency of invariant natural killer T cells. J. Allergy Clin. 35. Aref-Eshghi, E. et al. Evaluation of DNA methylation episignatures for Immunol. 133, 1420–1428 (2014). 1428.e1. diagnosis and phenotype correlations in 42 Mendelian Neurodevelopmental 6. Perreau, M. et al. Exhaustion of bacteria-specific CD4 T cells and microbial disorders. Am. J. Hum. Genet. 106, 356–370 (2020). translocation in common variable immunodeficiency disorders. J. Exp. Med. 36. Chater-Diehl, E. et al. Anatomy of DNA methylation signatures: Emerging 211, 2033–2045 (2014). insights and applications. Am. J. Hum. Genet. 108, 1359–1366 (2021). 7. Yu, G. P. et al. Regulatory T cell dysfunction in subjects with common variable 37. de la Rica, L. et al. PU.1 target genes undergo Tet2-coupled demethylation and immunodeficiency complicated by autoimmune disease. Clin. Immunol. 131, DNMT3b-mediated methylation in monocyte-to-osteoclast differentiation. 240–253 (2009). Genome Biol. 14, 1–21 (2013). 8. Vořechovský, I. et al. Family and linkage study of selective IgA deficiency and 38. Wang, Y. et al. WT1 recruits TET2 to regulate its target gene expression and common variable immunodeficiency. Clin. Immunol. Immunopathol. 77, suppress leukemia cell proliferation. Mol. Cell 57, 662–673 (2015). 185–192 (1995). 39. Xiong, J. et al. Cooperative action between SALL4A and TET proteins in 9. Abolhassani, H., Hammarström, L. & Cunningham-Rundles, C. Current genetic stepwise oxidation of 5-methylcytosine. Mol. Cell 64, 913–925 (2016). landscape in common variable immune deficiency. Blood 135, 656–667 (2020). 40. Costa, Y. et al. NANOG-dependent function of TET1 and TET2 in 10. Rodríguez-Cortez, V. C. et al. Monozygotic twins discordant for common establishment of pluripotency. Nature 495, 370–374 (2013). variable immunodeficiency reveal impaired DNA demethylation during naïve- 41. Grötsch, B. et al. The AP-1 transcription factor Fra1 inhibits follicular B cell to-memory B-cell transition. Nat. Commun. 6, 7335 (2015). differentiation into plasma cells. J. Exp. Med 211, 2199–2212 (2014). 11. Campos-Sanchez, E., Martínez-Cano, J., Del Pino Molina, L., López-Granados, 42. Wang, H. et al. Transcription factors IRF8 and PU.1 are required for follicular E. & Cobaleda, C. Epigenetic Deregulation in Human Primary B cell development and BCL6-driven germinal center responses. Proc. Natl. Immunodeficiencies. Trends Immunol 40, 49–65 (2019). Acad. Sci. U.S.A. 116, 9511–9520 (2019). 12. Kulis, M. et al. Whole-genome fingerprint of the DNA methylome during 43. Carotta, S. et al. The transcription factors IRF8 and PU.1 negatively regulate human B cell differentiation. Nat. Genet. 47, 746–756 (2015). plasma cell differentiation. J. Exp. Med. 211, 2169–2181 (2014). 13. Oakes, C. C. et al. DNA methylation dynamics during B cell maturation 44. Pérez-García, A. et al. CTCF orchestrates the germinal centre transcriptional underlie a continuum of disease phenotypes in chronic lymphocytic leukemia. program and prevents premature plasma cell differentiation. Nat. Commun. 8, Nat. Genet. 48, 253–264 (2016). 1–12 (2017). 14. Efremova, M., Vento-Tormo, R., Park, J.-E., Teichmann, S. A. & James, K. R. 45. Wang, H. et al. Widespread plasticity in CTCF occupancy linked to DNA Immunology in the era of single-cell technologies. Annu. Rev. Immunol. 38, methylation. Genome Res. 22, 1680–1688 (2012). https://doi.org/10.1146/annurev-immunol-090419-020340 (2020). 46. Hong, S. et al. Methyl-dependent and spatial-specific DNA recognition by the 15. Smallwood, S. A. et al. Single-cell genome-wide bisulfite sequencing for orthologous transcription factors human AP-1 and Epstein-Barr virus Zta. assessing epigenetic heterogeneity. Nat. Methods 11, 817–820 (2014). Nucleic Acids Res. 45, 2503–2515 (2017). 16. Angermueller, C., Lee, H. J., Reik, W. & Stegle, O. DeepCpG: accurate 47. Zuo, Z., Roy, B., Chang, Y. K., Granas, D. & Stormo, G. D. Measuring prediction of single-cell DNA methylation states using deep learning. Genome quantitative effects of methylation on transcription factor-DNA binding Biol. 18, 1–13 (2017). affinity. Sci. Adv. 3, eaao1799 (2017). 17. Souren, N. Y. et al. DNA methylation signatures of monozygotic twins 48. Liu, B. et al. Affinity-coupled CCL22 promotes positive selection in germinal clinically discordant for multiple sclerosis. Nat. Commun. 10, 1–12 (2019). centres. Nature 592, 133–137 (2021). 18. Vadgama, N. et al. De novo single-nucleotide and copy number variation in 49. Arens, R. et al. Signaling through CD70 regulates B cell activation and IgG discordant monozygotic twins reveals disease-related genes. Eur. J. Hum. production. J. Immunol. 173, 3901–3908 (2004). Genet. 27, 1121–1133 (2019). 50. Kobata, T. et al. CD27-CD70 interactions regulate B-cell activation by T cells. 19. Landrum, M. J. et al. ClinVar: public archive of relationships among sequence Proc. Natl. Acad. Sci. U.S.A. 92, 11249–11253 (1995). variation and human phenotype. Nucleic Acids Res. 42, D980–D985 (2014). 51. Carrasco, Y. R., Fleire, S. J., Cameron, T., Dustin, M. L. & Batista, F. D. LFA-1/ 20. Cingolani, P. et al. A program for annotating and predicting the effects of ICAM-1 interaction lowers the threshold of B cell activation by facilitating B single nucleotide polymorphisms, SnpEff: SNPs in the genome of Drosophila cell adhesion and synapse formation. Immunity 20, 589–599 (2004). melanogaster strain w1118; iso-2; iso-3. Fly 6, 80–92 (2012). 52. Giovannone, N. et al. Galectin-9 suppresses B cell receptor signaling and is 21. Jühling, F. et al. Metilene: fast and sensitive calling of differentially methylated regulated by I-branching of N-glycans. Nat. Commun. 9, 3287 (2018). regions from bisulfite sequencing data. Genome Res. 26, 256–262 (2016). 53. Cao, A. et al. Galectin-9 binds IgM-BCR to regulate B cell signaling. Nat. 22. Zhou, W. et al. DNA methylation loss in late-replicating domains is linked to Commun. 9, 3288 (2018). mitotic cell division. Nat. Genet. 50, 591–602 (2018). 54. Adachi, T., Wakabayashi, C., Nakayama, T., Yakura, H. & Tsubata, T. CD72 23. Berman, B. P. et al. Regions of focal DNA hypermethylation and long-range negatively regulates signaling through the antigen receptor of B cells. J. hypomethylation in colorectal cancer coincide with nuclear lamina-associated Immunol. 164, 1223–1229 (2000). domains. Nat. Genet. 44, 40–46 (2011). 55. Li, D. H. et al. CD72 down-modulates BCR-induced signal transduction and 24. Driessen, G. J. et al. B-cell replication history and somatic hypermutation diminishes survival in primary mature B lymphocytes. J. Immunol. 176, status identify distinct pathophysiologic backgrounds in common variable 5321–5328 (2006). immunodeficiency. Blood 118, 6814–6823 (2011). 56. Thon, V. et al. Defective integration of activating signals derived from the T 25. Del Pino-Molina, L. et al. Impaired CpG demethylation in common variable cell receptor (TCR) and costimulatory molecules in both CD4+ and CD8+ T immunodeficiency associates with B cell phenotype and proliferation rate. lymphocytes of common variable immunodeficiency (CVID) patients. Clin. Front. Immunol. 10, 1–11 (2019). Exp. Immunol. 110, 174–181 (1997). 26. Bravo González-Blas, C. et al. cisTopic: cis-regulatory topic modeling on 57. Brugnoni, D. et al. CD4+ cells from patients with Common Variable single-cell ATAC-seq data. Nat. Methods 16, 397–400 (2019). Immunodeficiency have a reduced ability of CD40 ligand membrane expression 27. Iwata, A. et al. Quality of TCR signaling encoded by differential enhancer after in vitro stimulation. Pediatr. Allergy Immunol. 7, 176–179 (1996). affinities for BATF/IRF4. Nat. Immunol. 18, 563–572 (2017). 58. Pons, J., Ferrer, J. M., Martínez-Pomar, N., Iglesias-Alzueta, J. & Matamoros, 28. Bao, K. et al. BATF modulates the Th2 locus control region and regulates N. Costimulatory molecules and cytokine production by T lymphocytes in CD4+ T cell fate during antihelminth immunity. J. Immunol. 197, 4371–4381 common variable immunodeficiency disease. Scand. J. Immunol. 63, 383–389 (2016). (2006). 29. Ren, G. et al. CTCF-mediated enhancer-promoter interaction is a critical 59. Seidel, M. G. et al. The European Society for Immunodeficiencies (ESID) regulator of cell-to-cell variation of gene expression. Mol. Cell 67, Registry Working Definitions for the Clinical Diagnosis of Inborn Errors of 1049–1058.e6 (2017). Immunity. J. Allergy Clin. Immunol. Pract. 7, 1763–1770 (2019). 30. Picelli, S. et al. Full-length RNA-seq from single cells using Smart-seq2. Nat. 60. Andrés-León, E., Núñez-Torres, R. & Rojas, A. M. miARma-Seq: a comprehensive Protoc. 9, 171–181 (2014). tool for miRNA, mRNA and circRNA analysis. Sci. Rep. 6, 1–7 (2016). 31. Lai, A. Y. et al. Dna methylation profiling in human b cells reveals immune 61. Babraham Bioinformatics—FastQC a quality control tool for high throughput regulatory elements and epigenetic plasticity at alu elements during b-cell sequence data. https://www.bioinformatics.babraham.ac.uk/projects/fastqc/. activation. Genome Res. 23, 2030–2041 (2013). Version 0.11.5 released. (08-03-2016). 16 NATURE COMMUNICATIONS | (2022)13:1779 | https://doi.org/10.1038/s41467-022-29450-x | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-022-29450-x ARTICLE 62. Li, H. & Durbin, R. Fast and accurate long-read alignment with Burrows- Health (PI16/00759) and European Regional Development Fund-European Social Fund Wheeler transform. Bioinformatics 26, 589–595 (2010). —FEDER-FSE) (C.R-G.), Grupo DISA (OA18/017) (C.R.-G.), the UK Biotechnology and 63. Van der Auwera, G. A. et al. From fastQ data to high-confidence variant calls: Biological Sciences Research Council (BBS/E/B/000C0426) (G.K.) and Medical Research the genome analysis toolkit best practices pipeline. Curr. Protoc. Bioinforma. Council (MR/S000437/1) (G.K.). We are indebted to the donors for participating in this 43, 11.10.1–11.10.33 (2013). research. We thank Antonio Garcia-Gomez for graphical design support, Sarah Teich- 64. Quinlan, A. R. BEDTools: The Swiss-Army tool for genome feature analysis. mann for her useful feedback, Hamish King for helping with single-cell germinal center Curr. Protoc. Bioinform. 2014, 11.12.1–11.12.34 (2014). dataset availability, Xi Chen for performing scATAC-seq analysis, Kirsty Ambridge and 65. Karczewski, K. J. et al. The mutational constraint spectrum quantified from Elena Prigmore for their involvement in single-cell RNA library generation, Martin Prete variation in 141,456 humans. Nature 581, 434–443 (2020). for creating online visualizations for our cell atlas and Esther Castaño and Beatriz 66. Auton, A. et al. A global reference for human genetic variation. Nature 526, Barroso from CCiTUB Cytometry Unit for their support with single-cell sorting and Dr. 68–74 (2015). Carla Gianelli and Dr. Rebeca Rodríguez Pena for the patient follow-up in the CVID 67. Sherry, S. T. et al. DbSNP: the NCBI database of genetic variation. Nucleic cohort. Acids Res. 29, 308–311 (2001). 68. Chen, X., Miragaia, R. J., Natarajan, K. N. & Teichmann, S. A. A rapid and robust method for single cell chromatin accessibility profiling. Nat. Commun. Author contributions J.R.-U., R.V.-T. and E.B. conceived the study; J.-R.U., L.P.-M. and L.C. designed and 9, 1–9 (2018). performed cell sorter, cell activation and KREC experiments; J.R.-U., L.C. and H.H. 69. Villani, A. C. et al. Single-cell RNA-seq reveals new types of human blood performed 10× Chromium experiments; J.R-U., S.J.C. and L.M. performed sample and dendritic cells, monocytes, and progenitors. Science. 356, 1–12 (2017). library preparation; J.R-U., A.A., M.J.B, C.C.-F. L.G.-A., L.-F. H., E.A.L. F.K., K.P., S.V.D. 70. Clark, S. J. et al. Genome-wide base-resolution mapping of DNA methylation and V.K. performed computational analysis; J.R.-U., A.A., J.M., K.W., R.V.-T. and E.B. in single cells using single-cell bisulfite sequencing (scBS-seq). Nat. Protoc. 12, analyzed and interpreted the data with contributions from L.P.-M, F.C.-M., V.R.-C., 534–547 (2017). E.L.-G., C.R.-G., O.S. and G.K.; L.P.-M., M.M-S., E.L.-G. and C.R-G. contributed with 71. Martin, M. Cutadapt removes adapter sequences from high-throughput clinical material and clinical interpretation of the results; J.R-U. and E.B. wrote the sequencing reads. EMBnet.journal 17, 10 (2011). manuscript with contributions from A.A. and R.V-T.; R.V.-T. and E.B. co-directed 72. Krueger, F. & Andrews, S. R. Bismark: a flexible aligner and methylation caller the study. All authors read and accepted the paper. for Bisulfite-Seq applications. Bioinformatics 27, 1571–1572 (2011). 73. Ritchie, M. E. et al. Limma powers differential expression analyses for RNA- sequencing and microarray studies. Nucleic Acids Res. 43, e47 (2015). Competing interests 74. McLean, C. Y. et al. GREAT improves functional interpretation of cis- The authors declare no competing interests. regulatory regions. Nat. Biotechnol. 28, 495–501 (2010). 75. Garcia-Alonso, L. et al. Transcription factor activities enhance markers of drug sensitivity in cancer. Cancer Res. 78, 769–780 (2018). Additional information 76. Dobin, A. et al. STAR: Ultrafast universal RNA-seq aligner. Bioinformatics 29, Supplementary information The online version contains supplementary material 15–21 (2013). available at https://doi.org/10.1038/s41467-022-29450-x. 77. Wolf, F. A., Angerer, P. & Theis, F. J. SCANPY: large-scale single-cell gene expression data analysis. Genome Biol. 19, 1–5 (2018). Correspondence and requests for materials should be addressed to Javier Rodríguez- 78. Heaton, H. et al. Souporcell: robust clustering of single-cell RNA-seq data by Ubreva, Roser Vento-Tormo or Esteban Ballestar. genotype without reference genotypes. Nat. Methods 17, 615–620 (2020). 79. Purcell, S. et al. PLINK: A tool set for whole-genome association and Peer review information Nature Communications thanks Rhys Allan, César Cobaleda, population-based linkage analyses. Am. J. Hum. Genet. 81, 559–575 (2007). and Pierre Milpied for their contribution to the peer review of this work. 80. Pijuan-Sala, B. et al. A single-cell molecular map of mouse gastrulation and early organogenesis. Nature 566, 490–495 (2019). Reprints and permission information is available at http://www.nature.com/reprints 81. Wolock, S. L., Lopez, R. & Klein, A. M. Scrublet: computational identification of Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in cell doublets in single-cell transcriptomic data. Cell Syst. 8, 281–291.e9 (2019). published maps and institutional affiliations. 82. Vento-Tormo, R. et al. Single-cell reconstruction of the early maternal–fetal interface in humans. Nature 563, 347–353 (2018). 83. Arutyunyan, A. & Bonder, M. J. ventolab/CVID: December-2021. https:// doi.org/10.5281/ZENODO.5792221 (2021). Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give Acknowledgements appropriate credit to the original author(s) and the source, provide a link to the Creative We thank the CERCA Program/Generalitat de Catalunya and the Josep Carreras Commons license, and indicate if changes were made. The images or other third party Foundation for institutional support. This publication is part of the Human Cell Atlas; material in this article are included in the article’s Creative Commons license, unless www.humancellatlas.org/publications. This study was funded by: Spanish Ministry of indicated otherwise in a credit line to the material. If material is not included in the Science and Innovation (grant number PID2020-117212RB-I00/AEI/10.13038/ article’s Creative Commons license and your intended use is not permitted by statutory 501100011033) (E.B.), Instituto de Salud Carlos III (ISCIII), Ref. AC18/00057, associated regulation or exceeds the permitted use, you will need to obtain permission directly from with i-PAD project (ERARE European Union program) (E.B.), the Jeffrey Modell the copyright holder. To view a copy of this license, visit http://creativecommons.org/ Foundation (E.B.), Wellcome Sanger core funding (grant no. WT206194) (R.V.-T.), the licenses/by/4.0/. Chan Zuckerberg Initiative (grant 2020-216799) (R.V.-T. and E.B.), an EMBO short- term fellowship (J.R.U.), Fondo de Investigación Sanitaria Instituto de Salud Carlos III (FIS PI16/01605) (L.P.-M.), the Spanish Ministry of Science, Innovation and Universities © The Author(s) 2022 (SAF2017-89109-P; AEI/FEDER, UE) (H.H.), Instituto de Salud Carlos III, Ministry of 1 Epigenetics and Immune Disease Group, Josep Carreras Research Institute (IJC), Badalona, Barcelona, Spain. 2Chromatin and Disease Group, Cancer Epigenetics and Biology Programme (PEBC), Bellvitge Biomedical Research Institute (IDIBELL), L’Hospitalet de Llobregat, Barcelona, Spain. 3Cellular Genetics Department, Wellcome Sanger Institute, Wellcome Genome Campus, Hinxton, UK. 4Centre for Trophoblast Research, University of Cambridge, Cambridge, UK. 5European Molecular Biology Laboratory, European Bioinformatics Institute (EMBL-EBI), Wellcome Genome Campus, Hinxton, Cambridgeshire, UK. 6European Molecular Biology Laboratory, Genome Biology Unit, Heidelberg, Germany. 7 Division of Computational Genomics and Systems Genetics, German Cancer Research Center (DKFZ), Heidelberg, Germany. 8Clinical Immunology Department, La Paz University Hospital and Lymphocyte Pathophysiology in Immunodeficiencies Group, La Paz Institute for Health Research (IdiPAZ), Center for Biomedical Network Research on Rare Diseases (CIBERER U767), Madrid, Spain. 9Epigenetics Programme, Babraham Institute, Cambridge, UK. 10Instituto de Parasitología y Biomedicina López-Neyra, Consejo Superior de Investigaciones Científicas (IPBLN-CSIC), Granada, Spain. 11Bioinformatics Unit, Babraham Institute, Cambridge, UK. 12Department of Immunology, Hospital Universitario de NATURE COMMUNICATIONS | (2022)13:1779 | https://doi.org/10.1038/s41467-022-29450-x | www.nature.com/naturecommunications 17 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-022-29450-x Gran Canaria Dr Negrín, Las Palmas de Gran Canaria, Spain. 13CNAG-CRG, Centre for Genomic Regulation (CRG), Barcelona Institute of Science and Technology (BIST), Barcelona, Spain. 14Universitat Pompeu Fabra (UPF), Barcelona, Spain. 15Faculty of Medicine, Department of Rheumatology and Clinical Immunology, Medical Center, University of Freiburg, Freiburg, Germany. 16Faculty of Medicine, Center for Chronic Immunodeficiency (CCI), Medical Center, University of Freiburg, Freiburg, Germany. 17Department of Clinical Sciences, University Fernando Pessoa, Las Palmas de Gran Canaria, Spain. 18Wellcome Sanger Institute, Wellcome Trust Genome Campus, Hinxton, UK. 19Epigenetics in Inflammatory and Metabolic Diseases Laboratory, Health Science Center (HSC), East China Normal University (ECNU), Shanghai, China. 20 These authors contributed equally: Javier Rodríguez-Ubreva, Anna Arutyunyan. 21These authors jointly supervised this work: Roser Vento- Tormo, Esteban Ballestar. ✉email: jrodriguez@carrerasresearch.org; rv4@sanger.ac.uk; eballestar@carrerasresearch.org 18 NATURE COMMUNICATIONS | (2022)13:1779 | https://doi.org/10.1038/s41467-022-29450-x | www.nature.com/naturecommunications