As previously reported, our genotype data were generated on genotyping arrays from Illumina, using standard quality assurance cut-offs (detailed in the supplementary methods description) (14). Ancestry was assigned based on principal component analysis. Our previous publication outlines how we defined a PM_2.5 sensitive pediatric asthma phenotype, assessed which SNPs are associated with this PM_2.5 sensitivity in children with asthma, and used this data to create a sPRS (14).

For the current study, we use the Open Targets Genetics Variant to Gene (V2G) tool to assess which genes are implicated by the variants included in the sPRS (18–20). We then use the Metascape and Enrichr tools to assess if any biological pathways are significantly enriched for members of this gene set (21–25). Q-values were adjusted via the Benjamini-Hochberg procedure to account for multiple testing (26, 27). Enriched terms were hierarchically clustered based on Kappa scores, with clusters defined by a similarity threshold of > 0.3 (28). Transcription factor target enrichment was further assessed using the ChEA 2022 dataset via Enrichr, and Metascape was used to perform complete pathway analysis for the target genes of enriched transcription factors.

Patients with asthma and at least one year of follow-up data after their initial asthma diagnosis were selected from the biobank at [Institution censored for review] as described previously (14). All subjects have provided informed consent to both genomic analysis and EMR mining as approved by our institutional IRB. Asthma cases were identified using a previously validated electronic medical record (EMR)-derived phenotype (29). To avoid confounding our results by population structure, only patients from our largest ancestry group (AFR) were selected for this study. Within this group, eight patients with an ancestry-specific sPRS z-score above one were matched with eight patients with an ancestry-specific sPRS z-score below negative one. Patients were matched by birth year, sex, sample collection month, and by living in a non-high-income (median household income as determined using census data <75^th percentile) area with frequent air pollution exposure (incidence of air quality index >50 >75^th percentile). This matching strategy was designed to ensure similar long-term air pollution exposure while minimizing confounding by holding key environmental and demographic variables constant.

Please refer to the supplementary methods section for details regarding sample processing, sequencing, and data processing approach.

Uniform Manifold Approximation and Projection (UMAP) was used to define and visualize cell clusters (30). For cell annotation, SingleR was applied using the celldex::DatabaseImmuneCellExpression reference (31). To assess the robustness of cell-type assignments, complementary validation analyses were performed using two independent immune reference atlases (MonacoImmuneData and Blueprint/ENCODE) (32, 33).

For differential expression analyses, gene-level counts were aggregated at the donor level within each of the 15 SingleR-defined cell types to generate pseudobulk expression matrices. Donor-level differential expression was tested using edgeR’s negative binomial quasi-likelihood framework, with matched pair included as a blocking factor to account for the paired study design (34). Genes were considered differentially expressed if they passed a Benjamini-Hochberg false discovery rate (FDR) threshold of 0.05.

Given the modest sample size and corresponding limitations in power at the donor level, we additionally conducted secondary, exploratory analyses to provide biological context for the observed expression patterns. Within each annotated cell type, genes meeting a nominal significance threshold of P < 0.05 in the pseudobulk analysis were assembled into cell-type–specific gene sets for downstream enrichment analyses. These gene sets were evaluated for enrichment of previously reported PM_2.5-associated transcriptional signatures using the Rummagene tool on the Enrichr platform (query term “PM_2.5”), with enrichment significance evaluated using adjusted P values. Pathway and biological process enrichment analyses were performed using Metascape, with cell type-specific gene lists submitted together, allowing enrichment to be computed per list and recurrent pathway themes to be identified across immune compartments.

We also performed perturbagen enrichment analyses using the LINCS L1000 Signature Search (L2S2) platform accessed through Enrichr. For each immune cell type, nominally significant genes were partitioned into up-regulated and down-regulated sets and submitted separately to identify perturbagens whose transcriptional signatures were either concordant (“Up”) or inversely aligned (“Down”) with the observed expression patterns. All secondary analyses were prespecified as exploratory and hypothesis-generating.

To explore the potential pulmonary relevance of sPRS-associated variants, we queried the GTEx database to identify expression quantitative trait loci (eQTL) in lung tissue. (See acknowledgement section) We then used the LungMAP portal to identify pulmonary cell types known to express genes affected by these pulmonary eQTLs (35, 36). Similarly, as analyses results implicated transforming growth factor beta (TGF-β1), we used LungMAP to identify pulmonary cell types known to express TGF-β1.

Of the 52 variants comprising the sPRS, 44 were resolvable to rsIDs, and 43 of these were successfully linked to candidate genes using the V2G tool (Table 1, Supplementary Table 1). Notably, a disproportionate number of these genes were identified as targets of the transcription factors RUNX2, SMAD2/3, or PAX3-FKHR (q-value 3.19 x 10^-6, 2.18 x 10^-5, and 7.5 x 10^-5, respectively, Figure 1) (37, 38). Although the target networks for RUNX2 and PAX3-FKHR were originally defined in cancer contexts and thus may not be directly relevant to pediatric asthma, these associations may suggest the involvement of similar downstream cascades.

To further investigate this, we conducted a Metascape enrichment analysis on the gene sets corresponding to RUNX2 and PAX3-FKHR target genes. This analysis revealed significant enrichment of the term “GO:0043408 regulation of MAPK cascade” for the target genes of both RUNX2 and PAX3-FKHR (P = 4.56 x 10–¹⁷ and P = 2.72 x 10^-11, respectively, Supplementary Table 2). Taken together, the genes implicated by the sPRS converge around the SMAD and MAPK pathways.

Cell-type annotations were validated across independent reference datasets after harmonization to major peripheral blood immune lineages, demonstrating consistent classification of the major immune populations included in downstream analyses. Donor-level pseudobulk differential expression analysis did not identify statistically significant gene-level differences that passed false discovery rate correction in this small, matched sample. To provide biological context for cell-level transcriptional variation, we therefore examined transcriptional patterns across immune populations in a prespecified secondary analysis.

In this secondary analysis, 14 out of 15 cell types had genes meeting the nominal significance threshold (Figure 2, Supplementary Table 3). These gene sets showed overlap with gene sets previously found to have differential expression after PM_2.5 exposure (Table 2, Supplementary Table 4). Because cases and controls were matched on ambient PM_2.5 exposure but selected based on sPRS, these findings are consistent with the hypothesis that individuals with a high sPRS may exhibit greater transcriptional perturbation for a given exposure.

To characterize recurrent biological themes across immune populations, we next evaluated pathway and process enrichment using Metascape, with multiple-testing correction applied. Enriched terms clustered into three broad functional superfamilies (Supplementary Table 5, Table 3).

The first superfamily encompassed pathways related to transcription and protein synthesis, including “intracellular protein transport” (P = 3.55 x 10-15), “regulation of protein modification process” (P = 1.37 x 10-14), “Selenocysteine synthesis” (P = 3.26 x 10-11), “Nervous system development” (P = 7.17 x 10-12), “Vesicle-mediated transport” (P = 8.71 x 10-12), “negative regulation of protein metabolic process” (P = 2.42 x 10-11).

The second superfamily centered on immune system activation, with significant enrichment observed for terms such as “leukocyte activation” (P = 3.42 x 10-22), “Cellular responses to stress” (P = 4.79 x 10-21), “positive regulation of immune response” (P = 1.71 x 10-17), “positive regulation of cytokine production” (P = 3.04 x 10-17), “regulation of leukocyte activation” (P = 6.75 x 10-17), “regulation of cellular response to stress” (P = 3.11 x 10-16), “Adaptive Immune System” (P = 6.18 x 10-15), “SARS-CoV Infections” (P = 2.01 x 10-12).

Finally, the third superfamily was comprised of various signaling pathways. Specifically, “Interleukin-4 and Interleukin-13 signaling” (P = 0.0001), “MAPK family signaling cascades” (P = 1.27 x 10-6), “Signaling by Receptor Tyrosine Kinases” (P = 2.38 x 10-14), “cell surface receptor protein tyrosine kinase signaling pathway” (P = 3.99 x 10-11), “cellular response to cytokine stimulus” (P = 3.93 x 10-12), “regulation of MAPK cascade” (P = 2.11 x 10-11).

To further contextualize these immune transcriptional signatures, we queried the LINCS L1000 Signature Search (L2S2) database to identify small-molecule perturbagens whose induced transcriptional profiles were either concordant (“UP-UP” or “DOWN-DOWN”) or inversely aligned (“UP-DOWN” or “DOWN-UP”) with the observed expression patterns across immune cell types. Perturbagens were prioritized based on replication of enrichment across multiple immune populations, with statistical strength used as a secondary criterion.

Several compounds demonstrated recurrent inverse alignment with the observed gene expression signatures (UP-DOWN quadrant), most notably tretinoin, BMS-387032 (a CDK inhibitor), and mitoxantrone. Each of these agents showed significant enrichment across 7 of the 10 immune cell types in which an L2S2 hit was observed in the UP–DOWN quadrant, with minimum adjusted P values ranging from 4.4 × 10–⁶ to 4.4 × 10^-16. Compounds demonstrating concordant alignment (UP-UP quadrant) across the largest number of immune populations included BRD-K08177763, niclosamide, and tozasertib, each enriched in 7 of 10 cell types with a significant UP-UP signature (minimum adjusted P values between 4.5 × 10–⁹ and 5.3 × 10^-5). In contrast, relatively few perturbagens showed consistent alignment in the DOWN-UP or DOWN-DOWN quadrants, and these signals were limited to two contributing cell types. Consistent with the exploratory nature of this analysis, these results are interpreted as hypothesis-generating pathway concordance signals rather than evidence of therapeutic relevance.

We next examined the impact of the sPRS variants by assessing which of them are expression quantitative trait loci (eQTL) in lung tissue. We found that rs9836522 is an eQTL for PFN2 with the allele conferring higher PM_2.5 sensitivity (the ‘risk allele’) associated with higher PFN2 expression (P = 0.01, Supplementary Figure 1). Furthermore, the risk alleles for rs13221603 and rs10155910 were associated with lower expression of AGMO and TMEM140 respectively (P = 0.021 and P = 0.023 respectively) in lung tissue. To pinpoint the cellular context of these associations, we queried the Human Lung Reference Cell Atlas (version 1.0). PFN2 expression was predominantly observed in ionocytes, pulmonary neuroendocrine cells (PNEC), and Deuterosomal cells (P = 3.18 x 10^-76, P = 1.20 x 10^-65, and P = 3.09 x 10–⁵⁷ respectively). Similarly, AGMO expression was enriched in PNEC and systemic venous endothelial cells (SVEC) (P = 2.29 x 10–⁵ and P = 0.0001), while TMEM140 expression was most prominent in lymphatic endothelial cells and CAP2 cells (P = 3.80 x 10–²⁹ and 3.63 x 10^-18).

Since signaling pathways implicated by the sPRS (SMAD2/3 and MAPK) and by the gene expression analysis (, MAPK, Receptor Tyrosine Kinase receptors, and IL4/IL13) converge around TGF-β1, we also assessed its expression in lung tissue. This analysis revealed that TGF-β1 is preferentially expressed by interstitial macrophages, patrolling monocytes, inflammatory monocytes, cDC1 cells, NK cells, and pDC cells (P = 3.98 x 10^-25, P = 7.94 x 10^-24, P = 7.08 x 10^-12, P = 2.09 x 10^-11, P = 3.39 x 10^-11, and P = 4.27 x 10^-10, respectively). Additionally, we assessed if differentially expressed genes overlapped with gene-sets previously linked to TGF-β1. This was found to be the case for Naïve B-Cells (adjusted P = 5.14 x 10^-20), CD14⁺ monocytes (adjusted P = 7.61 x 10^-9), CD16⁺ monocytes (adjusted P = 4.46 x 10^-8), CD4⁺ naïve T cells (adjusted P = 4.75 x 10^-112), CD8⁺ naïve T cells (adjusted P = 2.80 x 10^-18), and NK cells (adjusted P = 2.19 x 10^-37) (39, 40).