A comparative study of in vitro air–liquid interface culture models of the human airway epithelium evaluating cellular heterogeneity and gene expression at single cell resolution
Respiratory Research volume 24, Article number: 213 (2023)
The airway epithelium is composed of diverse cell types with specialized functions that mediate homeostasis and protect against respiratory pathogens. Human airway epithelial (HAE) cultures at air–liquid interface are a physiologically relevant in vitro model of this heterogeneous tissue and have enabled numerous studies of airway disease. HAE cultures are classically derived from primary epithelial cells, the relatively limited passage capacity of which can limit experimental methods and study designs. BCi-NS1.1, a previously described and widely used basal cell line engineered to express hTERT, exhibits extended passage lifespan while retaining the capacity for differentiation to HAE. However, gene expression and innate immune function in BCi-NS1.1-derived versus primary-derived HAE cultures have not been fully characterized.
BCi-NS1.1-derived HAE cultures (n = 3 independent differentiations) and primary-derived HAE cultures (n = 3 distinct donors) were characterized by immunofluorescence and single cell RNA-Seq (scRNA-Seq). Innate immune functions were evaluated in response to interferon stimulation and to infection with viral and bacterial respiratory pathogens.
We confirm at high resolution that BCi-NS1.1- and primary-derived HAE cultures are largely similar in morphology, cell type composition, and overall gene expression patterns. While we observed cell-type specific expression differences of several interferon stimulated genes in BCi-NS1.1-derived HAE cultures, we did not observe significant differences in susceptibility to infection with influenza A virus and Staphylococcus aureus.
Taken together, our results further support BCi-NS1.1-derived HAE cultures as a valuable tool for the study of airway infectious disease.
The human airway epithelium is composed of diverse cell types which collectively function to maintain airway integrity, execute mucociliary clearance, and regulate airway immune responses. Component cell types include basal cells, the multipotent stem cells that serve as progenitors for other airway epithelial populations , secretory cells (encompassing both goblet cells and club cells), which secrete mucus and other anti-microbial and immunomodulatory peptides , and ciliated cells, which propel directional transport of mucus through coordinated ciliary activity . In addition to these abundant cell types, airway basal cell precursors also give rise to rare cell types including neuroendocrine cells, tuft-like cells, and ionocytes, among others. Neuroendocrine cells sense diverse chemical and mechanical stimuli and signal these changes to the central nervous system [3, 4]. Tuft-like (or brush) cells are a chemosensory cell population that have been linked to inflammatory processes through their production of interleukin-25 and leukotrienes [5, 6]. Ionocytes, present at approximately 1–2% of airway epithelial cell frequency, express high levels of ion transporters, including CFTR, the chloride transporter dysfunctional in patients with cystic fibrosis [5, 7].
Accounting for this complexity and cellular specialization is important for physiologically relevant study of the airway epithelium in the context of infectious disease. Many respiratory viruses, such as influenza A virus (IAV)  and SARS-CoV-2  preferentially infect ciliated cells. Human strains of IAV enter cells using α2,6-linked sialic acids (SA), present in a graded fashion along the upper airway, where they are predominantly expressed by ciliated cells and to a lesser extent, goblet (secretory) cells [10, 11]. In contrast, avian influenza strains (H5N1) utilize α2,3-linked SA which are expressed largely by alveolar type II cells in the lower respiratory tract . For SARS-CoV-2, ciliated cells express both the viral receptor ACE2 and protease TMPRSS2 required for fusion at the cell membrane, and motile cilia have been shown to be required for efficient replication in in vitro airway epithelial models [2, 9]. Airway mucus, produced by secretory cells, is also critical for innate defense against respiratory infections. Staphylococcus aureus (S. aureus), Pseudomonas aeruginosa (P. aeruginosa) and Haemophilus influenzae (H. influenzae) can all elicit increased productions of mucins, which are important for airway defense but can also be detrimental in patients with underlying respiratory disease [12,13,14]. Experimental model systems that account for this cellular and corresponding functional diversity are important for studies of host defense to, and disease caused by, respiratory pathogens.
Primary differentiated polarized human airway epithelium (HAE) cultures at air–liquid interface are an established, physiologically relevant in vitro system used to study mechanisms of airway homeostasis and disease [15,16,17,18]. HAE cultures are classically derived from primary human bronchial epithelial cells acquired via airway biopsy, from which precursor cells can be differentiated into pseudostratified epithelia when grown on permeable supports at air–liquid interface . HAE cultures recapitulate many features of the human airway epithelium in vivo, including characteristic airway cell types and secreted extracellular environment [19, 20]. As such, HAE cultures have been widely used in airway research [8, 14,15,16,17,18, 20,21,22]. An important limitation of HAE cultures derived from primary precursor cells is their restricted passage capacity before losing differentiation potential and/or becoming senescent . This complicates the design of larger experiments, as sufficient material may require multiple donor sources, thereby introducing confounding genetic variation and increasing cost. Although achieved in some studies [16, 21], the limited passage capacity of primary-derived HAE cultures also poses feasibility challenges to any potential genetic manipulation, which can require extended, multi-passage selection procedures.
BCi-NS1.1 is an airway basal cell line isolated from a healthy, non-smoking male donor which has been engineered to express hTERT, thereby extending the practical lifetime of these cells up to 40 passages without losing differentiation potential . In contrast, unmodified primary precursor cells can be passaged only 3–4 times before losing differentiation potential or entering senescence . Previous studies have demonstrated that BCi-NS1.1 cells effectively differentiate into HAE cultures with similar morphology and cell type composition to primary HAE cultures . The BCi-NS1.1 system has therefore been leveraged for a variety of applications in airway research [23,24,25,26]. However, the effects of hTERT expression on rare airway cell composition and associated gene expression programs in differentiated HAE cultures have been incompletely characterized, and any potential functional consequences for infection models with respiratory pathogens have not been directly explored. Furthermore, while BCi-NS1.1 HAE cultures should exhibit less culture-to-culture variation than HAE cultures derived from different donors, this has not been directly assessed.
Here, we perform an in-depth comparison of HAE cultures derived from BCi-NS1.1 cells (n = 3 independent differentiations) and HAE cultures derived from independent, commercially-available primary bronchial epithelial cells (n = 3 donors). Initially, we evaluate culture microarchitecture and canonical human airway epithelial cell type markers by histology and immunofluorescence labeling. We then use single cell RNA-Sequencing (scRNA-Seq) for comprehensive and unbiased assessment of constituent cell populations, differentiation states, and gene expression programs. To assess potential functional differences relevant to airway infectious disease, we evaluate HAE culture responses to infection with two common pathogens known to infect and/or colonize the airway, IAV and S. aureus [27, 28]. Lastly, as HAE culture systems are employed extensively in cystic fibrosis research, we evaluate the expression and function of CFTR channels in BCi-NS1.1-derived HAE cultures. Overall, we sought to thoroughly define the similarities and differences of hTERT-expressing BCi-NS1.1-derived HAE cultures and primary-derived HAE cultures to inform appropriate applications of these powerful in vitro experimental models.
BCi-NS1.1 cells (passage 17–25) (kind gift of Dr. Ronald Crystal) and Normal Human Bronchial Epithelial cells (Lonza, cat. no. CC-2541) were seeded in Pneumacult™-Ex Plus Medium (StemCell) and passaged at least two times before plating (7 × 104 cells/well) on rat-tail collagen type-I coated permeable Transwell membrane supports (Corning Inc; 6.5 mm diameter, 0.4 µm pore size) to generate HAE cultures. HAE cultures were maintained in Pneumacult™-Ex Plus Medium until confluent, then grown at air–liquid interface with Pneumacult™-ALI Medium in the basal chamber for approximately 3 weeks to form well-differentiated, polarized cultures.
Histology, immunofluorescence labeling and imaging
Embedding, sectioning and Hematoxylin and Eosin (H&E)/Periodic acid-Schiff (PAS) Alcian blue staining of HAE cultures were performed as described previously by our laboratory. For immunofluorescence labeling and imaging of HAE cultures, prepared slides were deparaffinized in Xylene (Crystalgen, cat. no. 301-038-4000) for 30 min, then rehydrated in a graded series of ethanol dilutions in water. Antigen unmasking was performed by heating slides submerged in 1X Citrate buffer (pH 6.0) (Abcam, cat. no. ab64214) in a 1.3 kW microwave for 1 min at 100% power, then for 8 min at 10% power. Slides were cooled to room temperature, washed twice with BupH Tris Buffered Saline (ThermoFisher, cat. no. 28376,) 0.05% Tween-20 (ThermoFisher, cat. no. BP337-500), and blocked with Pierce™ SuperBlock™ T20 (TBS) Blocking Buffer (ThermoFisher, cat. no. 37581) 0.5% Triton X-100 (ThermoFisher, cat. no. 85111) for 30 min. Blocking buffer was discarded and slides were incubated with primary antibody diluted in blocking buffer overnight at 4 °C. Primary antibodies included α-villin-1 diluted 1:500 (Novus Biologicals, cat. no. NBP1-85335-25ul), α-cytokeratin 5 diluted 1:500 (Abcam, cat. no. Ab75869), α-CC10 diluted 1:400 (Santa Cruz Biotechnology, cat. No. 365992), α-MUC5B diluted 1:2000 (ThermoFisher, cat. no. PA5-82342), α-MUC5AC diluted 1:1000 (ThermoFisher, cat. no MA5-12178), α-p63 diluted 1:1250 (Abcam, cat. no. Ab124762), α-FoxJ1 diluted 1:300 (ThermoFisher, cat no. 14-9965-82), and α-CFTR diluted 1:40 (ThermoFisher, cat. No. 1080-MSM9-P1). Slides were washed three times with wash buffer, then incubated with fluorescent secondary antibodies. Secondary antibodies (all diluted 1:500 in blocking buffer) were donkey anti-rabbit 647 (ThermoFisher, cat. no. A31573) or donkey anti-mouse 647 (ThermoFisher, cat. no. A31571). Slides were incubated in secondary for one hour at room temperature. Slides were washed twice with wash buffer and three times with water, then mounted with coverslips using mounting medium from Pierce™ Peroxidase IHC Detection Kit (ThermoFisher, cat. no. 36000). Slides were imaged using a Keyence BZ-X810 Fluorescent Microscope (Keyence, cat. no. BZ-X810) and images were analyzed using BZ-X800 Analyzer.
Cell type quantification by immunofluorescence labeling
Immunofluorescence images of the villin-1, KRT5, CC10, MUC5AC and MUC5B-labeled HAE cultures were quantified using Imaris microscopy image analysis software (v9.8.2). The stitched image depicting the entire length of each HAE culture transection for each cell type stain was used for image quantification. Two technical replicate images of each biological replicate BCi-NS1.1 or primary HAE culture were used for image quantification. For each marker, signal area was determined using the surfaces creation tool and setting the intensity threshold to exclude any background as determined by HAE culture controls labeled without primary antibody. Total area of phalloidin signal was used as the total area of each HAE culture transection, and the percentage occupancy of cell type marker is represented as a percentage of this phalloidin labeled area.
Isolation and processing of single cell suspensions from HAE cultures for scRNA-Seq
HAE cultures maintained for 3 weeks post-airlift were processed for scRNA-Seq. Single cell suspensions of HAE cultures were prepared as described previously . To enable multiplexing and doublet detection, cells were labeled with “cell hashing” antibodies as described previously . Briefly, approximately 200,000 cells per sample were resuspended in staining buffer (PBS, 2% BSA, and 0.01% Tween) and incubated for 10 min with Fc block (TruStain FcX, BioLegend) and FcR blocking reagent (Miltenyi Biotec). Cells were then incubated with oligonucleotide-conjugated hashing antibodies (generated in-house by the New York Genome Center as described ) for 30 min at 4 °C. After labelling, cells were washed three times in staining buffer. After the final wash, cells were resuspended in PBS supplemented with 0.04% BSA, filtered, and counted. Cells were pooled (6 samples per pool in equal proportions), super-loaded to the 10X Genomics Chromium Controller (approximately 25,000 cells loaded to 1 Chromium lane) and processed with the Chromium Next GEM single-cell 5′ library and gel bead kit according to manufacturer’s protocols. Hashtag-oligonucleotide (HTO) additive oligonucleotide primer was spiked into the cDNA amplification PCR, and the HTO library was prepared as described previously . Gene expression and HTO libraries were pooled and sequenced in multiplex on the Illumina Novaseq 6000 system with the following read configuration: read1 28 cycles, Index read1 8 cycles, read2 91 cycles. Gene expression libraries were sequenced to an estimated depth of 32,000 read pairs per cell, and HTO libraries were sequenced to an estimated depth of 2,400 read pairs per cell.
Read mapping and quantification
CellRanger count (v6.1.2, 10X Genomics Inc.) was used to assign cell barcodes and map reads to transcriptome reference (GRCh38-2020-A, GENCODE v32/Ensembl 98) with default parameters. HTO quantification was also performed using CellRanger count via feature barcode detection.
Data preprocessing and integration
scRNA-Seq analyses were performed in R (v4.0.3) using Seurat (v4.0) . Samples were demultiplexed by HTO counts using the HTODemux function with default parameters. Demultiplexing quality control (QC) was assessed by inspection of a dimensionality reduction plot of HTO counts by t-SNE. Cell barcode × feature (gene) count matrices were filtered to exclude cells with > 18.7% mitochondrial gene content (dead or dying cells) or those with feature (gene) counts outside of the interval from 1505 to 5857 or unique molecular identifier (UMI) counts outside of the interval from 4569 to 24992, which correspond to 1.2 and 1.5 standard deviations from the mean number of counts (log), respectively. Intra-hash heterotypic doublet detection and removal were performed with scDoubletFinder (v1.4.0) . Gene expression data were normalized with scTransform (v.0.3.3)  using top 3,000 variable genes and including a regression factor for “cell cycle difference,” determined by subtracting the G2/M score from the S score for each cell calculated using the gene sets described by Tirosh and colleagues . To minimize the effect of donor-specific variation on cluster assignment, data integration was performed across all replicates using 2,000 gene features and Seurat commands PrepSCTIntegration, FindIntegrationAnchors and IntegrateData.
Data clustering and annotation
Principal component analysis (PCA) was performed on normalized, integrated data, and the first 30 principal components (PCs) were used for clustering and nearest neighbor detection. An iterative approach was used to remove low quality clusters of cells: initial clustering was performed using the graph-based Louvain algorithm with multilevel refinement  and a resolution parameter of 1.5. Two clusters (333 cells total) were excluded from downstream analysis due to low UMI counts and high mitochondrial content, respectively. Data were then re-clustered and assessed for stability at a range of resolutions using Clustree (v.0.5.0) . A final resolution parameter of 1.4 was selected from these results. Clusters were annotated to canonical airway cell populations based on expression of established marker genes [4, 5, 7, 20]. Cell annotations were verified by comparing assignments to predictions obtained by label transfer from a recently published scRNA-seq dataset of similar HAE cultures .
Pseudobulk differential gene expression analyses
Differential gene expression analyses were conducted on “pseudobulk” gene expression profiles using edgeR (v3.32.1)  and scran (v1.18.7) , a strategy found to provide better type I error control than individual cell-based differential expression workflows . Briefly, per cell gene counts were summed to aggregate “pseudobulk” profiles for each cell population replicate, with a minimum threshold for inclusion of 25 cells. Gene expression linear models were constructed for each cell population with at least two samples per condition (BCi-NS1.1, primary) passing inclusion thresholds. Genes expressed in fewer than 10% of cells in any pseudobulk sample were excluded from differential expression tests. Differential gene expression was defined as an adjusted p-value < 0.05 (Benjamini-Hochberg) and an estimated log2 fold-change > 1.
Pseudobulk gene set enrichment analyses
Gene set enrichment analyses were performed using CAMERA  for the molecular signatures database (MSigDB)  C2 canonical pathways gene set (n = 3050 sets) and HALLMARK gene set (n = 50 sets) collections. Significant set enrichment was defined as an adjusted p-value < 0.05.
Pseudotime trajectory inference
Cell differentiation trajectories were inferred using Slingshot . Analyses included only cells annotated along the basal-secretory axis. Principal curves were fit in integrated UMAP space and rooted to the basal cell cluster. Differential progression testing was assessed using a Kolmogorov–Smirnov Goodness of Fit Test between conditions using the progressionTest function with default parameters.
Differential abundance testing
Cell abundance was modeled using edgeR as described , without normalization for total cell count. Significant differential abundance was defined as any change in cell population frequency at a p-value < 0.05.
Infection of HAE cultures with influenza A virus
HAE cultures were washed apically twice with room temperature PBS, after which 5 × 106 plaque-forming units (PFU) of influenza A/California/07/2009 virus (BEI resources, NR-13663) was added apically in 50 µl of PBS and incubated for one hour at 37 °C. HAE cultures were then washed apically once with room-temperature PBS and incubated at air–liquid interface at 37 °C for 24 h. Mock-infected cultures were treated with the same protocol as IAV-infected cultures but using only PBS. After 24 h, cultures were harvested for quantification of viral PFU, RT-qPCR, LDH quantification, imaging, or western blot analyses. For viral PFU, a 200 µl apical wash was collected from the cultures, frozen at − 80 °C and later quantified by plaque assay on MDCK cells. For RT-qPCR, cultures were washed once, and the Transwell membrane was excised and submerged in buffer RLT (Qiagen) and frozen at − 80 °C. RNA extraction was performed using the RNeasy Kit (Qiagen, cat no. 74104), cDNA synthesis using SuperScript™ III First-Strand Synthesis System (ThermoFisher, cat no. 18080051) and RT-qPCR using PowerUp™ SYBR™ Green Master Mix (ThermoFisher, cat no. A25741) with RPS11 as a housekeeping gene. For LDH quantification a 200 µl apical wash was collected from the cultures and LDH release was measured in 50 µl of apical wash in triplicate using CyQUANT™ LDH Cytotoxicity Assay (ThermoFisher, cat. no. C20300). For western blot analysis, cultures were washed once and the cells were lysed by adding 200 µl of 1 × Invitrogen™ 4X Bolt™ LDS Sample Buffer (cat. no. B0007, Fisher Scientific). Lysate was frozen at -20 °C for western blot analysis. For imaging, cultures were submerged in 4% PFA, incubated overnight at 4 °C, then transferred to PBS. Staining and imaging of infected cultures was performed as described previously by our laboratory .
Infection of HAE cultures with S. aureus
S. aureus USA300 strain AH-LAC  was routinely grown at 37 °C on tryptic soy agar (TSA). For infection of HAE cultures, overnight bacterial cultures were grown in 5 mL tryptic soy broth (TSB) in 15-mL conical tubes under shaking conditions at 180 rpm and 37 °C with a 45° angle. A 1:100 dilution of overnight culture was subcultured into TSB and incubated for another 3 h as above. Bacteria were washed once with 5 mL PBS then diluted using Optical Density 600 (OD600) measurements to an OD of 0.8, corresponding to a concentration of 5 × 108 CFU/mL in PBS. HAE cultures were washed apically twice with room temperature PBS, then 1 × 106 CFU of AH-LAC was added apically in 50 µl of PBS and incubated for one hour at 37 °C. PBS was aspirated after 1 h of infection and HAE cultures were washed apically three times with room-temperature PBS, then incubated at 37 °C for 18 h. After 18 h, cultures were harvested for either quantification of bacterial CFU, LDH quantification or imaging. For quantification of bacterial CFU, the Transwell membrane was excised and submerged in 1% Saponin to lyse cells, and lysate was serially diluted and plated on TSA. Samples were harvested and prepared for LDH quantification and imaging as described above.
Treatment of HAE cultures with IFN-β
For treatment with IFN-β, basolateral media was removed from the cultures and media containing 250 U (500 U/mL) human IFN-β (Millipore Sigma, cat. no. IF014, lot 3481402) was added to the basolateral chamber. IFN-β-treated cultures were incubated for 24 h at 37 °C. After 24 h, IFN-β-treated cultures were harvested and prepared for either RT-qPCR or western blot analysis as described above.
CFTR functional assays
To assess the function of CFTR channels in hTERT-derived HAE cultures, transepithelial electrical resistance (TEER) of fully differentiated cultures was measured at baseline. Cultures were then treated with either untreated Pneumacult™ ALI basal media as a control, 0.2% DMSO as a vehicle control, or 20µM forskolin and 50µM 3-Isobutyl-1-methylxanthine (IBMX) to activate CFTR channels. TEER of cultures treated with each condition were subsequently measured at 10 min, 30 min, and 60 min post-treatment.
To assess similarities and differences between HAE cultures derived from the hTERT-expressing bronchial basal epithelial cell line BCi-NS1.1 and unmodified, primary-derived precursor cells, we generated HAE cultures from three donors of primary normal human bronchial epithelial cells, as well as three biological replicates of HAE cultures from BCi-NS1.1 (Fig. 1a).
Both primary and BCi-NS1.1 precursor cells differentiate into pseudostratified epithelia that include basal, secretory, and ciliated cells
We first characterized the histological morphology of HAE cultures derived from BCi-NS1.1 cells and primary cells. H&E staining of HAE culture sections demonstrated that cultures derived from both BCi-NS1.1 and primary precursors adopt a polarized pseudostratified epithelial morphology with ciliated cells facing the apical side (Fig. 1b). PAS-Alcian blue staining revealed the presence of intracellular mucins (dark blue), consistent with the presence of secretory cells (Fig. 1b). Overall, individual cultures exhibited subtle differences in epithelium thickness, secretory cell frequency, and secretory cell hyperplasia, but these differences were not associated with precursor cell source.
Next, we determined the presence and intra-epithelial organization of ciliated, basal, and secretory cells by immunofluorescence microscopy. Cultures were immunolabeled for canonical cell type-specific markers: villin-1 for ciliated cells, cytokeratin-5 (KRT5) for basal cells, and CC10 (encoded by the SCGB1A1 gene), MUC5B and MUC5AC for secretory cells [19, 20, 23, 47]. We observed KRT5 expression in small, brick-shaped cells lining the basolateral surface consistent with basal cells (Fig. 1c, Additional file 4: Fig. S1a). Image quantification of this marker identified KRT5-positive cells to account for approximately 31% of total cell area in both BCi-NS1.1- and primary-derived HAE cultures (Fig. 1d). We observed villin-1 in columnar cells lining the apical surface, consistent with ciliated cells (Fig. 1c, Additional file 4: Fig. S1b). These cells accounted for approximately 20% of total cell area in both BCi-NS1.1- and primary-derived HAE cultures (Fig. 1e). Lastly, secretory cell makers were present in morphologically heterogeneous cells, including large globular cells facing the apical surface (Fig. 1c, Additional file 4: Fig. S1c). CC10-positive cells made up on average 28% of total cell area in both BCi-NS1.1- and primary-derived HAE cultures, although the number of cells positive for this marker was more variable between cultures than either KRT5 or villin-1 (Fig. 1f). Cells positive for either MUC5AC or MUC5B generally made up less than 20% of the total cell area (Fig. 1c, g, Additional file 4: Fig. S1d). While technical constraints prohibited multiplex staining for CC10 and mucins, dual staining for MUC5B and MUC5AC revealed that most but not all MUC-positive cells expressed one of these two mucins but not both (Fig. 1c, g, Additional file 4: Fig. S1d). Although variable across cultures, the proportions of MUC5B vs MUC5AC positive cells (Fig. 1g, Additional file 4: Fig. S1d) were not associated with precursor cell source. Taken together, immunofluorescence labeling of cell-type-specific markers revealed that ciliated, basal, and secretory cells are present in both BCi-NS1.1 and primary-derived cultures at similar frequencies.
scRNA-Seq resolves constituent cell populations of BCi-NS1.1- and primary-derived HAE cultures
To assess the cellular composition at high resolution along with associated gene expression patterns of BCi-NS1.1- and primary-derived HAE cultures, we performed scRNA-Seq. In total, we analyzed 12,801 single cells (post quality-control filtering) from three samples of each HAE culture source (n = 3 independent differentiations of BCi-NS1.1, n = 3 different donors of primary cells). Unsupervised graph-based clustering identified 12 airway epithelial cell populations, which were annotated based on RNA expression of canonical marker genes (Fig. 2a) and label transfer predictions from a recently reported scRNA-Seq analysis of HAE cultures. scRNA-Seq data included basal cells (KRT5+, TP63hi), proliferating cells (MKI67+, TOP2A+), suprabasal cells (KRT6A+, KRT16+),”intermediate “ (i.e. along the basal-to-secretory differentiation axis) cells, three populations of secretory cells (I: MUC5AClo MUC5Blo, II: MUC5AClo MUC5Bhi, and III: MUC5AChi MUC5Blo), deuterosomal cells (DEUP1+, FOXN4+), ciliated cells (FOXJ1+,
CAPS+), tuft-like cells (LRMP+, ASCL2+)[4, 7], ionocytes (FOXI1+, CFTR+)[5, 7], and neuroendocrine cells (ASCL1+, HOXB5+)(Fig. 2b). Cell representation and annotation was consistent across samples (Additional file 4: Fig. S2a). The relative frequencies of most cell populations were similar between primary and BCi-NS1.1-derived HAE cultures, with the exception of MUC5AChi secretory III cells, deuterosomal cells, and ciliated cells, which were more highly represented in BCi-NS1.1-derived HAE cultures (Fig. 2c; Additional file 1). Notably, both culture sources exhibited rare ionocytes, tuft-like cells, and neuroendocrine cells at similar frequency. Unexpectedly, we detected far fewer ciliated cells in scRNA-Seq data than by immunofluorescence (~ 2% of cells vs ~ 20% of total cell area). Based on our prior experience analyzing HAE by scRNA-Seq, we attribute this discrepancy to technical issues during sample processing resulting in the selective loss of ciliated cells prior to 10X Chromium loading. With severely reduced ciliated cell numbers, we focused the remainder of our scRNA-Seq analyses on robustly sampled basal and secretory cells. Despite these limitations, scRNA-Seq identified the expected constituent cell populations of the airway epithelium with few differences in frequency across culture types.
Pseudotime analysis of scRNA-Seq data infers a common basal-to-secretory cell differentiation trajectory but differences in distribution of BCi-NS1.1-derived HAE cultures
As we observed subtle differences in the frequency of cell populations between primary- and BCi-NS1.1-derived HAE cultures, we investigated the possibility that these differences in cellular composition could be related to the airway cell differentiation process. To this end, we inferred differentiation trajectories by pseudotime analysis. We restricted our analysis to basal, suprabasal, intermediate, and secretory populations based on the established basal-to-secretory differentiation path  and to eliminate potential confounding factors from actively proliferating cells and rare cell types for which we lacked robust sampling. A single lineage was identified for both culture types: basal-to-suprabasal-to-intermediate-to-secretory (with component secretory populations ordered I, III, II) (Fig. 2d). Differential distribution testing indicated significant differences in progression along this trajectory between the culture conditions, with primary-derived HAE cultures exhibiting relatively more cells in basal and “early” secretory states in comparison to the BCi-NS1.1-derived cultures (Fig. 2e). A closer inspection of the pseudotime distribution for each cell population illustrated several shifts in the median pseudotime values: in comparison to primary-derived HAE cultures, the median basal and secretory I pseudotime values were higher in BCi-NS1.1-derived HAE cultures, while the median pseudotime value for intermediate cells was lower (Fig. 2f). This suggests that, while general cell population abundances are consistent across culture types, there exist modest but significant differences in the distribution of BCi-NS1.1-derived HAE cultures along the basal-to-secretory differentiation axis.
HAE culture cell populations exhibit generally similar gene expression patterns across BCi-NS1.1- and primary-derived cultures
We next assessed gene expression patterns in each cell population to compare HAE cultures from the primary- and BCi-NS1.1-derived HAE cultures. “Pseudobulk” gene expression profiles were aggregated (by cell population and replicate) and clustered by Spearman correlation distance (Fig. 3a). Data clustered primarily by major cell population (i.e. basal, secretory, ionocyte, ciliated/deuterosomal), and pairwise comparisons further emphasize similarities across the different precursor cell groups by cell population. Interestingly, BCi-NS1.1-derived ‘intermediate’ cell populations clustered with BCi-NS1.1-derived suprabasal cell populations, a result consistent with our pseudotime analysis. Similarly, PCA revealed a strong segregation by cell population in the first three PCs (accounting for 72.6% total variation), with basal/ciliated, basal/secretory, and ionocyte divergence captured by the first three components, respectively (Fig. 3b). PC4 appeared to segregate data by culture type but described a relatively small amount (6.1%) of total variation. Taken together, these results indicate that by gene expression, cell populations are largely similar across BCi-NS1.1 and primary HAE cultures, with a relatively small amount of variation associated with differences in the precursor cell groups.
Pseudobulk differential gene expression reveals higher interferon-stimulated gene (ISG) expression in some BCi-NS1.1-derived HAE culture cell populations
Despite the high degree of similarity in HAE culture cell populations across precursor cell groups, the relatively minor gene expression differences observed between BCi-NS1.1 and primary-derived HAE culture cell populations in PCA could be biologically meaningful. Therefore, for each cell population, we conducted differential gene expression analysis directly comparing BCi-NS1.1 and primary HAE cultures. In total, we detected 1,052 genes differentially expressed in at least one cell population of those with sufficient cell numbers for testing (basal, suprabasal, intermediate, secretory I, secretory II, Fig. 3c, Additional file 2). A “core” set of 30 genes were detected as differentially expressed (same direction) in all cell populations tested (Fig. 3d). The majority of differentially expressed genes (DEGs) were detected in only one cell population (Fig. 3d). As expected, TERT expression was higher in all BCi-NS1.1-derived HAE culture cell populations, though it failed to meet differential expression thresholds in all but secretory I and II cells due to sparse detection typical of droplet scRNA-Seq (Fig. 3c, Fig. S2a). We also compared expression levels of the cell type markers at single cell resolution (Fig. 2b) to the results of pseudobulk contrasts. Notably, although BCi-NS1.1 basal cells exhibit modestly lower expression of TP63 in the single cell marker gene dot plot (Fig. 2b), this difference (log2FC = − 0.53) did not clear differential expression thresholds in pseudobulk contrasts. MUC5AC expression was found to be significantly higher in secretory I and II populations from BCi-NS1.1-derived HAE cultures in pseudobulk analyses (consistent with pattern in marker gene dot plot, Fig. 2b), but overall expression was considerably lower than that observed in secretory III cells. Insufficient secretory III cell numbers across biological replicates precluded formal pseudobulk differential expression testing between culture sources.
To place DEGs in biological context, we performed gene set enrichment testing on the C2 canonical pathway gene sets from the Molecular Signatures Database. Several gene sets related to interferon (IFN) signaling were found to be enriched in BCi-NS1.1 cultures, particularly in secretory II cells (Fig. 3e, Additional file 4: Fig. S2a). Relatedly, we observed differential expression of ISGs MX1, IFIT1, and ISG15 in secretory II cells and ISG15 and MX1 expression in basal cells (indicated in Fig. 3c). Thus, while gene expression is largely similar, a relatively small number of genes are differentially expressed between BCi-NS1.1 and primary-derived HAE culture cell populations. Among these, we observed an elevated ISG signature in some BCi-NS1.1-derived cell populations.
JAK-STAT pathway activity and expression of select ISGs are higher at steady state in BCi-NS1.1 relative to primary HAE cultures, but are induced to similar levels upon IAV infection
As HAE cultures are frequently employed for infection studies by our group and others [8, 15,16,17,18, 49] we sought to further characterize potential differences in ISG expression and corresponding signaling between these two culture types at steady state, in response to IFN stimulation, and in response infection with respiratory pathogens. We apically infected both types of cultures with two human respiratory pathogens: IAV and S. aureus. As a positive control for ISG induction, we treated cultures basolaterally with human IFN-β. Mock-infected cultures served as controls. We measured the expression of select ISGs (identified in scRNA-Seq analyses) by RT-qPCR and by western blot (Fig. 4a–f). We also measured phosphorylated STAT1 to assess JAK-STAT pathway activity. Mock-treated BCi-NS1.1-derived HAE cultures showed higher expression of ISGs examined (MX1, IFIT1, ISG15) than primary-derived HAE cultures, as observed in our scRNA-Seq analysis (Fig. 3c, e). This difference was statistically significant for MX1 at both RNA and protein levels (Fig. 4a, b). Concomitantly, phosphorylated STAT1 was elevated in BCi-NS1.1-derived cultures at steady-state (Fig. 4b).
We next sought to determine whether these differences in steady-state ISG expression result in functional differences upon infection with common respiratory pathogens. As expected, we observed mRNA upregulation of ISGs (MX1, IFIT1 and ISG15) upon IAV infection (Fig. 4c, d) and IFN-ß treatment (Fig. 4e, f). Interestingly, despite the higher baseline ISG expression in BCi-NS1.1-derived cultures at steady-state, ISG RNA expression was upregulated to similar levels in in primary cell-derived HAE cultures and BCi-NS1.1-derived cultures. MX1 protein levels, while different between culture types at baseline (Fig. 4b), were also induced to similar levels upon IAV infection (Fig. 4d) and IFN-ß treatment (Fig. 4f). Finally, phosphorylated STAT1 levels, while higher in BCi-NS1.1-derived cultures at steady state (Fig. 4b), exhibited an inverse pattern upon IAV infection (Fig. 4c). No difference in phosphorylated STAT1 levels was observed upon IFN-ß stimulation (Fig. 4f). We did observe more variable ISG upregulation (RNA) in primary-derived HAE cultures than in BCi-NS1.1 cell-derived HAE cultures (Fig. 4c, e), possibly due to donor-associated variation.
Primary and BCi-NS1.1-derived HAE cultures exhibit similar pathogen loads and tissue damage upon infection
We next determined whether the observed increased baseline ISG expression in BCi-NS1.1-derived HAE cultures affects pathogen loads upon infection. Thus, at 24- or 18-h post-infection we measured IAV plaque-forming units (PFU) or S. aureus colony-forming units (CFU), respectively, in apical washes from BCi-NS1.1- and primary cell-derived HAE cultures. In addition, we measured pathogen-elicited tissue damage by lactate dehydrogenase (LDH) release and fluorescence microscopy. We found that S. aureus loads were equal between both precursor cell types (Fig. 4g). IAV loads trended slightly higher in primary cell-derived HAE cultures (Fig. 4h, l, k, Additional file 4: Fig. S4a), but these differences did not reach statistical significance (p = 0.1 and p = 0.2 for plaque assay and fluorescence data respectively, Mann–Whitney U test). IAV infection resulted in low-level LDH release with no detectable differences across precursor cell types (Fig. 4j). S. aureus infection resulted in considerable LDH release, with higher levels in BCi-NS1.1-derived cultures. These patterns were consistent with corresponding fluorescence microscopy of infected cultures to assess cellular damage (Fig. 4l, Additional file 4: Fig. S4b). Taken together, these results demonstrate that, despite some minor differences, BCi-NS1.1- and primary-derived HAE cultures behave similarly upon infection with two major respiratory pathogens, IAV and S. aureus.
BCi-NS1.1-derived HAE cultures express functional CFTR channels
As with any cell culture derived from a single donor, BCi-NS1.1 could exhibit donor-specific properties compatible or incompatible with different research applications. As HAE cultures are widely utilized for studies of airway function in cystic fibrosis, in which dysfunctional CFTR channels result in chronic airway disease, we assessed the expression of CFTR in HAE cultures. In scRNA-Seq data, as expected, mRNA expression of CFTR was highest in ionocytes in both BCi-NS1.1-derived and primary derived cultures (Fig. 2b, Additional file 3). Ionocyte CFTR mRNA expression levels appeared generally similar across cultures, though insufficient ionocyte numbers precluded a formal statistical comparison. Immunofluorescence labeling (Fig. 5a) confirmed similar CFTR protein expression patterns in both BCi-NS1.1-derived and primary derived cultures. To assess CFTR function in BCi-NS1.1-derived HAE cultures, we measured transepithelial electrical resistance (TEER) of cultures treated with activators of CFTR channels (forskolin and IBMX). Cultures treated with media or vehicle controls exhibited steady TEER over the course of one hour. However, cultures treated with activators of CFTR channels exhibited a decrease in TEER corresponding to CFTR channel opening beginning at 10 min post-activation and extending for up to one hour (endpoint of the experiment). These results indicate that CFTR channels in BCi-NS1.1-derived HAE cultures function similarly to those previously described in primary HAE cultures. 
Combining scRNA-Seq, immunofluorescence, and functional analyses, we have shown that HAE cultures derived from the hTERT-expressing airway basal cell line BCi-NS1.1 exhibit similar cell type composition and similar cell type-associated gene expression to HAE cultures derived from primary cells. While we observed evidence of heightened JAK/STAT pathway activity and corresponding increased baseline expression of several ISGs in BCi-NS1.1-derived HAE cultures at steady-state, both precursor cell types induced ISG expression to similar levels upon infection with IAV or upon IFN treatment. IAV infections resulted in apparent but statistically insignificant higher viral loads in primary-derived cultures as compared to BCi-NS1.1-derived cultures. These modest differences are consistent with the elevated steady-state ISG expression observed in BCi-NS1.1-derived cultures, but could also be a consequence of general variation across donor sources. We do not believe that this potential minor dissimilarity should discount the use of BCi-NS1.1 in respiratory virus studies, but should perhaps be considered in infection experiments incorporating multiple culture sources. Finally, we found that BCi-NS1.1-derived HAE cultures express functional CFTR channels, a characteristic important for the study of cystic fibrosis in this model.
Of note, a comprehensive comparison of component cell type gene expression patterns in BCi-NS1.1-derived and primary-derived HAE cultures in this study was somewhat limited due to the unexpectedly low frequency of ciliated cells in our scRNA-Seq data. The underrepresentation of ciliated cells in scRNA-Seq does not reflect the true composition of BCi-NS1.1-derived and primary-derived HAE cultures, but is likely a technical artifact related to cell dissociation and processing prior to loading the 10X Genomics Chromium controller. The proportion of ciliated cells as quantified by villin-1 labeling in our samples is similar across precursor cell types (Fig. 1b, c, e) and is within the expected range for HAE [17, 18]. Indeed, our group and others [17, 18, 22] have effectively detected ciliated cells at expected frequencies by 10X Genomics scRNA-Seq methods. For the unexpectedly few ciliated cells detected in scRNA-Seq data, most of which are from BCi-NS1.1 samples, we observed robust expression of ciliated cell marker genes. In addition, ciliated cells from BCi-NS1.1-derived cultures exhibit highly concordant scRNA-Seq gene expression patterns with ciliated cells from primary-derived HAE cultures in a recently published scRNA-Seq dataset . Although insufficient scRNA-Seq ciliated cell counts precluded a statistically robust differential expression analysis, these observations, along with our additional histological and functional comparisons, are generally consistent with similar ciliated cell features in BCi-NS1.1-derived and primary-derived HAE cultures.
Overall, our results further support BCi-NS1.1-derived HAE cultures as a relevant model system that recapitulates airway epithelial biology similarly to primary-derived HAE cultures. These cultures are of particular benefit in experiments that demand large quantities of HAE cultures/replicates beyond those available from primary cell sources. Additionally, the extended passage capacity of BCi-NS1.1 cells has enabled sophisticated genetic engineering techniques (with corresponding antibiotic selection) for mechanistic studies of airway cell function and disease that can be technically challenging or infeasible in primary cells. Despite overall similarities, we did observe modest differences between culture types. These differences may stem from either donor-to-donor variation, as all BCi-NS1.1-derived cultures were necessarily derived from the same donor while primary-derived cultures were derived from three distinct donors, and/or directly from engineered expression of hTERT in BCi-NS1.1 cells. Future studies in which hTERT expression is engineered in multiple independent primary cell sources may disentangle these mechanisms. Further, determining the specific effects of hTERT in HAE cultures will be important in establishing and interpreting hTERT-engineered culture models of specific demographic groups or disease states for which primary-derived HAE culture systems have been proven highly informative, such as cystic fibrosis, COPD, and asthma.
Availability of data and materials
scRNA-Seq data generated for this study are available via the NCBI Gene Expression Omnibus repository with accession number GSE225765. All analysis code is available on GitHub at https://github.com/BradRosenbergLab/airwayepithelialculturecomparison.
Knight DA, Holgate ST. The airway epithelium: Structural and functional properties in health and disease. Respirology. 2003;8:432–46. https://doi.org/10.1046/j.1440-1843.2003.00493.x.
Eon Kuek L, Lee RJ. First contact: the role of respiratory cilia in host-pathogen interactions in the airways. Am J Physiol Lung Cell Mol Physiol. 2020;319:603–19. https://doi.org/10.1152/ajplung.00283.2020.-Respi.
Branchfield K, Nantie L, Verheyden JM, Sui P, Wienhold MD, Sun X. Pulmonary neuroendocrine cells function as airway sensors to control lung immune response. Science. 1979;2016(351):707–10. https://doi.org/10.1126/science.aad7154.
Hewitt RJ, Lloyd CM. Regulation of immune responses by the airway epithelial cell landscape. Nat Rev Immunol. 2021;21:347–62. https://doi.org/10.1038/s41577-020-00477-9.
Plasschaert LW, Žilionis R, Choo-Wing R, Savova V, Knehr J, Roma G, et al. A single-cell atlas of the airway epithelium reveals the CFTR-rich pulmonary ionocyte. Nature. 2018;560:377–81. https://doi.org/10.1038/s41586-018-0394-6.
Ualiyeva S, Hallen N, Kanaoka Y, Ledderose C, Matsumoto I, Junger WG, et al. Airway brush cells generate cysteinyl leukotrienes through the ATP sensor P2Y2. vol. 5. 2020.
Montoro DT, Haber AL, Biton M, Vinarsky V, Lin B, Birket SE, et al. A revised airway epithelial hierarchy includes CFTR-expressing ionocytes. Nature. 2018;560:319–24. https://doi.org/10.1038/s41586-018-0393-7.
Ibricevic A, Pekosz A, Walter MJ, Newby C, Battaile JT, Brown EG, et al. Influenza Virus Receptor Specificity and Cell Tropism in Mouse and Human Airway Epithelial Cells. J Virol. 2006;80:7469–80. https://doi.org/10.1128/jvi.02677-05.
Wu CT, Lidsky P v., Xiao Y, Cheng R, Lee IT, Nakayama T, et al. SARS-CoV-2 replication in airway epithelia requires motile cilia and microvillar reprogramming. Cell 2023;186:112–130.e20. https://doi.org/10.1016/j.cell.2022.11.030.
Gagneux P, Cheriyan M, Hurtado-Ziola N, Brinkman Van Der Linden ECM, Anderson D, McClure H, et al. Human-specific Regulation of α2–6-linked Sialic Acids. Journal of Biological Chemistry 2003;278:48245–50. https://doi.org/10.1074/jbc.M309813200.
Shinya K, Ebina M, Yamada S, Ono M, Kasai N, Kawaoka Y. Influenza virus receptors in the human airway. 2006.
Dohrman A, Miyata S, Gallup M, Li J-D, Chapelin C, Coste A, et al. Mucin gene MUC 2 and MUC 5AC upregulation by Gram-positive and Gram-negative bacteria. vol. 1406. 1998.
Chen R, Lim JH, Jono H, Gu XX, Kim YS, Basbaum CB, et al. Nontypeable Haemophilus influenzae lipoprotein P6 induces MUC5AC mucin transcription via TLR2-TAK1-dependent p38 MAPK-AP1 and IKKβ- IκBα-NF-κB signaling pathways. Biochem Biophys Res Commun. 2004;324:1087–94. https://doi.org/10.1016/j.bbrc.2004.09.157.
Kim YO, Jung MJ, Choi JK, Ahn DW, Song KS. Peptidoglycan from staphylococcus aureus increases MUC5AC gene expression via RSK1- CREB pathway in human airway epithelial cells. Mol Cells. 2011;32:359–65. https://doi.org/10.1007/s10059-011-0118-3.
Davis AS, Chertow DS, Moyer JE, Suzich J, Sandouk A, Dorward DW, et al. Validation of Normal Human Bronchial Epithelial Cells as a Model for Influenza A Infections in Human Distal Trachea. J Histochem Cytochem. 2015;63:312–28. https://doi.org/10.1369/0022155415570968.
Vaidyanathan S, Salahudeen AA, Sellers ZM, Bravo DT, Choi SS, Batish A, et al. High-Efficiency, Selection-free Gene Repair in Airway Stem Cells from Cystic Fibrosis Patients Rescues CFTR Function in Differentiated Epithelia. Cell Stem Cell. 2020;26:161-171.e4. https://doi.org/10.1016/j.stem.2019.11.002.
Kelly JN, Laloli L, V’kovski P, Holwerda M, Portmann J, Thiel V, et al. Comprehensive single cell analysis of pandemic influenza A virus infection in the human airways uncovers cell-type specific host transcriptional signatures relevant for disease progression and pathogenesis. Front Immunol 2022;13. https://doi.org/10.3389/fimmu.2022.978824.
Ravindra NG, Alfajaro MM, Gasque V, Huston NC, Wan H, Szigeti-Buck K, et al. Single-cell longitudinal analysis of SARS-CoV-2 infection in human airway epithelium identifies target cells, alterations in gene expression, and cell state changes. PLoS Biol 2021;19. https://doi.org/10.1371/JOURNAL.PBIO.3001143.
Gray TE, Guzman K, William Davis C, Abdullah LH, Nettesheim P. Mucociliary Differentiation of Serially Passaged Normal Human Tracheobronchial Epithelial Cells. Am J Respir Cell Mol Biol. 1996;14:104–12.
Garcıá SR, Deprez M, Lebrigand K, Cavard A, Paquet A, Arguel MJ, et al. Novel dynamics of human mucociliary differentiation revealed by single-cell RNA sequencing of nasal epithelial cultures. Development (Cambridge) 2019;146. https://doi.org/10.1242/dev.177428.
Chu HW, Rios C, Huang C, Wesolowska-Andersen A, Burchard EG, O’Connor BP, et al. CRISPR-Cas9-mediated gene knockout in primary human airway epithelial cells reveals a proinflammatory role for MUC18. Gene Ther. 2015;22:822–9. https://doi.org/10.1038/gt.2015.53.
Carraro G, Langerman J, Sabri S, Lorenzana Z, Purkayastha A, Zhang G, et al. Transcriptional analysis of cystic fibrosis airways at single-cell resolution reveals altered epithelial cell states and composition. Nat Med. 2021;27:806–14. https://doi.org/10.1038/s41591-021-01332-7.
Walters MS, Gomi K, Ashbridge B, Moore MAS, Arbelaez V, Heldrich J, et al. Generation of a human airway epithelium derived basal cell line with multipotent differentiation capacity. Respir Res. 2013;14:26–30. https://doi.org/10.1186/1465-9921-14-135.
Zhou H, Brekman A, Zuo W-L, Ou X, Shaykhiev R, Agosto-Perez FJ, et al. POU2AF1 Functions in the Human Airway Epithelium To Regulate Expression of Host Defense Genes. J Immunol. 2016;196:3159–67. https://doi.org/10.4049/jimmunol.1502400.
Iverson E, Griswold K, Song D, Gagliardi TB, Hamidzadeh K, Kesimer M, et al. Membrane-Tethered Mucin 1 Is Stimulated by Interferon and Virus Infection in Multiple Cell Types and Inhibits Influenza A Virus Infection in Human Airway Epithelium. MBio 2022;13. https://doi.org/10.1128/mbio.01055-22.
Song D, Iverson E, Kaler L, Boboltz A, Scull MA, Duncan GA. MUC5B mobilizes and MUC5AC spatially aligns mucociliary transport on human airway epithelium. 2022.
Espersen F, Gabrielsen J. Pneumonia Due to Staphylococcus aureus During Mechanical Ventilation. vol. 144. 1981.
Gillet Y, Vanhems P, Lina G, Bes M, Vandenesch F, Floret D, et al. Factors predicting mortality in necrotizing community-acquired pneumonia caused by Staphylococcus aureus containing panton-valentine leukocidin. Clin Infect Dis. 2007;45:315–21. https://doi.org/10.1086/519263.
de Vries M, Mohamed AS, Prescott RA, Valero-Jimenez AM, Desvignes L, O’connor R, et al. A Comparative Analysis of SARS-CoV-2 Antivirals Characterizes 3CL pro Inhibitor PF-00835231 as a Potential New Treatment for COVID-19 2021. https://doi.org/10.1128/JVI.
Mimitou EP, Cheng A, Montalbano A, Hao S, Stoeckius M, Legut M, et al. Multiplexed detection of proteins, transcriptomes, clonotypes and CRISPR perturbations in single cells. Nat Methods. 2019;16:409–12. https://doi.org/10.1038/s41592-019-0392-0.
Stoeckius M, Zheng S, Houck-Loomis B, Hao S, Yeung BZ, Mauck WM, et al. Cell Hashing with barcoded antibodies enables multiplexing and doublet detection for single cell genomics. Genome Biol 2018;19. https://doi.org/10.1186/s13059-018-1603-1.
Hao Y, Hao S, Andersen-Nissen E, Mauck WM, Zheng S, Butler A, et al. Integrated analysis of multimodal single-cell data. Cell. 2021;184:3573-3587.e29. https://doi.org/10.1016/j.cell.2021.04.048.
McGinnis CS, Murrow LM, Gartner ZJ. DoubletFinder: Doublet Detection in Single-Cell RNA Sequencing Data Using Artificial Nearest Neighbors. Cell Syst. 2019;8:329-337.e4. https://doi.org/10.1016/j.cels.2019.03.003.
Hafemeister C, Satija R. Normalization and variance stabilization of single-cell RNA-seq data using regularized negative binomial regression. Genome Biol 2019;20. https://doi.org/10.1186/s13059-019-1874-1.
Tirosh I, Izar B, Prakadan SM, Wadsworth II MH, Treacy D, Trombetta JJ, et al. Dissecting the multicellular ecosystem of metastatic melanoma by single-cell RNA-seq. n.d.
Rotta R, Noack A. Multilevel local search algorithms for modularity clustering. ACM Journal of Experimental Algorithmics, vol. 16, Association for Computing Machinery; 2011. https://doi.org/10.1145/1963190.1970376.
Zappia L, Oshlack A. Clustering trees: a visualization for evaluating clusterings at multiple resolutions. Gigascience 2018;7. https://doi.org/10.1093/gigascience/giy083.
Robinson MD, McCarthy DJ, Smyth GK. edgeR: A Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics. 2009;26:139–40. https://doi.org/10.1093/bioinformatics/btp616.
Lun ATL, McCarthy DJ, Marioni JC. A step-by-step workflow for low-level analysis of single-cell RNA-seq data [version 1; referees: 5 approved with reservations]. F1000Res 2016;5. https://doi.org/10.12688/F1000RESEARCH.9501.1.
Squair JW, Gautier M, Kathe C, Anderson MA, James ND, Hutson TH, et al. Confronting false discoveries in single-cell differential expression. Nat Commun 2021;12. https://doi.org/10.1038/s41467-021-25960-2.
Wu D, Smyth GK. Camera: A competitive gene set test accounting for inter-gene correlation. Nucleic Acids Res 2012;40. https://doi.org/10.1093/nar/gks461.
Liberzon A, Subramanian A, Pinchback R, Thorvaldsdóttir H, Tamayo P, Mesirov JP. Molecular signatures database (MSigDB) 3.0. Bioinformatics 2011;27:1739–40. https://doi.org/10.1093/bioinformatics/btr260.
Street K, Risso D, Fletcher RB, Das D, Ngai J, Yosef N, et al. Slingshot: Cell lineage and pseudotime inference for single-cell transcriptomics. BMC Genomics 2018;19. https://doi.org/10.1186/s12864-018-4772-0.
Lun ATL, Richard AC, Marioni JC. Testing for differential abundance in mass cytometry data. Nat Methods. 2017;14:707–9. https://doi.org/10.1038/nmeth.4295.
Boles BR, Thoende M, Roth AJ, Horswill AR. Identification of genes involved in polysaccharide- independent Staphylococcus aureus biofilm formation. PLoS One 2010;5. https://doi.org/10.1371/journal.pone.0010146.
Başak K, Kumbul Doguç D, Aylak F, Karadayı N, Gültekin Lütfi Kırdar Kartal F. Effects of Maternally Exposed Food Coloring Additives on Laryngeal Histology in Rats. vol. 33. 2014.
Wang G, Lou HH, Salit J, Leopold PL, Driscoll S, Schymeinsky J, et al. Characterization of an immortalized human small airway basal stem/progenitor cell line with airway region-specific differentiation capacity. Respir Res. 2019;20:1–14. https://doi.org/10.1186/s12931-019-1140-9.
Liberzon A, Birger C, Thorvaldsdóttir H, Ghandi M, Mesirov JP, Tamayo P. The Molecular Signatures Database Hallmark Gene Set Collection. Cell Syst. 2015;1:417–25. https://doi.org/10.1016/j.cels.2015.12.004.
V’kovski P, Gultom M, Kelly JN, Steiner S, Russeil J, Mangeat B, et al. Disparate temperature-dependent virus–host dynamics for SARS-CoV-2 and SARS-CoV in the human respiratory epithelium. PLoS Biol 2021;19. https://doi.org/10.1371/journal.pbio.3001158.
Nilsson HE, Dragomir A, Lazorova L, Johannesson M, Roomans GM. CFTR and tight junctions in cultured bronchial epithelial cells. Exp Mol Pathol. 2010;88:118–27. https://doi.org/10.1016/j.yexmp.2009.09.018.
Moreau-Marquis S, Bomberger JM, Anderson GG, Swiatecka-Urban A, Ye S, O GA, et al. The F508-CFTR mutation results in increased biofilm formation by Pseudomonas aeruginosa by increasing iron availability. Am J Physiol Lung Cell Mol Physiol 2008;295:25–37. https://doi.org/10.1152/ajplung.00391.2007.-Enhanced.
Horndahl J, Svärd R, Berntsson P, Wingren C, Li J, Abdillahi SM, et al. HDAC6 inhibitor ACY-1083 shows lung epithelial protective features in COPD. PLoS One 2022;17. https://doi.org/10.1371/journal.pone.0266310.
Schindler VEM, Alhamdan F, Preußer C, Hintz L, Alhamwe BA, Nist A, et al. Side-Directed Release of Differential Extracellular Vesicle-associated microRNA Profiles from Bronchial Epithelial Cells of Healthy and Asthmatic Subjects. Biomedicines 2022;10. https://doi.org/10.3390/biomedicines10030622.
We thank Dr. Ronald Crystal for his kind gift of BCi-NS1.1 cells. We thank Dr. Margaret Scull for critical reading of the manuscript and helpful suggestions. We also thank all Dittmann and Rosenberg lab members for their thoughts and input on this project. Figs. 1a and 5b were created using Biorender.com.
Research was partially supported by the following grants from NIH/NIAID: R01AI143639 to M.D., AI137336 to V.J.T., R01 AI151029 to B.R.R. and U01 AI150748 to B.R.R. Work was further supported by The Vilcek Institute of Graduate Biomedical Sciences, and by NYU Grossman School of Medicine Startup funds. Both the Experimental Pathology Research Laboratory and the NYU Genome Technology Core are supported by NYU Cancer Center support grant P30CA016087 and by NYU Langone’s Laura and Isaac Perlmutter Cancer Center. Research reported in this paper was supported by the Office of Research Infrastructure of the National Institutes of Health under award number S10OD026880. The content is solely the responsibility of the authors and does not necessarily represent the official views of the National Institutes of Health.
Ethics approval and consent to participate
Consent for publication
The authors declare no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Cell population frequencies by precursor cell source. Cell frequency data for constituent HAE culture cell populations by culture progenitor group.
Pseudobulk differentially expressed genes across culture progenitor groups. Differentially expressed genes calculated from per-cell type pseudobulk profile contrasts for BCi-NS1.1 vs. primary-derived HAE cultures for.
Marker genes by cell population. Positively expressed marker genes for each cell population calculation from the edgeR GLM.
Additional immunofluorescence labelling of fixed and cross-sectioned HAE cultures from BCi-NS1.1-derived and primary-derived donors. Figure S2. scRNA-Seq cell population assignment by replicate and marker gene expression patterns. Figure S3. a Violin plots of log-normalized expression values for ISGs IFIT1, ISG15, and MX1 along with TERT by cell population. Plots are annotated with the FDR from pseudobulk differential expression testing across BCi-NS1.1-derived and primary HAE cultures where significant. b Violin plots of “Hallmark Interferon Alpha Response” gene set scores at single cell resolution by cell population. Low sampling of secretory III, deuterosomal, ciliated cells and ionocytes required their exclusion from pseudobulk contrasts, but they are included here for completeness. Figure S4. Top-down images of HAE cultures with individual channels stained for phalloidin (red) and DAPI (blue) and infected with a. IAV (green) and b. S. aureus (green). Scale bars = 1 mm.
About this article
Cite this article
Prescott, R.A., Pankow, A.P., de Vries, M. et al. A comparative study of in vitro air–liquid interface culture models of the human airway epithelium evaluating cellular heterogeneity and gene expression at single cell resolution. Respir Res 24, 213 (2023). https://doi.org/10.1186/s12931-023-02514-2
- Airway epithelium
- Innate immunity
- Respiratory infection
- Air–liquid interface epithelial culture
- Single cell transcriptomics