The study was approved by the National Medical Ethics Committee of the Republic of Slovenia (approval 0120-187/2019/10) and the Swedish Ethical Review Authority (approval 2019-05260). Written informed consent was obtained from all the donors (their legal guardians) who participated in the study.

Three separate cohorts of patients were recruited. For the discovery and validation cohorts we recruited Slovenian patients with allergic asthma, GINA (Global Initiative for Asthma) steps 1-4, not (yet) eligible for the biological therapy (Supplementary Table S1). The age range of the patients was 6–20 years, and included children above 6 years of age (the minimum age to be eligible for omalizumab treatment) (14) and adolescents/young adults (within the age range of 10–24 years) (15). A confirmed allergy to either grass pollen or house dust mites was an additional inclusion criterion. Patients were excluded with current asthma exacerbation or symptoms of acute respiratory infection. The clinical cohort was comprised of Swedish patients with allergic asthma, multiple confirmed allergies, aged 9–16 years, GINA steps 2-5, who were selected for omalizumab therapy (Supplementary Table S2).

From the discovery and validation cohort patients 12 mL of venous blood was collected into K2EDTA tubes (BD Vacutainer, BD, New Jersey, USA). The blood plasma was separated by centrifugation at 300×g, 20 min, room temperature (RT); next, the white blood cells (WBC) were isolated by gradient centrifugation (1.200×g, 10 min, RT) using a density medium (1.094 g/mL, mixed from High Density Spin Medium and Density Diluent Medium, both pluriSelect, Leipzig, Germany) and SepMate50 tubes (Stem Cell Technologies, Vancouver, Canada).

From the clinical cohort patients, 2.5 mL of venous blood was collected into the PAXgene RNA tubes (BD Vacutainer), and stored at -20 °C within 2h of collection, and then at -80 °C within 24h after freezing at -20 °C. The samples were stored at -80 °C until RNA isolation.

The freshly isolated WBC from the discovery and validation cohorts were seeded in RPMI-1640 media (Gibco, ThermoFisher Scientific, Waltham, MA, USA) supplemented with 10% of the respective donor’s own blood plasma and 1% penicillin/streptomycin (Sigma-Aldrich, St. Louis, MO, USA) at a density of 2.5×10⁶ cells/mL. One half of the cells were treated with 3.32 µg/mL omalizumab (Xolair, Novartis, Basel, Switzerland) for 24 hours at standard cell culture conditions - 37 °C, 5%CO₂, >95% relative humidity (Figure 1). The specific concentration of omalizumab used was calculated from the manufacturer’s recommended dose for subcutaneous administration to patients (14).

The basophils were isolated using fluorescence activated cell sorting (FACS Aria Fusion, BD Biosciences, New Jersey, USA). First, the single cells (singlets) were gated using the FSC-A/FSC-H gating. Next, the CD45⁺ leukocytes (cells) were defined using the CD45/SSC plot, and, finally, a CD123/HLA-DR plot was used to gate for the CD123⁺HLA-DR^neg population, corresponding to the basophils. The representative plots are shown in Figure 1. The following conjugated primary antibodies were used: Mouse IgG anti-human CD45-PE (#1P-160-T100), Mouse IgG anti-human CD123-APC (#1A-700-T100) and Mouse IgG anti-human HLA-DR-FITC (#1F-474-T100, all EXBIO, Prague, Czechia). Antibody concentrations were used recommended by the manufacturer. On a subset of samples, the sorting efficacy was confirmed via flow cytometry of the sorted samples (Supplementary Figure S1A).

Following the 24-hour incubation, parts of the omalizumab pre-treated and control (no pre-treatment) WBC were used to perform BAT using a patient-specific allergen- either 1.0 µg house dust mite allergen (#ALYOST46) or a 10.0 µg grass pollen allergen mix composed of 5 different types of grass species- i.e., 2.0 µg of each allergen (#ALYOS106, both Stallergenes Greer, Baar, Switzerland) for 30 minutes. The concentrations correspond to the ones used for the clinical testing with BAT at the local medical center (University Medical Center Maribor, Slovenia). 180 IU/mg heparin and 20 µM CaCl₂ (both Sigma-Aldrich) were added to the media together with allergens, to mitigate cell aggregation and enable basophil activation, respectively.

At the same time, the cells were also stained for the basophils-specific markers (see Isolation of the basophils) and CD63 (Mouse anti-Human CD63-V450, #561984, BD Biosciences), to determine the % of the total and of the CD63-expressing (CD45⁺CD123⁺HLA-DR^negCD63⁺⁺) basophils on the flow cytometer (Figure 1). The BAT results were calculated as the % of the basophils expressing CD63 on their surface (CD63⁺⁺ basophils) in the total basophil population.

A corrected activity readout was used to assess the in vitro patients’ response to omalizumab. For this, the proportion of CD63-expressing basophils was normalized to the total WBC as a % of the CD63⁺⁺ basophils among the CD45⁺ leukocytes. Next, the basophil loss (missing basophils) was quantified as the decrease in the % of total basophils in the allergen challenged samples relative to the corresponding samples without the allergen challenge. These calculations were performed separately for omalizumab treated and untreated conditions. Then, the proportion of missing basophils was added to both the total and the CD63⁺⁺ basophil count in the allergen challenged samples. The final readout was defined as the proportion of activated basophils (CD63⁺⁺ basophils + missing basophils) in (total basophils + missing basophils). All the fold changes of basophil activity (allergen vs. no treatment, allergen + omalizumab vs. allergen) were calculated from these corrected activity readouts.

The RNA was isolated from the sorted basophils (the in vitro testing) using the RNeasy Plus Mini Kit (Qiagen, Hilden, Germany) and from PAX RNA tubes (the clinical cohort) using the PAXgene Blood miRNA Kit (PreAnalytiX, Hombrechtikon, Switzerland), following the manufacturers’ recommendations. The RNA isolated from the PAX RNA tubes was depleted using the GLOBINclear Human Globin mRNA Removal Kit (ThermoFisher Scientific, Waltham, MA, USA). The concentration was determined using the Synergy2 microplate reader (Biotek, Winooski, VT, USA), and the RNA`s integrity was assessed using the Bioanalyzer 2100 (Agilent Technologies, Waldbronn, Germany) or by gel electrophoresis.

For the RNA sequencing (RNAseq), cDNA libraries were prepared using the Stranded mRNA Prep Ligation Kit (Illumina, San Diego, CA, USA) according to the manufacturer’s instructions. Paired-end RNA sequencing (2×74 bp) was performed on a NextSeq 550 system (Illumina) using a NextSeq 500/550 Mid Output Kit v2.5 (Illumina). The raw FASTQ data were processed as previously described (16). An RNAseq data analysis was performed as described previously (17) in the R programming environment (R Core Team, Vienna, Austria). Importantly, mean–variance modeling at the observation level transformation (VOOM) (18), followed by linear models and empirical bayes in the limma R package (19), were applied to obtain the differentially expressed genes between the groups. The differential gene expression was calculated as a log₂-fold change in expression between the in vitro better and poor/non-responders to omalizumab treatment.

For the qRT-PCR, the cDNA was generated using the High-Capacity cDNA Reverse Transcription Kit (Applied Biosystems, Thermo Fisher Scientific, Waltham, MA, USA). The qRT-PCR was performed with LightCycler480 SYBR Green I Master (Roche, Mannheim, Germany) on a QuantStudio 12 K Flex Real-Time PCR System (Applied Biosystems). The forward and reverse primers used to detect gene-specific mRNA were designed with the Primer3 software, version 4.1.0, and produced by Sigma-Aldrich. The primer sequences are listed in the Supplementary Table S3. The gene expression was calculated using the 2^−ΔΔCt method. Ct values above 40 were considered negative. The geometric mean value of the beta-2-microglobulin (B2M), glyceraldehyde 3-phosphate dehydrogenase (GAPDH) and 18S ribosomal N5 RNA (RNA18SN5) was used as the internal control.

The interactome of the proteins coded by the differentially expressed genes found with the RNAseq analysis was built with the CytoScape platform (version 3.10.3) (20), based on the protein–protein interactions downloaded from the BioGrid database (build version 5.0.250; https://thebiogrid.org/). The non-protein coding genes were inserted into the interactome manually.

A gene ontology (GO) analysis was performed using the ShinyGO 0.85 gene-set enrichment tool (21), with the following parameters and selected options: comparison of the selected gene set against a background of all the genes detected by the above RNAseq, GO-terms from the GO Biological Process database with at least 3, and no more than 500, associated genes considered, enrichment calculated based on hypergeometric distribution followed by the false discovery rate (FDR) correction. Genes of interest for downstream testing were selected based on the above GO analysis.

All the data are presented as the mean ± standard error. Statistical analyses of the qRT-PCR and the clinical data were performed using the GraphPad Prism software (GraphPad Software Inc., San Diego, CA, USA). A D’Agostino-Pearson omnibus normality test showed that most of the data sets had abnormal distributions. Hence, to assess the statistical significance of the differences between the two groups, a Wilcoxon matched pairs test, or the nonparametric Mann-Whitney test were used in cases of paired or unpaired data, respectively. The nonparametric Spearman r was calculated to assess the correlations between the datasets. A receiver operating characteristic (ROC) analysis was calculated using the Wilson/Brown method and a 95% confidence interval. The significance limit was set at p<0.05 for all the tests. A Benjamini-Hochberg step-up procedure (FDR) was used to adjust the p-values for multiple testing where applicable. Where mentioned specifically, the data outliers were determined using Tukey’s IQR rule and excluded from the analysis.

To establish an in vitro model for testing the patients’ response to omalizumab, we used BAT. In clinics, BAT is used typically to confirm a patient’s allergic response to a specific allergen. The patient’s WBC are exposed to the tested allergens, and a positive result is defined as >5% of basophils expressing CD63 on their surface (CD63⁺⁺ basophils) (22). For our model, we used allergens (grass pollen mix, house dust mites) to which the patients were confirmed to be allergic, i.e., patient-specific allergens that guaranteed a positive BAT result, whereby the intensity of the basophil activation is dependent on the allergen’s concentration (Supplementary Figure S1B). It turns out that the 1× concentration (i.e., the concentration optimized for clinical testing at the local medical center, see Materials and Methods section) is indeed the one that gives the maximum activation. On the other hand, when an allergen to which the patient is not allergic is used, the percentage of CD63⁺⁺ basophils remained similar as in no treatment control, even when the concentration of the allergen was doubled (Supplementary Figure S1B).

Next, we made sure that the WBC samples from allergic asthmatics contained cells (basophils) with free FcϵRI not occupied by IgE (Supplementary Figure S1D), confirming that omalizumab could have a measurable effect in our experimental setting. Namely, omalizumab blocks free IgE in a patient’s blood/serum, and prevents it from binding to the FcϵRI receptors on the surface of the basophils, thus preventing their activation when exposed to the corresponding allergen. It cannot, however, block IgE already bound to the FcϵRI (23). Indeed, in our model, pre-treatment of the WBC samples with omalizumab prior to the challenge with the patient-specific allergen resulted in a reduced % of CD63⁺⁺ basophils upon the challenge, but had no major effect on the % of CD63⁺⁺ basophils in the control (no allergen challenge) sample (Supplementary Figure S1C).

After confirming the feasibility of our in vitro setting, we used samples from 6 juvenile asthmatics (3 males, 3 females, mean age 14.2 years, GINA steps 1-3, Supplementary Table S1) with a confirmed allergy to either grass pollen or house dust mites to test it (the discovery cohort). In all 6 cases, exposure of the WBC to the patient-specific allergen increased the share of CD63⁺⁺ basophils markedly (Supplementary Figure S2A). While analyzing the numbers of basophils though, we found that the allergen exposure not only increased the proportion of CD63⁺⁺ basophils, but, in many cases, also decreased the total proportion of basophils in the WBC (Supplementary Figure S2B). Such a phenomenon was described previously, and may be attributed to strong activation and subsequent extensive degranulation, cell death or structural collapse, and, thus, changed scatter characteristics on the flow cytometer (24). We used a reverse gating strategy to test whether such a deleterious over-activation effect had occurred in our setting. We gated the total events (WBC) for the basophil-specific markers (CD45⁺CD123⁺HLA-DR^neg) in two steps, and then used the FSC-A/FSC-H dot plot to discriminate between the intact cells (i.e., cells with the expected scatter characteristics) and the cells with altered scatter characteristics (Supplementary Figure S2C). It turned out that, on average, 6.57 ± 2.23% of the CD45⁺CD123⁺HLA-DR^neg events (putative basophils) were present in the altered scatter characteristics gate under control (no treatment) conditions, which doubled to 14.73 ± 4.32% after exposure to the patient-specific allergen (Supplementary Figure S2D), while a handful of events was scattered outside the defined gates. The omalizumab pre-treated samples had, on average, a similar percentage of the putative basophils with altered scatter characteristics as those without omalizumab pre-treatment (7.10 ± 2.02% without allergen challenge, 14.10 ± 3.99% with allergen challenge). Looking at the individual donors, the allergen challenge increased the % of putative basophils with altered scatter characteristics in 5/6 samples (mean Fc 2.73 ± 0.60%), while the omalizumab pre-treatment counteracted this partially in 4/6 samples (mean Fc 1.08 ± 0.27%, Supplementary Figure S2E). As this result confirmed that some of the basophils in our experimental setting may indeed alter the scatter characteristics after the allergen challenge, we included the information on the percentage of “missing” basophils in our assessment of the omalizumab (non)-response.

Thus, the allergen exposure increased the % of activated (CD63⁺⁺ + missing) basophils from ~1.3 – 16.7% to 10.3 – 76.3% (a 2.07 – 26.3-fold increase, mean increase 7.89 ± 3.75-fold; Figure 2A). The omalizumab pre-treatment had a varied effect, sometimes decreasing the share of the activated basophils after the allergen challenge, and sometimes not (Fc 0.28 - 1.15, mean Fc 0.69 ± 0.12, Figure 2B). This effect of omalizumab was independent of the patient’s age, sex, GINA step, ACT score, therapy, allergen used, or the fold-change in the activated basophils after the allergen challenge (Supplementary Figure S3A-H). These results thus confirm the concept that the variable response to omalizumab treatment can be measured in an in vitro experimental setting independent of the actual therapy in a patient (25). Based on the fold change in the percentage of activated basophils from allergen to allergen + omalizumab treated samples (BAT Fc), the 6 patients were ranked into two equal-sized groups (3 better and 3 poor/non-responders to omalizumab), in order to enable downstream comparative transcriptomic analysis.

In parallel to the allergen challenge, we also isolated the basophils from the same WBC samples via the FACS (mean purity 97.2 ± 1.7%, Figure 1; Supplementary Figure S1A), and used them for the cell type-specific whole transcriptome analysis (RNAseq). The expression profiles were compared of the basophils from the better and the poor/non-responder group of patients. 71 genes were expressed differentially (i.e., |log2FC|>2.0, p<0.05; Figure 3A; Supplementary Table S4). 18 genes were upregulated, most prominently HBB, HBA1 and HBA2 (log₂FC>5.0) and 53 genes were downregulated. Among the latter, for 9 genes- APOBEC3B, CCDC163, CMPK2, DUSP2, EPSTI1, HERC5, HLA-C, OAS3 and RSAD2 the p-value also remained significant after the adjustment for multiple testing (Supplementary Table S4). No upregulated genes remained significant after the adjustment. The interactome of the 71 genes revealed a single major network with 16 nodes representing individual genes, while the majority of the nodes (46 nodes, of which four represent non-protein coding genes) had no direct connections to each other (Figure 3B). Among the nodes with no direct connections in the interactome were also 7/9 genes with the adjusted p<0.05 (APOBEC3B, CCDC163, CMPK2, EPSTI1, HERC5, OAS3 and RSAD2), while HLA-C was part of a small two-node network, and only DUSP2 belonged to the single major network.

Next, we performed a Gene Ontology (GO) analysis of the 9 differentially expressed genes with adjusted p<0.05. The analysis revealed 96 over-represented GO terms (Supplementary Table S5), however, associated mostly with only one or two of the nine genes (69 terms and 19 terms, respectively). Only 9 GO terms were associated with more (3, 4) genes, of which 8 had an FDR p<0.01. All of these 8 GO terms are connected to a (negative) regulation of the viral life cycle/(defense) response to a virus (Figure 3C; Supplementary Table S5), and associated with the same four of the differentially expressed genes with the adjusted p<0.05- APOBEC3B, HERC5, OAS3 and RSAD2, so these genes were identified as the most interesting targets for further study.

To test our in vitro model further and analyze the above most interesting target genes, an additional 22 juvenile asthmatics (16 males, 6 females, mean age 13.8 years, GINA steps 1-4, Supplementary Table S1) with confirmed allergy to either grass pollen or house dust mites were recruited as the validation cohort. The same as in the discovery cohort, the WBC were isolated, pretreated with omalizumab/no treatment, and then exposed to patient-specific allergens. Again, exposure of the WBC to the patient-specific allergens increased the proportion of the CD63⁺⁺ basophils markedly, and in practically all cases decreased the total number of basophils (Supplementary Figures S2A, B, respectively). In the samples from patients PT09 and PT21 the allergen challenge did not increase the proportion of the CD63⁺⁺ basophils above 5% (i.e., no positive reaction), hence, the two samples were excluded from further analysis. As with the discovery cohort, pretreatment with omalizumab had a varied effect on the fold change in the % of activated basophils between omalizumab pre-treated + allergen challenged samples vs. allergen challenged samples with no omalizumab pre-treatment (Fc 0.32 – 6.00, mean Fc 1.25 ± 0.27, Figure 4A). The results were independent of the patient’s age at inclusion, sex, GINA step, ACT score, FeNO level, comorbidities (allergic conjunctivitis, atopic dermatitis), regular therapy, application of immunotherapy, allergen used in the test, or the fold-change in activated basophils after the allergen challenge (Supplementary Figure S4C-K). Of note, the BAT Fc result in patient 19 (PT 19) deviated markedly from all the others (Figure 4A). Since Tukey’s IQR rule identified this result as an outlier, we excluded this sample from the further analysis.

As before, in parallel to the allergen challenge and omalizumab response assessment, basophils were isolated from the same WBC samples, and used for the qPCR validation of the four most interesting target genes. As the separation of the patients/samples into the better and poor/non-responders is necessarily at least partially arbitrary, given the lack of standardized cut-offs for response to omalizumab in vitro, we focused primarily on correlation analyses, while using better/poor responder categorization only for the complementary analyses. We checked whether the expression of the four genes correlated with the response of the cells to the omalizumab pre-treatment after the allergen challenge (the BAT Fc results). Indeed, the expression of RSAD2 in the basophils without omalizumab pretreatment correlated with the BAT Fc results significantly negatively, also after adjusting the p-value for multiple testing (Spearman r = -0.629, p = 0.005, adjusted p = 0.021, Figure 4B), which is consistent with the RNAseq results. No significant correlation between HERC5, APOBEC3B or OAS3 expression and the BAT Fc results were found (Figures 4C-E). Additionally, based on a margin derived from the rank-based separation of the discovery cohort into two equal size groups (Fc ~0.6), we explored the group-wise differences in the validation cohort, which yielded only 3 better responders (PT 11, 16 and 28). In this small group, the mean baseline RSAD2 expression in the basophils showed a trend toward higher expression compared to the poor/non-responder group (p=0.053, Supplementary Figure S5A). Further, the exploratory ROC analysis suggested a discriminatory capacity of the pre-treatment RSAD2 expression in the basophils (AUC = 0.867, 95% CI 0.667 – 1.000), although the statistical significance was borderline (p=0.051, Supplementary Figure S5B). The detailed threshold-based performance metrics are provided in Supplementary Table S6.

To assess the clinical relevance of our results on RSAD2, we examined its expression in the WBC of a small clinical cohort- 6 juvenile asthma patients (2 males, 4 females, mean age 13.8 years, GINA steps 2-5, Supplementary Table S2) selected to start the omalizumab therapy. Based on a clinical evaluation and improvement as perceived by the patients after 3 and 6 months of therapy, 4 of the patients were assessed as responders and 2 as non-responders to omalizumab. Of the two, one was a clinical non-responder, while the other reported no perceived improvement, and discontinued the omalizumab treatment due to side-effects after the 3 months follow-up. In general, the RSAD2 expression increased after 3 months of omalizumab therapy compared to the expression before the start of therapy, and stayed increased in the non-responsive group, but not in the responders after 6 months, although the differences were not significant (Figure 5A). Of note, the pre-therapy RSAD2 expression was higher in the responders group compared to the non-responders (not significant).

As an additional exploratory clinical comparison, we also assessed whether the RSAD2 expression tracked with the available clinical measures. Indeed, the pre-therapy RSAD2 expression showed a positive, but non-significant, association with the baseline FEV1/FVC ratio (forced expiratory volume in 1 s/forced vital capacity- an indicator of airways obstruction (26)) after 3 months of therapy (Spearman r = 0.900, p = 0.083, Figure 5B). No correlation was found between the pre-therapy RSAD2 expression and the FeNO (fractional exhaled NO) levels, FEV1/FVC ratios or ACT (asthma control test) score at therapy induction (Supplementary Figures 6A–D). Similarly, the patients’ GINA step, age at inclusion, sex, comorbidities, and regular therapies besides omalizumab did not influence the pre-therapy RSAD2 expression (Supplementary Figures 6E–M). Taken together, the clinical cohort exhibited RSAD2 expression trends compatible with those observed in our in vitro analyses.