This study recruited a total of 10 individuals between 3 and 6 years of age who were referred for elective dental procedures under general anesthesia ( Table 1 ). Families were approached before the day of the procedure to discuss the study protocol and explain how it would deviate from the standard anesthetic protocol. On the day of the procedure, prilocaine cream was applied to both antecubital fossae for at least 15 min to facilitate the insertion of a 22-G cannula into the cephalic or median cubital veins. Analgesic cryogenic ethyl chloride spray was applied after skin preparation, and ultrasound guidance was used if the vein could not be visualized or felt. The first blood sample was collected before the administration of any anesthetic gases (pre-SEVO). At this point, parents were asked to confirm consent for inhalational induction. A mixture of nitrous oxide and oxygen was then administered, with the concentration of sevoflurane gradually increased to 8%. Once the child was asleep, nitrous oxide was turned off, and a second blood sample was collected via the cannula (post-SEVO). Hemodynamic and oxygen stability were maintained throughout this step. Sevoflurane was then reduced to 3%, and additional anesthetic drugs were administered to facilitate airway management before proceeding with the planned dental procedures. This study was carried out with the approval of the UK National Research Ethics Service (EXE-T1D, REC ref 17/EM/0255).

Fresh heparinized blood was diluted 1:1 with complete medium (RPMI-1640 GlutaMAX supplemented with penicillin–streptomycin [Thermo Fisher Scientific, Waltham, MA, USA]). The diluted blood was carefully overlaid into Leucosep tubes containing Lymphoprep (Axis-Shield Diagnostics Ltd, Dundalk, Ireland) and centrifuged at 1,000 × g for 15 min at room temperature (RT). The plasma layer was carefully removed and discarded, and the Peripheral Blood Mononuclear Cells (PBMC) layer was washed twice with complete media at 300 and 200 × g for 10 min each. After determining the total cell number, Peripheral Blood Mononuclear Cells (PBMC) layer were centrifuged at 400 × g for 5 min, resuspended in CryoStor freezing medium (Sigma-Aldrich, St. Louis, MO, US), and transferred into cryogenic vials (Corning, Corning, NY, USA). Vials were placed in a controlled-rate freezing vessel (CoolCell, Biocision, San Rafael, CA, USA) at – 80°C overnight before being transferred to temperature-monitored LN2 tanks for long-term storage.

Leukocyte populations were characterized from fresh whole blood (collected in 2.7 mL Ethylenediaminetetraacetic Acid (EDTA) tubes) using five multiparameter flow cytometry panels previously validated and assessed for technical reproducibility across multiple laboratories (12, 13). In brief, surface marker staining ( Supplementary Tables S1 - S4 ) was performed on 100–200 μL of well-mixed fresh whole blood for 45 min at RT, followed by red blood cell lysis for 8 min at RT (10× BD Fluorescence-Activated Cell Sorting (FACS) lysing solution diluted in Double Distilled Water (ddH₂O), BD Biosciences, USA). Staining for the lineage flow cytometry panel was performed using BD Trucount tube (BD Biosciences, Franklin Lakes, NJ, USA), allowing for the calculation of absolute cell numbers in addition to cell frequencies. The lineage tube was vortexed to ensure homogenization and kept on ice until acquisition (lyse, no wash). All other staining panel tubes were centrifuged at 500 × g for 5 min, washed in 2 mL of FACS buffer (1 × phosphate-buffered saline [PBS; Invitrogen, USA] containing 0.2% BSA [Sigma-Aldrich, USA] and 2 mM EDTA [Sigma-Aldrich, USA]) by centrifugation at 500 × g for 5 min, resuspended in 200 μL of FACS buffer, and kept on ice until acquisition. Intracellular marker staining ( Supplementary Table 5 ) of 100 μL of fresh whole blood consisted of two 15-min incubations at RT: first with CD45RA-BV785, followed by 10 μL of fixative reagent (buffer 1 from PerFix-nc kit, Beckman Coulter, Brea, CA, USA), respectively. Permeabilizing reagent (buffer 2 from PerFix-nc kit, Beckman Coulter, USA) was added to the intracellular master mix of fluorescently labeled antibodies prepared in advance, and cells were stained intracellularly for 60 min at RT. After incubation, 3 mL of plain 1 × PBS (Invitrogen, USA) was added for 5 min, followed by centrifugation at 500 × g for 5 min. Cells were washed in 3 mL of 1 × R3 reagent (10 × buffer 3 diluted in ddH₂O from PerFix-nc kit, Beckman Coulter, USA), resuspended in 200 μL of 1 × R3 reagent, and kept on ice until acquisition.

RNA-seq was performed on venous whole blood samples collected from participants pre- and post-SEVO. Samples were stabilized using the Tempus Spin RNA Isolation Kit (Applied Biosystems, Waltham, MA, USA), as previously described (14). Base calls were processed to FASTQ files on BaseSpace (Illumina, San Diego, CA, USA), and quality trimming was applied to remove low-confidence base calls from the ends of reads. FASTQ files were aligned to the University of California Santa Cruz (UCSC) Human genome assembly version 38.91 using STAR v.2.4.2a. Duplicate filtering, quality control, and metrics analysis were performed using tools from the Picard 1.134 software suite (1). For quality control, samples were required to have > 1 million human-aligned mapped reads and a median coefficient of variation (CV) coverage < 0.6. Genomic variants called from the RNA-seq reads were used to verify RNA-seq sample identity using kinship comparisons (14). Sex-chromosome-associated gene expression was used to confirm that participant sex matched metadata annotations. Genes were filtered to include those with a trimmed mean of M-values (TMM) normalization count of at least 0.8 in at least two samples and further restricted to protein-coding genes. The final dataset included 19 samples comprising 14,590 genes.

Cryopreserved PBMCs from five individuals (pre- and post-SEVO paired time points) were cultured in RPMI-1640 GlutaMAX medium (Thermo Fisher Scientific, USA) supplemented with penicillin–streptomycin (Thermo Fisher Scientific, USA), Fungizone/amphotericin B (Thermo Fisher Scientific, USA), 10% fetal bovine serum (FBS, Thermo Fisher Scientific, USA), and InVivoMAb antihuman CD40 (Bio X Cell, Lebanon, NH, USA). Cells were either left unstimulated or stimulated overnight at 37°C with 50 ng/mL superantigen enterotoxin B (SEB; Sigma-Aldrich, USA). Following cell stimulation, cells were harvested and stained at 4°C for 30 min with fluorescently labeled antibodies, TotalSeq-C reagents ( Supplementary Table S6 ), and a combination of TotalSeq-Hashtag antibodies from BioLegend and Abcam. After staining, cells were washed twice with FACS buffer (1 × PBS (BioLegend, San Diego, CA, USA) containing 1% human AB serum (Sigma-Aldrich, US) and 2 mM EDTA (Sigma-Aldrich, USA]) by centrifugation at 400 × g for 5 min, followed by staining with live/dead (L/D) marker (7-AAD viability staining solution, BioLegend, USA). Samples were sorted on a BD FACSAria II flow cytometer (gating strategy shown in Supplementary Figure S6 ) and processed using the 10 × Chromium Next GEM Single Cell 5′ sequencing platform with Feature Barcode technology (10× Genomics, USA). Libraries were sequenced on a NextSeq 2000 sequencing system (Illumina, USA).

Principal component analyses (PCA) were performed on centered and scaled frequency data and log2-transformed expression data using the prcomp function from the built-in stats package in RStudio (version 4.2.2). Agglomerative hierarchical clustering analysis based on Euclidean distances was performed using the pheatmap package (version 1.0.12) in RStudio. Clustering was performed between immune parameters/samples on centered and scaled frequency data using Ward’s method, and between the 100 most variable genes/samples on centered and scaled log2-transformed expression data using the complete linkage method. The 100 most variable genes were identified by calculating the variance across all libraries using base R functions for each gene. Genes related to sex were excluded from this list. Comparison of immune cell frequencies between pre- and post-SEVO samples was performed using multiple paired t-tests with or without correction for multiple comparisons using the Bonferroni–Dunn method, in GraphPad Prism version 9.0.0. To assess the effect of gas induction relative to interindividual variation, a linear mixed model was applied. The ratio of variance explained by inhaled anesthetic induction versus variance explained by interindividual differences was calculated (15). Volcano plots were generated using the R package ggplot2 (version 3.4.4) on differentially expressed gene results from limma (16). Single-cell RNA-sequencing FASTQ files were generated using Cell Ranger (v7.2.0) mkfastq pipeline from the bcl files. Samples were then processed through the Cell Ranger multipipeline to generate feature-barcode matrices. Dehashing and doublet marking were performed using the Seurat package. Data were imported into Seurat (v4.2.0) to generate a category file compatible with the cloupe file format for downstream analysis using the Loupe Browser (10× Genomics, USA). Differences in the proportion of cells within each of the four clusters identified by graph-based clustering were assessed by Chi-squared analysis. A p-value of < 0.05 was considered statistically significant. Comparisons without labels indicate nonsignificant differences.

Flow cytometric staining was conducted using five panels of monoclonal antibodies to measure the frequency of granulocyte, T cell, B cell, dendritic cell (DC), monocyte, and natural killer (NK) cell subsets, and assess T-cell activation and proliferation (example Boolean gating schemes are shown in Supplementary Figures S1 - S5 ). A total of 100 leukocyte subpopulations were enumerated. To evaluate the effect of anesthetic induction, unbiased hierarchical clustering analysis was conducted on standardized population frequencies from pre- and post-SEVO samples, with the resulting clustered heatmap shown in Figure 1A . In all cases, paired blood samples from each individual clustered together, indicating that samples collected after inhaled anesthetic induction were highly representative of those collected prior to induction. These findings were supported by principal component analysis (PCA), which showed that paired samples were located in close proximity on the PCA plot ( Figure 1B ).

In addition to assessing whether sevoflurane affected the overall composition of the immune system, we evaluated its impact on individual immune cell populations using multiple paired t-tests (with or without correction for multiple comparisons). After applying the Bonferroni–Dunn correction, the frequency of a single population (CD4⁺ T cells) was significantly altered following induction ( Figure 1C ). Without correction for multiple comparisons, this number increased to 22/100 immune cell subsets, including NK cells, monocyte subpopulations, and CD4⁺ and CD8⁺ T cells ( Supplementary Table S7 ; Supplementary Figure S7 ). To contextualize these differences relative to interindividual variation, we applied a linear mixed model, which revealed that for 92/100 cell populations, the variance explained by gas induction was at least 10 times smaller than the variance explained by individual differences ( Supplementary Table S8 ). Four cell populations (% CD56hi NK, % CD56loCD16⁺ NK, % CD69⁺ CD8, and % Total NK) exhibited an effect where gas induction accounted for more than 30% of the variance. Overall, these results suggest that inhaled anesthetic induction has a minimal effect on a limited number of cell subsets, especially when compared to the variation observed between individuals.

When assessing absolute cell counts using the BD Trucount tube (allowing calculation of cells/µL of blood), we observed a small but significant increase in total leukocyte number following sevoflurane induction (mean fold change ± SD: 1.26 ± 0.43). This significant increase was observed across all 17 subpopulations assessed, with the greatest effects seen in NK cells, monocytes, including their subsets (monocyte mean fold change: 3.03 ± 1.68; NK cell mean fold change: 1.95 ± 2.05). These changes also showed the highest interindividual variability ( Supplementary Figures S8, S9 ).

Consistent with the increased fold change in cell number observed by flow cytometry, we also detected a significant increase in PBMC yield following sevoflurane administration (mean ± SD fold change: 1.53 ± 0.71) ( Supplementary Figure S10 ).

Whole blood transcriptional profiling identified a total of 14,590 protein-coding genes, of which only one—BNC2—was significantly altered following sevoflurane induction (fold change = 3.71, FDR = 0.02) ( Figure 2A ). Unbiased hierarchical clustering analysis of the 100 most variable genes was performed to investigate gene expression differences between pre- and post-SEVO samples, with the resulting clustered heatmap shown in Figure 2B . Paired samples from the same individual clustered together, indicating that interindividual variation in gene expression was more pronounced than changes induced by inhaled anesthetic induction. These findings were further supported by PCA, which showed that pre- and post-SEVO samples from each subject were closely clustered on the PCA plot ( Figure 2C ).

To investigate whether sevoflurane affected the functional response of CD4⁺ T cells, we isolated activated T cells using an activation-induced marker (AIM) assay following incubation with the oligoclonal activator SEB. These cells were then profiled using single-cell RNA sequencing (scRNA-seq). As expected, cells from activated and nonactivated culture conditions clustered separately ( Figure 3A ). Analysis by time point revealed no discernible differences in clustering between pre- and post-SEVO samples ( Figure 3B ). To investigate the phenotype of the responding cells, graph-based clustering was performed on SEB-stimulated conditions, which separated the cells into four distinct clusters ( Figure 3C ). Enumeration of the proportion of cells within each cluster revealed no differences between samples collected before and after inhaled anesthetic induction ( Figure 3D ).