Integrated study of systemic and local airway transcriptomes in asthma reveals causal mediation of systemic effects by airway key drivers
Genome Medicine volume 15, Article number: 71 (2023)
Systemic and local profiles have each been associated with asthma, but parsing causal relationships between system-wide and airway-specific processes can be challenging. We sought to investigate systemic and airway processes in asthma and their causal relationships.
Three hundred forty-one participants with persistent asthma and non-asthmatic controls were recruited and underwent peripheral blood mononuclear cell (PBMC) collection and nasal brushing. Transcriptome-wide RNA sequencing of the PBMC and nasal samples and a series of analyses were then performed using a discovery and independent test set approach at each step to ensure rigor. Analytic steps included differential expression analyses, coexpression and probabilistic causal (Bayesian) network constructions, key driver analyses, and causal mediation models.
Among the 341 participants, the median age was 13 years (IQR = 10–16), 164 (48%) were female, and 200 (58.7%) had persistent asthma with mean Asthma Control Test (ACT) score 16.6 (SD = 4.2). PBMC genes associated with asthma were enriched in co-expression modules for NK cell-mediated cytotoxicity (fold enrichment = 4.5, FDR = 6.47 × 10−32) and interleukin production (fold enrichment = 2.0, FDR = 1.01 × 10−15). Probabilistic causal network and key driver analyses identified NK cell granule protein (NKG7, fold change = 22.7, FDR = 1.02 × 10−31) and perforin (PRF1, fold change = 14.9, FDR = 1.31 × 10−22) as key drivers predicted to causally regulate PBMC asthma modules. Nasal genes associated with asthma were enriched in the tricarboxylic acid (TCA) cycle module (fold enrichment = 7.5 FDR = 5.09 × 10−107), with network analyses identifying G3BP stress granule assembly factor 1 (G3BP1, fold change = 9.1 FDR = 2.77 × 10−5) and InaD-like protein (INADL, fold change = 5.3 FDR = 2.98 × 10−9) as nasal key drivers. Causal mediation analyses revealed that associations between PBMC key drivers and asthma are causally mediated by nasal key drivers (FDR = 0.0076 to 0.015).
Integrated study of the systemic and airway transcriptomes in a well-phenotyped asthma cohort identified causal key drivers of asthma among PBMC and nasal transcripts. Associations between PBMC key drivers and asthma are causally mediated by nasal key drivers.
Asthma is a chronic respiratory disease that affects millions of people of all ages worldwide [1, 2]. Individuals with asthma experience wheezing, cough, chest tightness, and/or shortness of breath that can lead to impaired quality of life, emergency department visits, hospitalizations, and mortality . Although important progress has been made in asthma research, many gaps remain in our mechanistic understanding of this common disease .
Asthma is a heterogeneous disorder of airway hyperresponsiveness  where the presence of inflammatory cells including type 2 cells, eosinophils, and basophils in the local airway is common . Transcriptome studies of upper airway samples from individuals with asthma have helped to characterize local biology associated with asthma and asthma-related phenotypes [6,7,8,9,10,11,12]. However, measures of systemic inflammation, such as eosinophilia and neutrophilia, have also been associated with asthma and asthma phenotypes [13,14,15]. Local inflammation can be sensed by hematopoietic progenitor cells in the bone marrow, which leads to increased programming of myeloid cells that enter the circulation . The recruitment of circulating immune cells to the airway during inflammation  is a cross-talk that bridges systemic and local inflammation. These prior studies suggest a role for systemic inflammatory processes in asthma, which have been further examined via blood transcriptome-based investigations [17,18,19]. Interestingly, although airway and blood gene expression associations with asthma have been respectively reported, prior investigations of asthma phenotypes that concurrently examined blood and airway samples found differential gene expression in the airway only with no detectable differences in blood [20, 21]. The details and causal relationships between systemic and local transcriptomics in asthma merit further examination.
In this study, we hypothesized that (1) both systemic and local gene expression are associated with asthma, and (2) systemic inflammation associated with asthma is causally mediated by airway gene expression. Leveraging parallel peripheral blood mononuclear cell (PBMC) and nasal transcriptome data generated from 341 individuals and using causal network approaches with discovery and tests sets to ensure rigor, we identified, validated, and characterized systemic and local gene signatures of asthma as well as causal relationships between their key drivers.
Figure 1 provides an overview of the study flow. The study encompassed recruitment and sample collection from 341 participants followed by transcriptome data generation, discovery and test set assignment, differential gene expression and weighted gene coexpression network analyses, cellular deconvolution, and probabilistic causal (Bayesian) network and key driver analyses. These steps were carried out for PBMC and nasal samples in parallel, and key drivers identified from the PBMC and nasal data were then jointly analyzed in causal mediation models (Fig. 1). The methods for each of these steps are detailed below.
Study population and sample collection
Three hundred forty-one participants with and without asthma were recruited from the Mount Sinai Health System, New York, NY. All participants or parents of minors provided written informed consent for study participation. The study conformed to the principles of the Helsinki Declaration and was approved by the Mount Sinai Institutional Review Board (Study 15–00202). Inclusion criteria included those with persistent asthma (based on physician diagnosis, asthma symptoms ≥ 2x/week, and demonstration of a bronchodilator response or positive methacholine challenge) and controls without asthma. Non-asthmatic controls had normal spirometry with no bronchodilator response and no personal nor family history of asthma. Questionnaires addressing asthma-related symptoms and history were completed by all participants. Allergen sensitization by serum specific IgE measurement to 10 environmental allergens, the asthma control test (ACT) , pre- and post-bronchodilator spirometry lung function testing following American Thoracic Society guidlines , and peripheral blood collection were also performed on all subjects. All participants were off asthma medications for at least 4 weeks and nasal medications for at least 2 weeks at the time of sampling. PBMCs were immediately isolated from whole blood samples by Ficoll-Paque density gradient centrifugation and then cryopreserved. Additionally, all subjects were invited to undergo nasal brushing with a sterile cytology brush, with 292 subjects (Additional file 1: Table S1) agreeing to this additional procedure. Nasal brushings were immediately placed in RNALater (Qiagen, Valencia, CA) and stored at – 80 °C.
Transcriptome data generation
RNA from PBMC and nasal samples was extracted using Qiagen RNeasy Mini Kit (Valencia, CA). RNA quality and quantity were measured using 2100 Bioanalyzer and Qubit fluorometer. Sequencing libraries for PBMC RNA were prepared using the Ribo-Zero Gold kit (Illumina) and sequenced on Illumina NovaSeq (Illumina) with paired-end 150 bp reads generated with 40–50 million reads per sample. Sequencing libraries for nasal RNA were prepared using TruSeq RNA Sample Prep Kit v2 protocol (Illumina) and sequenced on Illumina HiSeq 2500. Paired-end 100 bp reads were generated with 40–50 million reads per sample. The reads were mapped to GRCH38 using STAR v2.4.0g1 aligner, and the number of reads mapped per gene was estimated using HTseq . The gene expression profiles were normalized to counts per million (cpm) using edgeR r package . Genes with ≤ 5 counts per million in > 10% of samples were removed to reduce noise on low counts and low abundance genes. After quality control and filtering, 12484 PBMC genes and 13996 nasal genes remained for analysis.
Discovery and test set assignment
The 341 participants were randomly split 2:1 into a discovery set (n = 228) and test set (n = 113). At each major step, the discovery set was used for initial models, and significant findings (FDR ≤ 0.05) that could also be found in the test set of independent participants were considered validated and carried forward.
Differential gene expression analysis and WGCNA
DESeq2  was used to identify PBMC and nasal genes associated with asthma. Age, sex, and race were included as covariates in the models. Genes with Benjamini–Hochberg corrected p values ≤ 0.05 in the discovery set and associated with asthma with the same direction of effect in the test set were deemed validated PBMC or nasal asthma genes, respectively.
Weighted gene coexpression network analysis (WGCNA)  was used to build PBMC and nasal gene co-expression networks separately. PBMC and nasal gene co-expression modules enriched for PBMC or nasal asthma genes were then identified by Fisher’s exact test (enrichment score > 1 and FDR ≤ 0.05). Gene ontology (GO) analysis was performed on each module’s gene set using DAVID , and biological processes ranked by fold enrichment and FDR ≤ 0.05 were selected as module names.
Probabilistic causal network and key driver analysis
Probabilistic causal (Bayesian) networks [31, 32] for the PBMC and nasal transcriptome data were built separately for each’s respective discovery and test sets using RIMBAnet [31, 32], a software to construct probabilistic causal networks using Markov chain Monte Carlo simulations. Inputs included the corresponding PBMC or nasal transcriptome data and eQTL data . To map eQTLs for the PBMC and nasal transcriptomes respectively, two-path mode from STAR aligner  was used to align RNAseq reads to the reference genome (GRCh38), and variants were called using the “HaplotypeCaller” GATK  tool after filtering against known RNA editing sites, intronic variants within four nucleotides from splice donor, and acceptor site variants called in repeated sequences. We performed genotype phasing and imputation using BEAGLE 5.1  against the 1000 Genome reference haplotypes for hg38 . After filtering out variants with low imputation accuracy (dosage R-squared < 0.8) and low minor allele frequency (MAF < 0.01), we hard called the variants using PLINK2 . Genetic ancestry of the study population was calculated based on the variant calls , and this genetic ancestry together with age and gender were used as covariates when mapping cis eQTLs [33, 40]. We identified the eQTLs where cis SNPs (within 1 Mb of the transcription start or end site of the gene) were associated with gene expression at FDR ≤ 0.05 using matrixeqtl  (Additional file 1: Tables S2-S4).
To construct the probabilistic causal networks, genes from the corresponding transcriptomic datasets were discretized into 3 states, high expression, low expression, and no expression by K-means clustering (K = 3) [31, 32, 41]. Markov chain Monte Carlo simulations were used to reconstruct 1000 networks, and the fit of each network was assessed by the Bayesian Information Criterion [31, 32, 41]. The final causal network was built by retaining edges present in ≥ 30% of the 1000 networks.
Key driver analysis using the KDA package  was then performed to identify each asthma module’s key drivers using the causal network and module members as targets. Subnetworks used as background for the enrichment analysis were identified by selecting nodes K-steps away from the nodes in the module member list. The enrichment of module member genes was assessed in k-step (k varies from 1 to K) downstream neighborhood stemming from each node. We used K = 7 in the analyses.
Causal mediation analysis
PBMC and nasal key drivers associated with asthma were tested in causal mediation analysis using the robust structural modeling equation implemented in the lavaan R package . We tested (a) the degree to which PBMC key drivers (mediator) mediated the effects of nasal key drivers (independent variable) on asthma (outcome) and (b) the degree to which nasal key drivers (mediator) mediated the effects of PBMC key drivers (independent variable) on asthma (outcome). For each causal mediation analysis, three regression models were implemented to estimate the effects of (1) independent variables on outcomes, (2) independent variables on mediators, and (3) mediators on the outcomes. Model 1 estimates the direct effect of independent variables on the outcome, while models (2) and (3) estimate indirect effects (mediation). FDR \(\le\) 0.05 was used as the threshold for significance.
The cohort included 341 participants recruited from the Mount Sinai Health System, New York, USA, of whom 200 (58.7%) individuals had persistent asthma based on physician diagnosis, symptoms ≥ 2x/week, and demonstration of bronchodilator response on lung function testing and/or positive methacholine challenge. One hundred forty-one (41.3%) individuals had no asthma based on no personal or family history of asthma and demonstration of normal spirometry without a bronchodilator response (Table 1). Participants were primarily children with median age of 13 years (IQR 10–16). Those with asthma were younger, more likely to be Black or Latino, and had lower ACT scores, FEV1%, and FEV1/FVC values. Many participants with asthma had high rates of emergency department visits and hospitalizations for asthma (Table 1). Peripheral blood for PBMC isolation was collected from all 341 participants. Nasal brushings were collected from a subset (n = 292, 85.6%) who agreed to nasal brushing; the characteristics of this subset did not significantly differ from the whole cohort (Additional file 1: Table S1).
Systemic processes in asthma: PBMC asthma genes are enriched in co-expression modules for NK cell-mediated cytotoxicity and interleukin production
To investigate systemic processes in asthma, we isolated PBMCs from participants’ peripheral blood and generated PBMC transcriptome profiles using RNA sequencing . The 341 participants were randomly split into a discovery set (n = 228) and test set (n = 113). To identify PBMC transcripts significantly associated (FDR ≤ 0.05) with asthma, differential gene expression analysis on the discovery set was performed with age, sex, and race/ethnicity as covariates. PBMC transcripts associated with asthma in the discovery set (Fig. 2) were then tested in the independent test set, and those that were also associated with the same direction of effect in the test set were deemed validated “PBMC asthma genes” (Additional file 1: Table S5).
We then sought to uncover functional biological context for the identified PBMC asthma genes. Weighted gene co-expression network analysis (WGCNA)  using the discovery set identified 24 PBMC co-expression modules representing broader constructs of biological function. Nine of these 24 modules were significantly enriched (Fisher’s exact test FDR ≤ 0.05) for PBMC asthma genes. Gene Ontology (GO) enrichment analyses revealed that the modules most strongly enriched for PBMC asthma genes (ranked by effect size and then lowest FDR) were the modules for “NK cell-mediated cytotoxicity” (fold enrichment = 4.5, FDR = 6.47 × 10−32) and “interleukin production” (fold enrichment = 2.0, FDR = 1.01 × 10−15). We will henceforth refer to these two modules as the NK cell and interleukin PBMC asthma modules, and their respective member genes are shown in Additional file 1: Tables S6, S7). The remaining seven modules largely represented cellular maintenance processes. Notably, genes in causal paths for Th2 cell differentiation regulation gene ontology process were enriched in the interleukin PBMC asthma module (enrichment score 19.0, Fisher’s exact p value = 7.24e − 23).
Causal relationships and key drivers of the PBMC asthma modules
To characterize causal relationships among genes within the identified PBMC asthma modules, we next performed probabilistic causal (Bayesian) network analysis [31, 32, 45]. In this network construction, causal relationships were statistically inferred between PBMC transcripts using expression quantitative trait loci mapped in this cohort as priors [31,32,33] (Additional file 1: Tables S3-S4). We then used key driver analysis (KDA) , a method for applying dynamic and statistical neighborhood searches on constructed probabilistic causal networks, to identify key drivers predicted to exert the greatest causal impact on downstream genes in each PBMC asthma module. Here, we also used a discovery and test set approach, where a probabilistic causal network was first built with the discovery set and key drivers were identified for each PBMC asthma module (Additional file 1: Table S8, Table S9). Independent probabilistic causal network construction with the test set and key driver analysis were then separately performed for each module. Key drivers identified in the discovery set and also found in the test set were deemed validated “PBMC key drivers” of each module.
The PBMC key drivers for the NK cell module are shown in Fig. 3A and for the interleukin module in Fig. 4A. As key drivers predicted to exert the greatest causal impact on downstream genes in each module, these key drivers appear at the top of each Eiffel tower plot, with causality flowing from top to bottom. In addition to the key drivers, we label a few downstream genes of interest for each module. Biological relationships between the key drivers for each module, based on the probabilistic causal network and prior experimental work, are shown in Figs. 3B and 4B [46,47,48,49,50,51,52,53,54,55,56,57,58,59,60,61].
Local airway processes in asthma: nasal asthma genes, nasal asthma modules, and nasal key drivers
To characterize local airway processes in asthma, we next performed the same analytic flow with the nasal transcriptome data  that had been generated in parallel with the PBMC data (Fig. 1). The same discovery and test set assignments were used for the available nasal transcriptome data. Nasally expressed genes found to be associated with asthma by differential expression analysis in the discovery set that were also significantly associated with asthma in the independent test set were deemed “nasal asthma genes” (Additional file 1: Table S10).
To explore how immune cells in the nasal samples might influence the nasal asthma genes identified, we performed cellular deconvolution  and found small cell fractions of myeloid, mast cells, and plasma cells in the nasal samples that did not differ between subjects with and without asthma (Additional file 1: Table S11). Adding these cell fractions as additional model covariates did not significantly change the nasal genes identified; nasal asthma genes from this sensitivity analysis strongly overlapped with those found from the original model (enrichment score 12.7, Fisher’s exact p value = 6.3e − 306).
WGCNA using the nasal transcriptome discovery set identified 21 nasal modules, of which 8 were significantly enriched with nasal asthma genes (FDR ≤ 0.05). Of these 8 modules, 2 had module eigenvalues significantly associated with asthma at FDR ≤ 0.05: the module for tricarboxylic acid (TCA) cycle and the module for metabolic pathways (Additional file 2: Fig. S1). The TCA module eigenvalue was associated with asthma with greatest effect size and lowest FDR (Additional file 2: Fig. S1) and was therefore considered the top ranked “nasal asthma module.” The metabolic pathways module, whose eigenvalue was also nominally associated with asthma but with lower effect size and higher FDR, was considered a broad module representing cellular maintenance.
Construction of a probabilistic causal network with the nasal transcriptome discovery set and key driver analysis for the TCA nasal asthma module identified nasal key drivers. An independent probabilistic causal network was then built with the nasal transcriptome test set with key driver analysis also performed. Key drivers identified for the TCA nasal asthma module in the discovery set (Additional file 1: Table S12) that were also validated in the test set were considered “nasal key drivers.” These nasal key drivers included G3BP1 and INADL, two members of the TCA module that were also nasal asthma genes themselves (Additional file 1: Table S10).
Relationship between systemic and local processes in asthma: nasal key drivers causally mediate associations between PBMC key drivers and asthma
Our next goal was to investigate the relationship between systemic and local processes in asthma. Specifically, we sought to examine causal mediation between PBMC key drivers associated with asthma (PRF1 and NKG7 in the NK cell module; CAPZA2, ATP6AP2, CTSS, and RAB3D from the interleukin module) and nasal key drivers associated with asthma (G3BP1 and INADL from the nasal TCA module). These causal mediation models revealed that the nasal key drivers of asthma significantly mediate the association between PBMC key drivers in the NK cell module (PRF1 and NKG7) and asthma (FDR = 0.0076 to 0.01) (Fig. 5). No causal mediation was observed in the interleukin module, and there was no statistically significant finding for the converse of PBMC key drivers mediating associations between nasal key drivers and asthma. Permutation testing  with 1,000,000 iterations confirmed our finding of nasal key drivers significantly mediating associations between PBMC key drivers and asthma in the NK cell module and lack of significance for mediation by PBMC key drivers.
In this integrated study of systemic and local processes in asthma, we generated and leveraged parallel PBMC and nasal transcriptome profiles  from 341 individuals in a well-characterized asthma cohort to identify systemic and local key drivers of asthma, as well as causal mediation relationships underlying systemic-local cross-talk in this complex disease. By capturing both systemic and local transcriptomics and employing innovative network analyses, this study moves beyond prior transcriptomic studies of asthma by (1) finding evidence for both systemic and local transcriptional signatures of asthma in the same cohort (Additional file 1: Table S5, Additional file 1: Table S10); (2) identifying key drivers of these asthma signatures in the PBMC and airway compartments (Figs. 3 and 4, Additional file 1: Table S8, Table S9, Table S12); and (3) characterizing causal relationships between the detected systemic and local key drivers (Fig. 5). The discovery and test set approach that we employed, where each discovery step required validation in a test set of independent participants, brings rigor to the findings reported.
Through differential gene expression, WGCNA, and module enrichment analyses, we found that NK cell-mediated cytotoxicity and interleukin production were the biological processes for the top PBMC modules enriched for asthma-associated PBMC genes, while TCA cycle was the top nasal module enriched for asthma-associated nasal genes. These findings from unbiased whole transcriptome analyses of the PBMC and nasal compartments highlight mechanistic pathways that have to date been overshadowed by focus on commonly targeted pathways (e.g. type 2 inflammation) . The importance of these biological processes is supported by prior findings. For example, increased NK cell cytotoxicity activity is present in the peripheral blood of asthmatics vs. controls, and children with acute exacerbated asthma have higher numbers of circulating NK cells compared to those with stable asthma . NK cells activate in response to well-known allergy triggers including allergens, RSV, and other respiratory viruses . Furthermore, compared to controls, patients with asthma have increased IL-4-producing NK cells in their blood . Our finding of interleukin production as a top PBMC asthma module is not surprising, given the broad representation of this module and the pervasive roles that interleukins play in orchestrating inflammatory processes in asthma, including IL-4, IL-5, and IL-13 in type 2 inflammation [65, 66]. Indeed, the pivotal role of circulating interleukins has been the basis for the explosion of biologic medications as effective systemic treatments for asthma . Moreover, we found significant enrichment of Th2 cell differentiation genes in the interleukin PBMC asthma module. In contrast, finding that the TCA module was the top nasal asthma module was less expected, although evidence from the cellular and molecular biology field supports that members of the TCA cycle act mechanistically to shape smooth muscle airway bronchoconstriction, hyperresponsiveness, and airway remodeling in asthma [68, 69].
While illuminating to identify genes and modules associated with asthma in both the PBMC and nasal compartments, we thought it would be important to move beyond association toward new insight on key drivers of systemic and local processes in asthma. The identification of key drivers prioritizes gene signature lists by highlighting the transcripts at root levels of causality for more efficient therapeutic targeting and mechanistic parsing. Here, we mapped and leveraged expression quantitative trait loci as priors to build probabilistic causal networks with directional connectivity, which then enabled us to identify key drivers of the PBMC and nasal asthma modules (Figs. 3 and 4, Additional file 1: Table S8, Table S9, Table S12). These key drivers are statistically predicted to drive the regulatory state of each asthma module. For the NK cell PBMC asthma module, we identified perforin (PRF1) and NK cell granule protein (NKG7) as key drivers that were each also associated with asthma (Fig. 3A). Increased expression of PRF1 in peripheral blood lymphocytes of both allergic and intrinsic asthmatics has been previously reported,  and NKG7 was found in a study of publicly available gene expression data to be a marker of severe asthma . Both PRF1 and NKG7 demonstrated prominent upstream roles in NK cell-mediated cytotoxicity when causal relationships derived from the constructed network and evidence from the literature were combined to provide functional biologic context (Fig. 3B). Given the broader scope of the interleukin production PBMC asthma module, more key drivers were identified for this module (Fig. 4A), which were predicted by the causal network and prior literature to regulate a wide variety of cellular activities including those involved in prevalent type 2 asthma, such as allergen presentation and eosinophil chemotaxis (Fig. 4B). For the nasal asthma module, G3BP stress granule assembly factor 1 (G3BP1), which promotes innate immune transcriptional responses via NF-kB and c-Jun N-terminal kinase pathways,  and the protein-coding gene InaD-like protein (INADL) involved in epithelial migration  were both identified as nasal key drivers that were also associated with asthma.
In contrast to prior investigations of transcriptomic associations with asthma and asthma-related phenotypes in blood only [17,18,19], upper airway samples only [6,7,8,9,10,11], and both where only airway but no systemic associations were found [20, 21], we found evidence for both systemic and airway transcriptomic associations with asthma in the same cohort. Our sample size of 341 participants was relatively large for a dual transcriptome study [12, 20, 21], augmenting our power for detection and also enabling the network and discovery/test set approaches we used to ensure rigor and advance insight beyond associations. Our paired findings from the systems and local domains provided us with a unique opportunity to begin to address the fundamental question of the degree to which asthma is a local process with systemic findings or a systemic process with local findings. While this is a somewhat existential question without an absolute answer, our causal mediation analyses of the relationships between the identified PBMC asthma key drivers, nasal asthma key drivers, and asthma status in this population begin to offer some insight. Here, we found that associations between PBMC key drivers and asthma were causally mediated by nasal key drivers, but not the converse (Fig. 5). Specifically, for the PBMC module key drivers PRF1 and NKG7, their association with asthma was causally mediated by the nasal asthma key drivers G3BP1 and INADL. There was no evidence for the converse of PBMC key drivers causally mediating associations between nasal key drivers and asthma. This suggests that asthma is a systemic process causally mediated by airway key drivers. Given circulating immune cells can be recruited to the airway during inflammation [13, 16], our findings speak to the cross-talk that bridges systemic and local inflammation and immunity. While we found significant causal mediation relationships for PBMC key drivers of the NK cell module (PRF1 and NKG7), there were no detectable causal mediation effects for PBMC key drivers of the interleukin module, which may have been due to the broader, less specific scope of that module.
We recognize the limitations to our study. While our investigation of parallel PBMC and nasal transcriptome profiles in asthma stands apart in successfully finding both systemic and airway gene signatures, modules, and key drivers of asthma as well as causal relationships between them, our study did not endotype asthma. It is possible that different findings would be detected if analyses were limited to particular endotypes of asthma, which was outside the scope of this study but a future direction we will consider. Additionally, while we employed a discovery and test set approach where each was analyzed independently at every step to ensure rigor in our reported findings, the discovery and test sets were sampled from the same population. It is possible that findings would be distinct for asthma populations in other parts of the world. Finally, the key drivers were identified via mathematical and statistical models, and we acknowledge that further work is needed to replicate and experimentally validate these findings to deepen our understanding of their role in asthma.
In this network study of systemic and local processes in asthma, we generated and examined dual PBMC and nasal transcriptome profiles from 341 individuals to identify systemic and local key drivers of asthma, as well as causal mediation relationships between them. In addition to identifying both systemic and local transcriptional signatures of asthma, we found that associations between systemic key drivers of asthma and asthma are causally mediated by nasal key drivers of asthma. Our study moves beyond anatomically isolated transcriptome-wide association analyses to elucidate causal relationships between systemic and airway key drivers of asthma.
Asher MI, Rutter CE, Bissell K, Chiang C-Y, El Sony A, Ellwood E, et al. Worldwide trends in the burden of asthma symptoms in school-aged children: Global Asthma Network Phase I cross-sectional study. Lancet. 2021;398:1569–80.
Asher MI, Garcia-Marcos L, Pearce NE, Strachan DP. Trends in worldwide asthma prevalence. Eur Respir J. 2020;56:2002094.
Porsbjerg C, Melen E, Lehtimaki L, Shaw D. Asthma. Lancet. 2023;401:858–73.
Edwards MR, Saglani S, Schwarze J, Skevaki C, Smith JA, Ainsworth B, et al. Addressing unmet needs in understanding asthma mechanisms: From the European Asthma Research and Innovation Partnership (EARIP) Work Package (WP)2 collaborators. Eur Respir J. 2017;49:1602448.
Gauvreau GM, El-Gammal AI, O’Byrne PM. Allergen-induced airway responses. Eur Respir J. 2015;46:819–31.
Sajuthi SP, Everman JL, Jackson ND, Saef B, Rios CL, Moore CM, et al. Nasal airway transcriptome-wide association study of asthma reveals genetically driven mucus pathobiology. Nat Commun. 2022;13:1632.
Zhu Z, Camargo CA Jr, Raita Y, Freishtat RJ, Fujiogi M, Hahn A, et al. Nasopharyngeal airway dual-transcriptome of infants with severe bronchiolitis and risk of childhood asthma: a multicenter prospective study. J Allergy Clin Immunol. 2022;150:806–16.
Do AN, Chun Y, Grishina G, Grishin A, Rogers AJ, Raby BA, et al. Network study of nasal transcriptome profiles reveals master regulator genes of asthma. J Allergy Clin Immunol. 2021;147:879–93.
Chun Y, Do A, Grishina G, Grishin A, Fang G, Rose S, et al. Integrative study of the upper and lower airway microbiome and transcriptome in asthma. JCI Insight. 2020;5:e133707.
Forno E, Zhang R, Jiang Y, Kim S, Yan Q, Ren Z, et al. Transcriptome-wide and differential expression network analyses of childhood asthma in nasal epithelium. J Allergy Clin Immunol. 2020;146:671–5.
Hekking PP, Loza MJ, Pavlidis S, de Meulder B, Lefaudeux D, Baribaud F, et al. Pathway discovery using transcriptomic profiles in adult-onset severe asthma. J Allergy Clin Immunol. 2018;141:1280–90.
Park HW, Weiss ST. Understanding the molecular mechanisms of asthma through transcriptomics. Allergy Asthma Immunol Res. 2020;12:399–411.
Crisford H, Sapey E, Rogers GB, Taylor S, Nagakumar P, Lokwani R, et al. Neutrophils in asthma: the good, the bad and the bacteria. Thorax. 2021;76:835–44.
Carr TF, Zeki AA, Kraft M. Eosinophilic and noneosinophilic asthma. Am J Respir Crit Care Med. 2018;197:22–37.
Jackson DJ, Akuthota P, Roufosse F. Eosinophils and eosinophilic immune dysfunction in health and disease. Eur Respir Rev. 2022;31:210150.
Chavakis T, Mitroulis I, Hajishengallis G. Hematopoietic progenitor cells as integrative hubs for adaptation to and fine-tuning of inflammation. Nat Immunol. 2019;20:802–11.
Bigler J, Boedigheimer MJ, Schofield JPR, Skipp PJ, Corfield J, Rowe A, et al. A severe asthma disease signature from gene expression profiling of peripheral blood from U-BIOPRED cohorts. Am J Respir Crit Care Med. 2017;195:1311–20.
Lemonnier N, Melen E, Jiang Y, Joly S, Menard C, Aguilar D, et al. A novel whole blood gene expression signature for asthma, dermatitis, and rhinitis multimorbidity in children and adolescents. Allergy. 2020;75:3248–60.
Jiang Y, Gruzieva O, Wang T, Forno E, Boutaoui N, Sun T, et al. Transcriptomics of atopy and atopic asthma in white blood cells from children and adolescents. Eur Respir J. 2019;53:1900102.
Hoda U, Pavlidis S, Bansal AT, Takahashi K, Hu S, Ng Kee Kwong F, et al. Clinical and transcriptomic features of persistent exacerbation-prone severe asthma in U-BIOPRED cohort. Clin Transl Med. 2022;12:e816.
Altman MC, Gill MA, Whalen E, Babineau DC, Shao B, Liu AH, et al. Transcriptome networks identify mechanisms of viral and nonviral asthma exacerbations in children. Nat Immunol. 2019;20:637–51.
Schatz M, Kosinski M, Yarlas AS, Hanlon J, Watson ME, Jhingran P. The minimally important difference of the Asthma Control Test. J Allergy Clin Immunol. 2009;124(719–23): e1.
Graham BL, Steenbruggen I, Miller MR, Barjaktarevic IZ, Cooper BG, Hall GL, et al. Standardization of Spirometry 2019 Update. An Official American Thoracic Society and European Respiratory Society Technical Statement. Am J Respir Crit Care Med. 2019;200:e70-e88.
Anders S, Pyl PT, Huber W. HTSeq–a Python framework to work with high-throughput sequencing data. Bioinformatics. 2015;31:166–9.
Robinson MD, McCarthy DJ, Smyth GK. edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics. 2010;26:139–40.
Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15:550.
Langfelder P, Horvath S. WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics. 2008;9:559.
Sherman BT, Hao M, Qiu J, Jiao X, Baseler MW, Lane HC, et al. DAVID: a web server for functional enrichment analysis and functional annotation of gene lists (2021 update). Nucleic Acids Res. 2022;50:W216–W221.
Newman AM, Steen CB, Liu CL, Gentles AJ, Chaudhuri AA, Scherer F, et al. Determining cell type abundance and expression from bulk tissues with digital cytometry. Nat Biotechnol. 2019;37:773–82.
Ordovas-Montanes J, Dwyer DF, Nyquist SK, Buchheit KM, Vukovic M, Deb C, et al. Allergic inflammatory memory in human respiratory epithelial progenitor cells. Nature. 2018;560:649–54.
Zhu J, Sova P, Xu Q, Dombek KM, Xu EY, Vu H, et al. Stitching together multiple data dimensions reveals interacting metabolomic and transcriptomic networks that modulate cell regulation. PLoS Biol. 2012;10: e1001301.
Zhu J, Zhang B, Smith EN, Drees B, Brem RB, Kruglyak L, et al. Integrating large-scale functional genomic data to dissect the complexity of yeast regulatory networks. Nat Genet. 2008;40:854–61.
Irizar H, Chun Y, Arditi Z, Do A, Grishina G, Grishin A, et al. Examination of host genetic effects on nasal microbiome composition. J Allergy Clin Immunol. 2022;150:1232–6.
Dobin A, Davis CA, Schlesinger F, Drenkow J, Zaleski C, Jha S, et al. STAR: ultrafast universal RNA-seq aligner. Bioinformatics. 2013;29:15–21.
van der Auwera GA, O’Connor BD. Genomics in the Cloud: Using Docker, GATK, and WDL in Terra. 1st ed. O’Reilly Media; 2020.
Browning SR, Browning BL. Rapid and accurate haplotype phasing and missing-data inference for whole-genome association studies by use of localized haplotype clustering. Am J Hum Genet. 2007;81:1084–97.
Browning BL, Zhou Y, Browning SR. A one-penny imputed genome from next-generation reference panels. Am J Hum Genet. 2018;103:338–48.
Chang CC, Chow CC, Tellier LC, Vattikuti S, Purcell SM, Lee JJ. Second-generation PLINK: rising to the challenge of larger and richer datasets. GigaScience. 2015;4:7.
Hart T, Moffat J. BAGEL: a computational framework for identifying essential genes from pooled library screens. BMC Bioinformatics. 2016;17:164.
Shabalin AA. Matrix eQTL: ultra fast eQTL analysis via large matrix operations. Bioinformatics. 2012;28:1353–8.
Zhu J, Wiener MC, Zhang C, Fridman A, Minch E, Lum PY, et al. Increasing the power to detect causal associations by combining genotypic and expression data in segregating populations. PLoS Comput Biol. 2007;3: e69.
Zhang B, Zhu J. Identification of key causal regulators in gene networks. Proceedings of the World Congress on Engineering. London, UK. 2013;II:5–8.
Rosseel Y. lavaan: An R Package for Structural Equation Modeling. J Stat Softw. 2012;48:1–36.
Zhang L, Chun Y, Irizar H, Arditi Z, Grishina G, Grishin A, et al. Integrated study of systemic and local airway transcriptomes in asthma reveals causal mediation of systemic effects by airway key drivers. Synapse. 2023. https://doi.org/10.7303/syn51061925.
Watson CT, Cohain AT, Griffin RS, Chun Y, Grishin A, Hacyznska H, et al. Integrative transcriptomic analysis reveals key drivers of acute peanut allergic reactions. Nat Commun. 2017;8:1943.
Ogawa K, Tanaka K, Ishii A, Nakamura Y, Kondo S, Sugamura K, et al. A novel serum protein that is selectively produced by cytotoxic lymphocytes. J Immunol. 2001;166:6404–12.
Voskoboinik I, Whisstock JC, Trapani JA. Perforin and granzymes: function, dysfunction and human pathology. Nat Rev Immunol. 2015;15:388–400.
Ng SS, De Labastida RF, Yan J, Corvino D, Das I, Zhang P, et al. The NK cell granule protein NKG7 regulates cytotoxic granule exocytosis and inflammation. Nat Immunol. 2020;21:1205–18.
Dotiwala F, Mulik S, Polidoro RB, Ansara JA, Burleigh BA, Walch M, et al. Killer lymphocytes use granulysin, perforin and granzymes to kill intracellular parasites. Nat Med. 2016;22:210–6.
Kishimoto Y, Hiraiwa M, O’Brien JS. Saposins: structure, function, distribution, and molecular genetics. J Lipid Res. 1992;33:1255–67.
Pazdrak K, Young TW, Straub C, Stafford S, Kurosky A. Priming of eosinophils by GM-CSF is mediated by protein kinase CbetaII-phosphorylated L-plastin. J Immunol. 2011;186:6485–96.
Ciaffoni F, Tatti M, Boe A, Salvioli R, Fluharty A, Sonnino S, et al. Saposin B binds and transfers phospholipids. J Lipid Res. 2006;47:1045–53.
Ciciro Y, Sala A. MYB oncoproteins: emerging players and potential therapeutic targets in human cancer. Oncogenesis. 2021;10:19.
Rodriguez-Cruz A, Vesin D, Ramon-Luing L, Zuniga J, Quesniaux VFJ, Ryffel B, et al. CD3(+) macrophages deliver proinflammatory cytokines by a CD3- and transmembrane TNF-dependent pathway and are increased at the BCG-infection site. Front Immunol. 2019;10:2550.
Song Q, Meng B, Xu H, Mao Z. The emerging roles of vacuolar-type ATPase-dependent Lysosomal acidification in neurodegenerative diseases. Transl Neurodegener. 2020;9:17.
Tian X, Jin RU, Bredemeyer AJ, Oates EJ, Blazewska KM, McKenna CE, et al. RAB26 and RAB3D are direct transcriptional targets of MIST1 that regulate exocrine granule maturation. Mol Cell Biol. 2010;30:1269–84.
Huang Y, Mao X, van Jaarsveld RH, Shu L, Terhal PA, Jia Z, et al. Variants in CAPZA2, a member of an F-actin capping complex, cause intellectual disability and developmental delay. Hum Mol Genet. 2020;29:1537–46.
Bhagat G, Naiyer AJ, Shah JG, Harper J, Jabri B, Wang TC, et al. Small intestinal CD8+TCRgammadelta+NKG2A+ intraepithelial lymphocytes have attributes of regulatory cells in patients with celiac disease. J Clin Invest. 2008;118:281–93.
Braud VM, Allan DSJ, O’Callaghan CA, derstro¨m KS, D’Andrea A, Ogg GS, et al. HLA-E binds to natural killer cell receptors CD94:NKG2A, B and C. Nature. 1998;391:795–9.
Brown R, Nath S, Lora A, Samaha G, Elgamal Z, Kaiser R, et al. Cathepsin S: investigating an old player in lung disease pathogenesis, comorbidities, and potential therapeutics. Respir Res. 2020;21:111.
Yeh TW, Okano T, Naruto T, Yamashita M, Okamura M, Tanita K, et al. APRIL-dependent lifelong plasmacyte maintenance and immunoglobulin production in humans. J Allergy Clin Immunol. 2020;146(1109–20): e4.
Schadt EE, Molony C, Chudin E, Hao K, Yang X, Lum PY, et al. Mapping the genetic architecture of gene expression in human liver. PLoS Biol. 2008;6: e107.
Maspero J, Adir Y, Al-Ahmad M, Celis-Preciado CA, Colodenco FD, Giavina-Bianchi P, et al. Type 2 inflammation in asthma and other airway diseases. ERJ Open Res. 2022;8:00576–2021.
Gorska MM. Natural killer cells in asthma. Curr Opin Allergy Clin Immunol. 2017;17:50–4.
Chung KF. Targeting the interleukin pathway in the treatment of asthma. Lancet. 2015;386:1086–96.
Lambrecht BN, Hammad H, Fahy JV. The cytokines of asthma. Immunity. 2019;50:975–91.
Brusselle GG, Koppelman GH. Biologic therapies for severe asthma. N Engl J Med. 2022;386:157–71.
Yeung BHY, Huang J, An SS, Solway J, Mitzner W, Tang WY. Role of isocitrate dehydrogenase 2 on DNA hydroxymethylation in human airway smooth muscle cells. Am J Respir Cell Mol Biol. 2020;63:36–45.
Huang SK. A fresh take on the “TCA” cycle: TETs, citrate, and asthma. Am J Respir Cell Mol Biol. 2020;63:1–3.
Arnold V, Balkow S, Staats R, Matthys H, Luttmann W, Virchow JC Jr. Increase in perforin-positive peripheral blood lymphocytes in extrinsic and intrinsic asthma. Am J Respir Crit Care Med. 2000;161:182–6.
Jiang Y, Deng S, Hu X, Luo L, Zhang Y, Zhang D, et al. Identification of potential biomarkers and immune infiltration characteristics in severe asthma. Int J Immunopathol Pharmacol. 2022;36:3946320221114194.
Reineke LC, Lloyd RE. The stress granule protein G3BP1 recruits protein kinase R to promote multiple innate immune antiviral responses. J Virol. 2015;89:2575–89.
Wells CD, Fawcett JP, Traweger A, Yamanaka Y, Goudreault M, Elder K, et al. A Rich1/Amot complex regulates the Cdc42 GTPase and apical-polarity proteins in epithelial cells. Cell. 2006;125:535–48.
We thank the families who participated in this study.
This study was funded by the National Institutes of Health R01 AI118833.
Ethics approval and consent to participate
All participants or parents of minors provided written informed consent for study participation. The study conformed to the principles of the Helsinki Declaration and was approved by the Mount Sinai Institutional Review Board (Study 15–00202).
Consent for publication
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Table S1. Characteristics of the participants who provided nasal samples. Table S2. NK cell-mediated cytotoxicity PBMC module eQTLs. Table S3. Interleukin production PBMC module eQTLs. Table S4. Nasal TCA cycle module eQTLs. Table S5. PBMC asthma genes. Table S6. PBMC NK cell-mediated cytotoxicity module genes. Table S7. PBMC interleukin production module genes. Table S8. Key driver analysis results for the PBMC NK cell-mediated cytotoxicity module. Table S9. Key driver analysis results for the PBMC interleukin production module. Table S10. Nasal asthma genes. Table S11. Nasal Myeloid, Mast cell and Plasma cell fractions. Table S12. Key driver analysis results for the nasal TCA cycle module.
Fig. S1. Associations between nasal module eigenvalues and asthma.
About this article
Cite this article
Zhang, L., Chun, Y., Irizar, H. et al. Integrated study of systemic and local airway transcriptomes in asthma reveals causal mediation of systemic effects by airway key drivers. Genome Med 15, 71 (2023). https://doi.org/10.1186/s13073-023-01222-2