Comprehensive sequencing of the lung neuroimmune landscape in response to asthmatic induction

TL;DR

This study shows how asthma affects the lung's neuroimmune system by changing gene activity in specific cell clusters and locations.

Contribution

The study reveals spatial and cluster-specific changes in the neuroimmune transcriptome in response to asthma.

Findings

Asthmatic induction increases allergen-sensing genes in nodose/jugular ganglia, with a caudal-to-rostral spatial pattern.

Lung transcriptome changes show upregulated neural excitation genes in neuroendocrine-like cells and neurons.

Distinct nodose neuron clusters align with previously identified cell types, highlighting spatial organization.

Abstract

Evidence demonstrates that sensory neurons respond to pathogenic/allergic infiltration and mediate immune responses, forming an integral part of host defense that becomes hypersensitized during allergy. Our objective was to investigate how asthmatic induction alters the pulmonary neuroimmune transcriptome. We hypothesized that asthmatic induction would upregulate genes in the vagal ganglia (nodose/jugular ganglia), which would be associated with asthmatic immunity, and that these would be clustered, primarily in nodose neurons. Furthermore, lungs would increase transcripts associated with nerve activation, and these would be centered in neural and neuroendocrine-like cells.

Standard RNA sequencing, single nucleus-RNA sequencing, and spatial RNA sequencing of vagal ganglia. Standard RNA-sequencing and spatial RNA-sequencing of lungs in naïve and mice that have undergone asthmatic induction with Alternaria alternata.

Bulk RNA-seq revealed that genes related to allergen sensing were increased in asthmatic ganglia nodose/jugular ganglia compared to control ganglia. These genes were associated with nodose clusters as shown by single-nucleus RNA sequencing, and a distinct caudal-to-rostral spatial arrangement was presented as delineated by spatial transcriptomics. The distinct clusters closely match previous identification of nodose neuron clusters. Correspondingly, the lung transcriptome was altered with asthmatic induction such that transcripts associated with neural excitation were upregulated. The spatial distribution of these transcripts was revealed by spatial transcriptomics to illustrate that these were expressed in neuroendocrine-like cells/club cells, and neurons.

These results show that the neuroimmune transcriptome is altered in response to asthmatic induction in a cell cluster and spatially distinct manner.

Figures

Alternaria alternata alters the neuroimmune transcriptome. (A) Model of A alternata asthmatic induction (14, 22, 23). Image created with BioRender and used with a CC-BY-NC-ND license. (B) Serum levels of IgE were elevated after induction of asthma with A alternata (two-sided t-test). (C) Goblet cell metaplasia was induced with A alternata (two-sided t-test). N = 4 per group. (D) Representative images, ×20, scale bar = 20 µm. (E) PCA shows distinct clustering patterns based on allergic induction in lungs. (F) Volcano plot displaying differentially expressed genes from lungs (x-axis represents log2FC in gene expression, and y-axis represents −log10 adjusted p-value). (G) Heatmap of differentially expressed lung transcripts between control and A alternata groups. (H) PCA showing difference between gender and allergic induction in vagus. (I) Volcano plot displaying differentially expressed vagus genes (x-axis represents log2FC in gene expression, and y-axis represents −log10 adjusted p-value). (J) Heatmap of differentially expressed vagus transcripts between control and A alternata groups. PCA, principal component analysis.

Single-nucleus RNA-seq analysis shows Alternaria alternata produces increased sensory transcripts in the vagal ganglia. (A) Unbiased clustering reveals 11 distinct clusters in vagal ganglia in control and asthmatic ganglia. (B) UMAP projection of vagal ganglion cells from asthmatic and control conditions. (C) General cellular annotation using canonical gene markers. (D) Specific cellular annotation using RCTD method with Kupari et al. dataset as the reference. (E) Canonical gene markers used to identify general cell types. (F) Prdm12 and Phox2b delineate jugular and nodose cell subsets. (G) Key nodose ganglion gene markers identified by Kupari et al. (18) are prominently shown in nodose clusters. (H) Key jugular ganglion gene markers identified by Kupari et al. (18) are prominently shown in jugular clusters. UMAP, uniform manifold approximation and projection; RCTD, Robust Cell Type Deconvolution.

CellChat analysis. (A) Number and strength of interactions identified by CellChat. (B) Information flow of each signaling pathway between Alternaria alternata and control vagal neurons. (C) Significant outgoing signals from neuron-related pathways between A alternata and control vagal neurons. (D) Changes in signaling patterns among cell populations, with a focus on how A alternata alters communication dynamics. (E) Network of cell populations involved in NRG signaling, indicating which cells are primarily sending, receiving, mediating, and influencing after A alternata exposure. NRG, neuregulin pathway. (F) Comparison of NRG signaling pathway activity between A alternata exposed and control vagal neurons, emphasizing shifts in cellular communication patterns.

Spatial distribution of vagal transcriptome in the ganglia and lungs. (A) Cell types identified by single-nucleus sequencing, projected onto spatial map of vagus. (B) Location of jugular and nodose cell types in vagus. (C) Spatial lung data were deconvoluted using LungMAP single-cell dataset to identify all cell types and projected onto UMAP coordinates. (D) Lund spatial dataset deconvoluted with our vagal single-nucleus RNA sequencing and Kupari et al. (18) single-cell RNA sequencing datasets. (D) Spatial distribution of jugular and nodose clusters. (F) Enlarged area outlined in red from (C) showing spatial distribution of ciliated, club, and MANCs. (G) Enlarged H&E stain of lung section. (H) Enlarged area outlined in red from (E) showing the spatial distribution of nodose and jugular cells. UMAP, uniform manifold approximation and projection; MANCs, mesenchymal alveolar niche cells.

Full text

Introduction

Asthma is characterized by variable airflow obstruction where exaggerated airway narrowing occurs due to aberrant inflammation in response to innocuous substances (1). Thus, inflammation is central to the pathogenesis of asthma (2). The surveillance of the lung by vagal afferent nerves occurs to protect the lung from foreign substances in the form of chemosensitive and mechanosensitive reflexes (3–5). Vagal sensory neuron hypersensitization is known to induce quick airway reflexes such as cough and bronchoconstriction (3, 4) in response to their stimulation, which is both dependent on and independent of immune cell-mediated stimulation (6). Vagal sensory afferent nerves also play an important role in driving asthmatic airway hyperresponsiveness and inflammation (4). Specifically, sensory neurons are sensitized by cytokines (7, 8), pathogens (9, 10), and allergens (11, 12), which release neuropeptides that activate T cells to induce inflammation (8, 13), B cells to release immunoglobulins (14), and mast cells to degranulate (15, 16) or traffic eosinophils (8, 13). Therefore, the interplay between sensory neurons, stromal cells, and immune cells dictates asthmatic pathogenesis.

Neuron subtypes can be classified by genomic features that directly reflect both ontogeny and function, with a taxonomy containing four principal types, namely cold, mechano-heat, A-low threshold mechanoreceptors, and mechano-heat-itch and C-low threshold mechanoreceptors, where each type branches into several neuron subtypes with unique response properties (17). This heterogeneity among neuron types is the reason for the cellular basis for discrimination among somatic sensory modalities and forms the basis for the ability to perceive, explore, and interpret the internal organ function (18, 19). Recent transcriptomic analyses have allowed the identification of distinct neural subsets involved in interoception, chemosensation, mechanosensation, and a host of homeostatic processes (18, 19) including single-cell RNA-seq analysis, which demonstrates molecularly unique subtypes of vagal sensory neurons (3, 18). Furthermore, molecularly distinct subsets of vagal sensory neurons have specific nerve terminals depending on their sensory modalities and functions (19, 20), and their functional properties are largely determined by their molecular identity (3, 20). This provides important information to predict possible phenotypic changes to vagal immune-sensing ganglia in response to disease. To date, however, few data are available to elucidate how vagal ganglia and respective target organs are transcriptionally altered in response to disease states. This information would be valuable for deterministic investigation of the mechanistic underpinnings of asthmatic inflammation.

In this study, we were motivated to investigate the simultaneous phenotypic changes induced by asthmatic sensitization, with Alternaria alternata, to the vagal ganglia and lungs. The A. alternata fungus significantly contributes to atopic human allergy and asthma (21), and in mouse models, it reproducibly induces heightened immunoglobulin E (IgE), eosinophilia, and mucus production (14, 22–24). Driven by previous investigations that provided distinct delineation of vagal sensory neuron subsets, we wanted to provide a complete transcriptomic profile of vagal sensory neurons in response to asthmatic induction along with lung-associated transcriptomic changes. We show that, globally, the lungs appear to upregulate neural transcription factors, whereas the vagal sensory afferent upregulates immune transcription factors. We used “bulk” RNA sequencing of the lungs and vagal ganglia, single-nucleus RNA sequencing of paired vagal ganglia, and spatial transcriptomic profiling of asthmatic lungs and vagal ganglia to showcase how asthmatic induction alters the neuro-immune system of the lungs.

Methods

Mice

All animal experiments were approved by The Lundquist Institute at Harbor UCLA Institutional Animal Care and Use Committee protocol #32183. Mice were housed in a specific pathogen-free animal facility at The Lundquist Institute. C57BL/6J mice were purchased from the Jackson Laboratory (Bar Harbor, ME, USA). At the start of the experiments, the mice were 5–8 weeks old. Age-matched male and female mice were used for experiments. In accordance with our animal protocol, all mice were anesthetized with 5% isoflurane and exsanguinated, followed by immediate dissection of tissues.

A. alternata asthmatic induction

A. alternata extract was purchased from CiteQ Biologics (Groningen, The Netherlands) (09.01.26) and dissolved to yield a final concentration of 25 µg/mL in phosphate-buffered saline (PBS). Mice were inoculated intranasally with 50 µL on days 0, 3, 6, and 9, where experiments took place 16 hours later (22). All control mice received PBS in the same volume intranasally.

Immunoglobulin analysis

IgE levels in serum were measured using enzyme-linked immunosorbent assay (ELISA; 501128838; Thermo Fisher Scientific, Waltham, MA, USA) according to the manufacturer’s instructions and analyzed using a BioTek Synergy H1 plate analyzer (Agilent Technologies, Santa Clara, CA, USA). Samples were compared using a two-sided unpaired t-test.

Histology

The lungs or vagal ganglia were embedded in paraffin and cut on a microtome (Leica, Wetzlar, Germany) at 5 µm and stained with hematoxylin (Epredia 6765001) and eosin (Epredia 6766007 or periodic acid–Schiff reagent) (Sigma395B-1kit). Histological slides were mounted with Permount toluene solution. The stained area was assessed using ImageJ; cell counts were normalized to tissue length. Samples were compared using a two-sided unpaired t-test.

Bulk RNA sequencing

Vagal (nodose/jugular) ganglia were dissected from A. alternata and PBS mice. Ten vagal ganglia from five mice were pooled together per sample. A total of n = 3 samples per group, per sex, were obtained. Three sets of the lungs (n = 3) were used per group. RNA was isolated using the inTron Easy spin Total RNA extraction Kit (Boca Scientific 17221, Dedham, MA, USA). RNA purity was assessed as >2.00, 260/280 ratio with spectrophotometry. RNA integrity was assessed prior to sequencing using the Agilent Bioanalyzer, and samples with RNA integrity number (RIN) >7 were used for library construction. RNA sequencing was performed by Novogene Corporation Inc. (Sacramento, CA, USA). mRNA was purified from total RNA using poly-T oligo-attached magnetic beads. To generate the cDNA library, the first cDNA strand was synthesized using a random hexamer primer and M-MuLV Reverse Transcriptase (RNase H−). Second-strand cDNA synthesis was subsequently performed using DNA Polymerase I and RNase H. Double-stranded cDNA was purified using AMPure XP beads, and the remaining overhangs of the purified double-stranded cDNA were converted into blunt ends via exonuclease/polymerase. After 3′ end adenylation, a NEBNext Adaptor with a hairpin loop structure was ligated to prepare for hybridization. To select cDNA fragments of 150–200 bp in length, the library fragments were purified using the AMPure XP system (Beckman Coulter, Beverly, MA, USA). Finally, PCR amplification was performed, and PCR products were purified using AMPure XP beads. The samples were sequenced on an Illumina NovaSeq 6000 with ≥20 million read pairs per sample. A full RNA-seq analysis pipeline was conducted by The Novogene Corporation, which included alignment with HISAT2 (25), principal component analysis (PCA) on the gene expression values [Fragments Per Kilobase of transcript per Million mapped reads (FPKM)], differential gene expression (DEG) analysis using DESeq2 (26), and finally a functional analysis including Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) ontology.

Single-nucleus RNA sequencing

Vagal ganglia were dissected from A. alternata and PBS mice. Ganglia were immediately snap-frozen in -80°C. Eight-to-ten ganglia (from four to five mice) were pooled together per group; n = 2 samples per group. Nuclei were isolated with the Minute™ single-nucleus isolation kit (Invent Biotechnologies Inc., Plymouth, MN, USA; #BN-020). The isolated nuclei were purified, centrifuged, resuspended in PBS with RNase Inhibitor, and diluted to 700 nuclei/µL for standardized 10x capture and library preparation protocol using 10x Genomics Chromium Next GEM 3′ Single Cell Reagent kits v3.1 (10x Genomics, Pleasanton, CA, USA). Libraries were sequenced using an Illumina NovaSeq 6000 (Illumina, San Diego, CA, USA). The libraries were sequenced with ~150 million PE50 reads per sample on Illumina NovaSeq. The raw sequencing reads were analyzed with the mouse reference genome (mm10) using Cell Ranger v7.1.0. To further clean the data, CellBender was used to remove background noise (27). A rudimentary cell quality control was employed after CellBender by removing cells using percentile-based thresholds on unique molecular identifier (UMI) count, gene count, and mitochondrial read fraction. Uniform manifold approximation and projections (UMAPs) were created after 1) finding highly variable genes using the seurat_v3 algorithm, 2) normalizing counts per cell, 3) log scaling counts, 4) scaling counts of 2,000 highly variable genes, and 5) performing principal component analysis on those scaled values for the highly variable genes. A nearest-neighbor graph was constructed with 20 neighbors based on cosine distance in principal component space (top 25 principal components).

Spatial RNA sequencing

The lungs were dissected, flushed with PBS, perfused with 4% paraformaldehyde with 20cmH_2_O pressure for 10 minutes, and then further fixed at 4°C for 24 hours. Vagal ganglia were fixed in 4% paraformaldehyde for 2 hours. Tissues were processed, and 5-μm paraffin slices were placed on the Xenium slide. DNA probes from a custom-designed panel (Supplementary Table S4), informed by bulk RNA-seq and single-nucleus RNA-seq, were used for spatial sequencing. The selected gene sets were hybridized through a process involving ligation hybridization and amplification, followed by rolling circle amplification to generate fluorescently labeled targets. This process provides a robust spatial resolution and minimizes false positives, as ligation requires complete probe binding for signal generations. A series of iterative hybridization, high-resolution imaging, and probe removal cycles were conducted to visualize all targeted probes. Spatial sequencing and initial imaging-based data acquisition were conducted using the Xenium imaging analyzer. Downstream data analysis and visualization of spatial sequencing were completed with Xenium Ranger (v1.7.1) using the default settings and Seurat R package (v5.10).

Bioinformatics analysis

The R-package Seurat (version 5.0.3) was used for most of the single-nucleus data analysis. The four independent biological replicates generated a total of 7,670 sequenced cells with an average of approximately 117,958 reads, 613 genes, and 910 UMIs detected per nucleus. All genes expressed in less than three nuclei were removed, and the filtered gene barcode matrices from the individual runs were merged. The cell–gene matrix was further filtered for mitochondrial DNA-derived gene expression (20% was set as the high cut-off), the number of genes detected per cell (>200 as the low cut-off and <5,000 as the high cut-off), and the number of UMIs detected per gene (>400 as the low cut-off and <10,000 as the high cut-off), removing 1,322 nuclei. Counts were normalized using the SCTransform function, which applies a regularized negative binomial regression model to account for technical noise and nucleus-specific biases. This normalization step adjusts for differences in sequencing depth and captures biological variability across cells. Normalized data were used for PCA to reduce the dimensionality of the dataset. UMAP and t-distributed stochastic neighbor embedding (t-SNE) were applied to visualize the cells in a low-dimensional space, capturing similarities between cells based on their gene expression profiles. The clustering of nuclei was performed using the original Louvain algorithm on the PCA-reduced data, defining clusters based on shared transcriptional profiles. Differential expression analysis was conducted to identify genes that were differentially expressed between nuclei clusters using the standard Wilcoxon rank sum test. Genes showing significant differential expression (adjusted p-value < 0.05) and high fold change (logFC > 0.25 or <-0.25) were considered cluster-specific markers and used to characterize distinct cellular populations based on their biological functions and expression profiles. Clusters were annotated based on known marker genes and enriched biological functions using Gene Ontology analysis. Cell type identities were assigned by comparing marker genes with reference datasets (18). Clusters and gene expression patterns were visualized using feature plots, violin plots, and dot plots to illustrate differential expression and cluster-specific marker genes.

Correlation analysis of bulk and single nuclei

For our analysis, we utilized Transcript per Million (TPM) values from bulk RNA sequencing data and Counts Per Million (CPM) values from single-cell RNA sequencing data. To assess the relationship between these datasets, we calculated the correlation between TPM values from the bulk RNA sequencing and Cluster PM values from the single-nucleus RNA sequencing. This approach allowed us to evaluate the consistency and coherence of gene expression profiles across different types of RNA sequencing data. To evaluate the correlation between our dataset and the previously published Kupari dataset, we first calculated the CPM for each dataset separately, followed by a log2 transformation of the CPM values. We then performed a correlation analysis on the log2-transformed CPM values.

CellChat

We used CellChatv2.1.1 to elucidate the intercellular communication networks in our single-nucleus RNA-seq dataset (28). We employed the full CellChat database, excluding the “Non-protein Signaling” category, for our cell–cell communication analysis. This tool utilizes known ligand–receptor interactions to infer communication patterns between different cell types. We created CellChat objects from our Seurat objects using the createCellChat function. We then followed the standard CellChat pipeline, which involved computing the communication probability using the computeCommunProb and computeCommunProbPathway functions. Subsequently, we calculated network centrality scores and the contribution of known ligand–receptor pairs using the netAnalysis_computeCentrality and netAnalysis_contribution functions, respectively. We visualized the cell communication results using the standard functions provided by the CellChat package.

Deconvolution of Xenium lung and neuron dataset

To deconvolute our Xenium data from the lung and ganglia, we utilized the Robust Cell Type Deconvolution (RCTD) method inside of the spacexr package (v2.2.1) in conjunction with appropriate reference datasets (29). For the lung Xenium data, we employed a dataset on LungMap (30, 31), which provided a detailed reference for lung cell type-specific expression profiles. We integrated the LungMap dataset with the lung Xenium data using the create.RCTD function with the LungMap dataset as the reference and our lung Xenium data as the query. RCTD’s deconvolution algorithm was then applied using the run.RCTD function to estimate the proportion of each lung cell type in the Xenium samples, allowing us to accurately infer cell type compositions. Subsequently, the ImageDimPlot function was used with the group.by parameter set to “predicted.celltype” to visualize the predicted cell types overlayed on our lung Xenium section.

For the vagal ganglia Xenium data, we integrated our single-nucleus RNA sequencing dataset with the published Kupari et al. (18) dataset. This integrated dataset served as a reference for cell type-specific expression profiles relevant to ganglia. Using this combined reference, we applied RCTD to the ganglia Xenium data to deconvolute the cell type proportions. The deconvolution process involved comparing the expression patterns observed in the Xenium data with the reference matrix to identify the cellular composition. To ensure the robustness of our findings, we performed cross-validation and sensitivity analyses to validate the stability and reliability of the deconvolution results for both tissue types. Additionally, we used our ganglia single-nucleus dataset integrated with the Kupari et al. (18) dataset as a reference for the RCTD algorithm to see if we could identify ganglion cells innervating the lung.

Results

Confirmation of asthmatic induction

Mice subjected to A. alternata sensitization displayed enhanced goblet cell metaplasia and heightened IgE concentrations in serum (Figures 1A–D) compared to PBS controls, demonstrating the induction of an allergic asthmatic phenotype consistent with our previous results, which demonstrated increased IgE, eosinophilic infiltration, and airway hyperresponsiveness to methacholine challenge (14, 22–24).

Alternaria alternata alters the neuroimmune transcriptome. (A) Model of A alternata asthmatic induction (14, 22, 23). Image created with BioRender and used with a CC-BY-NC-ND license. (B) Serum levels of IgE were elevated after induction of asthma with A alternata (two-sided t-test). (C) Goblet cell metaplasia was induced with A alternata (two-sided t-test). N = 4 per group. (D) Representative images, ×20, scale bar = 20 µm. (E) PCA shows distinct clustering patterns based on allergic induction in lungs. (F) Volcano plot displaying differentially expressed genes from lungs (x-axis represents log2FC in gene expression, and y-axis represents −log10 adjusted p-value). (G) Heatmap of differentially expressed lung transcripts between control and A alternata groups. (H) PCA showing difference between gender and allergic induction in vagus. (I) Volcano plot displaying differentially expressed vagus genes (x-axis represents log2FC in gene expression, and y-axis represents −log10 adjusted p-value). (J) Heatmap of differentially expressed vagus transcripts between control and A alternata groups. PCA, principal component analysis.

Lung transcriptome

Bulk sequencing of whole lungs showed an average of 54,922,129 reads of 29,897 genes across samples. Principal component analysis revealed that 56% and 21% were sufficient to describe the differences between samples (Figure 1E). We then performed a two-group comparison of each ganglion group to the other two groups and identified the genes that were enriched or depleted, DESeq2, p ≤ 0.05. Applying these criteria, we identified 1,783 unique genes that were differentially regulated in A. alternata mice compared to PBS, which included a downregulation of 433 genes and an upregulation of 1,350 genes (Figure 1F). Gene expression differences are displayed as heatmaps in Figure 1G. A full list of these genes with their transcript LogF2C numbers is shown in Supplementary Table S1.

Next, we used GO and KEGG orthology analysis to identify likely gene sets, which were differentially expressed in A. alternata lungs compared to PBS control lungs. GO analysis showed a differential upregulation of eosinophil migration (GO:0072677, p_adj_ = 0.000337), eosinophil chemotaxis (GO:0048245, p_adj_ = 0.000458), and mucosal immune response (GO:0002385, p_adj_ = 0.00152), which demonstrated a robust asthmatic phenotype. Interestingly, upregulation of response to pain (GO:0048265, p_adj_ = 0.001067) and behavioral response to pain (GO:0048266, p_adj_ = 0.009631) showed an augmented upregulation of sensory neuron genes (Supplementary Figure S1A). KEGG orthology revealed a differential regulation of neuroactive ligand–receptor interaction (mmu04080, p_adj_ = 6.08 × 10^−7^) along with calcium signaling pathway (mmu04020, p_adj_ = 0.002555), asthma (mmu05310, p_adj_ = 0.000791), and cytokine–cytokine receptor interaction (mmu04060, p_adj_ = 0.0000359; Supplementary Figure S1B).

Vagal transcriptome

Bulk sequencing of whole vagal ganglia showed an average of 55,046,103 reads of 33,423 genes across all samples. The two principal components (sex vs. disease state) revealed that 50% and 12% were sufficient to describe the differences between samples (Figure 1H). We then performed a two-group comparison of each ganglion group to the other two groups and identified the genes that were enriched or depleted. We filtered transcripts using the p-value ≤ 0.001. Applying these criteria, we identified 2,347 unique differentially regulated genes in A. alternata mice compared to PBS with an upregulation of 726 and a downregulation of 1621 genes compared to PBS mice (Figure 1I). Gene expression differences are displayed as heatmaps in Figure 1J.

Next, we used Gene Ontology and KEGG orthology to identify likely gene sets that were differentially expressed in A. alternata ganglia compared to PBS control ganglia. GO analysis revealed a single significant upregulation of G protein-coupled receptor (GPCR) activity (GO:0004930, p_adj_ = 0.043113), ion transmembrane transporter activity (GO:0015075, p_adj_ = 0.048606), and ion transport (GO:0006811, p_adj_ = 0.048224), which were upregulated (Supplementary Figure S1C). Specific GO analysis revealed that transcripts associated with immunogenic processes were upregulated in the vagal ganglia (Supplementary Figure S1C). KEGG orthology revealed a single significant upregulation of the neuroactive ligand–receptor (mmu04080, p_adj_ = 0.009449; Supplementary Figure S1D). A full list of these genes with their transcript LogF2C numbers is shown in Supplementary Data S2.

Single-nucleus RNA sequencing of vagal ganglia

A total of 2,640 nuclei from A. alternata and 5,030 nuclei from PBS-treated mouse vagal ganglia were initially analyzed. After quality control to remove low-quality nuclei and potential doublet nuclei, 1,627 from A. alternata and 3,416 nuclei from PBS-treated mice passed the quality control (QC) requirements. After feature selection, removal of confounding sources of variation, and principal component analysis, the nuclei were clustered using the Seurat Louvain-based algorithm. The first round of analysis at a resolution of 0.5 produced 13 clusters, and after merging highly similar groups, 11 clearly distinctive clusters remained (Figure 2A, Supplementary Figure S2). All neuron clusters included nuclei from each biological replicate, demonstrating that the transcriptional profiles were robust and reproducible over the >5,000 number of nuclei collected between two separate experiments and groups (Figures 2B–D).

Single-nucleus RNA-seq analysis shows Alternaria alternata produces increased sensory transcripts in the vagal ganglia. (A) Unbiased clustering reveals 11 distinct clusters in vagal ganglia in control and asthmatic ganglia. (B) UMAP projection of vagal ganglion cells from asthmatic and control conditions. (C) General cellular annotation using canonical gene markers. (D) Specific cellular annotation using RCTD method with Kupari et al. dataset as the reference. (E) Canonical gene markers used to identify general cell types. (F) Prdm12 and Phox2b delineate jugular and nodose cell subsets. (G) Key nodose ganglion gene markers identified by Kupari et al. (18) are prominently shown in nodose clusters. (H) Key jugular ganglion gene markers identified by Kupari et al. (18) are prominently shown in jugular clusters. UMAP, uniform manifold approximation and projection; RCTD, Robust Cell Type Deconvolution.

Of the 11 clusters, five clusters expressed neuronal marker genes, including tubulin Beta 3 (Tubb3), synaptosome-associated protein 25 (Snap25), ELAV-like binding protein 4 (Elavl4), and RNA binging fox-1 homolog 3 (Rbfox3). These clusters also expressed either the nodose marker paired-like homeobox 2b (Phox2b) or the jugular marker PR domain zinc finger protein 12 (Prdm12, Figures 2E, F), which are established neuronal transcripts to separate jugular from nodose neurons (18, 32). Importantly, *Phox2b^+^

• and *Prdrm12^+^
• neurons were distinctly separated and contained unique gene sets (Figures 2E–H), which closely corroborated previous data (Supplementary Figure S3, 4D). For instance, cluster 4, which showed higher expression of Prdm12, also showed higher expression of paired-related homeobox protein-like 1 (Prrxl1), SAM Domain, SH3 Domain And Nuclear Localization Signals 1 (Samsn1), and synaptotagmin 7 (Syt7). Clusters 1 and 6, which are also *Phox2b^+^ *, showed expression of other nodose marker genes such as homeobox B3 (Hoxb3), neurotrophic tyrosine kinase, receptor, type 2 (Ntrk2), purinergic receptor P2X, ligand-gated ion channel, 2 (P2rx2), formin homology 2 domain containing 3 (Fhod3), EYA transcriptional coactivator and phosphatase 1 (Eya1), and neuronal growth regulator 1 (Negr1). Clusters 0 and 2 have many nodose neuron markers but also showed specific expression of homeobox 5 (Hoxb5). Cluster 8 showed higher expression of genes marking Schwann cells such as peripheral myelin protein 22 (Pmp22), sex-determining region Y-box 10 (Sox10), periaxin (Prx), myelin protein zero (Mpz), and myelin-associated protein (Ncmap). Cluster 7 showed higher expression of glial cell marker genes: aspartoacylase (Aspa), which plays a role in the metabolism of myelin, in addition to laminin alpha 2 (Lama2), cadherin 19 (Cdh19), and ATP-binding cassette A8a (Abca8a, Figures 2G, H), indicating that these may be satellite glial cells (SGCs). Cluster 3 may represent a precursor or transitional state leading to cluster 7 (SGCs), as they share many marker genes, but cluster 3 showed expression of marker genes that are typically associated with cell differentiation such as adhesion G protein-coupled receptor B3 (Adgrb3) and lysophosphatidic acid receptor 1 (Lpar1). Additionally, clusters 5 and 9 expressed endothelial cell (EC) marker genes: endomucin (Emcn), platelet and endothelial cell adhesion molecule 1 (Pecam1), and friend leukemia integration 1 (Fli1). Cluster 5 expressed genes involved in angiogenesis and vascular integrity, such as solute carrier organic anion transporter family, member 1a4 (Slco1a4), ATP binding cassette subfamily G member 2 (Abcg2), lymphoid enhancer-binding factor 1 (Lef1), friend leukemia integration 1 (Fli1), and FMS-like tyrosine kinase 1 (Flt1), which suggest that these may be nuclei from vascular endothelial cells. Cluster 9 may represent a more specialized subset of metabolically active endothelial cells, potentially involved in lipid metabolism or related to the vasculature of adipose tissues due to the higher expression of genes such as platelet glycoprotein 4 (Cd36), fatty acid binding protein 4 (Fabp4), and peroxisome proliferator-activated receptor gamma (Pparg). To further annotate the nucleus types, we utilized a highly detailed single-cell RNA-seq dataset published by Kupari et al. (18), which resulted in the annotation shown in Figure 2D.

A total of 222 genes identified as differentially expressed in the bulk RNA-seq data were also found to be differentially expressed between control and asthma conditions in the nodose and jugular of the single-nucleus RNA-seq data. Many of these genes have previously been linked to inflammation and may play roles in the response to asthma. The key differences between clusters are shown in Supplementary Table S3. For example, transforming growth factor–beta receptor II (Tgfbr2) and transforming growth factor–beta receptor III (Tgfbr3) were both downregulated in both datasets and were shown to be expressed in the jugular neurons in the single-nucleus dataset. Additionally, beta-2 microglobulin (B2m), Cd36, endothelial PAS domain protein 1 (HIF2a/Epas1), heparan sulfate proteoglycan 2 (Hspg2), integrin alpha 1 (Itga1), and adhesion G protein-coupled receptor F5 (Adgrf5) were downregulated in both datasets and specifically in jugular and nodose neurons, all of which are known to be involved in immune responses. Other genes previously shown to play roles in inflammation or asthma conditions were upregulated in both datasets, such as activity-dependent neuroprotective protein (Adnp), ArfGAP with RhoGAP domain, ankyrin repeat and PH domain 2 (Arap2), jun D proto-oncogene (Jund), matrix remodeling associated protein 7 (Mxra7), and solute carrier family 39 member 11 (Slc39a11).

Using CellChat, 1,260 and 1,416 interactions were identified between control and asthma, respectively, (Figure 3A), which are grouped into 56 pathways (Figure 3B). CellChat further revealed 156 additional interactions with a higher association that were grouped into 12 different pathways (Figures 3C, D). Those pathways include those with key roles in tissue repair, fibrosis, and airway remodeling, such as BMPa, FGF, CDH, and COLLAGEN, which are crucial processes in chronic asthma (33–36). Other interactions were linked to neuregulin (NRG) and SEMA3, which are involved in neuronal signaling, inflammation, and immune regulation, all of which are integral to asthma pathophysiology (37–40). TENASCIN, AGRN, and CADM contribute to tissue and synapse remodeling, likely affecting airway structure and neurogenic inflammation (41–44). Glutamate and DHEA are linked to neuronal signaling and immune modulation, influencing bronchoconstriction and inflammation in asthma (45–49). Lastly, ADGRG is involved in immune cell migration and activation, impacting the immune response during asthma hypersensitivity (50–52). Together, these pathways were shown to have more interactions in the asthma group, which contribute to the inflammation, immune response, and structural changes in asthma, particularly in response to allergens like A. alternata. Cell types with the most significant changes in sending or receiving signals after asthma induction were predominant in jugular and nodose cells (Figure 3D). One of the prominent CellChat pathways upregulated with asthma induction was NRG (Figures 3B, C). This pathway is increased, with nodose and jugular neurons being the primary senders to Schwann cells and glial cells as influencers (Figures 3E, F). This is consistent with previous literature that shows NRG signaling not only supports the maintenance and repair of neural structures but also modulates immune responses by influencing the differentiation and activation of immune cells, particularly through its interaction with Schwann cells (53). Differentially regulated ligand–receptor pairs between conditions, with jugular and nodose neurons as the source and other neurons are the target, are shown in Supplementary Figure S5.

CellChat analysis. (A) Number and strength of interactions identified by CellChat. (B) Information flow of each signaling pathway between Alternaria alternata and control vagal neurons. (C) Significant outgoing signals from neuron-related pathways between A alternata and control vagal neurons. (D) Changes in signaling patterns among cell populations, with a focus on how A alternata alters communication dynamics. (E) Network of cell populations involved in NRG signaling, indicating which cells are primarily sending, receiving, mediating, and influencing after A alternata exposure. NRG, neuregulin pathway. (F) Comparison of NRG signaling pathway activity between A alternata exposed and control vagal neurons, emphasizing shifts in cellular communication patterns.

Spatial sequencing

The selected 493 genes included 300 genes provided by the Xenium Mouse Tissue Atlassing panel and 92 custom genes selected from genes of interest from bulk sequencing and single-nucleus sequencing (Supplementary Table S4). Of the analyzed genes, we identified vagal cell clusters from those identified in single-nucleus sequencing (Figure 4A). The genes identified are consistent with upregulation in both snRNA-seq and bulk RNA-seq data where the agreement is (Supplementary Figure S4C). Then, we mapped the clusters identified in single-nucleus sequencing spatially to identify their tissue distribution, which confirms the spatial distribution and the clustering of cell types as presented previously (3, 18) (Figure 4B). These clusters are reminiscent of all clusters identified from snRNA-seq (Figure 4B).

Spatial distribution of vagal transcriptome in the ganglia and lungs. (A) Cell types identified by single-nucleus sequencing, projected onto spatial map of vagus. (B) Location of jugular and nodose cell types in vagus. (C) Spatial lung data were deconvoluted using LungMAP single-cell dataset to identify all cell types and projected onto UMAP coordinates. (D) Lund spatial dataset deconvoluted with our vagal single-nucleus RNA sequencing and Kupari et al. (18) single-cell RNA sequencing datasets. (D) Spatial distribution of jugular and nodose clusters. (F) Enlarged area outlined in red from (C) showing spatial distribution of ciliated, club, and MANCs. (G) Enlarged H&E stain of lung section. (H) Enlarged area outlined in red from (E) showing the spatial distribution of nodose and jugular cells. UMAP, uniform manifold approximation and projection; MANCs, mesenchymal alveolar niche cells.

We then used a published lung single-cell dataset to deconvolute the spatial lung dataset (Figure 4C). When we superimposed the clusters from single-nucleus vagal sequencing onto the lung spatial analysis, we identified that all jugular and nodose clusters were spatially distributed in the airways and specific alveolar regions (Figures 4C–E). Specifically, the jugular and nodose clusters were predominantly found in areas surrounding the airways, where they overlayed with ciliated/club cells (Figures 4E–H). The nodose clusters were highly prevalent in alveolar regions, overlapping with alveolar cell types, including mesenchymal alveolar niche cells (MANCs) (Figures 4E–H). The innervated cell types by the vagal clusters were identified as club cells and neuroendocrine cells. This is consistent with previous data characterizing these cells with immunohistochemical analyses (32, 54).

Discussion

Vagal ganglia are key multimodal sensors of the lungs. As a multitude of the neuron types housed in the paired nodose and jugular ganglia innervate the lung, as demonstrated both physiologically and transcriptomically (3, 18, 55), their susceptibility to changes in response to asthma should be profound. Indeed, single-cell RNA sequencing has identified that specific cytokines reprogram lung-innervating vagal neurons (56). Our study provides additional transcriptomic analysis and is one of the most comprehensive transcriptomic analyses of vagal ganglia in response to asthmatic induction and correspondingly identifies cell types associated with neural patterns in the lung. Additionally, our data show the changes between control and asthmatic vagal ganglia. Furthermore, our single-nucleus and spatial mapping identify vagal clusters that innervate the lung and provide a map to identify the spatial locations of these clusters in the lung. We also show changes in specific vagal clusters between asthma and control congruent with bulk RNA-seq. Our results allow for a number of functional predictions to be made but, importantly, also provide tools for direct experimental strategies to address their morphology, physiology, connectivity, and function in response to allergic sensitization.

Our transcriptomic analysis of the lungs and vagal ganglia, in tandem, further underscores the importance of vagal surveillance of the lungs. Bulk sequencing of the lungs revealed GO and KEGG pathway upregulation in A. alternata lungs compared to PBS controls, which are typically associated with a neural phenotype. Of note, Mrgprg (MAS-related GPR, member G, 6.99 Log2FC), Runx1 (runt-related transcription factor 1, 1.06 Log2FC), P2rx2 (purinergic receptor P2X, ligand-gated ion channel, 2, 1.99 Log2FC), Cacna1 (calcium channel, voltage-dependent, P/Q type, alpha 1A subunit, 2.48 Log2FC), and Ntrk1 (neurotrophic tyrosine kinase, receptor, type 1, 2.00 Log2FC) were highly and significantly upregulated in the lungs in response to A. alternata induction. That nerve-associated receptors, ion channels, and nerve injury markers were increased in A. alternata lungs shows that vagal hypersensitization is a significant component of the asthmatic genotype/phenotype. This is consistent with an increase in sensory neuron hypersensitization and likely their influence in innervating neuroendocrine cells and the epithelium. In particular, Mas-related G protein-coupled receptors are significantly tied to allergy where they contribute to asthmatic inflammation when stimulated either on sensory neurons (57) or on mast cells (58). In addition to this novel finding, the lungs also upregulate Tgfb, which is consistent with an allergic phenotype. Genes that were downregulated also showed how allergic sensitization affects the lung transcriptome, as GO:0031012 extracellular matrix was downregulated in the lungs of asthmatic mice, demonstrating a disruption of the epithelial barrier.

The upregulation of particular genes in the vagal ganglia as a result of A. alternata induction supports the importance of the neuroimmune sensing capabilities of the vagal ganglia. To establish these data, single-nucleus RNA sequencing confirmed bulk sequencing findings and showed how the clustering of specific genes occurs. Single-nucleus sequencing resulted in seven unbiased clusters, which were consistent with previous analyses by others (18). Interestingly, when jugular and nodose clusters are combined, A. alternata increases gene expression in jugular and nodose clusters, which is consistent with the predictions from Kupari et al. to be lung-innervating ganglia (18). The agreement between bulk and single-nucleus sequencing and differentially expressed genes between A. alternata and control mouse ganglia demonstrate the upregulation of immune features and specifically immune features associated with immune memory. The single-nucleus sequencing is consistent with spatial transcriptomics and shows the rostral–caudal distribution of these clusters consistent with a nodose/jugular distribution (32). Furthermore, our single-nucleus sequencing of vagal ganglia in response to A. alternata induction of asthma showed an upregulation of immune response genes such as Bcl6 (B-cell lymphoma 6), Il6 (interleukin-6), Itk (IL2-inducible T-cell kinase), and Cyp26a1 (cytochrome P450 26 subfamily A), which was similar to a recent single-cell RNA sequencing analysis of lung-innervating vagal neurons in response to ovalbumin induction of asthma (56). The upregulation of immune response genes in the vagal ganglia demonstrates the importance of vagal sensory neurons in the immune regulation of the lung.

The transcriptomic distribution within the lungs showed the diversity of vagal cell clusters within the lungs. The vagal cell clusters are consistent with the physiologic discernment of lung cells and the single-cell lung atlases publicly available (59, 60). Interestingly, when we superimpose the information from single-nucleus sequencing of the vagal ganglia onto the lung spatial transcriptomic map, we see the distinct distribution of these clusters within a finite space. The vagal cell clusters are distributed within the airways and distinct alveolar regions and were seen to “innervate” club cells/neuroendocrine-like cells. To our knowledge, this is one of the first transcriptomic analyses to corroborate the functional descriptions from previous investigations. The most robust changes in response to A. alternata induction occurred in both jugular and nodose cells. This is consistent with the upregulation of gene types in the single-nucleus and spatial vagal ganglion sequencing and the majority of changes occurring with immune transcripts.

Based on the cross-activation of gene sets within the vagal ganglia and lungs and the close spatial profiles of vagal cluster sets within the lungs, our data show how sensory neurons are transcriptomically affected by allergic sensitization. Importantly, CellChat analysis identified increased cell-to-cell communication, including increased NRG signaling, which may reflect enhanced communication between neurons and supporting cells within the ganglia and lungs, contributing to hypersensitization in the asthmatic phenotype. As our data show consistency to previously published gene sets in the vagus, the changes seen in response to asthmatic induction are confined to certain cell types that are specific to lung-innervating cell types (3, 18). As our spatial sequencing dataset only involved 470 genes, only 80 overlapped with the CellChat database. Of the 80 overlapping genes, only 17 were identified as highly variable ligand–receptor pairs suitable for use in the inference process. This limited gene coverage significantly restricted our ability to infer robust and meaningful cellular interactions, as the dataset does not capture the breadth of signaling pathways necessary for meaningful analysis.

Our data are among the first to corroborate the functional purpose of these cell clusters and give a genetic basis to potential sensitization in response to asthmatic induction. Our data provide a rich dataset to test the functional consequences of allergic sensitization of vagal cell subsets.