We designed an observational longitudinal 12-month study to analyze the upper respiratory tract microbiome in adults with CVID ( Figure 1 ). The study involved CVID patients aged over 18 who were regularly followed by the Referral Care Centre for Primary Immunodeficiencies at Sapienza University of Rome, Italy. Patients were diagnosed according to the ESID criteria for CVID (20). Twenty-six healthy donors (HD) recruited among administrative employees of Sapienza University were also invited to participate in the study. At baseline, we collected demographics, IgG trough levels, IgA and IgM serum levels, and peripheral immune phenotype, including frequencies of B cell and Switched Memory B cells (MBC), to group patients according to the EUROCLASS classification (21). We stratified patients as having undetectable (<0.01 g/L) or detectable (≥0.01 g/L) IgA serum levels and as having IgM serum levels above or upper 0.20 g/L (2 DS lower than the reference). We also collected CVID-related health issues, including the presence of bronchiectasis (by CT scan), chronic obstructive pulmonary disease (COPD) (22), systemic autoimmunity and autoimmune cytopenia, enteropathy, Granulomatous and Lymphocytic Interstitial Lung Diseases (GLILD) (23), and concomitant treatments to allow comparison of outcome measures. Based on the data collected, we classified patients according to the prevalent manifestations into the infective or complicated phenotype (24). After enrolment, clinicians recorded monthly respiratory infections and antibiotic courses for six months (T0-T6). Six months after enrollment (T6), oropharyngeal swabs were collected. Participants were excluded if they had experienced symptoms of acute respiratory infections or had taken antibiotics in the month before the sampling. After swab collection, respiratory acute infections were monthly recorded for additional six months (T6-T12) through interviews and clinical evaluations. During the study time, patients were allowed to continue their therapies, including immunoglobulin replacement therapy (IgRT) and antibiotic prophylaxis. The study followed the Good Clinical Practice guidelines and the Declaration of Helsinki. The Ethical Board of the Policlinico Umberto I of Rome approved this study. The patients/participants provided their written informed consent to participate in this study.

Oropharyngeal sampling was carried out using nylon flocked swabs (Copan, Brescia, Italy) placed in milk-tryptone-glucose-glycerol (STGG) medium (25). One swab was collected from each patient, refrigerated, and transferred to the laboratory at Istituto Superiore di Sanità within four hours from collection, then stored at -80°C until further processing. Total DNA was extracted from each swab using the NucleoBond^® AXG DNA kit (MACHEREY-NAGEL GmbH & Co., Duren, Germany), following the manufacturer’s instructions. Briefly, samples were vortexed, and 800 μl of STGG medium was used to extract total DNA. Libraries were prepared by PCR enrichment of 16S rRNA encoding sequences using universal primers for the V3-V4 hypervariable region of the 16S rRNA gene and sequenced by MiSeq technology (Illumina, Harvard, California).

IgG, IgA and IgM serum levels were assessed by immunoturbidimetric assay according to the manufacturer’s instructions (BIO plastics, Rome, Italy).

PBMCs were stained with the appropriate combination of fluorochrome-conjugated antibodies to identify B and T cell subsets. B cell subsets were identified based on the expression of CD19, CD27, CD24, and CD38 markers by flow-cytometry. T cells were identified based on the expression of CD3+, CD4+, CD8+ by flow-cytometry. FCS files were analyzed using the FlowJo ™ v10.8.1 software (BD, Biosciences). Prior to analysis, all samples were normalized to equal cell numbers.

Continuous variables were described using median and interquartile ranges (IQR), and categorical variables using frequencies and percentages. Comparisons of continuous parameters between treatment groups were calculated with a t-test if normally distributed and with a Mann-Whitney U test if not normally distributed (as tested by the Kruskal-Wallis’s test). Differences in frequencies between groups were calculated by using the χ2 exact test. Secondary analyses were performed by multiple logistic and a binomial regression model to ascertain risk factors related to alpha diversity (Chao), dysbiosis index, and acute respiratory infections recorded in the T6-T12. Odds ratios (OD) and 95% confidence intervals (CIs) were calculated. Differences were considered significant when p values<0.05. Analyses were carried out using SPSS software, version 18 (SPSS, Chicago, IL) and GraphPad Prism version 8.0.0 for Windows, GraphPad Software, San Diego, CA USA.

Sequence data were processed using the Quantitative Insights Into Microbial Ecology (QIIME) pipeline (v.1.9.1.) (26). Data pre-processing included sample demultiplexing and trimming of Illumina adapters and primers. Paired-end reads were first assembled into longer joined sequences, and then sequences shorter than 200 bp and greater than 1000 bp were discarded. Chimeric sequences were removed before downstream processing. Sequences were clustered de novo into Operational Taxonomic Units (OTUs) using UCLUST at 97% similarity (As 97% identity in rRNA genes identifies the same taxa, we clustered sequences from the same taxa together). OTUs were assigned through a comparison with the SILVA database using Mothur. Alpha and beta diversity analyses were performed at all taxonomic levels (Phylum, Class, Order, Family, and Genus) (27). Samples were first characterized in terms of alpha diversity: sample richness was explored in terms of the observed number of OTUs and Chao1 index (Chao, 1987). The t-test and the paired t-test were used to determine differences between samples. Sample-to-sample differences in microbiome composition were quantified using beta-diversity and visualized in low-dimensional space through Principal Component Analysis (PCA). PERMANOVA was used to detect global community differences in PCA. Differences in microbiome composition between CVID patients and HDs were studied by Linear discriminant analysis effect size (LEfSe) (28) to identify taxa that differed consistently between sample types. LEfSe employs the non-parametric factorial Kruskal-Wallis’s sum-rank test (a = 0.05) to identify taxa with significantly different abundances between categories, followed by Linear discriminant analysis (LDA) to estimate the effect size of each feature of the differential abundance. The differences in abundance were regarded as statistically significant when the logarithmic LDA score was >2.0. If multiple taxonomic levels with different ranks showed significance in the same taxon, the lowest-ranked taxa were regarded as responsible. All sequencing data associated with this study were uploaded to the NCBI bio project database: PRJNA747877.

A CVID-specific microbial dysbiosis index was calculated according to Gevers et al. (29). We first defined CVID- and HD-taxa as the taxa enriched in CVID and HD groups, respectively. Then, we computed CVID- and HD scores by summing the relative abundance of CVID and HD taxa in each sample, respectively. Finally, the dysbiosis index was computed by subtracting the HD- and the CVID score and multiplying the result by 100, as shown in the following formula:

The equations should be inserted in editable format from the equation editor.

Using a linear classification model, we also investigated if taxonomic global microbial community composition alone would be sufficient to distinguish CVID patients from HD (CVID vs. HD). We used SIAMCAT version 1.8.1 (Statistical Inference of Associations between Microbial Communities And host phenoTypes, https://siamcat.embl.del/), a pipeline for analyzing microbiome data (30). SIAMCAT provides a full pipeline supporting data pre-processing, statistical association testing, and statistical modeling. Briefly, we first filtered microbial species based on their maximal abundance in any of the samples using 3% as a threshold using the filter.features function (parameters: filter.method = ‘abundance’, cutoff = 0.03). We then log-normalized microbial abundances using the normalise.features function (parameters: norm.method = “log.unit”, norm.param = list (log.n0 = 1e-06, n.p = 2, norm.margin = 1)). After splitting the data for cross-validation using the create.data.split function (parameters num.folds = 5, num.resample = 2), we trained a Random Forest classifier using the train.model function (parameters: method = “randomForest”). To perform and evaluate model predictions on the test set, we used the make.predictions and the evaluate.predictions functions. Ultimately, we generated the ROC curves and the model interpretation plot using the model.evaluation.plot and model.interpretation.plot functions. All statistical analyses were performed in R.

Oropharyngeal samples were obtained from 80 patients. After sequencing, three samples were excluded owing to the low number of reads (below 9,000). Five patients were excluded as they had an infection treated by antibiotics in the month preceding the sampling (T5-T6). The demographic and clinical data of the 72 CVID patients included in the analysis are reported in Table 1 . As controls, we enrolled 26 asymptomatic HD, age- and sex-matched (age, median 47.7 years, range, females 13, 50%). Frequencies of smokers were comparable among groups ( Supplementary Table 1 ). All CVID patients were treated by IgRT. Forty-seven percent of participants had a diagnosis of COPD, and 50% had bronchiectasis identified at a CT scan. Twenty-eight percent of CVID patients had undetectable serum IgA (<0.01 g/L) and 28% IgM levels<0.20 g/L. Fifty-seven percent of patients were classified as having a complicated phenotype: 30.6% had one or more autoimmune diseases 16.7% GLILD, 5.6% enteritis, and 15.2% malignancy ( Table 1 ). Seven patients (9.7%) received antibiotic prophylaxis with azithromycin (500 mg thrice weekly). Sixty-three antibiotic courses (range 1-5 per patient) were prescribed to 34 patients in the period T0-T5. Details on antibiotic treatments are provided in Table 1 .

According to the Chao1 index, the alpha diversity, a measure of the richness and evenness of a sample, was lower in CVID compared to HD (p=0.003, Figure 2A ). Beta-diversity analysis, which measures differences between communities, showed a separation between HD and CVID groups (PERMANOVA p=0.042, Figure 2B ). Microbiome profiles were also investigated in terms of taxonomic assignment and relative abundance ( Supplementary Figure 1 ). We found 24 taxa at different taxonomic levels with significant differential abundance in CVID patients compared to HD, with the order Lactobacillales, including the family Streptococcaceae being more abundant in CVID, and the order Bacteroidales, Flavobacteriales, Fusobacteriales, Selenomonadales, Campylobacterales, including the families Prevotellaceae, Fusobacteriaceae, Veillonellaceae, Campylobacteraceae, and Flavobacteriaceae being more abundant in HD ( Figure 2C ; Supplementary Table 3 ). To quantify the strength of the association between CVID status (CVID vs. HD) and global oropharyngeal microbiome observed by metagenomic analysis, we also used a classification concept based on LASSO regression (30). Cross-validation accuracy was high for CVID status (mean AUC=82%), confirming a different global oropharyngeal microbiome composition in CVID and HD ( Figure 3 ).

Intra-CVID differences in species richness and diversity (Chao1) were analyzed by univariate analysis, revealing a decreased microbial alpha diversity associated with age ≥45 years, low IgM serum levels (<0.20 g/dL), undetectable IgA serum levels (<0.1 g/dL), low B cells (≤1% of lymphocytes), and severe/moderate COPD status ( Figures 4A–E ; Supplementary Table 2 ). A multiple logistic regression model using a stepwise selection confirmed that the Chao1 index was related to IgM defect, age, and COPD moderate/severe status ( Table 2 ). A reduction trend in Chao1 was observed in those treated with antibiotic prophylaxis, even if not statistically significant ( Figure 4F ). Differently, bronchiectasis, days of antibiotic treatments (T0-T5), GLILD coexistence, having a CVID complicated phenotype, and the concomitant treatment with steroid inhalers were not associated with a reduction in alpha diversity in this cohort ( Supplementary Table 2 ). At the study time, all participants were receiving immunoglobulin replacement treatment (median IgG TL 6.7 g/L, Table 1 ). No contribution of IgG TL to the Chao1 index was observed in this cohort ( Supplementary Table 2 ). Beta-diversity analyses revealed a significant separation between CVID lacking IgA and CVID with detectable IgA serum levels (PERMANOVA p = 0.014, Supplementary Figure 2 ). We further analyzed the impact of the identified variables on the bacterial community composition. Compared to those with IgM ≥0.20 g/L, CVID patients with IgM< 0.20 g/L had lower levels of Prevotella and Veillonella, two genera usually related to a healthy status ( Figure 5A ). Furthermore, patients with undetectable IgA serum levels had a higher relative abundance of the genera Haemophilus and Streptococcus compared to those with serum IgA levels ≥0.01 g/L ( Figure 5B ). A similar pattern was observed in CVID patients with moderate/severe COPD, who showed higher relative abundance for Streptococcus, Haemophilus, and Rothia genera compared with patients with mild or without COPD ( Figure 5C ). Compared to controls, patients >45 years had an even higher abundance of Streptococcus genera and lower levels of Prevotella and Veillonella than younger patients ( Supplementary Figure 3 ). No differences in bacterial community composition were identified among the seven patients treated with azithromycin prophylaxis. We further analyzed the effect of antibiotic therapy to treat acute respiratory infection, only identifying a reduced relative abundance of Prevotella in those treated with one or more courses of quinolones (p=0,008).

To explore the possible impact of disease-related variables on the overall severity of dysbiosis, we calculated a CVID-specific dysbiosis index. We used the relative abundance of the taxa differentiating CVID patients and controls to calculate the index according to the method as previously described (12, 29). The dysbiosis index was higher in CVID patients than in HD. It was inversely associated with the alpha diversity (R 0.6 p<0.0001, Supplementary Figure 4 ), confirming that the abundance of selected taxa primarily captures the dysbiosis in CVID patients. In a linear regression model, the dysbiosis index was directly related to be treated by prophylaxis with macrolides (OR 31.9, p=0.009, 95%CI 8,6 to 55.3), and negatively related to being female (OR -21.4, p=0.006, 95%CI from -36.2 to -6.7) ( Table 3 ). As only few patients were treated with antibiotic prophylaxis, we also tested the model after removing its contribution and we further identified the moderate/severe COPD status as having a direct relation with dysbiosis index (OR 18.7, p=0.036, 95%CI from 1.3 to 36.1) ( Supplementary Table 4 ).

Next, we investigated whether the observed oropharyngeal microbiota differences in CVID patients were associated with infections in the T6-T12 period. In binomial regression model including all CVID participants, Haemophilus (p = 0.029, OR 5.52, 95%CI 0.606- 10.427, Supplementary Table 5 ) positively correlated with having a respiratory acute infection within 60 days from sampling, whereas Anoxybacillus (p=0.026, OR 1.47, CI95% 0.18-2.74) positively correlated with having a respiratory acute infection within six months ( Supplementary Table 6 ).