Single-Cell Transcriptomic Analysis Links Nonmyelinating Schwann Cells to Proinflammatory Response in the Lung

The peripheral nervous system (PNS) is in close contact with our environment, especially at barrier tissues such as the skin, gut, and lung (1). Schwann cells (SCs) are the principal glial cells of the PNS and have been shown to play a role in barrier defenses in the gut and skin (2, 3). The neurons at these sites are protected and supported by two types of SCs, myelinating SCs (mSCs) and nonmyelinating SCs (NMSCs) (4). mSCs wrap around large-caliber axons (≥1 µm in diameter) and produce myelin for axon insulation, whereas NMSCs ensheath small-diameter axons (<1 µm in diameter). In the lung, most neurons are nonmyelinated and are associated with NMSCs (5).

Studies in sciatic nerve injury, digit tip regeneration, and skin and gut injury animal models have shown that peripheral glia play important roles in tissue repair (2, 3, 6–8). The response of SCs (NMSCs and mSCs) to tissue injury, including injury to the PNS, relies on their unique ability to dedifferentiate into an immature SC state, subsequently expand to restore tissue homeostasis, and then redifferentiate to promote tissue regeneration and functional recovery (9). Upon injury of the PNS, axonal degeneration leads to the activation of SCs to produce chemokines, cytokines, and growth factors that lead to the recruitment of immune cells, promote the expansion of mesenchymal stem cells, and facilitate repair of injured axons (10–12).

In the lung, the crosstalk between sensory neurons and the immune system plays a critical role in the maintenance of barrier function and host defense; this neuroimmune interaction has been implicated in chronic inflammatory conditions of the airways such as chronic obstructive pulmonary disease and asthma (13). Despite their intimate association with neurons in the lung, the molecular characteristics and role of NMSCs in tissue homeostasis are not well understood. Thus, a systematic characterization of NMSCs both in vitro and in vivo is essential to begin to understand the functional relevance of this cell type in airway physiology and pathophysiology.

Lung NMSCs were first characterized by electron microscopy (5) and were subsequently studied by immunohistochemistry using markers such as glial fibrillary acidic protein (GFAP) (14). These imaging studies demonstrated a close association of NMSCs and axons of pulmonary nerves. In addition, lung NMSCs in mice can be genetically labeled by GFAP promoter–driven expression of GFP (14). Such GFAP reporter mice have been widely used to mark and characterize astrocytes in the CNS and SCs in the PNS (15–18).

In this study, we used GFAP reporter mice to determine the spatial association of SCs with neurons in lung tissue. Then, we isolated GFAP-labeled SCs, including both mSCs and NMSCs, from mouse lungs and characterized their transcriptome using single-cell RNA sequencing (scRNA-seq). Furthermore, we interrogated the LPS-induced proinflammatory response of lung SCs both in vitro and in vivo. Taken together, our results suggest that lung NMSCs, but not mSCs, may be actively involved in the innate immune response in the lung.

Human GFAP (hGFAP) promoter cre (15) was crossed with ROSA26.Loxp-Stop-Loxp.tdTomato mice to generate hGFAP-tdTomato reporter mice. Mice used in these studies were females, with an age range of 12–16 wk. For all animal takedowns, mice were euthanized by CO₂, and lungs were perfused with PBS/heparin (5 U/ml). For acute lung inflammation study, oropharyngeal aspiration (OA) of 50 μl of PBS for 2 h or 50 μl of LPS at 0.2 mg/kg for 2 h was used. All animal activities were performed as required by the Institutional Animal Care and Use Committee of Genentech, the Animal Welfare Act, and in accordance with the Guide for the Care and Use of Laboratory Animals.

SCs were sorted from hGFAP-tdTomato mice by selecting for CD45⁻EPCAM⁻CD31⁻tdTomato^hi cells, cultured for 8 d in DMEM plus 10% FBS (Invitrogen) and 1× penicillin/gentamicin (Invitrogen), and then stimulated with PBS or 20 ng/ml LPS for either 2, 4, or 16 h. Bone marrow–derived macrophages (BMDMs) were generated using male wild-type C57BL/6J mice (The Jackson Laboratory) by flushing out bone marrow from the tibia and femurs using 1× PBS with a 25G needle. After a 1500 rpm spin, the cells were resuspended in 4 ml of 1× PBS and layered over 4 ml of Ficoll. The cells were centrifuged at 1500 rpm, without break, for 15 min at room temperature. The cells at the interphase were collected and resuspended in 1× PBS up to 10 ml and centrifuged for 5 min at 1500 rpm. The cell pellet was resuspended in BMDM media, DMEM plus 10% FBS (Invitrogen), 1× penicillin/gentamicin (Invitrogen), and 50 ng/ml M-CSF (PeproTech) for 7 d, with fresh media introduced on days 3 and 5. On day 7, BMDMs were split into 12-well plates using BMDM media and stimulated on day 8 with PBS or 20 ng/ml LPS for 4 h.

hGFAP-tdTomato mice were used to sort out glial cells by briefly selecting for CD45⁻EPCAM⁻CD31⁻tdTomato^hi cells. Glial cells were cultured for 8 d in 10% FBS (Life Technologies), DMEM (Life Technologies), and penicillin/streptomycin (Life Technologies). Cells were stimulated with PBS or 20 ng/ml LPS for 2 or 4 h. After stimulations, RNA was harvested using the RNeasy micro purification column kit (Qiagen). RNA quantification and purity were analyzed with NanoDrop 2000 (Thermo Scientific).

The fastq sequence files for all RNA-seq samples were filtered for read quality and rRNA contamination. The remaining reads were then aligned to the mouse reference genome (GRCm38) using the GSNAP alignment tool (19). Alignments were produced using the following GSNAP parameters: “-M 2 n 10 -B 2 -i 1 N 1 w 200000 -E 1 –pairmax-rna = 200000 –clip-overlap”. These steps and the downstream processing of the resulting alignments to obtain read counts and normalized reads per kilobase transcript per million mapped reads per gene (over all exons of RefSeq gene models) were implemented in the Bioconductor package HTSeqGenie (v3.12.0). Only uniquely mapped reads were used for further analysis.

Lungs from hGFAP-tdTomato mice were used for all imaging studies. Briefly, animals were anesthetized using isoflurane and perfused with PBS/heparin (5 U/ml) followed by tissue fixation using 4% paraformaldehyde for lung section staining; after fixation, lungs were incubated in 30% sucrose in 1× PBS overnight, followed by dry ice freezing in OCT (Tissue-Tek). Cryosections of 25 μm were placed on glass slides and processed by incubating for 10 min in 4% paraformaldehyde, washed three times with PBS, permeabilized in 0.3% TX-100 in 1× PBS for 30 min, blocked in blocking buffer (6% donkey serum [Jackson ImmunoResearch] + 1% BSA [Invitrogen] + 0.2% TX-100 in 1× PBS) for 1 h, primary Ab in blocking buffer overnight at 4°C, and washed three times for 5 min, followed by secondary Ab (Jackson ImmunoResearch) incubation at room temperature for 1 h, washed three times for 5 min, stained with 5 μg/ml DAPI (Invitrogen) for 5 min, and washed in 1× PBS, concluding with coverslip mounting in ProLong Gold mounting media (Invitrogen). Primary Abs used were anti–red fluorescent protein (Rockland Immunochemicals), anti-Tuj1 (Abcam), anti–myelin basic protein (MBP; Abcam), and anti–ZO-1 (Thermo Fisher Scientific). Secondary donkey Abs were conjugated with Alexa Fluor 488, Alexa Fluor 546, and Alexa Fluor 647 (Jackson ImmunoResearch). Conjugation of primary Abs was performed using a fluorophore labeling kit (Invitrogen).

Briefly, animals were anesthetized using isoflurane and perfused with PBS/heparin (5 U/ml) followed by tissue fixation using paraformaldehyde 4% (19). Lung tissue was harvested and subjected to immunolabeling using primary conjugated Abs and a previously published staining protocol avoiding methanol treatment steps (20). Tissue clearing was performed using the FluoClearBABB approach (21), and whole-mount images were then acquired using a Leica SP8 microscope equipped with a white light laser and a Leica BABB immersion lens (HCX APO L ×20/numerical aperture 0.95 immersion media). To visualize hGFAP-tdTomato endogenous fluorescence, imaging was performed without the tissue clearing steps. Acquired data were visualized on a power workstation using Imaris (Bitplane).

Lungs from wild-type mice were perfused with 1× PBS, isolated, and rapidly frozen in OCT (Tissue-Tek) Cryomolds. Lung cryosections of 14 μm were used for the RNAscope assay. A RNAscope multiplex fluorescent kit v2 (Advanced Cell Diagnostics was used following the manual for fresh-frozen sample preparation and pretreatment. Briefly, for pretreatment conditions, sections were fixed for 45 min with 4% paraformaldehyde, dehydrated in 50, 75, and 100% ethanol, treated with hydrogen peroxide for 10 min, and protease 3 was used for 30 min; all other steps were followed using the protocol by Advanced Cell Diagnostics. Mouse probes used were as follows: Sox10 (catalog no. 435931-C2), Apod (catalog no. 580181-C3), positive three-plex control (catalog no. 320881), and negative three-plex control (catalog no. 320871), all probes from Advanced Cell Diagnostics. Fluorophore dyes used were TSA Vivid fluorophore 570 (Bio-Techne, catalog no. 323272) and TSA Vivid fluorophore 650 (Bio-Techne, catalog no. 323273); dyes were diluted 1:1000 in Advanced Cell Diagnostics probe dilution buffer. Imaging was performed on a Leica SP8 confocal microscope.

Upon CO₂ euthanasia, serum was collected and lungs were perfused with 1× PBS-EDTA via the right ventricle. Lungs were collected for either histology or tissue processing for single-cell dissociation. All lung samples per cohort were processed at the same time by transferring the lungs in C tubes (Miltenyi Biotec) and adding 5 ml of digestion media consisting of 0.05% Liberase, 0.1% DNase I (Roche), 1% BSA (Sigma-Aldrich), and 0.1% Dispase II (Sigma-Aldrich) in RPMI 1640. Samples were digested using a gentleMACS tissue dissociator (Miltenyi Biotec) using lung dissociation settings, and samples were incubated at 37°C for 30 min in a shaker at 120 rpm. After digestion and incubation, samples were further homogenized using the lung dissociation settings, and the supernatants were obtained and centrifuged at 1600 rpm for 5 min. The pelleted supernatants were treated with ACK (ammonium-chloride-potassium) lysing buffer, followed by resuspension in 1× PBS and passing through a 40-μm cell strainer. Samples were pelleted and resuspended in 1× PBS, stained with Live/Dead fixable dye (Invitrogen) at a 1:1000 dilution, incubated on ice for 15 min, washed in 1× PBS, resuspended in FACS buffer (PBS + 2.5 mM EDTA + 5% BSA) with FcR blocking reagent (Miltenyi Biotec), and stained with the appropriate conjugated fluorescent Abs. Abs used were as follows: anti–CD45-BV650 (BD Biosciences, catalog no. 563410), anti–CD31-BV421 (BD Biosciences, catalog no. 563356), and anti–EPCAM-BUV396 (BD Biosciences, catalog no. 740281). Samples were run and analyzed on a Symphony analyzer (BD Biosciences).

In vitro analytes from day 8 sorted hGFAP-tdTomato glial cells were stimulated with PBS or 20 ng/ml LPS for 16 h, and day 8 BMDMs and day 8 culture sorted hGFAP-tdTomato cells stimulated with PBS or 20 ng/ml LPS for 4 h were measured using a Luminex bead assay (Millipore).

The Chromium Next GEM (gel beads in emulsion) Single Cell 3′ GEM, Library & Gel Bead Kit v3.1 (10x Genomics PN-1000121) was used for library preparation according to the manufacturer’s guidelines. The cell-reverse transcription mix was prepared to aim for 10,000 cells per sample and applied to the Chromium Controller for GEM generation and barcoding. Then, samples were subjected to post–GEM-reverse transcription cleanup, cDNA amplification (11 cycles with v3.1), and library construction according to the user’s manual. Sample index PCR was done with 12 cycles. Libraries were then quantified by a Qubit dsDNA high-sensitivity assay kit (Thermo Fisher Scientific, catalog no. Q33230) and profiled by a Bioanalyzer high-sensitivity DNA kit (Agilent Technologies, catalog no. 5067-4626). Libraries were sequenced by HiSeq 4000 (Illumina) following the 10x Genomics sequencing specifications.

scRNA-seq data were processed with a Cell Ranger analysis pipeline. Briefly, reads were demultiplexed based on perfect matches to expected cell barcodes. Transcript reads were aligned to the appropriate species genome using GSNAP (2013-10-10) (22). Only uniquely mapping reads were considered for downstream analysis. Transcript counts for a given gene were based on the number of unique molecular identifiers (UMIs) for reads (up to one mismatch). Exonic reads were used to determine transcript count. Cell barcodes from empty droplets were filtered by requiring a minimum number of detected transcripts. Data quality for individual libraries was assessed based on total read depth, percentage of reads with valid barcodes, percentage of demultiplexed reads in detected cells, number of detected cells, and number of analyzed cells. Sample quality was further assessed based on the distribution of per-cell statistics, such as total number of reads, percentage of reads mapping uniquely to the reference genome, percentage of mapped reads overlapping exons, number of detected transcripts (UMIs), number of detected genes, and percentage of mitochondrial transcripts.

After this initial quality control, cells with <500 unique detected genes and >30% mitochondrial UMIs were discarded. After the filtering step, the gene × cell matrix of raw UMI counts was log normalized using NormalizeData() in Seurat v4 (23). All libraries were integrated using FindIntegrationAnchors() and IntegrateData() functions in Seurat v4. Then, we scaled the integrated data, performed dimensionality reduction by principal component analysis, and calculated uniform manifold approximation and projection (UMAP) coordinates and Louvain clustering for all cells using Seurat v4. For PBS- and LPS-treated samples, we implemented an additional filtering step to ensure a high-resolution lung glia dataset by removing all non–glial cell clusters. The resulting digital gene expression matrix was carried forward for clustering.

Mouse lung scRNA-seq data from Travaglini et al. (24) was downloaded from the National Institutes of Health’s Sequence Read Archive under BioProject accession number PRJNA632939 (https://www.ncbi.nlm.nih.gov/bioproject/?term=PRJNA632939) in FASTQ format and processed with the Cell Ranger analysis pipeline 5.0. Cell-type annotations were imported from Travaglini et al. (24), and the digital gene expression matrix with cell-type annotation was used for Seurat’s canonical correlation analysis with sorted GFAP-tdTomato cells from lung tissues.

To compare transcriptional programs induced by LPS stimulation on lung glial cells, we performed differential expression (DE) analysis between LPS- and PBS-treated samples using edgeR and limma R packages (25–28). For single-cell datasets, we generated pseudo-bulk expression profiles by summing all transcript counts of each cell-type cluster by each sample. We first constructed a design matrix that reflected the experimental design and filtered out low-abundance genes. Then, normalization factors were computed using the calcNormFactors() function from edgeR, and the pseudo-bulk count matrix was log-CPM (counts per million) normalized using the normalization factors. Then, we ran the voomWithQualityWeights() function from limma to model the mean-variance trend in the log-CPM across genes and then limma::lmFit() to fit linear models to the log-CPM for each gene. Next, we set up contrasts for comparing log fold changes between PBS and LPS groups and ran the limma::contrasts.fit() function to reorganize the fitted linear model so that contrast is represented by its own coefficients. Finally, limma::topTable() was used to retrieve summary statistic DE analysis. Threshold for significance at a false discovery rate of 5% and log fold changes of 2 were used for DE analysis on single-cell datasets, whereas a false discovery rate of 5% and log fold changes at 3 significance thresholds were used for bulk RNA-seq datasets.

Gene Ontology (GO) enrichment analysis was performed using the enrichGO() function from the clusterProfiler R package in which p values are calculated based on the hypergeometric distribution and corrected for testing of multiple biological process GO terms using the Benjamini–Hochberg procedure (29). GO terms were accessed using the AnnotationHub R package.

GraphPad Prism 6 was used for statistical analysis using the unpaired Student t test. Statistical details are provided in the figure legends.

RNA-seq data generated in this paper have been deposited in NCBI Gene Expression Omnibus under accession number GSM7511743 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSM7511743). This study does not report any original code. Any additional information required to reanalyze the data reported in this study is available from the corresponding authors upon request.

Lung SCs have been morphologically characterized both by immunohistochemistry and through inspection of cells expressing a fluorescent reporter under the control of the GFAP promoter (14). We first attempted to validate the location of SCs in the mouse lung using cross-sectional analysis of tissues as well as three-dimensional (3D) imaging. For these analyses, we used a well-studied GFAP transgenic mouse in which Cre recombinase is driven from the human GFAP promoter (16). This line was crossed to a ROSA26.LSL.tdTomato line to express tdTomato in GFAP-expressing cells (this line is referred to as GFAP-tdTomato). Lung cross-sections revealed that SCs were localized near the bronchi (Fig. 1A). Additionally, whole-mount sections costained with the neuronal marker TUJ1 demonstrated that SCs were localized along neuronal tracts (Fig. 1B). We further visualized SCs through high-resolution 3D imaging of a whole-lung major lobe from GFAP-tdTomato mice and compared their localization, patterning, and abundance to mSCs that are marked by MBP expression. The imaging from this study suggests that although SCs, including both NMSCs and mSCs, are colocalized with pulmonary nerves of large bronchi and organized in large neuronal bundles exhibiting elongated patterning (Fig. 1C, Supplemental Fig. 1), NMSCs appear to be a more abundant SC type ensheathing the pulmonary nerves than mSCs (Fig. 1C–E). Specifically, we observed that NMSCs are colocalized with single fibers of neurons that extend from large to small bronchi down to the alveoli where they manifest a set of gradually elongated cell processes (from an intricate mesh to a serpentine patterning) (Fig. 1C–E). In contrast, but consistent with previous findings (5, 14), our high-resolution 3D analysis showed that MBP-labeled mSCs are mainly associated with the larger-caliber pulmonary axons (Fig. 1C, Supplemental Fig. 1). Moreover, our 3D imaging revealed that NMSC processes extended into the alveoli space in close proximity to lung stromal and alveolar epithelial cells (Supplemental Video 1), suggesting that these cells might communicate with each other.

We then used flow cytometry to quantify SCs in GFAP-tdTomato mouse lungs. We gated out immune (anti-CD45), endothelial (anti-CD31), and epithelial (anti-EPCAM) markers and positively identified SCs by tdTomato fluorescence (Fig. 2A). Using this strategy, we were able to consistently isolate ∼3000–4000 SCs per mouse lung (Fig. 2A, 2B). The same strategy was also used to successfully sort and culture lung SCs (Fig. 2C).

We next characterized the molecular profile of CD45⁻CD31⁻EPCAM⁻tdTomato⁺ sorted lung SCs by annotating the cell populations resulting from our scRNA-seq analysis (Fig. 3A, Supplemental Fig. 2A). Because established lung single-cell atlases lack identifiable SCs, we integrated our sorted SC dataset with a previously published mouse lung atlas to ensure unbiased annotation of our data (24). Our integration analysis highlighted that the bulk of our sorted cells formed a distinct cluster that included both NMSCs and mSCs, among a total of 22 total cell-type clusters in the lung (Fig. 3B, 3C). Glia-specific markers (Plp1, Sox10, Gfap, Gpr37l1, and Mpz) were distinctly and differentially expressed in the SC (NMSC/mSC) cluster (Fig. 3D, 3E, Supplemental Fig. 2B). A subsequent analysis of the NMSC/mSC cluster identified three subclusters including two predominant populations of NMSCs (NMCS1 and NMSC2) and one minor population of mSCs (Fig. 3E). Differential gene expression analysis identified several markers specific to each one of the three populations: myelinating genes including Mbp, Mpz, Mag, and Pmp22 were expressed in mSCs; Slitrk6, Apod, and Gfra2 were expressed in NMSC1s; and Ntrk2 and Ntrk3 were expressed in NMSC2s (Fig. 3E, 3F). Furthermore, using RNAscope technology, we confirmed that NMSC1 (Sox10^hiApod^hi) can be differentiated based on Apod expression in the lungs, and most SCs (both NMSCs and mSCs) appear to be localized along the airways (Fig. 3G). GO analysis linked these three SC populations with distinct functions: mSCs expressed genes associated with myelination and axon ensheathment biological processes, whereas NMSC1s and NMSC2s expressed genes involved in axon development and nerve processes, as well as myeloid leukocyte recruitment and the response to bacterium and LPS (Supplemental Fig. 2C). Collectively, these data suggest that lung NMSCs may not only participate in conventional neuronal processes such as axon support and protection, but they also may serve as an active modulator of lung inflammation.

To explore the lung SC response to proinflammatory cues, we purified CD45⁻CD31⁻EPCAM⁻tdTomato⁺ cells from the lung and cultured them for 8 d to allow their expansion. We then stimulated them with PBS or LPS on day 8 for either 2 or 16 h and performed either RNA-seq or Luminex on the samples (Fig. 4A). RNA-seq analysis of cultured SCs validated the molecular identity of both mSCs and NMSCs through the identification of canonical markers (e.g., Plp1, Sox10, Gfap, Ncam1, Mbp, Mpz) known to be expressed in these cells (Supplemental Fig. 2D). The analysis of bulk RNA-seq showed that LPS robustly induces chemokine and cytokine gene expression in these lung SCs (Fig. 4B, 4C). Consistently, after a 16-h stimulation, we observed that lung SCs release elevated levels of chemokines and cytokines such as LIF, KC (CXCL1), MIP2, and MCP1 (Fig. 4D–G). Correspondingly, the transcript levels of these genes (Lif, Cxcl1, Cxcl2, Ccl2) were also increased after 2 h of LPS treatment (Fig. 4H–K). We further compared LPS-induced proinflammatory responses in lung SCs and BMDMs. Both Luminex and RNA-seq analyses showed that LPS-induced cytokine and chemokine expression levels were comparable between these two cell types, although we still noticed a cell type–specific response to LPS stimulation on a small number of genes (Supplemental Fig. 3). Collectively, our findings support a role of SCs in lung inflammation.

We next sought to determine the transcriptional response of lung SCs (mSCs and NMSCs) to LPS in vivo. We administered LPS to the lung via the OA route for 2 h to elicit a proinflammatory response in the lung, which can be confirmed by dramatic myeloid cell infiltration (Supplemental Fig. 4A). After OA administration of PBS or LPS for 2 h, lungs were immediately isolated and processed for the sorting of tdTomato⁺ cells for scRNA-seq analysis (Fig. 5A). Consistent with our previous results (Fig. 3B, 3E), the integrative analysis from this in vivo experiment verified the existence of two distinct NMSC populations (NMSC1 and NMSC2) and one mSC population in the lung (Fig. 5B). In addition, these three types of SCs were found in both LPS- and PBS-treated lungs (Fig. 5C, Supplemental Fig. 4B). Furthermore, we evaluated the LPS-induced transcriptional response in NMSCs and mSCs by performing differential gene expression analysis. The results highlighted that the NMSC1s exhibits a more robust expression of Tlr4 and a greater induction of proinflammatory cytokine and chemokine genes than do NMSC2s (Fig. 5D, 5E, 5G, Supplemental Fig. 4C). GO analyses also link NMSC subtypes to the innate immune response (Fig. 5F, 5H), and this observation appears to be associated with Tlr4 expression (Supplemental Fig. 4C). Contrary to our observations in NMSCs, the transcriptomic profiles of mSCs showed little change upon LPS challenge. Thus, these data suggest a potential role for NMSCs, but not mSCs, in lung inflammation.

In the current study, we used a cell type–specific reporter mouse (GFAP-tdTomato) to characterize the structure and spatial distribution of lung SCs (NMSCs and mSCs), and to profile transcriptional signatures unique to both cell types. We identified two distinct NMSC subpopulations that are molecularly differentiated from one another by the expression of Slitrk6 (NMSC1s) and Ntrk2 (NMSC2s). In unperturbed conditions, both NMSC1s and NMSC2s expressed genes that are required for maintenance of the PNS. After lung insults, NMSCs expressed genes associated with a proinflammatory response whereas mSCs exhibited no statistically significant transcriptional change. Taken together, our analyses support a notion that NMSCs may play a previously underestimated role in lung inflammation.

3D imaging confirmed that most lung SCs are NMSCs, and that NMSCs exhibit different patterning depending on their anatomical location. Staining of NMSCs revealed a pattern of cells that tightly colocalized with bundles and single fibers of neurons extending from large to small bronchi. Additionally, the staining revealed a pattern of elongated cell processes with a serpentine structure that extend to the alveoli. Glial cells from other visceral organs such as the gut have been shown to regulate intestinal barrier function (30). Thus, the intimate interaction with peripheral neurons together with the proximity of NMSCs to the alveoli could be indicative of a role of these cells in the front line of mucosal defense in the alveoli.

A role of glial cells in modulating activity of the immune response has been described in the gut, where enteric glial cells regulate gut inflammation and tissue repair after pathogen-induced intestinal damage (2, 31). In the lung, Tlr4 expression appears to be higher in NMSC1s than in NMSC2s. When stimulated with LPS in vivo, NMSC1s showed a greater proinflammatory response than did NMSC2s, as illustrated by elevated gene expression of chemokines and cytokines. This response was observed as early as 2 h after LPS treatment, at which point an initial recruitment of neutrophils and monocytes takes place. Preliminary data generated from our laboratory (data not shown) highlighted that the proinflammatory activity of NMSC1s is also maintained several hours after treatment, suggesting a potential function of these cells not only in promoting, but also in maintaining the innate immune response to pulmonary infections. Although enteric glial cells developmentally differ from NMSCs (32), our data provide initial evidence that supports the presence of common functional pathways in these two populations. The inflammatory responses in these cells may not depend on glial cell ontogenesis but are more likely dictated by their role in barrier tissues.

We found considerable differences at the level of spatial distribution, cell abundance, and gene expression between mSCs and NMSCs in the lung. mSCs show a uniform elongated serpentine structure that travels down the major bronchi along thick bundles of nerves. Our 3D imaging and single-cell analysis revealed that the mSC population represents the minority of lung SCs, corroborating a previously reported electron microscopy study (5). The top expressed genes in NMSCs were primarily involved in regulating neuronal and stromal homeostasis. These genes include Ntrk2, Ntrk3, and Col18a1, which are involved in neuron survival, proliferation, migration, and synapse formation and plasticity, and S1pr3, Matn4, and Wnt6, which are involved in the formation of extracellular matrices and angiogenesis. Further transcriptomic analysis showed that genes associated with myelination (Pmp22, Mbp, Mpz) were also expressed in most NMSCs, although the average expression was minimal as compared with the expression of these genes in mSCs. Thus, our findings suggest that in addition to regulating the activity of peripheral neurons, lung NMSCs may also play various nonneuronal functions and are therefore emerging as important local regulators of diverse pulmonary functions.

It would be important to assess and establish the functional role of NMSCs in the maintenance of lung homeostasis, airway inflammation and lung repair, and the impact that spatial localization plays on NMSC function. For example, single-cell ligand–receptor network analysis and RNAscope could help characterize NMSC interactions with both neuronal and nonneuronal cells. Additionally, coupling this analysis with loss-of-function studies would deepen our understanding of the function of lung NMSCs. In addition, multiple studies have shown that mSCs play a critical role in the repair of peripheral nerves through a cellular reprogramming process that generates cells that promote nerve regeneration. However, the role of NMSCs following neuronal damage has not been characterized yet. Thus, the use of conditional knockout mouse models could inform the function of genes expressed by NMSCs in the repair and regeneration of damaged peripheral axons, expanding our understanding of the roles of SCs (mSCs and NMSCs) in the development of peripheral neuronal disorders.

In summary, our work uncovered the spatial and molecular heterogeneity of lung NMSCs and provided an initial characterization of these cells as a new responder to proinflammatory cues. Further investigation will be required to elucidate the detailed functional role of NMSCs in pulmonary innate immune responses, which will improve our understanding of the pathogenesis of lung disorders and may inform therapeutic hypotheses targeting NMSCs in the lung.

All authors except A.M. are or were employees of Genentech/Roche. A.M. has no financial conflicts of interest.

We thank C.K. Poon, Alice Tan, and Terence Ho for FACS support, Laszlo Komuves for immunohistochemistry image analysis; Yuxin Liang, Qixin Bei, and Zora Modrusan for sequencing support; and the OMNI Bioinformatics group for helpful discussion on computational analysis. Figs. 3A, 4A, and 5A were created with BioRender.com.