IgD from atypical-like memory B cells and plasma cells targets commensal and environmental antigens

The human nasopharyngeal mucosa includes palatine, pharyngeal, lingual, and tubal tonsils, which are organized lymphoepithelial structures representing the human equivalent of the murine nasal-associated lymphoid tissue (Cerutti et al., 2011). Fissure-like openings on the tonsillar surface, termed crypts, convey commensal, airborne, and food antigens from the lumen of the crypt to the underlying stratified epithelium, which is inhabited by abundant immune cells, including IgD⁺IgM⁻ plasma cells (IgD-PCs) (Cerutti et al., 2011). These plasma cells (PCs) are thought to differentiate from IgD⁺IgM⁻ germinal center (IgD-GC) B cell precursors that have undergone IgM-to-IgD class-switch recombination (CSR) and somatic hypermutation (SHM) in local lymphoid follicles (Chen et al., 2020).

In both humans and mice, CSR from IgM to IgD occurs through an unconventional pathway that initiates IgD responses to both autologous and microbial antigens (Chen et al., 2009; Choi et al., 2017; Koelsch et al., 2007; Rouaud et al., 2014; Xu et al., 2022). In humans, the IgD-GC B cells involved in IgD responses express highly hypermutated Ig heavy (IgH) and light (IgL) chain genes (Chen et al., 2020) and differentiate into IgD-PCs characterized by biased usage of IgL chains of the type lambda (Igλ) (Chen et al., 2009; Arpin et al., 1998). Due to the pronounced autoreactivity of IgD class-switched cells and their putative inability to further switch to downstream IgG, IgA, or IgE isotypes, class switching to IgD has been proposed to have a tolerogenic function, serving as a “sink” for autoreactive B cells (Koelsch et al., 2007).

While tonsillar IgD-GC B cells are thought to constitute the major precursors of IgD-PCs (Chen et al., 2020; Arpin et al., 1998), direct evidence of the clonal relationship between IgD-GC B cells and IgD-PCs remains limited at best. It is also unclear whether tonsillar IgD-GC B cells generate IgD⁺IgM⁻ memory (IgD-ME) B cells and whether these IgD class-switched cells can further differentiate into IgD-PCs. In this regard, IgD-ME B cells have mostly been reported in the general circulation (Koelsch et al., 2007; Stensland et al., 2023) and their phenotypic, transcriptional, and molecular landscapes, including clonal relationship with IgD-PCs and canonical ME B cell subsets, remain elusive. Moreover, it is unknown whether IgD from tonsillar IgD-PCs targets environmental antigens commonly present in the nasopharyngeal mucosa. This possibility would echo recent findings, suggesting that IgD responses to some airborne or food allergens correlate with enhanced tolerance to these allergens (Shan et al., 2018; Boonpiyathad et al., 2020; Suprun et al., 2020; Itoh et al., 2021; Satitsuksanoa et al., 2025).

Here, we found that IgD⁺IgM⁻ B cells included a phenotypically and transcriptionally heterogeneous population of IgD-ME B cells that mostly inhabited epithelial and interfollicular areas of the tonsillar mucosa and shared some properties with atypical memory B cells (ABCs). In agreement with their extensive clonal affiliation with IgD-PCs, IgD-ME B cells differentiated into IgD-PCs via a mutation-intensive pathway that required a complex mix of innate and adaptive signals. This IgD response targeted multiple sinonasal and bronchial commensal bacteria in addition to common environmental antigens, including airborne allergens. By showing that germline IgD revertants exhibited attenuated reactivity to these antigens, our data suggest that IgD responses emerge from an IgD-ME B cell–centered mutation-dependent nasopharyngeal program presumably aimed at enhancing mucosal homeostasis, immunity, and tolerance.

Previous studies have identified a population of circulating IgD⁺IgM⁻CD27⁺ B cells, potentially representing IgD-ME B cells (Koelsch et al., 2007; Stensland et al., 2023), but their existence in aerodigestive tissues remains elusive. We implemented unbiased clustering of tonsillar IgD⁺IgM⁻ B cells by analyzing IgD, IgM, CD10, CD19, CD27, CD38, and other surface molecules through spectral flow cytometry. Cluster visualization onto a t-distributed stochastic neighbor embedding (t-SNE) projection revealed four major IgD⁺IgM⁻ B cell subsets, including CD10⁻CD27^hiCD38^hi IgD-PCs, CD10⁺CD38⁺ IgD-GC B cells, putative CD10⁻CD27⁺CD38⁻ IgD-ME B cells, and CD10⁻CD27⁻CD38^low anergic naive (IgD-AN) B cells (Fig. 1 A and Fig. S1 A). This last subset is autoreactive and downregulates IgM (Duty et al., 2009), thereby mimicking IgD class-switched IgD⁺IgM⁻ B cells.

Consistent with the known bias of IgD class-switched B cells for Igλ expression (Chen et al., 2009), IgD-ME B cells and IgD-PCs showed a high usage of Igλ (Fig. S1 A). In agreement with the presence of a distinct ABC subset within the ME B cell population (Winslow and Levack, 2025), IgD-ME but not IgD-GC B cells displayed less CD21 (Fig. S1 A), a phenotypic property of some ABCs (Sutton et al., 2021). Like ABCs (Sutton et al., 2021; Holla et al., 2021), IgD-ME B cells also expressed the activation-related molecules CD69, CD11c, CD95 (FAS), and CXCR3, but not CD43 (Fig. S1 A). Confocal microscopy confirmed the existence of IgD⁺IgM⁻ IgD-ME B cells, which were detected within interfollicular areas and were distinguishable from both IgD⁺IgM⁺ naive B cells and IgD⁺IgM⁻ IgD-PCs, the latter showing abundant intracellular IgD but not IgM (Fig. 1 B; and Fig. S1, B and C). As shown earlier (Quách et al., 2011), the lack of intracellular IgM also distinguishes IgD-ME from IgD-AN B cells.

To verify whether IgD-ME B cells derived from hypermutated GC precursors, we compared the mutational load of their IgH chain variable (IGHV) genes with that of unswitched naive and IgD-AN B cells, as well as class-switched IgG-ME and IgA-ME B cells. Each of these B cell subsets was identified by high-dimensional flow cytometry and then purified by sorting, IgG-ME and IgA-ME B cells being IgD⁻IgM⁻CD10⁻CD19⁺ CD38⁻CD27⁺ cells that either lacked IgA in addition to IgD and IgM or expressed IgA, respectively (Fig. S1 D). High-throughput next-generation sequencing of IGHV genes revealed that similar to IgD-GC B cells and IgD-PCs, as well as IgG-ME and IgA-ME B cells (Arpin et al., 1998; Liu et al., 1996), most IgD-ME B cells were hypermutated, whereas naive and IgD-AN B cells were not (Fig. 1 C).

Next, flow cytometry showed that tonsillar IgD-ME B cells expressed more pronounced ABC-like traits compared with tonsillar naive, IgG-ME, and, to a lesser extent, IgA-ME B cell subsets, including higher CD11c, CD95, CXCR3, FCRL4, and CD69 expression, but lower CD21 and comparable CD27 expression (Fig. S1 E). This phenotype was roughly similar to that of circulating and splenic IgD-ME B cells, although tonsillar IgD-ME B cells also evidenced clear differences (Fig. S1 E).

Bulk RNA sequencing (RNA-seq) was performed to further characterize differences and similarities between IgD-ME B cells and unswitched naive B cells or class-switched IgG-ME and IgA-ME B cells. Principal component analysis (PCA) showed that the transcriptome of IgD-ME B cells clustered away from the transcriptome of naive B cells, but next to the transcriptomes of IgG-ME and IgA-ME B cells (Fig. 1 D). Volcano plot and heatmap tools visualized over three thousand differentially expressed genes (DEGs) between IgD-ME and naive B cells (Fig. 1, E and F). Of these DEGs, CD70, CD80, FCRL4, AIM2, and MUC16 are known to mark ME B cells from the respiratory tract (Borst et al., 2005; Cancro and Tomayko, 2021; Yang et al., 2021; Moroney et al., 2020) and were enriched in IgD-ME B cells (Fig. 1 F). Of note, PC-associated syndecan-1 (SDC1), XBP1, and JCHAIN were also enriched (Fig. 1 F), which highlighted the activated state of IgD-ME B cells and the increased predisposition of at least some of these cells to undergo PC differentiation.

We then defined the topography of tonsillar IgD-ME B cells by using a recently published spatial transcriptomics dataset (Massoni-Badosa et al., 2024). A CD19⁺IGHD⁺IGHM⁻CD27⁺CD38⁻CD80⁺MME⁻ (CD10) TNFRSF13B⁺ (TACI) ITGAX⁺ (CD11c) FAS⁺ (CD95) CXCR3⁺NT5E⁺CR2⁻CD3G⁻CD4⁻CD5⁻CD7⁻CD8A⁻ gene signature designed on the basis of published data (Chen et al., 2009; Sutton et al., 2021; Holla et al., 2021; Massoni-Badosa et al., 2024; Gao et al., 2024) confirmed the interfollicular location of some IgD-ME B cells and mapped more of them within subepithelial and epithelial areas (Fig. 1, G–I). Finally, as shown by flow cytometry, the frequency of IgD-ME B cells within the class-switched IgM⁻CD10⁻CD19⁺CD38⁺ ME pool was higher in tonsils than in the spleen, but showed intermediate values in the peripheral blood (Fig. 1 J), which may reflect the induction of IgD responses in the aerodigestive mucosa, followed by the transition of some IgD-ME B cells from tonsillar inductive sites to distal effector sites via the general circulation (Chen et al., 2009). Thus, IgD-ME B cells are distinct from other tonsillar IgD⁺ B cell subsets and mostly inhabit epithelial and subepithelial areas, being nonetheless present also in interfollicular areas.

Next, we wondered about specific gene expression differences between IgD-ME B cells and class-switched IgG-ME or IgA-ME B cells in tonsils. Purified IgD-ME B cells were first confirmed to express a phenotypic hallmark of bona fide IgD class-switched cells (Chen et al., 2009; Arpin et al., 1998) by showing preferential Igλ expression compared with naive, IgG-ME, and IgA-ME B cells (Fig. S2 A). Bulk RNA-seq–derived volcano plots showed that IgD-ME B cells displayed hundreds of DEGs compared with IgG-ME or IgA-ME B cells (Fig. 2 A). As visualized by heatmaps (Fig. 2, B and C), IgD-ME B cells exhibited the upregulated expression of immune response gene sets typically associated with ABCs (Holla et al., 2021; Gao et al., 2024; Yang et al., 2022), including JCHAIN, SLAMF7, CD274 (PD-L1), ITGAX, FCGR2B, TNFRSF13B, CD86, CXCR3, TBX21 (T-BET), CD151, ITGB2 (integrin β2), ITGB7 (integrin β7), FAS, and IL10RA.

The enriched ABC-like gene signature of IgD-ME B cells further included the downregulated expression of lymphotoxin-β (LTB), CD55, CCR7, CD40, PAX5 (BSAP), and SELL (L-selectin) compared with IgG-ME or IgA-ME B cells (Holla et al., 2021; Gao et al., 2024; Yang et al., 2022) (Fig. 2 C). However, unlike ABCs (Gao et al., 2024), IgD-ME B cells showed downregulated thymocyte selection–associated HMG box (TOX) and ezrin (EZR) expression and no enrichment in zinc finger E-box–binding homeobox 2 (ZEB2) (Fig. 2 C), a transcription factor involved in ABC differentiation (Dai et al., 2024). Given its importance in ABC differentiation (Winslow and Levack, 2025; Yang et al., 2022), T-BET–encoding TBX21 was further quantified bioinformatically, which confirmed its enrichment in IgD-ME B cells compared with control naive, IgG-ME, and, although not significantly, IgA-ME B cells (Fig. S2 B).

The ABC-like properties of IgD-ME B cells were corroborated by gene set enrichment analysis (GSEA) (Fig. 2 D) and further substantiated when our bulk RNA-seq data were cross-referenced to recently published cellular indexing of transcriptomes and epitopes by sequencing (CITE-seq), as well as single-cell RNA-seq (scRNA-seq) data (Massoni-Badosa et al., 2024), as detailed in Materials and methods. Compared with classical IgG-ME and IgA-ME B cells, IgD-ME B cells showed more similarity to ABCs (Winslow and Levack, 2025; Yang et al., 2022) when ITGAX, IFNGR1, FGR, FCGR2B, SOX5, FCRL5, ZEB2, FCRL4, and CCR1 mRNAs were analyzed (Fig. 2 E). Of note, IFNGR1 encodes the α chain of the receptor for IFN-ɣ, a cytokine involved in T-BET–driven ABC differentiation (Obeng-Adjei et al., 2017; Jenks et al., 2018).

Some of the above gene expression differences reflected similar patterns detected earlier at the protein level, including differences in CD11c and FCRL4 expression. As shown by t-SNE projections, IgD-ME B cells exhibited more pronounced phenotypic differences with IgG-ME than IgA-ME B cells (Fig. 2 F), which could reflect the shared mucosal nature of IgD and IgA responses (Chen et al., 2020). Lastly, CD11c and CXCR3 expression, which is a hallmark of ABCs (Winslow and Levack, 2025; Yang et al., 2022), was further validated in IgD-ME B cells by confocal imaging (Fig. 2 G and Fig. S2 D). Thus, tonsillar IgD-ME B cells are distinct from class-switched IgG-ME and IgA-ME B cells and at least some of them express ABC-like properties.

Further analysis of a recently published ABC-derived scRNA-seq dataset (Massoni-Badosa et al., 2024) identified two distinct subclusters within the IgD-ME B cell cluster, IgD-ME 1 being larger than IgD-ME 2 (Fig. 3 A). These subclusters showed highly comparable expression of ME B cell– and ABC-associated mRNAs, including ITGAX, TBX21, MS4A1 (CD20), CD44, IFNGR1, and FGR (Fig. 3, A and B). However, IgD-ME 1 was enriched in FCGR2B, LTB, TNFRSF13B, SSR4, and gelsolin (GSN), whereas IgD-ME 2 was enriched in CD69, TOX, and SOX5 (Fig. 3, A and B). To gain some insights into the functional properties of these subclusters, all the genes differentially expressed by IgD-ME 1 and the IgD-ME 2 (Table S1) were processed through GSEA, followed by gene ontology (GO) analysis of the 20 most upregulated and 20 most downregulated GO terms (Fig. 3 C). Of the latter, 11 were selected based on their biological relevance (Fig. 3 C).

Compared with IgD-ME 2, IgD-ME 1 was enriched in gene sets associated with ribosomal biogenesis and assembly, protein translation, ribosomal RNA processing, adenosine triphosphate synthesis, and antigen presentation (Fig. 3 C). Compared with IgD-ME 1, IgD-ME 2 was enriched in gene sets linked to mitosis, reduced cell adhesion to the cell matrix, and guanosine triphosphate–targeting signaling (Fig. 3 C). These processes are associated with plasmablast expansion and cytoskeleton remodeling needed for organelle trafficking and changes in cell polarity and locomotion (Etienne-Manneville and Hall, 2002). Thus, IgD-ME B cells may encompass a larger IgD-ME 1 subcluster gearing up for increased antibody synthesis, and a smaller IgD-ME 2 subcluster getting ready for enhanced proliferation, as well as increased mobility, shape remodeling, and organelle trafficking.

We next wondered whether tonsillar IgD-ME B cells emerge from the GC, as at least some classical ME B cells do (Cancro and Tomayko, 2021). We found that IgD-ME B cells displayed a higher mean mutation count compared with canonical ME B cells (Fig. 4 A). Compared with these cells, IgD-ME B cells also showed more mutations in antigen-binding complementarity-determining regions (CDRs) and scaffold-providing framework regions (FWRs) (Fig. S3 A). Of note, the CDRs from IgD-ME B cells featured higher replacement-to-silent mutation ratio (R:S ratio) compared with other ME B cells and PC subsets (Fig. S3 B). As shown by the κ-deleting recombination excision assay (van Zelm et al., 2007), the higher mutational load of IgD-ME B cells was uncoupled from an increased in vivo proliferation history of these cells and their precursors compared with conventional IgM-ME, IgG-ME, and IgA-ME B cells (Fig. 4 B). Using data obtained from a published scRNA-seq dataset (Massoni-Badosa et al., 2024), the comparison of IgD-GC B cells with IgG-GC B cells evidenced that the former were significantly enriched in RAD52, which encodes the DNA repair RADiation sensitive 52 (RAD52) protein involved in the mismatch repair (MMR) pathway of SHM and essential for IgM-to-IgD CSR (Xu et al., 2022); POLQ, which encodes the error-prone polymerase θ (POLQ) required for both MMR and base excision repair (BER) pathways of SHM (Odegard and Schatz, 2006; Casellas et al., 2016); and, lastly, EXO1 and PMS1, which encode the exonuclease 1 and postmeiotic segregation increased 1 proteins involved in the MMR pathway of SHM (Odegard and Schatz, 2006; Casellas et al., 2016) (Fig. 4 C). Accordingly, the expression ratio of transcripts encoding MMR versus transcripts encoding BER proteins was increased in IgD-GC B cells compared with IgG-GC B cells (Fig. 4 D). Thus, the increased mutational load of IgD-ME B cells may stem from the unique SHM machinery of IgD-GC precursors rather than increased IgD-GC proliferation.

To gain additional insights into the molecular landscape of tonsillar IgD-ME B cells, we compared their Ig V_HDJ_H gene usage with that of other B cell and PC subsets expressing IgD, IgA, IgG, or IgM. As shown earlier (Duty et al., 2009), IgD-AN B cells expressed both IgD and IgM transcripts as naive B cells did. This co-expression and the possibility that rare CD27⁻ B cells from the IgD-ME subset could contaminate naive and IgD-AN B cells prompted us to focus on IgM rather than IgD mRNAs. We found that IgD-GC and IgD-ME B cells, as well as IgD-PCs, had a unique IGHV usage pattern, distinct from that of GC B cells, ME B cells, and PCs expressing IgM, IgG, or IgA and also distinct from that of naive and IgD-AN B cells (Fig. 4 E). Consistently, IgD class-switched cell subsets from all four donors largely clustered together and separately from all other B cell and PC subsets (Fig. 4 E).

Compared with ME B cells expressing IgM, IgG, or IgA, IgD-ME B cells used less IGHV3-7, IGHV3-13, and IGHV3-48, but more IGHV4-30-2 and IGHV4-31 (Fig. 4 F). Of note, IgD-GC B cells, IgD-ME B cells, and IgD-PCs also showed trends toward lower IGHV3-23 usage but higher IGHV3-30 and IGHV4-34 usage compared with their counterparts expressing IgM, IgG, or IgA (Fig. S3 C), which was also noted in previously published but smaller scale studies (Seifert et al., 2009; Zheng et al., 2004). Furthermore, IgD-GC B cells, IgD-ME B cells, and IgD-PCs used substantially more IGHJ6 than their counterparts expressing IgM, IgG, or IgA (Fig. 4 G). Of note, IGHJ6 is usually counterselected by GC B cells along with IGHV4-34 due to intrinsic autoreactivity (Zheng et al., 2004; Schickel et al., 2017). Finally, consistent with an earlier finding from tonsillar IgD-GC B cells (Koelsch et al., 2007), the length of H-CDR3 in IgD-ME B cells and IgD-PCs was higher than in their counterparts expressing IgM, IgG, or IgA (Fig. 4 H), which agrees with the enhanced autoreactivity of IgD clones (Duty et al., 2009; Zheng et al., 2004).

Next, we explored the degree of clonal overlap between IgD⁺IgM⁻ subsets. A clone was defined as a cluster of B cell receptor (BCR) sequences utilizing a unique combination of IGHV and IGHJ genes and an identical CDR3 amino acid sequence. IgD-ME B cells clonally overlapped with IgD-GC B cells and IgD-PCs (Fig. 4 I). This overlap was higher than that observed in other tonsillar B cell and PC subsets, except IgG-GC B cells and IgG-PCs, which also displayed a rather high overlap. In addition, IgD-ME B and IgD-GC B cells, as well as IgD-PCs, only rarely clonally overlapped with equivalent B cells or PCs expressing IgM, IgG, or IgA, or with naive or IgD-AN B cells (Fig. 4 I). Thus, tonsillar IgD-ME B cells likely emerge from a GC-driven response that generates clonally affiliated IgD-PCs. This IgD response is isotype-restricted and involves GC retention of IgD specificities encoded by intrinsically autoreactive IGHV and IGHJ genes.

Given that ME B cells can differentiate into PCs in response to extrafollicular cues (Inoue and Kurosaki, 2024; Eisenbarth et al., 2025), we wondered about the spatial relationship between IgD-ME B cells and IgD-PCs. We used a published spatial transcriptomics dataset (Massoni-Badosa et al., 2024) to map tonsillar IgD-PCs expressing an IGHD⁺IGHM⁻MME⁻CD19⁺CD27⁺CD38⁺TNFRSF17 (BCMA)⁺ gene signature designed on the basis of published data (Chen et al., 2009; Shan et al., 2018; Massoni-Badosa et al., 2024). Similar to IgD-ME B cells, IgD-PCs preferentially occupied epithelial and subepithelial areas of the tonsillar mucosa (Fig. 5, A and B).

Considering that some ME B cells, including a fraction of ABCs, exhibit a prominent PC differentiation–poised activation state (Gao and Cockburn, 2022), we next explored changes in key mRNAs of the gene expression program underpinning PC differentiation (De Silva et al., 2012; Ripperger and Bhattacharya, 2021). When IgD-ME B cells were compared with IgG-ME or IgA-ME B cells, bulk RNA-seq showed decreased expression of BACH2, FOXO1, PAX5, and SPIB genes involved in B cell identity preservation; increased expression of PRDM1 (BLIMP-1) and IRF4 genes required for PC differentiation; increased expression of XBP1 and ERN1 (IRE-1) genes essential for the unfolded protein response (UPR) preceding antibody secretion (Ripperger and Bhattacharya, 2021; Xiong et al., 2022); and increased expression of the JCHAIN gene (see Fig. 2 B), a hallmark of PC differentiation (Kawasaki et al., 2024). Of note, IgD-ME B cells further exhibited increased expression of the SOX5 gene (Fig. 5 C) implicated in ABC differentiation and in a preparatory phase of PC differentiation (Yang et al., 2022; Rakhmanov et al., 2014). Lastly, the IgD-ME 1 and IgD-ME 2 subclusters included scattered cells expressing low to moderate amounts of a complete PC gene signature (Fig. 5 D, top) designed based on differential expression (DE) analysis of a recently published tonsillar scRNA-seq dataset (Massoni-Badosa et al., 2024), including PRDM1, XBP1, and IRF4 (Fig. S4 A). When compared to IgD-ME 2, IgD-ME 1 from the same tonsils was enriched in the above PC gene signature (Fig. 5 D, middle), but expressed less abundant B cell identity genes, including PAX5, BACH2, and FOXO1 (Fig. S4 A).

Consistent with these data, bulk RNA-seq showed that IgD-ME B cells were also enriched in SDC1 (CD138), IL6R, IFNAR2, IL15RA, TLR8, and STAT3 mRNAs (Fig. 5 E), which are implicated in PC differentiation (Chen et al., 2009; Jego et al., 2003). When evaluated by spectral flow cytometry, BLIMP-1, IFNAR2, and TLR8 expression was indeed higher on IgD-ME compared with IgG-ME but not IgA-ME B cells along with BCMA expression (Fig. 5 F), a hallmark of PCs similar to BLIMP-1 expression (Schulz et al., 2025). Consistently, a fraction of IgD-ME but not IgG-ME or IgA-ME B cells expressed BLIMP-1 or BCMA (Fig. S4, B and C), despite displaying the ME-associated CD27⁺CD38^–/low signature shown earlier by our IgD-ME gating strategy.

We then explored whether IgD-ME B cells were enriched in gene programs linked to antibody secretion, including the UPR program (Ripperger and Bhattacharya, 2021; Xiong et al., 2022). GSEA showed that compared with IgG-ME or IgA-ME B cells, IgD-ME B cells were indeed enriched in gene sets linked to protein synthesis, posttranslational modification, and trafficking processes involving the endoplasmic reticulum and Golgi apparatus (Fig. 5 G). This gene signature echoed transmission electron microscopy (TEM) data showing that some IgD-ME B cells displayed PC-associated ultrastructural properties similar to those of IgD-PCs, including abundant rough endoplasmic reticulum (RER) and large perinuclear Golgi apparatus (Fig. 5 H). Confirming their ability to differentiate into IgD-PCs, IgD-ME B cells differentiated in vitro into CD27^highCD38^high plasmablasts/PCs (Fig. 5 I) and secreted IgD upon 6-day stimulation with CD40L and IL-21 (Fig. 5 J), which provide signals mimicking T follicular helper (T_FH) cell help in both GC and extra-GC contexts (Song and Craft, 2024). Thus, a large fraction of the IgD-ME cluster colocalizes with IgD-PCs within subepithelial and epithelial districts of the tonsillar mucosa and encompasses clones with enhanced PC but decreased B cell identity gene properties, which may underlie an increased propensity to differentiate into IgD-PCs.

Knowing that naive B cells undergo IgM-to-IgD CSR and secrete IgD in response to signals from BCR, CD40, TACI, TLRs, and cytokine receptors in vitro (Chen et al., 2009; Xu et al., 2022), we explored the impact of these signals in vivo by studying the IgD response in selected primary immunodeficiency (PID) patients with rare inborn errors of immunity. Compared with most controls, serum IgD was below the limit of detection in the following patients (Fig. 6 A and Table S2): one X-linked agammaglobulinemia (XLA) case with hypomorphic mutation of the BTK gene encoding the BCR-associated BTK kinase (Yang et al., 2022); two hyper-IgE syndrome (HIES) cases with DOCK8 mutations disrupting DOCK8-mediated TLR signaling via MyD88 (Jabara et al., 2012); 12 common variable immunodeficiency (CVID) cases with TNFRSF13B (10 cases) or NFKB1 (two cases) mutations disrupting TACI- and NF-κB p50-mediated CSR and PC differentiation (Peng et al., 2023); three hyper-IgM type 1 syndrome (HIGM1) cases with CD40L mutations disrupting CD40L-mediated CSR (Cunningham-Rundles and Ponda, 2005); two MyD88 deficiency cases with MYD88 mutations disrupting TLR-mediated MyD88-driven IgM production (Maglione et al., 2014); and one late-onset combined immunodeficiency (LOCID) case with hypomorphic RAG2 mutations disrupting RAG2-orchestrated BCR genesis and lymphopoiesis (Schmitt and Cooper, 2021; Dorna et al., 2019).

In the following patients, serum IgD was as abundant as in age-matched controls (Fig. 6 A and Table S2): 10 HIES cases with STAT3 mutations disrupting classical ME B cell differentiation (Speckmann et al., 2008; van de Veen et al., 2019); three congenital disorders of glycosylation cases with POMT1 (two cases) or GMPPB (one case) mutations disrupting protein O-mannosyltransferase 1– and GDP-mannose pyrophosphorylase B–mediated glycosylation (Jaeken and Matthijs, 2007); and one CVID case with an NFKB2 mutation disrupting NF-κB p52-mediated B cell survival (Peng et al., 2023).

Lastly, serum IgD was increased compared with age-matched healthy controls in the following patients (Fig. 6 A and Table S2): two HIES cases with PGM3 or CARD11 mutations disrupting phosphoglucomutase 3 (PGM3)– and CARD11-mediated glycosylation and BCR or T cell receptor (TCR) signaling (Minegishi, 2021); and twelve hyper-IgD syndrome (HIDS) cases with biallelic MVK mutations disrupting mevalonate kinase–regulated inflammasome signaling (Ozen and Bilginer, 2014). Pooled CVID cases or pooled XLA, HIGM1, MyD88 deficiency, and LOCID cases had less serum IgD than controls (Fig. 6 B), whereas pooled PGM3-HIES and CARD11-HIES cases or HIDS cases had more serum IgD than controls, with one PGM3-HIES case and two HIDS cases showing the highest IgD (Fig. 6 C).

Consistent with the presence of stronger IgD response early in life (Renz et al., 2018; Blanco et al., 2018), tonsillar IgD-PCs and, to a lesser extent, IgD-ME-B cells showed a nonstatistically significant but clear trend toward a reduction following 2–3 years of age (Fig. 6, D and E), which suggests a direct relationship between IgD-ME-B cells, IgD-PCs, and serum IgD. Similarly, the frequency of circulating IgD-ME-B cells was higher in PGM3-HIES and CARD11-HIES but not STAT3-HIES cases (Fig. 6, F–H), who also had a tendentially higher serum IgD than healthy controls (shown earlier). In the PGM3-HIES case, circulating IgD-PCs were as frequent as those from a healthy control with an unusually high frequency of IgD-ME B cells (Fig. 6, F–H).

Finally, bulk RNA-seq showed that compared with IgG-ME or IgA-ME B cells, IgD-ME B cells expressed more genes encoding POMT1, PGM3, and GMPPB (Fig. 6 I), which might be involved in the glycosylation of IgD-regulating molecules and/or IgD itself. In summary, IgD-ME B cell–centered IgD responses may be (1) positively regulated by NF-κB p50 signals from BCR, CD40, TACI, and/or TLRs; (2) constrained by PGM3-mediated glycosylation-dependent or CARD11-mediated BCR- or TLR-dependent signals; (3) independent of STAT-3 signals; and (4) more abundant in the first years of life, which are central to the generation of tolerance against innocuous environmental antigens.

We then explored whether IgD-ME B cells could serve as IgD-PC precursors. First, we mined published scRNA-seq data (Massoni-Badosa et al., 2024) to perform trajectory inference analysis. UMAP obtained as detailed in Methods showed two IgD-GC B cells, IgD-ME B cells, and IgD-PC subclusters, subcluster 1 being larger than subcluster 2 (Fig. 7 A, top). Of note, IgD-ME 1 was closely and sequentially connected to IgD-GC 1 and IgD-PC 1, whereas IgD-ME 2 was mostly linked to IgD-PC2, which in turn showed a dual link to IgD-GC 2 (Fig. 7 A, top). Compared with IgD-ME 2, IgD-ME 1 expressed more ABC-associated (i.e., FCRL4, TBX21) and PC-associated (i.e., TNFRSF13B) mRNAs, but less B cell identity (i.e., FOXO1, PAX5) and ABC differentiation–associated ZEB2 (Fig. 7 A, bottom) mRNAs. Compared with IgD-GC 2, IgD-GC 1 expressed less GC-associated (i.e., MME, CD38, BCL6), but more ME B cell– and PC-associated CD27 mRNAs (Fig. 7 A, bottom). Lastly, IgD-PC 1 expressed less ERN1 than IgD-PC 2 (Fig. 7 A, bottom). Thus, IgD-ME B cells, IgD-GC B cells, and IgD-PCs encompass distinct subpopulations seemingly organized in three developmental pathways that possibly underpin distinct IgD responses.

Given our earlier findings that IgD-ME B cells are more mutated that other class-switched ME B cells and considering that compared with IgG-GC B cells, IgD-GC B cells are enriched in transcripts encoding key components of the SHM machinery, we asked whether these unique properties were preferentially associated with IgD-GC 1 or IgD-GC 2. We found that compared with IgG-GC B cells, IgD-GC 2 was more enriched in RAD52, which is essential for IgM-to-IgD CSR (Xu et al., 2022), as well as EXO1, PMS1, XRCC4, MSH2, MSH6, and MLH1 (Fig. 7 B), which encode key components of the MMR machinery (Odegard and Schatz, 2006; Casellas et al., 2016). Thus, IgD-GC 2 may be more intensely engaged in SHM compared with IgD-GC 1.

To further understand the developmental history of IgD-ME B cells, we built phylogenetic trees from our BCR sequencing (BCR-seq) data by assigning V_HDJ_H clones into lineage clusters as detailed in Methods. Lineage trees were grouped according to their isotype composition, and all groups containing at least five trees per donor were considered to ensure statistical robustness. In this way, we generated groups of lineage trees encompassing IgD alone, IgM alone, combined IgG-IgA, or combined IgM-IgG-IgA. To characterize the developmental history of the cells from these lineage trees, we calculated the normalized tree height, i.e., the distance from the germline to the furthest tip in the tree divided by the number of non-germline sequences per tree. IgD lineage trees were higher than trees with IgM alone, both IgG and IgA, or a combination of IgM, IgG, and IgA (Fig. 7 C), suggesting that IgD class-switched B cells have a longer and more complex developmental history.

We also assessed the mutation activity in IgD class-switched lineage trees by calculating the average mutation count per tree edge. IgD lineage trees had a higher average mutation count per tree edge compared with trees including IgM alone, IgG and IgA, or IgM, IgG, and IgA (Fig. 7, C and D). As reflected by their smaller scale bar, IgD lineage trees accumulated mutations at a higher rate compared with IgM, IgA, or IgG trees (Fig. 8 A; and Fig. S4, D and E). Thus, IgD class-switched B cells form large clonal families entailing GC B cells, ME B cells, and PCs that do not further switch to downstream isotypes, but rather engage in intense SHM.

To identify IgD-PC precursors, we converted IgD trees into pairwise distance matrices, selected the nearest (minimal distance) non-PC neighboring sequence for each PC-originated sequence, and calculated the average nearest neighbor proportion (%) for each donor. IgD-ME B cells were the nearest neighbors to IgD-PCs (Fig. 8 B), representing their most frequent precursors. This trend could also be explained by the higher mutation rate of GC B cells, causing them to further diverge from the germline compared with IgD-ME B cells and IgD-PCs at the time of tonsillectomy. A switch proportion (SP) test performed on IgD lineage trees showed that the IgD-ME B cell-to-IgD-PC transitions had a higher SP statistic (δ) value in three of four donors compared with IgD-GC B cell-to-IgD-PC transitions (Fig. 8 C). Thus, IgD-PCs preferentially emerge from IgD-ME B cell subsets rather than IgD-GC B cells.

To verify whether class-switched IgD antibodies reacted to common nasopharyngeal antigens, we generated a set of recombinant monoclonal antibodies (rmAbs) from tonsillar IgD-ME and IgD-GC B cells. To establish reactivity thresholds, control IgG rmAbs were generated from circulating IgG-ME B cells specific to the receptor-binding domain of severe acute respiratory syndrome coronavirus 2 (De Campos-Mata et al., 2024). All these rmAbs encompassed Cγ1 to facilitate downstream analysis (Koelsch et al., 2007). As shown by ELISA, some IgD rmAbs from IgD-ME or IgD-GC B cells recognized multiple airborne and food antigens, including dust mite (Der p 1), pollen (Bet v 1, Pla a 3), milk (whole casein, α-s-casein, β-lactoglobulin, α-lactalbumin), and/or egg (ovalbumin) antigens (Fig. 9, A and B; and Fig. S5 A). Some IgD rmAbs from IgD-ME B cells also recognized viral (echovirus) and/or fungal (chitosan), but not algae-associated (laminarin), antigens (Fig. 9 B). In addition, some IgD rmAbs from IgD-ME B cells showed mild reactivity to single-stranded DNA (ssDNA) and/or double-stranded DNA (dsDNA) (Fig. 9 B). In addition to IgD rmAbs from tonsillar IgD-ME B cells, IgG rmAbs from IgG-ME cells inhabiting identical tonsils recognized similar environmental antigens (Fig. S5 B), showing that the humoral response to these antigens is not restricted to IgD, but rather involves multiple antibody classes.

Bacterial flow cytometry showed that IgD rmAbs from IgD-ME or IgD-GC B cells further recognized a cultured nasopharyngeal isolate, Streptococcus mutans, in addition to gut Bacillus cereus and Ruthenibacterium lactatiformans isolates. Some IgD rmAbs also recognized Bifidobacterium longum, but none recognized Escherichia coli and Roseburia intestinalis (Fig. 9, C and D; and Fig. S5, C–E). Of all IgD rmAbs from IgD-GC B cells, rmAbs 10, 47, and 48 showed polyreactivity, rmAbs 10 and 47 exhibiting detectable binding to most microbial antigens tested (Fig. S5, D and E). While IgD rmAbs 10 and 48 were encoded by IGHV3-30, rmAb 47 was encoded by IGHV4-34 (Table S3), two intrinsically autoreactive IGHV genes (Zheng et al., 2004; Schickel et al., 2017; Kobie et al., 2015; Ohm-Laursen et al., 2018; Pugh-Bernard et al., 2001).

Unlike IgD rmAb 10, which was unmutated, IgD rmAbs 47 and 48 were highly mutated (Fig. S5 F and Table S3), suggesting that SHM contributes to the polyreactivity and autoreactivity of some but not all IgD clones. Accordingly, germline revertants of rmAbs 47 and 48 generated by replacing mutated nucleotide residues with native residues from their putative germline sequence (Fig. S5 F) showed decreased binding to B. cereus and S. mutans, as well as airborne, food, viral, and self-antigens compared with mutated wild-type controls (Fig. 9, E–G). To evaluate the impact of Cδ and Cγ1 on the reactivity of IgD rmAbs, highly polyreactive IgD rmAbs 47 and 48 and mildly polyreactive IgD rmAbs 70 and 98 engineered to encompass Cδ were found to have grossly comparable reactivities to equivalent IgD rmAbs encompassing Cγ1 (Fig. 9 H). However, compared with Cγ1-expressing rmAbs, Cδ-expressing rmAbs displayed higher reactivity to Der p1 and α-s casein but lower reactivity to dsDNA (Fig. 9 H). Of note, these differences will need to be reevaluated in a larger group of rmAbs.

To confirm the above results in native secreted IgD, we enriched IgD polyclonal antibodies (pAbs) from serum as detailed in Methods. Similar to tonsillar IgD rmAbs, serum IgD pAbs from healthy controls or HIDS and HIES patients recognized airborne and food antigens in addition to ssDNA (Fig. 9 I and Fig. S5 G). Consistent with evidence suggesting the parallel induction of IgD and IgE responses by T helper type 2 (T_H2) cell signals and the constraining impact of IgD on allergen-driven IgE-mediated allergic responses (Shan et al., 2018), IgE pAbs from some HIES patients were found to bind to the environmental and food antigens Der p 1 and α-s-casein like IgD pAbs did (Fig. 9 I and Fig. S5, G and H). Thus, the IgD response may target some commensals in addition to innocuous environmental antigens, including potential allergens.

We then tested whether endogenous IgD pAbs could recognize endogenous commensal bacteria. Given the likely origin of circulating IgD-PCs from aerodigestive inductive sites (Chen et al., 2020; Ohm-Laursen et al., 2018; Ohm-Laursen et al., 2021), we first interrogated the reactivity of circulating IgD pAbs for commensal isolates. Serum IgD pAbs were found to recognize a cultured S. mutans nasopharyngeal isolate in addition to intestinal B. cereus and R. lactatiformans isolates (Fig. 10 A). Next, we studied sinonasal IgD responses to the microbiota from healthy donors and patients with various disorders (Table S4). Confocal imaging visualized sinonasal IgD⁺IgM⁻IgA⁻ IgD-PCs in addition to control IgD⁻IgM⁻IgA⁺ IgA-PCs and IgD⁻IgM⁺IgA⁻ IgM-PCs (Fig. 10 B).

Furthermore, bacterial flow cytometry determined that ∼3 and 48% of the total microbiota from sinonasal fluids were coated by IgD and IgA, respectively (Fig. 10 C). When only the IgD/IgA-coated bronchoalveolar microbiota was considered, about 5% of these bacteria were targeted by IgD (Fig. 10 C). The majority of these IgD-coated bacteria were cotargeted by IgA, only a minority being coated by IgD alone (Fig. 10 C). About 95% of the IgD/IgA-targeted sinonasal microbiota was coated by IgA alone (Fig. 10 C). Consistently, ELISAs detected sinonasal free IgD pAbs in addition to larger amounts of sinonasal free IgA pAbs (Fig. 10 D).

Knowing that the lung mucosa includes IgD clones from IgD-PCs (Ohm-Laursen et al., 2021), we also studied IgD responses in bronchial mucosa samples from cancer patients undergoing bronchoscopy (Table S5). Confocal imaging detected bronchial IgD⁺IgM⁻IgA⁻ IgD-PCs together with IgD⁻IgM⁻IgA⁺ IgA-PCs (Fig. 10 E). Moreover, bacterial flow cytometry determined that about 3 and 14% of the total bronchoalveolar microbiota were coated by IgD and IgA, respectively (Fig. 10 F). When only the IgD/IgA-targeted bronchoalveolar microbiota was considered, about 32% of this microbiota was coated by IgD (Fig. 10 F). The majority of these IgD-coated bacteria were targeted by IgD alone, with a smaller fraction being cotargeted by both IgD and IgA (Fig. 10 F). The remaining 68% of the IgD/IgA-targeted bronchoalveolar microbiota was coated by IgA alone (Fig. 10 F). Lastly, bacterial flow cytometry evidenced that the sinonasal and bronchoalveolar tracts included comparable frequencies of IgD- or IgA-coated bacteria, the frequency distribution being very broad (Fig. 10, G and H). As shown by ELISA, sinonasal IgD pAbs recognized allergens like Der p 1 and α-s-casein in addition to commensal bacteria (Fig. 10 I), which echoes the reactivity of tonsillar IgD rmAbs and serum IgD pAbs.

In summary, tonsils include a heterogeneous population of IgD-ME B cells that differentiate into IgD-PCs (Fig. 10 J). These IgD-PCs would migrate to both local and distal effector sites, including the sinonasal and bronchial mucosae, and secrete IgD, which binds to a fraction of the microbiota in addition to common nasopharyngeal antigens, including allergens. Aside from supporting a mutation-intensive pathway involved in the polyreactivity of secreted IgD, sustained exposure to these antigens may account for the ABC-like properties of at least some IgD-ME B cells.

Here, we have characterized a heterogeneous IgD-ME B cell population that inhabits epithelial, subepithelial, and interfollicular areas of the tonsillar mucosa and displays prominent ABC-like traits in addition to properties usually associated with PCs. IgD-ME B cells emerge from a SHM-intensive GC-driven program that gives rise to clonally affiliated IgD-PCs. The IgD secreted by these cells targets a fraction of the respiratory microbiota in addition to various environmental antigens, including some airborne and food allergens. By providing new insights into a hitherto elusive mucosal B cell population, our findings position IgD-ME B cells at the center of a largely neglected IgD response presumably aimed at enhancing homeostasis in the respiratory mucosa.

IgD-PCs were first identified in tonsils more than two decades ago (Arpin et al., 1998; Liu et al., 1996). By showing that IgD from nasopharyngeal IgD-PCs targets environmental antigens, including allergens, the present findings echo recent evidence showing that in allergic patients, increased allergen-specific IgD production correlates with a lower risk of anaphylaxis, higher likelihood of naturally outgrowing allergy, and increased tolerance in response to allergen immunotherapy (Shan et al., 2018; Boonpiyathad et al., 2020; Suprun et al., 2020; Itoh et al., 2021; Satitsuksanoa et al., 2025). These human findings along with additional mouse results (Shan et al., 2018) raise the possibility that IgD responses contribute to the establishment and/or maintenance of immune tolerance. Clarifying the biology of secreted IgD may help understand its impact on mucosal homeostasis, immunity, and inflammation (Chen et al., 2023; Itoh and Ohshima, 2023).

Secreted IgD may enhance tolerance in the context of a protective-regulatory brand of T_H2 cell–mediated humoral responses in which IgD antibodies from IgD-PCs interact with basophils and mast cells through a receptor complex encompassing galectin-9 and CD44 (Shan et al., 2018). By crosslinking this membrane-bound IgD, secondary allergen exposure would elicit basophil and mast cell release of IL-4, IL-5, and IL-13, which may promote allergen clearance by amplifying T_H2 cell–induced production of IgG and IgE by allergen-specific B cells (Shan et al., 2018). Concomitantly, signals from receptor-bound IgD would attenuate allergen-induced inflammation by interfering with IgE-mediated FcεRI-driven basophil and mast cell degranulation and by eliciting regulatory T cell expansion (Shan et al., 2018).

By characterizing the tissue topography, clonal history, molecular landscape, antigen reactivity, and in vivo requirements of tonsillar IgD-ME B cells, we have unveiled hitherto unappreciated facets of mucosal IgD responses. In the aerodigestive mucosa, IgD-ME B cells shared a preferential epithelial and subepithelial topography with IgD-PCs, suggesting a close relationship between these IgD class-switched cells. Consistently, clonotypic analysis and transcriptome trajectory and phylogenetic topology studies identified IgD-ME rather than IgD-GC B cells as the main direct precursors of IgD-PCs. These last showed a trend toward a reduction following early childhood, which is when tolerance to environmental antigens is acquired (Renz et al., 2018). Unlike IgD-PCs, nasopharyngeal IgG-PCs seemed to have a closer and possibly more direct affiliation to IgG-GC B cells than to IgG-ME B cells when the degree of clonotypic overlap was evaluated, which might reflect the more frequent involvement of IgD-ME B cells in anamnestic IgD responses to common mucosal antigens.

Besides exhibiting extensive clonal affiliation with IgD-PCs, the IgD-ME B cell population was enriched in some ABC-like phenotypic properties compared with IgG-ME or IgA-ME B cells, including enhanced CD11c, CD69, CD95, FCRL4, and CXCR3 expression. These ABC-like traits were not homogeneously expressed, but rather reflected the presence of two IgD-ME subclusters encompassing B cells likely lying at distinct stages of ABC activation and differentiation, as suggested by their differential CD69, FCGR2B, LTB, TNFRSF13B, TOX, SOX5, SSR4, and GSN expression. At the transcriptome level, the larger subcluster 1 was characterized by increased ribosomal biogenesis and enhanced protein translation, which likely reflects a stage preparatory to sustained antibody production. The transcriptome of the smaller subcluster 2 hinted at extensive GTPase signaling–induced cytoskeleton changes, which are usually required for increased locomotion, cell shape remodeling, and organelle trafficking.

These two IgD-ME subclusters further included some B cells at distinct stages of an early commitment to PC differentiation, as suggested by their DE of PRDM1, XBP1, IRF4, PAX5, BACH2, and FOXO1, subcluster 1 being more enriched in a PC-associated gene signature compared with subcluster 2. These properties correlated with the presence of minute BLIMP-1⁺ and BCMA⁺ fractions within the IgD-ME B cell population, which underwent PC differentiation and IgD secretion upon in vitro stimulation by CD40L and IL-21. As suggested by our analysis of PID patients with CD40L deficiency and by published properties of both T_FH and ME B cells (Inoue and Kurosaki, 2024; Song and Craft, 2024), T_FH cell–derived signals from CD40L may be implicated in both the emergence of IgD-ME B cells from IgD-GC B cells and the extra-GC differentiation of IgD-ME B cells into IgD-PCs. Consistent with more findings published earlier (Chen et al., 2009; Xu et al., 2022) and our present results from PID patients, additional BCR, TLR, and TACI signals from local commensals, as well as myeloid and/or epithelial cells, may be involved in the extrafollicular differentiation of at least some IgD-ME B cells into IgD-PCs (Inoue and Kurosaki, 2024; Song and Craft, 2024; Eisenbarth et al., 2025).

These intricate signaling requirements may underlie the complex transcriptional trajectories of the IgD-ME 1 and 2 subclusters, which engendered two distinct IgD-PC 1 and IgD-PC 2 developmental programs. While IgD-ME 1 resided along an IgD-GC 1 > IgD-ME 1 > IgD-PC trajectory, the smaller IgD-ME 2 subcluster followed an IgD-ME 2 > IgD-PC 2 trajectory. The first trajectory was characterized by a more robust commitment of GC-derived IgD-ME B cells into the PC differentiation program. Accordingly, compared with IgD-ME 2 cells, IgD-ME 1 expressed more PC-associated mRNAs. The second trajectory could reflect the recall of a distal pool of IgD-ME B cells that may include recirculating IgD-ME cells involved in an anamnestic IgD-PC response. A third IgD-GC 2 > IgD-PC 2 trajectory echoes published findings, showing that some GC B cells enter early stages of PC differentiation due to the attainment of a sufficiently high affinity for antigen (Kräutler et al., 2017; Ise et al., 2018; Inoue and Kurosaki, 2024). Of note, compared with IgD-GC 1, IgD-GC 2 was enriched in SHM-related, B cell identity (e.g., FOXO1, PAX5), and GC identity (e.g., MME, CD38, BCL6) mRNAs, which could implicate an earlier differentiation stage of IgD-GC 2 compared with IgD-GC 1 along the GC pathway.

In agreement with the general state of active readiness of ME B cells toward anamnestic antibody responses, some IgD-ME B cells displayed transcriptional and ultrastructural properties reminiscent of PCs, including abundant RER and Golgi network. Considering that both IgD-ME B cells and ABCs are enriched in PC-related gene signatures, such as JCHAIN (Lau et al., 2017; Wang et al., 2018; Nellore et al., 2023), the differentiation of IgD-PCs from IgD-ME B cells might involve intertwined ABC and PC developmental programs. The subsequent navigation of IgD-PCs to subepithelial and epithelial areas might implicate chemotactic fields orchestrated by CXCR3 via its ligand CXCL10. As suggested recently (Kawasaki et al., 2024), additional chemotactic clues might be provided by CXCL-like signals from the J-chain, which usually promotes the transcytosis of polymeric IgA and IgM across epithelial cells. The chemotactic properties of the J-chain could explain why IgD-PCs abundantly express this protein despite releasing monomeric rather than polymeric IgD (Chen et al., 2020; Brandtzaeg and Korsrud, 1984).

By pointing to the involvement of CD40, TACI, and TLRs in IgD responses in vivo, our study of PID patients confirmed earlier findings showing the positive impact of these signals in IgD CSR and production (Chen et al., 2009; Xu et al., 2022). In concert with the IFNAR (Chen et al., 2009; Xu et al., 2022; Jego et al., 2003), CD40, TACI, and TLRs could, in principle, trigger the differentiation of IgD-ME B cells into IgD-PCs through a pathway further involving CXCR3, an ABC-associated CXCL10 receptor known to contribute to PC differentiation (Xu et al., 2012).

By showing conserved IgD production in STAT3-HIES patients, our data suggest that similar to IgE responses, IgD responses develop independently of STAT-3, which instead is central to IgG and IgA production (Speckmann et al., 2008; van de Veen et al., 2019). Unlike STAT3-HIES and DOCK8-HIES cases, which exhibited normal or absent IgD responses, PGM3-HIES and, to a lesser extent, CARD11-HIES cases showed increased IgD responses, suggesting that protein glycosylation and components of BCR or TCR signaling may have a constraining impact on these responses. The unavailability of a sufficient number of pure B cells from these rare PID cases prevented us from establishing the B cell–intrinsic or B cell–extrinsic impact of STAT3, PGM3, or CARD11 on IgD production.

Although enriched in many ABC properties, including enriched ITGAX, ITGB2, ITGB7, TBX21, TNFRSF13B, SLAMF7, IL10RA, CD86, and CD151 expression (Yang et al., 2022; Winslow and Levack, 2025), IgD-ME B cells lacked other ABC traits, like enriched EZR expression. This and other differences (Sutton et al., 2021; Holla et al., 2021) may relate to the tissue and isotype specificity of IgD-ME B cells. Regardless, our identification of ABC-like cells within the IgD-ME B cell population echoes recent data, indicating that ABCs are a physiological component of the humoral response (Sutton et al., 2021; Gao and Cockburn, 2022). Similar to ABCs (Sutton et al., 2021; Gao et al., 2024; Gao and Cockburn, 2022), IgD-ME B cells expressed hypermutated IGHV genes. Their analysis revealed the complex developmental history of IgD-ME B cells, which transited through a GC reaction permitting the persistence of intrinsically autoreactive genes usually counterselected by GCs, including IGHV4-30, IGHV4-31, and, although not significantly, IGHV4-34 (Pugh-Bernard et al., 2001). Similarly, earlier small-scale studies show IGHV4-34 enriched in IgD class-switched B cells (Seifert et al., 2009; Zheng et al., 2004), raising the possibility that IgD responses retain some degree of autoreactivity to maximize their function. This residual autoreactivity may be instrumental for IgD to cross-target highly conserved epitopes shared by some environmental and commensal antigens with autologous antigens (Schickel et al., 2017; Leviatan et al., 2022).

Accordingly, the class-switched IgD repertoire was enriched not only in autoreactive IGHV4 family members, but also in some IGHV3 genes usually involved in antimicrobial antibody responses. Of note, two donors showed IgD class-switched B cells enriched in IGHV3-30, which encodes specificities targeting microbial carbohydrates (Bryson et al., 2016). Concomitantly, IgD class-switched B cells exhibited lesser usage of IGHV3-13, IGHV3-48, IGHV3-7, and (although not significantly) IGHV3-23, which encode specificities targeting a broad range of antigens (Wu et al., 2010). This peculiar IGHV composition was combined with the lack of IgD clonal affiliation with IgM, IgG, and IgA, which echoes studies published earlier (Zheng et al., 2004) showing the limited capacity of IgD class-switched B cells to further switch to downstream isotypes. The isotype-restricted clonal architecture of IgD responses suggests that such responses target a restricted and possibly unique set of antigens.

Most IgD antibodies coated a fraction of the sinonasal microbiota cotargeted by IgA, the latter being a dominant mucosal antibody (Cerutti, 2008). A larger fraction of IgD-only–coated microbiota was detected in the bronchial mucosa (Goeteyn et al., 2023), which might reflect the need by this mucosal district of controlling commensals through a noninflammatory, polyreactive, and flexible antibody like IgD (Chen et al., 2020; Übelhart et al., 2015). The flexibility of the antigen-binding Fab domain of IgD mostly stems from the uniquely long hinge region of its constant Fc domain, which may facilitate the heteroligation of biochemically distinct antigenic epitopes, including low-density epitopes (Chen et al., 2020; Übelhart et al., 2015).

Despite showing more mutations than IgG-ME B cells, IgD-ME B cells exhibited a comparable in vivo proliferation history, raising the possibility that their larger mutational load stems from the intrinsically hyperactive SHM machinery rather than more numerous GC transits. This possibility is consistent with the finding that compared with IgG-GC B cells, IgD-GC B cells, and particularly a subcluster of them, expressed more abundant transcripts encoding RAD52, which is essential for IgM-to-IgD CSR (Chen et al., 2009; Xu et al., 2022). Compared with IgG-GC B cells, IgD-GC B cells were also enriched in transcripts encoding additional components of the MMR and BER pathways of SHM (Odegard and Schatz, 2006; Casellas et al., 2016), the former pathway being more represented that the latter in IgD-GC compared with IgG-GC B cells. The reasons behind these distinctive features of IgD-GC B cells remain elusive, but might be linked to the unique signaling properties of IgD-BCRs. Indeed, in IgD-GC B cells, IgD-BCRs operate without any counterbalance from IgM-BCRs and, compared with these last, preferentially bind multivalent antigens with repetitive epitopes (Bellenot et al., 1988; Übelhart et al., 2015).

By showing that the reactivity of IgD antibodies to commensal and environmental antigens decreased in germline IgD revertants, our findings indicate that SHM is central to the GC-driven nasopharyngeal program underlying the induction of class-switched IgD. In this regard, some hypermutated IgD antibodies displayed a reactivity that spanned across different kingdoms of living organisms, with some IgD being able to recognize plant, animal, fungi, and protist antigens in addition to antigens from food, viruses, and airborne particles. Cross-kingdom poly- and cross-reactivity could be key to enhance noninflammatory immunity and tolerance by means of IgD-mediated targeting of highly conserved epitopes shared by environmental and commensal antigens from distant species (Schickel et al., 2017; Leviatan et al., 2022).

The widespread availability and common nature of IgD-reactive antigens are suggested by the clustering together of IGHV genes from IgD-GC B cells, IgD-ME B cells, and IgD-PCs from four different donors. The involvement in this process of superantigen-binding sites from Cδ (Chen et al., 2020; Seifert et al., 2009; Jendholm et al., 2008; Thurner et al., 2020) seems unlikely, as replacement mutations of IGHD genes from IgD-GC B cells, IgD-ME B cells, and IgD-PCs preferentially targeted H-CDR rather than H-FWRs, which suggests antigen-driven B cell selection in the GC (Neumeier et al., 2022; Yaari et al., 2015). The little or no role of Cδ superantigen-binding sites was confirmed by the largely comparable reactivity of tonsillar and serum IgD antibodies, including tonsillar IgD expressing either Cγ1 or Cδ.

A key strength of the present work relates to the integrated, multimodal analysis across tissue compartments—combining spatial localization, transcriptomics, repertoire phylogenetics, and recombinant antibody testing—to link IgD-ME B cell biology to IgD-PC output and antigen specificity. A limitation relates to the limited yield and viability of IgD-ME B cells, which hindered a systematic functional analysis of IgD-ME B cells, including their response to BCR, TACI, and TLR signals. An additional limitation relates to the limited power of studies involving rare PID patients, which was partly compensated by pooling together samples when possible. Furthermore, the composition of IgD-coated commensal communities remains unknown, although ongoing work will soon reveal it.

In summary, we have characterized a hitherto elusive GC-imprinted IgD-ME B cell subset inhabiting epithelial and interfollicular areas of the aerodigestive mucosa. This heterogeneous IgD-ME population shows massive clonal affiliation with IgD-PCs, which secreted IgD antibodies to both a fraction of the respiratory microbiota and common environmental antigens, including allergens. In light of recent findings (Shan et al., 2018; Boonpiyathad et al., 2020; Suprun et al., 2020; Itoh et al., 2021; Satitsuksanoa et al., 2025), our data raise the possibility that IgD responses contribute to mucosal homeostasis, including tolerance to common environmental antigens. Should that be the case, treatment strategies harnessing the IgD-ME B cell–IgD-PC axis might help attenuate inflammatory disorders of the respiratory tract, including allergy.

Tonsil samples were collected from donors tonsillectomized because of tonsillar hypertrophy. Heparinized blood samples were acquired from healthy adult donors. Splenic samples were obtained from deceased healthy donors or individuals subjected to posttraumatic splenectomy as described earlier (Puga et al., 2011). Sinonasal biopsies and mucosal swabs were obtained from otorhinolaryngological patients. BALs were obtained from patients undergoing bronchoscopy. The Ethical Committee for Clinical Investigation of the Hospital del Mar Research Institute approved the use of blood and tissue samples (CEIC-IMIM 2011/4494/I, 2014/5892/I and 2022/10464/I). Fresh samples and formalin-fixed or paraffin-embedded tissue sections were obtained from the MarBiobanc tissue repository with patient-signed informed consent. Fresh sinonasal biopsy samples, mucosal swabs, and BAL fluids (Tables S4 and S5) were obtained from the Otorhinolaryngology and Pulmonology Services at Hospital del Mar (Barcelona, Spain) during routine diagnostic procedures with patient-signed informed consent (CEIC/IMIM 2022/10464/I). Archived serum samples were obtained from biobanks at Hospital Sant Joan de Déu (Barcelona, Spain), Hospital Clinic (Barcelona, Spain), Mount Sinai Hospital (New York, NY, USA), Spedali Civili di Brescia (Brescia, Italy), and Freiburg University Hospital (Freiburg, Germany) with respective patient-signed informed consents.

Tonsillar mononuclear cells were obtained by perfusion of fresh tissue specimens with sterile phosphate buffer solution (PBS) (Biowest). Peripheral blood mononuclear cells (PBMCs) were isolated from heparinized blood samples by separation on a Ficoll-Hypaque gradient (GE Healthcare). Fresh spleen samples were enzymatically digested for 40 min at 37°C in Hank’s balanced salt solution (Lonza) supplemented with 1 mg/ml collagenase IV (Thermo Fisher Scientific), 50 ng/ml DNase (New England Biolabs), and 0.5% human serum (Sigma-Aldrich). Splenic mononuclear cells were separated on a Ficoll-Hypaque gradient. All mononuclear cell suspensions were aliquoted and cryopreserved in fetal bovine serum (FBS) (Gibco) supplemented with 10% dimethyl sulfoxide (Merck) until the time of analysis. Sinonasal biopsies were included in the paraffin.

Single-cell suspensions of tonsillar, splenic, or PBMCs were incubated at 4°C in sterile staining buffer comprised of PBS at pH 7.4 supplemented with 5 g/L bovine serum albumin (BSA) (Sigma-Aldrich), 2 mM ethylenediaminetetraacetic acid (Merck), and 2 μl Fc receptor blocking reagent (Miltenyi) for every 10⁷ cells. After 10-min incubation, cell suspensions were stained at 4°C for 30 min with a conventional or spectral B cell and antibody cocktail (Table S6). For conventional flow cytometry, 4′-6-diamindine-20-phenylindole (DAPI) (Sigma-Aldrich) was used to exclude dead cells, and cells were acquired with BD LSRFortessa Cell Analyzer (BD Biosciences). In spectral flow cytometry, cells were acquired with Aurora Cell Analyzer (Cytek). Dead cells were excluded using LIVE/DEAD Fixable Yellow Dead Cell Stain Kit (Invitrogen). All flow cytometry data were analyzed using FlowJo version 10 software (TreeStar). t-SNE projections from spectral flow cytometry data were generated using default settings from the FlowJo built-in t-SNE function. Clusters were defined based on the expression of CD10, CD38, CD27, Igλ, CD11c, CD43, CD21, CD95, CXCR3, and CD69 using the XShift clustering algorithm (Samusik et al., 2016) and projected onto t-SNE plots for visualization.

Tonsil mononuclear cell suspensions were incubated for 10 min at 4°C with 2 μl Fc-blocking reagent for every 10⁷ cells in sterile staining buffer. Cells were then stained with a cell sorting–specific antibody cocktail (Table S6). Cell suspensions were filtered through a 35-µm cell strainer, prior to sorting with a FACSAria II (BD Biosciences). DAPI was used to exclude dead cells. B cells were sorted into tubes containing complete RPMI 1640 (Biowest) and immediately used for downstream analysis. Discrete B cell or PC populations were selected according to forward scatter (FSC) and side scatter (SSC) parameters, as well as expression of specific surface molecules. Only IgA-ME and IgG-ME B cells expressing CD27 were included in functional, transcriptional, and TEM studies. Repertoire and phenotypic studies included both dominant CD27⁺ and small CD27⁻ fractions of IgA-ME and IgG-ME B cells (Ehrhardt et al., 2005; Ehrhardt et al., 2008).

FACSorted tonsillar naive, IgD-ME, IgM-ME, IgA-ME, and IgG-ME B cells were seeded (5 × 10⁴/well) in U-bottomed 96-well plates (Corning), and cultured for 6 days in complete RPMI 1640 supplemented with 10% FBS, 10 U/ml penicillin, 10 U/ml streptomycin (Biowest), and 2 mM L-glutamine (Biowest) with or without 100 ng/ml megaCD40L (Enzo Life Sciences) and 500 ng/ml IL-21 (PeproTech). Cells were microscopically monitored daily and counted using a Neubauer chamber. On day 6, the well contents were transferred into new tubes and centrifuged at 300 g for 5 min. Culture supernatants were saved for Ig quantification by ELISA, and differentiation was assessed by flow cytometry.

Formalin-fixed and paraffin-embedded 3-μm-thick human tonsil sections were dewaxed by overnight incubation at 60°C, followed by rehydration through treatment with a decreasing alcohol gradient. Heat-induced epitope retrieval was performed for 20 min in Tris-EDTA buffer (pH 9) and citrate buffer (pH 6). Sections were permeabilized by incubation with 0.2% Triton X (Merck) in PBS for 12 min and blocked with 5% BSA (Miltenyi) or 5% Fc receptor blocking reagent for 1 h at room temperature. Samples were stained with various combinations of primary antibodies (Table S7) for 2 h at room temperature. After washing, fluorochrome-conjugated secondary antibodies (Table S7) were added to the tissue together with DAPI and incubated for 1 h at room temperature. After washing, FluorSave Reagent (Millipore) was applied, and coverslips were fixed with CoverGrip Coverslip Sealant (Biotium). Confocal microscopy image acquisition was performed with a ZEISS LSM 980 microscope (ZEISS Group, Jena, Germany) equipped with a 40×/1.3 NA EC Plan-Neofluar oil objective. High-resolution images were taken through a tile scan. Each square region was 610 µm per side with a pixel size of 0.13 µm. Images were processed with ImageJ (FIJI) software (Schindelin et al., 2012) and QuPath 0.5.1 software (Bankhead et al., 2017).

To generate tonsillar IgD rmAbs from IgD-ME B cells, as well as IgG rmAbs from IgG-ME B cells, we took advantage of BCR single-cell data from a recently published human tonsil atlas (Massoni-Badosa et al., 2024). Seurat object of the CITE-seq dataset was used, as well as the associated Cell Ranger outputs of BCR-seq (Massoni-Badosa et al., 2024). Antibody sequences were analyzed using the R package Seurat (Hao et al., 2024). Firstly, ME B cells were selected according to the original data annotation. Their corresponding BCR-seq contigs were filtered, so that only cells containing a single heavy chain sequence were included. They were categorized as IgD-ME and IgG-ME B cells according to their IGH gene expression, i.e., IGHD for IgD-ME B cells and IGHG1, IGHG2, IGHG3, or IGHG4 for IgG-ME B cells. IgBLAST (version 1.22.0) (Ye et al., 2013) was used to compare each sequence of heavy and light chains with the putative germline and estimate the number of mutations. Variable region sequences were ordered as gBlocks from IDT. AbVec-hIGG1, AbVec-IGLC2, or AbVec-IGKC vectors were linearized with EcoRI-HF and Sal-HF, and EcoRI-HF and XhoI restriction enzymes (NEB) for the heavy and light chains, respectively. Linearized vectors were purified from gel using QIAquick Gel Extraction Kit (Qiagen), and respective variable region sequences were cloned using the Gibson assembly method (Gibson et al., 2009). Each new construct was verified by Sanger sequencing. 150 × 10⁶ Expi293 cells seeded in 50 ml of Expi293 Expression Medium (Thermo Fisher Scientific) were transiently transfected with 1 µg/ml of antibody-encoding plasmid DNA. After 5 days, cells were harvested by centrifugation for 5 min at 4,000 rpm. The resulting supernatant was ultracentrifuged for 30 min at 30,000 rpm and filtered with a 0.45-µm Nalgene nylon filter (Thermo Fisher Scientific). Next, the filtered sample was applied onto a PBS-equilibrated HiTrap Protein A MabSelect or HiTrap LambdaMabSelect SuRe purification column (Cytiva) to purify Cɣ1- or Cδ-expressing mAbs, respectively. The column was washed with 10 CV PBS and eluted with 100 mM glycine buffer, pH 3. The eluted sample was further subjected to glycine-to-PBS buffer exchange and concentrated using the Amicon Ultra 50 K filter. Purified mAbs and the corresponding flow-through controls were loaded onto an SDS-PAGE gel for IgH and IgL band visualization and quality control. Plasmids encoding human IgD rmAbs from IgD-GC cells (Koelsch et al., 2007) were provided by Patrick C. Wilson and processed according to an rmAb production and purification protocol as described above. Peripheral blood IgG, IgA, and IgM rmAbs specific for the receptor-binding domain of severe acute respiratory syndrome coronavirus 2 were expressed in parallel with tonsillar IgD mAbs as described earlier (De Campos-Mata et al., 2024) and used as controls. Of note, all rmAbs were synthesized to express the Fc domain from IgG1, i.e., Cγ1, which facilitated their purification, except for four Cδ-expressing rmAbs, as specified below. In any case, each rmAb expressed the native antigen-binding Fab domain.

Four rmAbs from IgD-GC B cells, i.e., rmAb 47, 48, 70, and 98, were engineered to encompass Cδ instead of Cγ1. Cγ1-encoding vectors of interest were linearized with SalI-HF and HindIII-HF (NEB) and gel-purified using QIAquick Gel Extraction Kit (Qiagen). The constant region of the IgD HC was amplified by PCR from the pFUSE-CHIg-hD vector (InvivoGen) and cloned into our vector of interest by the Gibson assembly (Gibson et al., 2009). Final constructs were sequence-verified by Sanger sequencing. Transfection, culture, harvest, and purification were performed as described above. Since all Cδ-expressing rmAbs contained λ light chain, HiTrap LambdaMabSelect SuRe purification column (Cytiva) was used to purify IgD mAbs.

Mutations in the germline sequences of IGHV and IgL chain variable genes from rmAbs 47 and 48 were determined using IgBLAST. Full-length IGHV genes from these IgD rmAbs, including CDR3 regions, were identified following published criteria (Corbett et al., 1997). Germline sequences were ordered as gBlocks from IDT. Vectors were linearized with EcoRI-HF and Sal-HF, and EcoRI-HF and XhoI restriction enzymes (NEB) for the heavy and light chains, respectively. Linearized vectors were purified from gel using the QIAquick Gel Extraction Kit (Qiagen), and respective germline sequences were cloned using the Gibson assembly method (Gibson et al., 2009). New constructs were verified by Sanger sequencing, and germline revertant mAbs were produced and purified as described above.

Total IgD pAbs from cell culture supernatants were measured using Human IgD ELISA Quantitation Set (Bethyl Laboratories) according to the manufacturer’s instructions. Briefly, 96-well flat-bottomed plates (Nunc) were coated for 1 h with 10 μg/ml goat anti-human IgD coating antibody in a coating buffer composed of 0.05 M carbonate–bicarbonate in PBS, pH 9.6. Wells were washed with a washing buffer composed of ultrapure water supplemented with 50 mM Tris, 0.14 M NaCl, and 0.05% Tween 20 at pH 8. Wells were blocked for 30 min at room temperature with a blocking buffer composed of ultrapure water supplemented with 50 mM Tris, 0.14 M NaCl, and 1% BSA at pH 8. Both samples and standards were diluted in a dilution buffer composed of ultrapure water supplemented with 50 mM Tris, 0.14 M NaCl, 1% BSA, and 0.05% Tween 20 and incubated for 1 h at room temperature. After washing, a horseradish peroxidase (HRP)–labeled anti-human IgD antibody was added at 13.3 ng/ml and incubated for 1 h at room temperature. Plates were washed and developed using TMB Substrate Reagent Set (BD Bioscience). Absorbance was read at 450 nm with an Infinite 200 PRO plate reader (Tecan).

Total IgD pAbs from serum or mucosal secretions were measured using IgD Human ELISA Kit (Abnova) according to the manufacturer’s instructions. In both assays, IgD concentrations were calculated by extrapolating sample absorbance values with values from the standard curve.

96-well flat-bottomed plates were coated with 1 μg/ml of goat anti-human Ig (SouthernBiotech) in a coating buffer and incubated overnight at 4°C. After washing with PBS supplemented with 0.05% Tween-20 (PBS-T), wells were blocked with 1% BSA in PBS for 2 h at room temperature. Samples were diluted in PBS with 1% BSA and 0.05% Tween-20. A standard curve was prepared using purified human IgA (Mpbio) at 250 ng/ml, and six 1:1 serial dilutions were performed. Both samples and standards were incubated for 2 h at room temperature. After washing, HRP-conjugated goat anti-human IgA (SouthernBiotech) antibodies were added and incubated for 45 min at room temperature. Plates were washed and developed using the TMB Substrate Reagent Set, and absorbance was read at 450 nm with an Infinite 200 PRO plate reader (Tecan). Concentrations were calculated by extrapolating sample absorbance values with values from the standard curve.

To evaluate the reactivity of IgD rmAbs and IgD or IgE pAbs from serum samples, 96-well half-area flat-bottom high-bind microplates (Corning) were coated overnight at 4°C with a given antigen (Table S8) diluted in PBS or with PBS alone for background subtraction. To measure Ig reactivity to insulin or chitosan, the antigen was diluted in 1% acetic acid and neutralized with an equimolar concentration of NaOH before coating. After washing with PBS-T, wells were blocked for 2 h at room temperature with PBS supplemented with 5% BSA and 0.1% Tween. To measure Ig reactivity to dsDNA, ssDNA, or echovirus antigen, plates were blocked with PBS supplemented with 2% BSA. Recombinant mAbs and IP-IgD were diluted to a concentration of 10 µg/ml in PBS supplemented with 1% BSA and 0.05% Tween-20. Serum samples for IgE reactivity were serially diluted 1:3, starting with a 1:9 dilution. mAbs and IP-IgD samples were serially diluted 1:4 six times, added to the antigen- or PBS-coated plate, and incubated for 2 h at room temperature. After washing, plates were incubated with HRP-conjugated anti-human IgG to detect mAbs; HRP-conjugated anti-human IgD or IgE to detect serum IgD or IgE, respectively; or HRP-conjugated anti-human Igλ (SouthernBiotech) to compare the reactivity of Cγ1-expressing vs Cδ-expressing mAbs, for 45 min at room temperature. Plates were washed with PBS-T and developed with the TMB substrate reagent set (BD Biosciences). The development reaction was stopped with 1M H₂SO₄. Absorbance was read at 450 nm with an Infinite 200 PRO plate reader (Tecan). To quantitate the level of each antigen, optical density (OD) values were calculated after subtraction of the background defined as OD₄₅₀ value of the corresponding sample dilution on plates coated with PBS but no antigen. All negative values were normalized to 0. Threshold values were determined by calculating the mean ± 2 SD of the highest concentration OD₄₅₀ value of all anti-RBD antibodies or control serum samples.

Bacterial isolates of E. coli (25992; ATCC), Bacteroides fragilis (25285; ATCC), B. longum (15707; ATCC), R. lactatiformans (100348; ATCC), B. cereus (11778; ATCC), S. mutans (25175; ATCC), Streptococcus agalactiae (13813; ATCC), and R. intestinalis (14610; DSM) were heat-inactivated at 65°C for 20 min. Bacterial isolates were then incubated with rmAbs or immunoprecipitated serum pIgD at a concentration of 10 μg/ml for 30 min at room temperature. After washing, bacterial pellets were incubated for 15 min at room temperature with polyclonal PE-conjugated anti-human IgG antibodies (Infrared Laboratories) or APC-conjugated anti-human IgD antibodies (BioLegend). Bacteria incubated only with secondary antibodies were used as negative control to set up the PE- or APC-positive gate. Finally, bacterial samples were washed and resuspended in PBS with SYTO BC (Thermo Fisher Scientific; 1:60,000). Contamination was minimized by passing all buffers and reagents through sterile 0.22-μm filters before use. Bacterial cells were analyzed using Cytek Aurora Cytometer (Cytek Bioscience) with low FSC and SSC thresholds to allow bacterial detection. FSC, SSC, and SYTO BC were set to a biexponential scale, and samples were gated first at B1 vs. B14 detectors for SYTO BC⁺ to define bacterial counts free from autofluorescence, then as SSC⁺, and assessed for PE- or APC-positive counts.

Nasal swabs were immersed in 500 μl sterile filtered solution of PBS containing a protease inhibitor cocktail (cOmplete, Roche) at the concentration recommended by the manufacturer, and mixed by vortexing for 5 min. After removing the swab, tubes were centrifuged at 400 × g for 5 min to pellet cells and debris. BALs were transferred to sterile tubes and centrifuged at 400 × g for 5 min to pellet cells and debris. The resulting supernatants were passed through 70-µM cell strainers and transferred into new tubes. Sinonasal or bronchoalveolar bacteria were then pelleted at 8,000 × g for 5 min, and the supernatants were stored for Ig concentration measurements by ELISA. Bacteria were washed once with filtered PBS to remove any remaining free Igs. Finally, bacterial samples were resuspended in filtered PBS, followed by a 30-min incubation at room temperature with APC-conjugated anti-human IgD (BioLegend), PE-conjugated anti-human IgA (SouthernBiotech), or a mixture of both along with their corresponding isotype controls. Samples were washed and analyzed as described earlier.

Low serum and sinonasal IgD concentrations limit the use of prediluted whole samples for ELISA testing. We therefore developed an experimental protocol for the immunoprecipitation of IgD pAbs. Dynabeads M-270 Epoxy were conjugated with an anti-human IgD antibody (SouthernBiotech), using the Dynabeads Co-Immunoprecipitation kit (Thermo Fisher Scientific), according to the manufacturer’s instructions. Then, 0.5 mg of prewashed anti-IgD–coupled beads was resuspended in a mixture of 50 μl of serum with 50 μl of a buffer comprised of 50 mM Tris, 0.14 M NaCl, 1 % BSA, and 0.05 % Tween 20, or 100 μl of nasal swab supernatant obtained as described above, and incubated for 20 min with tilting and rotation at 4°C. After incubation on the magnet, the nonbound fraction was saved, and the beads were washed three times according to the manufacturer’s protocol. The immunoprecipitated IgD was then eluted in 40 μl of elution buffer and stored at −20°C for further analyses.

Tonsillar B cells and PCs were sorted as described earlier. Cells were fixed with 2.5% glutaraldehyde in phosphate buffer of 0.1 M for 2 h at 4°C, postfixed with 1% osmium tetroxide with 0.8% potassium ferrocyanide for 2 h, and dehydrated with increasing concentrations of ethanol. Then, pellets were embedded in EPON resin (EMS, Hatfield) and polymerized at 60°C for 48 h. Sections of 70 nm in thickness were obtained with a Leica EM UC6 microtome (Wetzlar), stained with 2% uranyl acetate and Reynolds solution consisting of 0.2% sodium citrate and 0.2% lead nitrate, viewed by a JEM-1400 transmission electron microscope (JEOL), and imaged at 120 kV voltage.

Genomic DNA was isolated from sorted B cell subsets using a protocol for small cell numbers described previously (van der Burg et al., 2006). The replication history of B cell subsets was determined using the κ-deleting recombination excision circle assay as reported earlier (van Zelm et al., 2007). This assay is based on a quantification of coding joints and signal joints of an Igκ-deleting rearrangement (intron RSS-Kde) by qRT-PCR. The ΔCT between the signal joint and the coding joint exactly represents the number of cell divisions a B cell has undergone. The previously established control cell line U698 DB01 (van Zelm et al., 2007) contains one coding and one signal joint per genome and was used to correct for minor differences in efficiency of both real-time quantitative PCR (qPCR) assays.

Sorted tonsillar cells were centrifuged and homogenized using QIAshredder (Qiagen). RNA was isolated with RNeasy Micro Kit (Qiagen) according to the manufacturer’s instructions. RNA quality and quantity were assessed by Bioanalyzer using Agilent RNA 6000 Pico Kit (Agilent Technologies). Samples used for analysis had RNA integrity number ≥6.4. Next-generation sequencing libraries with polyA capture were prepared according to the protocol of NEBNext Ultra II Directional RNA Library Prep Kit for Illumina, version 4.0, 4/21 (New England Biolabs).

Raw sequencing reads in the fastq files were mapped with STAR version 2.7.8a (Dobin et al., 2013) Gencode release 41 based on the GRCh38.p13 reference genome and the corresponding GTF file. The table of counts was obtained with the featureCounts function in the package subread, version 2.0.3 (Liao et al., 2014). Filtering of lowly expressed genes was done by keeping genes having >1 CPM in at least four samples. Raw library size differences between samples were treated with the weighted trimmed mean method (Robinson and Oshlack, 2010) implemented in the edgeR package (version 3.36.0). The normalized counts were used to make unsupervised analysis, PCA, and clusters. Genes starting with IgH (except for IGHMBP2), IgK, or IgL (except for those starting with IgLL and IgLon5) were excluded. For the DE analysis, read counts were converted to log2 counts per million and the mean–variance relationship was modeled with precision weights using the voom approach in the limma package. SVA package (version 3.38.0) was used to compute surrogate variables (SVs) using the svaseq function, and all SVs were added to the design matrix with the primary variable. Donor variables were also included in the model. FDR was used to adjust for multiple comparisons. Genes were considered differentially expressed if adjusted P value (adj. P) <0.05 and |log2FC| >1. All analyses were done with R version 4.1.2 Pre-Ranked GSEA (Subramanian et al., 2005) implemented in the clusterProfiler (Yu et al., 2012) package version 4.0.0 was used to retrieve enriched functional pathways (Reactome from MSigDB, version 7.5.1). The ranked list of genes was generated using the −log(p.val)*signFC for each gene from the statistics obtained in the DE analysis. Gene lists for heatmaps were manually curated based on the Panther version 17.0 classification (Thomas et al., 2022) of the DEG lists.

The Visium Spatial Gene Expression data (10X Genomics) and the corresponding images of tonsillar tissue slides were adapted from a study published recently (Massoni-Badosa et al., 2024). The experimental approaches and bioinformatics pipeline for preprocessing and annotating the data are described in the original publication (Massoni-Badosa et al., 2024). The annotation of the tissue regions, initially described by the authors, was slightly modified in the present work. The Seurat R package (version 5.1.0) (Hao et al., 2024) was used to assess the specific RNA markers from each region, and the “Follicle” label was replaced by the “GC” label, since the differential markers displayed GC characteristics. A representative tissue slide was selected for the spatial visualization of a custom IgD-ME gene signature (IGHD⁺IGHM⁻CD27⁺MME⁻CD38⁻CD19⁺CD3G⁻CD4⁻CD5⁻CD7⁻ ​CD8A⁻ITGAX⁺CD80⁺NT5E⁺FAS⁺CXCR3⁺CR2⁻TNFRSF13B⁺) and IgD-PC gene signature (IGHD⁺IGHM⁻MME⁻CD19⁺CD38⁺XBP1⁺ ​TNFRSF17⁺CD3G⁻). The intensity of the signature in the different spots of the slide was estimated using the UCell R package (version 2.8.0) (Andreatta and Carmona, 2021) in combination with Seurat (version 5.1.0) (Hao et al., 2024).

For the comparison of ABCs and ME B cells, scRNA-seq and CITE-seq datasets were obtained from a study published recently (Massoni-Badosa et al., 2024). The experimental procedures and bioinformatics pipeline for preprocessing and annotating the data are described in the original publication (Massoni-Badosa et al., 2024). For the selection of IgD-ME or combined IgG-ME and IgA-ME B cells, the Seurat R package (Hao et al., 2024) was used to subset the mRNA cluster annotated as ME B cells, whereas the UCell R package (Andreatta and Carmona, 2021) was used to compute surface protein signatures specific to IgD-ME (IgD⁺IgM⁻), IgA-ME (IgM⁻IgA⁺), and IgG-ME (IgM⁻IgG⁺) B cells.

The Seurat R package was used to select IgD class-switched cells within a pre-established ME B cell cluster by subsetting cells with expression levels of IGHD >2. To discard cells expressing IgM from the analysis, a threshold of background IGHM expression of 1 was established based on IGHM expression detected in IgG-ME, IgG-GC, and IgG-PC clusters with IGHG1 >1 or IgA-ME, IgA-GC, and IgA-PC clusters with IGHA1 >1. This threshold was used to discard cells with expression levels of IGHM ≥1. After rescaling and normalization, two IgD-ME B cell subclusters were identified.

For the selection of ABCs, CD19⁺CD20⁺CD10⁻CD21⁻CD27⁻ mRNA signature described by Holla et al. (Holla et al., 2021) was computed using the UCell R package. The ABC- and PC-associated gene signatures were computed on subclusters 1 and 2 of IgD-ME B cells using the UCell R package. The gene signature composed of ITGAX, CD69, FCGR2B, TBX21, TOX, and LTB was assessed in the ABC subcluster, and genes included in the PC-associated signature were obtained by estimating the DEG (adj. P <0.05, Log2FC >1) of the total PC cluster with the FindAllMarkers function (Seurat). The DEGs between subclusters 1 and 2 of IgD-ME B cells were assessed by the FindMarkers function (Seurat), and GSEA of GO terms was used for a functional characterization of differential biological processes.

The IgD-GC and IgG-GC B cells were compared to assess DEGs using the FindMarkers function (Seurat) with the default filtering parameters (minimum X-fold difference of 0.25, detected in at least 0.1 of cells in either of the two populations). A custom panel of genes highly involved in SHM was generated to compare IgD-GC cells with IgG-GC cells or IgD-GC subclusters 1 and 2 with IgG-GC B cells. The relationship between the abundance of MMR and BER events in the IgD-GC B cells in comparison with IgG-GC B cells was assessed by estimating the ratio between the mean expression of genes involved in the GO term “mismatch repair” (accession GO:0006298) and the genes from the “base excision repair” term (accession GO:0006284).

The trajectory inference was based on the scRNA-seq human tonsil dataset provided by a study published recently (Massoni-Badosa et al., 2024). The experimental procedures and bioinformatics pipeline for preprocessing and annotating the data are described in the original publication (Massoni-Badosa et al., 2024). The Seurat R package (Hao et al., 2024) was used to identify IgD-GC B cell, IgD-ME B cell, and IgD-PC clusters, which were isolated by selecting IgD class-switched cells within pre-established ME, GC, and PC clusters. Selection was implemented by subsetting expression levels of IGHD >2. To remove any potential cells with IgM isotype, a threshold of background IGHM expression 1 was established based on the expression levels of IGHM in IgA and IgG class-switched cells, i.e., ME, GC, and PC clusters with IGHA1 >1 or IGHG1 >1. This threshold was used to discard cells with expression levels of IGHM ≥1 from the analysis. After rescaling and normalization, the annotation of GC B cell, ME B cell, and PC clusters was verified by analyzing differentially expressed markers (Seurat functions FindAllMarkers, and FindMarkers), and verifying the DE of MME, CD38, AICDA, and BCL6 in GC B cell subclusters, CD27, CD38, FCRL4, and TNFRSF13B in ME B cell subclusters, and CD27, CD38, MS4A1, JCHAIN, XBP1, IRF4, PRDM1, and ERN1 in PC subclusters. The trajectory of the clusters was inferred using the Monocle3 package (version 1.3.7) (Cao et al., 2019).

Sorted cells were lysed with QIAshredder columns (Qiagen) following the manufacturer’s protocol. Cellular RNA was isolated from sorted B cell subsets with RNeasy Micro Kit (Qiagen) following the manufacturer’s instructions, with a final elution in 14 μl RNase-free water. 10 μl RNA was used for reverse transcription into cDNA using TaqMan Reverse Transcription Reagents with random hexamers (Thermo Fisher Scientific). For IGHV gene PCR, 5 μl of cDNA was mixed with High-Fidelity Platinum PCR SuperMix (Thermo Fisher Scientific) containing 50 nM forward primers specific for the FWR 1 of VH1, VH2, VH3, VH4, VH5, or VH6 and 250 nM primers specific for Cδ, Cμ, Cγ, or Cα and encompassing corresponding Illumina Nextera sequencing tags along with unique molecular identifiers (UMIs) (Table S9). For this study, UMI-based correction was not performed (Gupta et al., 2015; Nouri and Kleinstein, 2020). Amplification was performed with an initial step at 95°C for 3 min, followed by 35 cycles at 95°C for 30 s, 58°C for 30 s, and 72°C for 30 s, supplemented with a final extension step of 72°C for 5 min. PCR products were purified with AMPure XP beads (Beckman Coulter), and Nextera XT indices were added by PCR under the following conditions: 98°C for 30 s, 5 cycles at 98°C for 10 s, 63°C for 30 s, and 72°C for 3 min. AMPure XP beads (Beckman Coulter) were used to purify PCR products, which were subsequently validated and pooled. The final pool was quantified by qPCR. Single-strand products were paired-end–sequenced twice on a MiSeq instrument (Illumina) with 600 Cycle version 2 Kit (2 × 300 bp). In total, 12,076,7500 IGHV gene sequences from four donors were obtained through next-generation sequencing.

Raw sequencing reads for each sample from the two MiSeq runs were merged. The mean base call Phred quality score of all samples was ≥30. Subsequently, VDJ alignment and clonotyping were performed using the MiXCR software package, version 3.0.12, considering the full VDJ region nucleotide sequence for clonotype assembly (Bolotin et al., 2015). In this way, a clone is defined as a unique VDJ region nucleotide sequence. We excluded clones containing stop codons and out-of-frame sequences. As described previously (Greiff et al., 2017), we retained only those clones with at least a read count of two and with a CDR3 amino acid sequence length of >3 amino acids. For IgD class-switched IgD⁺IgM⁻ B cells, we only included IgD-unique clones for each sample by excluding IgD clones that have an identical CDR3 amino acid sequence with IgM clones found in the same sample when sequenced with IgM-specific primers. Moreover, for each sample, we excluded those clones where the MiXCR-annotated constant gene did not match the reverse primer used to generate the sample (e.g., IgM-annotated sequences in a sample where IgG primers were used).

We used the Change-O toolkit to parse the preprocessed dataset in blast format using the command MakeDb.py (Gupta et al., 2015). Subsequently, we assigned our clones into clusters using the hierarchicalClones command from the SCOPer package (version 1.1) (Nouri and Kleinstein, 2020). Sequences from the same donor with identical V and J germline genes and 85% junction nucleotide sequence similarity were assigned to a single lineage cluster. Maximum parsimony phylogenetic trees were built from clonal lineages including 20–200 clonal sequences per lineage using the R package dowser (version 0.0.4) (Hoehn et al., 2022).

To perform nearest neighbor analysis, we included only lineage trees that exclusively contained IgD sequences, for a total of 32, 59, 82, and 101 clonal lineages (i.e., trees) for donors 1–4, respectively. IgD-only trees were converted into pairwise distance matrices using the R package ape (version 5.4-1) (Paradis and Schliep, 2019), and we selected the nearest non-PC neighboring sequence for each PC-derived sequence (i.e., minimal distance). We calculated the proportion of neighboring sequences to IgD-PC sequences (%). For the normalized tree height analysis, pairwise distance matrices were used to extract tree height (i.e., distance from the germline to the furthest tip in the tree) and size (i.e., number of non-germline sequences per tree) information. Subsequently, we divided the tree height by the tree size to calculate the normalized tree height. For the average SHM count per tree edge analysis, we extracted pairwise tree edge lengths for each node and its parent from phylogenetic tree objects and calculated the corresponding average. In both normalized height and average SHM count per tree edge analyses, we included phylogenetic trees belonging to one isotype group that contained >1 cell type and had an occurrence of at least five trees per donor to ensure statistical robustness.

For this analysis, we included clonal lineages that exclusively contained IgD sequences. We first examined the association significance between tree trait values (cell types) and tree topology using the parsimony score test (PS test), as described previously (Hoehn et al., 2022). We calculated the significance of the PS test statistic using a permutation test, where the test statistic was recomputed for 1,000 randomizations of trait values at the tips for each tree. We observed a significant PS test statistic in phylogenetic trees from all donors of the study. Following that, we used a restricted SP test as described previously (Hoehn et al., 2022) to examine the most likely ancestor to IgD-PC B cells across all trees for each donor. Similarly, using a permutation test, we calculated the significance of the SP test statistic for our phylogenetic trees within each donor by comparing observed test statistics with those obtained from 1,000 randomizations of trait values at the tips for each tree. Of note, PS and SP test statistics were considered significant when accompanied with P-values <0.05.

We used the SHazaM R package (version 1.0.2) to compute the count and frequency of V gene replacement (R) and silent (S) mutations (Gupta et al., 2015). We considered a VDJ clone to be mutated when harboring more than 2 mutations since the PCR step might be responsible for one to two mutations per clone (Hendricks et al., 2019). To compute the R:S ratio in mutated clones, we divided the number of R mutations by the number of S mutations. In sequences with R > 0 but no S mutations (S = 0), we followed an approach described previously (Ghraichy et al., 2020), where the number of S mutations was set to 1 to avoid mathematically undefined results.

Germline V gene repertoires for each sample were calculated by determining the frequency of VDJ clones that map to a particular V gene independently of the clone occurrence frequency in the repertoire (Greiff et al., 2017). We used the ComplexHeatmap R package (version 2.5.1) to generate our heatmap following a complete linkage hierarchical clustering algorithm, where pairwise distances among rows and among columns were calculated by the Euclidean distance (Gu et al., 2017).

The Ig gene repertoire analysis, RNA-seq statistical analysis, and graphics generation were performed in RStudio (R version 4.0.3) (Posit Team, 2026). We used packages ggpubr, version 0.4.0 (Kassambara, 2026), and rstatix, version 0.6.0 (Kassambara, 2023), for statistical testing, as well as packages ggplot2, version 3.3.3 (Wickham, 2016), and cowplot, version 1.1.1 (Wilke, 2025), for data visualization, unless stated otherwise. For phenotypic and functional studies, statistical analysis was performed using Prism 5.03 software (GraphPad). Further details regarding statistical testing are included in figure legends.

Figs. S1, S2, S3, S4, and S5 and Tables S1, S2, S3, S4, S5, S6, S7, S8, and S9 are available online and include extended data supporting the main findings of this study. These data encompass flow cytometry gating strategies, DEGs, detailed patient clinical information, reagent lists, and other supporting material that integrates the information provided in the main text. Fig. S1 shows phenotyping and gating strategies to identify tonsillar B cell subsets. Fig. S2 shows tonsillar IgD-ME B cells exhibit unique tissue-specific phenotypic properties. Fig. S3 shows tonsillar IgD-ME B cells and IgD-PCs display a mutation profile consistent with GC ontogeny. Fig. S4 shows tonsillar IgD-ME B cell clones form mutation-intensive and large lineage trees, and show more clonal relatedness to IgD-PC B cells rather than IgD-GC B cells. Fig. S5 shows IgD rmAbs and pAbs broadly react against common aerodigestive antigens, including allergens and commensal bacteria. Table S1 shows DEGs from IgD-GC 1 and 2, IgD-ME 1 and 2, and IgD-PC 1 and 2 subclusters. Table S2 shows selected characteristics of PID patient cohorts. Table S3 shows molecular properties of IgD rmAbs from tonsillar IgD-ME B cells, IgG-ME B cells, and IgD-GC B cells. Table S4 shows selected characteristics of sinonasal sample donors. Table S5 shows selected characteristics of BAL sample donors. Table S6 shows antibodies used for flow cytometry and cell sorting. Table S7 shows antibodies used for tissue immunofluorescence analysis. Table S8 shows antigens used to test IgD reactivity by ELISA. Table S9 shows primers used for Ig gene repertoire sequencing.

V(D)J-seq data are publicly available at NCBI’s Sequence Read Archive under accession code PRJNA1044580 at https://www.ncbi.nlm.nih.gov/bioproject/?term=PRJNA1044580. Bulk RNA-seq data are publicly available at Gene Expression Omnibus under the accession number GSE243121 at https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE243121. scRNA-seq, CITE-seq, and spatial transcriptomics (version 2.0) datasets were published by Massoni-Badosa et al., and are available at Zenodo under https://doi.org/10.5281/zenodo.8373756. Correspondence and requests for materials should be addressed to A. Cerutti ([email protected]) and V. Greiff ([email protected]).

We thank Patrick C. Wilson (Weill Cornell Medicine, New York, NY, USA) for providing plasmids encoding IgD rmAbs from tonsillar IgD-GC B cells; Menno van Zelm for helping with κ-deleting recombination excision circle assays; Kenneth Hoehn for assisting with the assembly of phylogenetic tree analysis; personnel from the flow cytometry facility at Parc de Recerca Biomèdica de Barcelona (PRBB) for providing technical support and helping with cell sorting; personnel from the genomics facility at PRBB for helping with experimental design and DNA sequencing; personnel from the Pathology Department of Hospital del Mar for providing tonsillar samples; personnel from MARGenomics service at PRBB for providing excellent technical support and guidance; and personnel from the electron microscopy service at Universitat Autònoma de Barcelona (UAB) for providing excellent technical support. We are indebted to the Biobanc de l’Hospital Infantil Sant Joan de Déu per a la Investigació integrated in the Spanish Biobank Network of ISCIII for the sample and data procurement. We also thank Javier Antunez-Sanchez (School of Life Sciences, the University of Warwick) for his help in BCR-seq data analysis. Finally, we are grateful to Daniel Garzón Segura for helping with cellular immunology techniques.

This work was funded and supported by Ministerio de Ciencia, Innovación y Universidades grant RTI2018-093894-B-I00 (to A. Cerutti), Ministerio de Ciencia, Innovación y Universidades grant SAF2014-52483R (to A. Cerutti), Ministerio de Ciencia, Innovación y Universidades grant PID2024-155263OB-I00 supported by MICIU/AEI /10.13039/501100011033 and by FEDER, UE (to A. Cerutti), UiO World-Leading Research Community, European Research Council, European Advanced Grant ERC-2011-ADG-20110310 (to A. Cerutti), Institute of Health Carlos III, Miguel Servet Research Program grant MS19/00002 (to G. Magri), Ministerio de Ciencia e Innovación, Ramon y Cajal grant RYC2021-031642-I (to G. Magri), Ministerio de Ciencia e Innovación, Ramon y Cajal grant RYC2021-032169-I (to M. Guzmán), EU Horizon 2020 iReceptorplus #825821 (to V. Greif), Research Council of Norway project #300740 (to V. Greif), and Research Council of Norway project #331890 (to V. Greif).

Author contributions: Roser Tachó-Piñot: data curation, formal analysis, investigation, methodology, project administration, visualization, and writing—review and editing. Habib Bashour: data curation, formal analysis, investigation, methodology, project administration, software, visualization, and writing—review and editing. Martyna Filipska: data curation, formal analysis, investigation, methodology, project administration, resources, validation, visualization, and writing—review and editing. Celia Corral-Vazquez: data curation, formal analysis, visualization, and writing—original draft, review, and editing. Mauricio Guzman: investigation, methodology, and visualization. Xavi Marcos-Fa: investigation, methodology, and visualization. Donata Martinuzzi: investigation and visualization. Hannah Honner: investigation. Pablo Canales Herrerias: investigation and writing—review and editing. Sonia Tejedor Vaquero: investigation and methodology. Alba Sáez Gordón: investigation. Julia Perera-Bel: formal analysis and resources. Jorge Domínguez Barragán: data curation, formal analysis, and software. Berta Arcos-Ribas: investigation. Leire de Campos-Mata: investigation. Andrei Slabodkin: conceptualization, methodology, validation, and writing—review and editing. Maria Chernigovskaya: formal analysis. Maria Luisa Rodríguez de la Concepción: investigation and writing—review and editing. Jose Gutierrez-Marcos: supervision. Ana García-García: resources and writing—review and editing. Andres Nascimento Osorio: data curation, investigation, visualization, and writing—review and editing. Mariona Pascal: conceptualization, resources, and writing—review and editing. Jordi Yagüe: investigation. Manel Juan: resources. Juan Ignacio Aróstegui: formal analysis, resources, and writing—review and editing. Rafael Hijano Esqué: investigation. Albert Sánchez Font: conceptualization and investigation. Stephan Ehl: resources and writing—review and editing. Bodo Grimbacher: resources and writing—review and editing. Marta Rizzi: resources and writing—review and editing. Laura Dotta: resources and writing—review and editing. Kang Chen: writing—review and editing. Raffaele Badolato: resources and writing—review and editing. Laia Alsina: conceptualization, investigation, resources, and writing—review and editing. Saurabh Mehandru: conceptualization and investigation. Charlotte Cunningham-Rundles: funding acquisition and writing—review and editing. Jorge Carrillo: investigation and writing—review and editing. Giuliana Magri: conceptualization, investigation, methodology, and writing—review and editing. Victor Greiff: funding acquisition, investigation, project administration, supervision, and writing—original draft, review, and editing. Andrea Cerutti: conceptualization, funding acquisition, methodology, project administration, resources, supervision, validation, and writing—original draft, review, and editing.