Introduction

Chronic obstructive diseases such as COPD and asthma are major public health challenges characterised by airflow limitation. There is a lack of specific curative treatments for these diseases, in part because of the paucity of integrative approaches for deciphering their pathophysiology, especially for the distal airways. On the one hand, animal and especially rodent models are widely used for preclinical studies, but their airways differ anatomically, histologically and physiologically from those of humans, and they do not fully recapitulate disease hallmarks [1]. To date, they have poorly predicted which drugs will help treat patients with chronic obstructive diseases. On the other hand, the lung has a complex structure leading to a peculiar stiffness and three-dimensional (3D) conformation that is difficult to model in vitro, highlighting the need to create adequate complex cellular models.

Various strategies have been proposed to model human bronchi in vitro [2]. Two-dimensional (2D) culture of bronchial epithelial cells (BECs) at the air–liquid interface (ALI) allows the reconstitution of a fully differentiated epithelium with functional cilia, with easy access to the apical and basal sides [3, 4], but it does not reproduce the 3D architecture essential for properly modelling the airways. 3D bioprinting is a promising method for building tissue with a relevant geometry [5, 6], especially with the recent development of suitable bioinks, such as those composed of alginate and decellularised tissue, enabling the printing of a tubular structure [7]. However, the length of the tubes obtained is limited, and 3D printing of the last branches of the bronchial tree remains below the resolution limit of current bioprinting techniques. The elegant microfluidic approach has been used to design so-called airway-on-a-chip [8, 9], providing unique advantages such as the ability to perfuse air, particles or nebulised drugs; to incorporate immune cells; and the potential application of cyclic stretching to reproduce breathing movements. While the microengineered scaffold restricts morphogenesis to a predefined geometry, it does not completely match the exact topology of a bronchus and does not offer a permissive and remodellable environment.

Other studies, starting from Benali et al. [10] to more recent ones [11–13], have favoured the third dimension with “bronchospheres”, in which BECs are embedded in a gel and self-organise into spheres with a liquid-filled lumen [14]. These self-renewing organoids exhibit ciliary activity and mucus production. In most of the bronchospheres, the apical side of the cells faces the inside of the organoid, limiting access to the lumen, although protocols to reverse airway organoid polarity have been developed [15, 16]. This closed and spherical architecture also impairs the removal of apical side-released dead cells inside and the introduction of an ALI. Finally, the formation of invasive tubular structures from lung organoids has been observed [17], but this process remains stochastic and rather limited.

Cellular Capsule Technology, adapted from a co-extrusion microfluidic technique, has been used to encapsulate cells and extracellular matrix proteins in a tubular scaffold made of alginate gel [18]. Cell-laden microfibres have been generated with various cell types, such as endothelial cells, epithelial cells, fibroblasts, skeletal muscle cells and nerve cells [18]. When co-encapsulated, endothelial and vascular smooth muscle cells self-organise into tubular structures that model vessels [19]. Although epithelial cell lines such as the human hepatoma G2 (HepG2) cell line, the mouse pancreatic β-cell (MIN6m9) cell line, the cervical cancer Henrietta Lacks (HeLa) cell line, the Madin-Darby Canine Kidney (MDCK) cell line and the EpH4-J3B1A mammary gland epithelial cell line have been previously successfully used [18, 20], the behaviour of primary BECs in this culture system has never been investigated.

Here, we aimed to construct a tubular organoid-like model of distal airways with human primary BECs. By using a tubular template made of alginate using Cellular Capsule Technology, we were able to drive the spontaneous self-organisation of lung epithelial cells to obtain a 3D “bronchioid” model with a physiologically relevant shape and size. Importantly, the tubular geometry allows lumen access to create an ALI.

MethodsPatients

Macroscopically normal lung resection material for the purification of basal epithelial cells was obtained from a cohort of patients requiring thoracic lobectomy surgery for nodules or cancer (pN0) (i.e. TUBE, sponsored by the University Hospital of Bordeaux, which includes its own local ethics committee (CHUBX 2020/54)). As per French law and MR004 regulations, patients received an information form, allowing them to refuse the use of their surgical samples for research. The research was conducted according to the principles of the World Medical Association Declaration of Helsinki. A total of 11 COPD patients with a clinical diagnosis of COPD according to the Global Initiative for Chronic Obstructive Lung Disease (GOLD) guidelines [21] and 18 non-COPD patients with normal lung function (i.e. forced expiratory volume in 1 s (FEV1) to forced vital capacity (FVC) ratio >0.70) were enrolled from the University Hospital of Bordeaux (table 1). Five other patients were also enrolled from the TUBE cohort for the comparison of bronchioids to classical ALI culture (supplementary table E1). All clinical data were collected at the Clinical Investigation Center (CIC1401) of the University Hospital of Bordeaux.

TABLE 1

Patient characteristics for bronchioid generation

Non-COPD patientsCOPD patientsp-valueSubjects (n)1811Age (years)67.3±8.667.0±5.20.92Sex>0.99 Male85 Female106Smoking history0.64 Current52 Former108 Nonsmoker31Pack years26.7±17.842.7±33.60.22PFT FEV1 (% predicted)90.9±29.6#73.2±22.8#0.08 FEV1/FVC ratio (%)75.7±10.5#61.8±13.8¶<0.0001Bronchial epithelial cell culture

Human basal BECs were derived from bronchial specimens as described previously [22]. Bronchial epithelial tissue was cultured in PneumaCult-Ex medium (STEMCELL Technologies, Vancouver, Canada) under a water-saturated 5% CO2 atmosphere at 37°C until basal BECs reached 80–90% confluence. For ALI culture, basal BECs were plated (1.105 cells per well) on tissue culture-treated nucleopore membranes (ref 3470, 6.5 mm diameter, 0.4 µm pore size; Transwell, Costar) in PneumaCult-ALI Medium (STEMCELL Technologies) applied at the basal side only to establish the ALI. The cells were maintained in culture for 9 or 21 days.

Microfluidic co-extrusion device and bronchioid formation

The design of a co-extrusion microfluidic device for bronchioid formation has been previously described for vessels [19]. Using this device with a 450 µm-diameter nozzle, we produced hollow alginate tubes filled with basal BECs. Cells were harvested using 0.05% trypsin-EDTA (Fisher Scientific, Illkirch, France) and resuspended in Matrigel (Corning, NY, USA) and DMEM (Sigma-Aldrich, Saint Quentin-Fallavier, France), both of which ranged from 33.5% to 35% (v/v). The cell loading of the core cell suspension (CCS) varied from 30% to 33% (v/v). Cells were then injected into the central zone isolated from a 2% (w/v) alginate solution (I3G80, Alliance Gums & Industries, Ile de France, France) by a 300 mM intermediate d-sorbitol solution (Sigma-Aldrich). The experiments were performed using the following injection flow rates: 2 mL·h−1 (alginate solution), 1 mL·h−1 (intermediate d-sorbitol solution) and 1 mL·h−1 (CCS). With these flow rates, the thickness of the alginate walls was ∼50 µm, resulting in a 350 µm-diameter cell tube (figure 1a). A 100 mM calcium bath at 37°C was used to trigger rapid gelation of the alginate gel while avoiding osmotic shock. Immediately after formation (within <10 min), the cell-laden tubes were transferred to Pneumacult Airway Organoid Seeding Medium (STEMCELL Technologies) at 37°C. After 24 h, the bronchioids were transferred to Pneumacult Airway Organoid Differentiation Medium (STEMCELL Technologies) to initiate differentiation. To inhibit Rho-associated protein kinase (ROCK), the airway organoid medium was supplemented with Y-27632 (Sigma-Aldrich) at a concentration of 10 µM. The media were changed daily.

FIGURE 1

The formation of an alginate tube filled with primary basal epithelial cells is called a “bronchioid”. a) Schematic of the device used to generate bronchioids. Created with BioRender.com. The three solutions were injected simultaneously by a computer-controlled pump inside a three-dimensional-printed device soaked in a 100 mM calcium bath. b) Left panel: brightfield images of bronchioids over time, with (+) or without (−) 10 µM Y-27632 (Y). The green and pink dotted lines indicate the external (alginate tube) and internal (epithelial tube) limits, respectively. Scale bars: 100 µm. Right panel: quantification of external and internal diameters (at least 10 measurements along the same tube for each condition, n=3 experiments). The medians are represented as horizontal lines. c) Immunoblots and analysis of myosin light chain (MLC) II and double-phospho-MLC (PP-MLC) (Thr18, Ser19) expression on day 2 in bronchioid tissue from three different donors (n=3). Total protein was used as a loading control. d) Upper panels: brightfield images of bronchioids cultured with 10 µM Y-27632. Scale bars: 100 µm. Lower panels: dot plots showing propidium iodide fluorescence (y-axis) versus annexin V-allophycocyanin fluorescence (x-axis) in cells dissociated from bronchioids and analysed by flow cytometry at the indicated time points. The percentages of propidium iodide− annexin V− cells are shown in pink. IS: intermediate d-sorbitol solution; AS: alginate solution.

Western blot

The bronchioids were assayed for myosin light chain (MLC) II and double-phospho-MLC (PP-MLC) (Thr18, Ser19) levels by Western blotting. Total protein extracts were obtained by transferring the bronchioids to Dulbecco's PBS (DPBS) (Sigma-Aldrich) for 30 min at 37°C with gentle agitation to dissolve the alginate. After centrifugation, the pellet was suspended in RIPA lysis buffer (Santa Cruz Biotechnology, Dallas, TX, USA) supplemented with protease inhibitor cocktail and calyculin A (Sigma-Aldrich). These extracts were loaded onto a 4–20% SDS‒PAGE gel (Bio-Rad, Hercules, CA, USA), transferred onto a nitrocellulose membrane and identified with a rabbit anti-target antibody (supplementary table E2). Horseradish peroxidase-coupled anti-rabbit secondary antibody was used for visualisation using a ChemiDoc imaging instrument (Bio-Rad). Protein expression was normalised to total protein loading (Stain-Free system, Bio-Rad).

Quantitative real-time PCR

Total RNA was extracted from bronchioids using a RNeasy Kit according to the manufacturer's instructions (Qiagen, Hilden, Germany). 500 ng of RNA were reverse-transcribed using the iScript Reverse Transcription Supermix (Bio-Rad). Complementary DNA (cDNA) samples were then analysed by quantitative PCR (qPCR) using Perfecta SYBR Green SuperMix (Quantabio, Beverly, MA, USA) through the CFX Connect Real-Time PCR Detection System (Bio-Rad). Primers (peptidylprolyl isomerase A (PPIA), forward 5′-CGGGTCCTGGCATCTTGT-3′ and reverse 5′-CAGTCTTGGCAGTGCAGATGA-3′; ribosomal protein L13 (RPL13), forward 5′-GGGAGCAAGGAAAGGGTCTTA-3′ and reverse 5′-CACCTGGACAATTCTCCGAGT-3′; glucuronidase β (GusB), forward 5′-CCATCTGGGTCTGGATCAAAA-3′ and reverse 5′-TGAAATCGGCAAAATTCCAAAT-3′; keratin 5 (KRT5), forward 5′-GCATCACCGTTCCTGGGTAA-3′ and reverse 5′-GACACACTTGACTGGCGAGA-3′; tumour protein 63 (TP63), forward 5′-CCTCAGGGAGCTGTTATCCG-3′ and reverse 5′-ATACTGGGCATGGCTGTTCC-3′; secretoglobin family 1A member 1 (SCGB1A1), forward 5′-TCCTCCACCATGAAACTCGC-3′ and reverse 5′-AGGAGGGTTTCGATGACACG-3′; dynein axonemal heavy chain 5 (DNAH5), forward 5′-AGAGGCCATTCGCAAACGTA-3′ and reverse 5′-CCCGGAAAATGGGCAAACTG-3′; anterior gradient 2 (AGR2), forward 5′-AAGGCAGGTGGGTGAGGAAATC-3′ and reverse 5′-GTCGAGAGTCCTTTGTGTCCTT-3′; mucin 5AC (MUC5AC), forward 5′-TACTCCACAGACTGCACCAACTG-3′ and reverse 5′-CGTGTATTGCTTCCCGTCAA-3′; epithelial cell adhesion molecule (EPCAM), forward 5′-CCATGTGCTGGTGTGTGAAC-3′ and reverse 5′-AACGCGTTGTGATCTCCTTCT-3′; and cadherin 1 (CDH1), forward 5′-TTACTGCCCCCAGAGGATGA-3′ and reverse 5′-TGCAACGTCGTTACGAGTCA-3′) were purchased from Sigma-Aldrich. mRNA expression was determined using the comparative 2-ΔΔCt method and normalised to the mRNA expression level of endogenous references (PPIA, RPL13 and GusB).

Generation of rhinovirus type 16-green fluorescent protein

The cDNA pA16-GFP plasmid encoding rhinovirus type 16 (RV16)-green fluorescent protein (GFP), which contains the rhinovirus sequence downstream of the T7 promoter, was a kind gift from Yury A. Bochkov (University of Wisconsin, Madison, WI, USA) [23, 24]. The reverse genetics technique developed by Lee et al. [24] was used to produce functional RV16-GFP virions. Briefly, for initial production, the pA16-GFP plasmid was linearised by incubation with the restriction enzyme EcI136II (EcoICRI, Thermo Scientific) at 10 U·μL−1. After purification using the QIAprep Spin Miniprep Kit (Qiagen), the DNA was transcribed in vitro into RNA by T7 polymerase and then transfected into HeLa cells using Lipofectamine (Invitrogen, Fisher Scientific). Finally, 24 h after transfection, the virus produced was purified and titrated by digital PCR after extraction and retrotranscription (Transcriptomic Platform, University of Bordeaux).

Infection with RV16-GFP

At days 4 to 5, an ∼5 cm-long bronchioid was transferred into a 32 mm-diameter DMEM-containing Petri dish. The bronchioid was cannulated by a 30 G needle (Microlance, BD Biosciences, Eysins, Switzerland) with a diameter of 300 μm, which is the order of magnitude of the inner diameter of the bronchioid. The needle was then connected to a 200 μL syringe containing DMEM (Fisher Scientific) with RV16-GFP at a concentration of 10 000 particles·μL−1 and supplemented with 0.1% eosin (Clinisciences, Nanterre, France). The solution was manually injected into the bronchioid lumen, with visual control via an eosin stain. The cultures were not washed after infection, but the opposite side of the bronchioid was “open”, so viral particles diffused out from the lumen into the medium. At 48 h post-infection, the supernatants were collected and frozen at −80°C for cytokine concentration measurements. The bronchioids were dissociated, and the cells were manually counted before being processed for fluorescence-activated cell sorting (FACS) analysis. The volume of injection corresponds to the volume of the lumen, which was estimated to be between 1 µL and 3 µL for 5 cm-long bronchioids containing between 12 000 and 40 000 cells, with measured lumen diameters of ∼200 µm. This allowed us to calculate a multiplicity of infection (MOI) between 0.5 and 1.5.

Cytokine production

The interferon-β (IFN-β), interferon-λ1/3 (IFN-λ1/3) and C-X-C motif chemokine ligand 8 (CXCL8) concentrations in the bronchioid medium were assayed 48 h post-infection by ELISA following the manufacturer's recommendations (Bio-Techne, Minneapolis, MN, USA). Values below the detection limit were considered zero.

Flow cytometry

Epithelial cells in ALI culture were washed in PBS containing 0.5% bovine serum albumin (BSA) (Sigma-Aldrich) and 2 mM EDTA (Invitrogen, Fisher Scientific) for 4 min at 37°C and harvested by incubation with 0.05% trypsin-200 mg·L−1 EDTA solution (Sigma-Aldrich) for 5 min, followed by the addition of trypsin-neutralising solution (Lonza, Basel, Switzerland), as described previously [25]. For bronchioids, alginate dissolution was achieved by incubating the epithelial cell-laden tubes in DPBS for 20 min under gentle agitation. The cells were then mechanically and enzymatically dissociated by pipetting and by incubating the suspension with 0.05% trypsin-200 mg·L−1 EDTA (Sigma-Aldrich) for 5 min. After centrifugation, the cell suspension was filtered through a 70 μm cell strainer (Corning) and immediately evaluated by flow cytometry to test for rhinovirus infection. Alternatively, annexin V-propidium iodide (Fisher Scientific) or DAPI (5 µg·mL−1, Sigma-Aldrich) was used for the detection of apoptotic and dead cells, respectively, according to the manufacturer's instructions. For differentiation, cells were stained with anti-EpCAM/peridinin chlorophyll protein complex/cyanine5.5 or its isotype control (supplementary table E2), fixed, and permeabilised using the Foxp3 buffer set (Fisher Scientific) or the IntraPrep Permeabilization Reagent Kit (Beckman Coulter, Villepinte, France). The cell pellet was then stained with anti-pan cytokeratin-fluorescein isothiocyanate (FITC), anti-cytokeratin 5-Alexa Fluor 488, anti-SCGB1A1-FITC, anti-acetylated-tubulin-phycoerythrin, anti-MUC5AC-Alexa Fluor 647 or their respective isotype controls (supplementary table E2) for 1 h. FACS data were acquired using a Canto II 4-Blue 2-Violet 2-Red laser configuration (BD Biosciences). Flow cytometry analysis was performed using Diva 8 (BD Biosciences).

Immunostaining

The bronchioids were fixed with 4% paraformaldehyde (Sigma-Aldrich) diluted in DPBS supplemented with calcium (Sigma-Aldrich) for 30 min at room temperature and immunostained as described previously [26]. Briefly, after fixation, permeabilisation and blocking steps were successively performed by incubating the bronchioids with cold DPBS supplemented with calcium containing 0.1% (v/v) Tween (P1379, Sigma-Aldrich) and with cold organoid washing buffer (OWB) composed of PBS containing 2 g·L−1 BSA and 0.1% (v/v) Triton X-100 (Biosolve, Dieuze, France) at 4°C for 15 min. The bronchioids were then incubated with primary antibodies diluted in OWB at 4°C overnight with mild rocking (supplementary table E2). The bronchioids were washed three times with OWB for 2 h each and then stained with secondary antibodies, DAPI and phalloidin diluted in OWB at 4°C overnight, with mild rocking (supplementary table E2), before three supplementary washes with OWB for 2 h each. Next, the bronchioids were transferred to a µ-Slide 4-Well chamber slide (Ibidi, Gräfelfing, Germany) and resuspended in a fructose–glycerol clearing solution composed of 60% (v/v) glycerol and 2.5 M fructose in distilled H2O or in DPBS supplemented with calcium.

Imaging

Fixed bronchioids were imaged using a confocal Leica SP8 WLL2 microscope on an inverted DMI6000 stand (Leica Microsystems, Mannheim, Germany) using a ×20 multi-immersion objective (0.75 numerical aperture HC Plan Apo CS2). White light laser 2 (WLL2) was used for the excitation of Alexa Fluor 647, 568 and 488 in the ranges of 650–750 nm, 590–640 nm and 500–550 nm, respectively. The confocal microscope was also equipped with a diode laser at 405 nm. A range of 410–480 nm was used to excite DAPI. The scanning was performed using a conventional scanner (400 Hz). The image size was 1024×1024 pixels, with a pixel size of 0.568 µm and a voxel depth of 1.040 µm. Detection was achieved with two internal photomultiplier tubes and two internal hybrid detectors. Fluorescence images were acquired with LAS X software (Leica Microsystems). Using Fiji software, the background value was subtracted from the original image, and the image contrast in the background-corrected image was adjusted. The images were then processed with Imaris software (version 9; Oxford Instruments, Abingdon, UK) for segmentation and 3D visualisation. The “clipping plane” function was used to truncate the image along a given plane in 3D space.

Live imaging to analyse ciliary beat frequency (CBF) was performed on an inverted Leica DMi8 microscope (Leica Microsystems) equipped with a Flash 4.0 sCMOS camera (Hamamatsu Photonics, Hamamatsu, Japan). The illumination system used was a Cool LED PE-4000 (CoolLED, Andover, MA, USA). The objective used was a HCX PL Fluotar L ×40 dry 0.6 numerical aperture PH2. A 37°C/5% CO2 atmosphere was created with an incubator box and a gas heating system (Pecon GmbH, Erbach, Germany). Brightfield images were acquired with MetaMorph software (Molecular Devices, Sunnyvale, CA, USA) for 2 s at a frame rate of 1000 images·s−1 in 9–10 different fields of the bronchioid region. The image size was 512×32 pixels, with a pixel size of 0.1625 µm. The number of images used for the analysis was 1024. This resulted in a temporal resolution of 0.97 Hz (1000 images·s−1 divided by 1024 images). If needed, recursive alignment of the stack of images was performed using the plugin “StackReg” of the Fiji software. Regions of interest were selected on the image generated by a maximum intensity projection using the plugin Stardist 2D of Fiji software (supplementary figure E6a). The mean grey intensity was measured in each of these regions over time, and fast Fourier transformation analysis was applied using Microsoft Excel to retrieve the different frequencies. These frequencies exhibit a discrete pattern, with an interval of 0.97 Hz between each of the possible values (supplementary figure E6b). The main frequency was taken as the CBF (supplementary figure E6b). Next, in each field, the median CBF was determined from the CBF values of the regions of interest (supplementary figure E6c). Nine to ten fields were evaluated per donor.

Perfusion

A 3–4 cm-long bronchioid was cannulated with a 30 G needle (Microlance, BD Biosciences) with a diameter of 300 μm. The needle was then connected to a microfluidic system (Fluigent, Le Kremlin-Bicêtre, France) composed of a microfluidic low-pressure generator (FLPG Plus), a microfluidic vacuum pressure-based controller (LineUp Push-Pull), a noninvasive flow sensor for liquid or air perfusion and a reservoir with a pressurisation cap. Air was injected into the lumen at one opened end of the bronchioid at a continuous flow rate of 0.2 mL·min−1. Air was released in the medium from the opposite tube end, generating visible air bubbles. Imaging was performed using an Olympus CKX53 microscope with an Olympus DP23 camera with objectives of ×10 or ×20. Perfusion by air was continued for 4 h. After the infusion was stopped, calcein (Invitrogen) was added to the culture medium at a final concentration of 1 μM, and the cells were incubated for 15 min at 37°C.

Single-cell RNA-sequencing data acquisition and analysis

A single-cell suspension was obtained from bronchioids at day 21 as described for flow cytometry analysis. Cell number and viability were determined by a Trypan blue exclusion assay before proceeding with the 10x Genomics (Pleasanton, CA, USA) protocol for 3′ transcript capture and single-cell library preparation. Approximately 18 000 cells were loaded on a 10x Chromium Next GEM Chip G. GEM generation, reverse transcription (RT), cleanup, cDNA amplification, fragmentation, end repair and A-tail preparation, and sample index tagging were performed using the Chromium Next GEM Single Cell 3ʹ GEM, Library & Gel Bead Kit v3.1 as per the manufacturer's instructions.

Quantification of the libraries was performed using a QubitTM 2.0 Fluorometer (Thermo Fisher), and the quality of the libraries was checked with a 2100 Bioanalyzer with a high-sensitivity DNA kit (Agilent Technologies, Les Ulis, France). Libraries were sequenced using a P3 flow cell on a NextSeq2000 in 100 bp paired-end mode (Illumina, Évry, France) to read depths of 111 543 mean reads·cell−1 (Patient 1) and 103 719 mean reads·cell−1 (Patient 2).

The raw sequencing data were aligned to the GRCh38 human genome assembly (Ensembl98), and the transcripts were quantified using the CellRanger software suite (v7.1.0; 10x Genomics). Subsequently, the count matrices were loaded into R (v4.2.3; www.r-project.org), and downstream analyses were conducted using the Seurat R package (v5.0.1, [27]). Cells were filtered based on mitochondrial gene counts (<10%) and the number of genes detected per cell (>2000). Gene counts were log-normalised (NormalizeData, default parameters), and highly variable genes (HVGs) were selected (FindVariableFeatures, default parameters). Batch correction was applied using the Seurat reciprocal principal component analysis (RPCA) workflow with “Patient” as a covariate. A nearest neighbour graph was constructed (FindNeighbors, dims=1:10), and a 2D embedding was computed with the uniform manifold approximation and projection (UMAP) algorithm (RunUMAP, dims=1:10). Unsupervised clustering was performed using the Louvain algorithm (FindClusters, resolution=0.5). Cell types were annotated based on the expression of known marker genes. Ionocytes and hillock-like cells did not cluster separately in the unsupervised clustering due to their small cell number and were selected based on UMAP coordinates, because they exhibited clearly distinct locations on the UMAP representation. Differential gene expression analysis was performed using the Wilcoxon rank-sum test implemented in the Seurat package (FindAllMarkers, only.pos=TRUE, min.pct=0.25), and the differential gene expression results can be found in supplementary table E3. Additionally, cell types were predicted with the CellTypist [28] Python package (v2.1.0) using Human_Lung_atlas.pkl as a model, with majority voting enabled. CellTypist-predicted annotations were grouped into broader categories (e.g. “Basal resting” and “Suprabasal” were categorised as “Basal”) to facilitate comparability with the manual annotations. The adjusted Rand index was computed to assess concordance between manual annotations and CellTypist predictions.

Primary lung tissue reference construction

Primary lung tissue single-cell RNA-sequencing (scRNA-seq) data along with associated cell type annotations, clinical metadata and technical metadata were obtained from the full Human Lung Cell Atlas (HLCA) as disclosed in March 2024 [29]. Cells annotated as airway epithelium were extracted from the dataset. To ensure the technical coherence of the primary tissue reference with the bronchioid data produced in the present study, only cells analysed with 10x Genomics 3′ v2 or v3 assays and derived from fresh tissue were retained for reference construction. Gene counts were log-normalised, and 3000 HVGs were selected and scaled for subsequent data integration using the Seurat RPCA workflow with “dataset” as a covariate. A nearest neighbour graph was constructed, and a 2D UMAP embedding was computed based on 30 RPCA dimensions.

Neighbourhood correlation analysis

Neighbourhood graph correlation-based similarity analysis [30] was employed to compare cell states in the bronchioid model to those in primary lung tissue from distinct anatomical locations based on the common coordinate framework (CCF) established in the HLCA [31]. Cells from the primary lung tissue reference (constructed as described above) were categorised into four CCF categories: CCF 0, corresponding to the inferior turbinate; CCF 0.36–0.5, corresponding to the trachea and the main bronchi; CCF 0.72–0.81, corresponding to lobular bronchi and distal lobular airways; and CCF 0.97, corresponding to the parenchyma [31]. For the bronchioid data as well as the primary lung tissue reference CCF categories, independent neighbourhood graphs based on transcriptional similarity were constructed with the miloR package (v1.6.0) using k=30 nearest neighbours and 30 RPCA dimensions of the respective data source [32]. The R package scrabbitr (v0.1.0) was used to calculate the mean expression profiles within each neighbourhood and to compute the Pearson correlations of expression of the 3000 HVGs from the primary lung tissue reference, excluding mitochondrial genes, between all pairs of bronchioid and CCF category neighbourhoods. The maximum correlation between each bronchioid neighbourhood and a neighbourhood in the CCF-stratified primary tissue reference was considered a measure of similarity.

Data and code availability

The scRNA-seq data generated in the current study are available from Zenodo under record number 10.5281/zenodo.10834205. The R code used for scRNA-seq data analysis has been deposited on GitHub: https://github.com/saliba-lab/Bronchioid_scRNAseq.

ResultsFormation of the tubular bronchial epithelial model from patient-derived basal cells

Basal BEC cultures were obtained from 29 patients (table 1, n=18 non-COPD patients, n=11 COPD patients). Using Cellular Capsule Technology, we encapsulated BECs in the presence of Matrigel into tubular structures made of alginate (figure 1a). Notably, cells were injected at a concentration sufficient to allow immediate coverage of the entire alginate wall as soon as the tube was generated. Of the 29 samples, only one culture from a non-COPD patient (3%) failed to generate bronchioids due to a technical problem (supplementary figure E1). Nine bronchioids (n=4 non-COPD patients, n=5 COPD patients) were fully used before 9 days for early characterisation (supplementary figure E1). Among the 19 bronchioids maintained for at least 11 days in culture, 6 (32%) were from patients with COPD, which is a similar proportion to the enrolled patient population. Taken together, these findings show that bronchioids are generated with a high success rate and that they can be cultured for several weeks in a similar way whether derived from COPD or non-COPD patients.

Within the first 24 h, epithelial cells formed a monolayer that covered the entire inner part of the tube (figure 1b). However, as previously shown with epithelial cell lines [20], the bronchioid collapsed 1–2 days after encapsulation (figure 1b). Inhibition of actomyosin contractibility using the ROCK inhibitor Y-27632 (10 µM) prevented epithelial detachment (figure 1b) while decreasing the level of PP-MLC (Thr18, Ser19) (figure 1c), suggesting that cell contractility was the main driving force for shrinkage.

The dissolution of the alginate gel was successfully achieved by calcium and magnesium cation depletion for 30 min, allowing the generation of a cell suspension. Analysis of this suspension by flow cytometry demonstrated that the tubular structures remained stable and viable for at least 9 days (figure 1d). Thus, in subsequent experiments, bronchioids were continuously treated with Y27632 to avoid collapse. We confirmed the epithelial nature of the cells by flow cytometry and showed that at day 1 the cells were almost all positive for EpCAM and cytokeratin, similar to differentiated BECs at the ALI (supplementary table E1, figure 2a). The mRNA levels of CDH1 (encoding E-cadherin) and EPCAM, which were measured by quantitative RT-PCR (RT-qPCR), remained stable over time (figure 2b). Confocal imaging and 3D reconstruction demonstrated the organisation of bronchial cells into a tubular structure with a typical epithelial morphology and a visible lumen, reminiscent of the structure of distal bronchi (figure 2c). The presence of zonula occludens-1 (ZO-1) and actin fibres at cell‒cell contacts confirmed the presence of tight junctions (figure 2d).

FIGURE 2

Characterisation of the epithelial nature of cells in the bronchioid model. a) Dot plots representing epithelial cell adhesion molecule (EpCAM)-peridinin chlorophyll protein complex (PerCP)-cyanine 5.5 fluorescence (y-axis) versus pan cytokeratin-fluorescein isothiocyanate fluorescence (x-axis) of cells dissociated from two-dimensional (2D) air–liquid interface (ALI) culture and bronchioids at day (D) 1. The gating strategy is shown in the left panels. The percentages of EpCAM+ cytokeratin+ cells are shown in pink. b) Expression of the genes CDH1 and EPCAM in bronchioids over time and in cells dissociated from 2D culture at day 9 after ALI introduction. The bronchioid and ALI culture samples were obtained from six and three different donors, respectively, and data are presented as mean±sd. Gene expression was normalised to that of the housekeeping genes PPIA, RPL13 and GusB and expressed relative to that of basal epithelial cells in 2D submerged culture. c, d) Longitudinal and transverse sections and three-dimensional (3D) views of 3D reconstructions obtained from Z-stack confocal images of a 2-day-old bronchioid stained for F-actin (phalloidin, magenta) and nuclei (DAPI, white) (c) and F-actin (magenta), zonula occludens-1 (ZO-1) (cyan) and nuclei (white) (d). Scale bars: 100 µm and 50 µm in the magnified lower panel.

Cell differentiation and organisation in the bronchioid model

RT-qPCR also revealed that the expression of the airway basal cell gene KRT5 remained stable, and that of TP63 progressively declined over time, whereas the expression of the club cell secretory marker SCGB1A1 was highly upregulated (figure 3a). The expression of the ciliated cell markers forkhead box J1 (FOXJ1) and DNAH5, the transcription factor AGR2 and the goblet cell marker MUC5AC progressively increased in culture (figure 3a). At day 9, the expression of these genes was comparable with BECs at the ALI for the same time, with the exception of EPCAM, SCGB1A1 and MUC5AC, whose expression was greater, and AGR2, whose expression was reduced under the ALI condition (figures 2b and 3a). The cellular composition and viability were very similar in three different portions sampled along the tube (supplementary figure E2), indicating a homogenous cellular structure throughout the bronchioid. Using confocal microscopy, we found that ciliated cells were located in the internal layer lining the bronchioid lumen at day 11 (figure 3b and supplementary video E1). Most goblet cells also faced the lumen, with some MUC5AC+ cells more externally localised (figure 3b). Positive staining for MUC5AC was also found inside the lumen, suggesting mucus secretion. Although all intermediate states of differentiation are probably not detected by simple staining, especially during the last 2 weeks in culture, flow cytometry confirmed the progressive loss of basal cells over time, with a concomitant increase in club, goblet and ciliated cells in non-COPD-derived bronchioids until day 21 (figure 3c, d). In comparison, COPD-derived bronchioids were characterised by a prominent mucous phenotype, with an already high proportion of goblet cells at day 9 (supplementary figure E3a, b). This phenotype persisted until day 21, at the expense of ciliated cells, whose percentage remained low. This resulted at day 21 in a five-fold increase in the percentage of goblet cells and an eight-fold decrease in the percentage of ciliated cells in COPD bronchioids compared with non-COPD ones (supplementary figure E3c).

FIGURE 3

Differentiation induction of primary bronchial epithelial cells into bronchioids. a) Expression of the genes TP63, KRT5, SCGB1A1, FOXJ1, DNAH5, AGR2 and MUC5AC in bronchioids over time and in cells dissociated from two-dimensional (2D) culture at day (D) 9 after air–liquid interface (ALI) introduction. The bronchioid and ALI culture samples were obtained from six and three different donors, respectively, and data are presented as mean±sd. Gene expression was normalised to that of the housekeeping genes PPIA, RPL13 and GusB and expressed relative to that of basal epithelial cells in 2D submerged culture. b) Longitudinal and transverse sections and three-dimensional (3D) views of 3D reconstructions obtained from Z-stack confocal images of 11-day-old bronchioids stained for acetylated α-tubulin (magenta), mucin-5AC (MUC5AC) (cyan) and nuclei (white). Arrows: ciliated cells; arrowheads: goblet cells. Scale bars: 50 µm. c) Histograms of representative cell counts (y-axis) versus DAPI fluorescence (x-axis) of cells dissociated from bronchioids at the indicated time points. Dot plots represent cytokeratin-5/secretoglobin family 1A member 1 (SCGB1A1)/MUC5AC/acetylated tubulin fluorescence (y-axis) versus forward side scatter (FSC) (x-axis). DAPI− cells, cytokeratin 5+, SCGB1A1+, MUC5AC+ and acetylated tubulin+ cells are shown in dark blue, green, orange, pink and purple, respectively. d) Percentages of DAPI−, cytokeratin 5+, SCGB1A1+, MUC5AC+ and acetylated tubulin+ cells over time. The samples are from four different non-COPD donors and data are presented as mean±sd.

Single-cell characterisation of non-COPD bronchioids

To further characterise the bronchioid model, scRNA-seq was conducted on bronchioids from two non-COPD patients after 21 days of culture (Methods). In total, 5087 high-quality single cells (Patient 1=2860 cells, Patient 2=2227 cells) were analysed (Methods and supplementary figure E4a) and categorised into nine distinct populations without noticeable differences between the donors (figure 4a and supplementary figure E4b). These populations could be annotated based on established markers, such as basal (KRT5, TP63), secretory (SCGB1A1, SCGB3A1, MUC5B), goblet (MUC5AC, CEACAM5), deuterosomal (FOXN4, CDC20B), ciliated (FOXJ1, TPPP3), ionocyte (FOXI1, CFTR) and hillock-like (KRT13, KRT14) cells (figure 4a, b and supplementary table E3). Additionally, populations of differentiating basal cells (KRT5 and AGR2) and proliferating cells (MKI67) were identified. To validate the manual cell type annotations, CellTypist [28] with the HLCA [31] model were used for cell type prediction and demonstrated a high level of agreement (adjusted Rand index=0.691, supplementary figure E4c). Notably, all cell types, except for hillock-like cells (14 cells, 0.24%), were consistently observed in bronchioids derived from both patients, even if the proportions were not strictly identical (figures 4c and figure E4b). Although differentiated secretory and multiciliated lineage cells, as well as rare ionocytes and hillock-like cells, were revealed in our bronchioid scRNA-seq data, a distinct population of SCGB3A2+SFTPB+ terminal airway-enriched secretory cells [33, 34] was not identified as a cluster despite the co-expression of marker genes in a few cells (supplementary figure E4d, e). Finally, to evaluate whether the bronchioid model accurately recapitulates small airways specifically at the transcriptional level, we used data from the HLCA and the CCF introduced therein [31]. A primary tissue reference across four distinct anatomical regions based on CCF was assembled, comprising cells from eight datasets marked broadly as airway epithelium in the HLCA (Methods, supplementary figure E5a–c). Bronchioid model cell states were compared to this primary tissue reference by performing a neighbourhood graph correlation-based similarity analysis (Methods) [30]. Overall, bronchioid cell states were more similar to cell states from the airways (CCF 0.32–0.5 and CCF 0.72–0.81) than to those from the nose (CCF 0) or parenchyma (CCF 0.97). Comparing the upper and lower airways, the similarity analysis results indicated that cell states from the distal airways (CCF 0.72–0.81) tended to more closely resemble bronchioid cell states (figure 4d, e). In summary, these data are consistent with the induction of distal airway fate in bronchioids, with identified cell states recapitulating distal airway epithelial cells at a global transcriptional level with high fidelity and differentiated cells localised to their corresponding in vivo positions.

FIGURE 4

Single-cell transcriptome profiling in the bronchioid model. a) Uniform Manifold Approximation and Projection (UMAP) representation of single-cell RNA-sequencing (scRNA-seq) transcriptomic data showing the different cell types detected in 21-day-old bronchioids from two different non-COPD donors (Patients 1 and 2). b) Dot plot showing scaled mean expression (colour) and percentage of expressing cells (dot size) of selected marker genes in the indicated cell groups. c) Relative abundance of cell types identified in bronchioids from Patients 1 and 2. d) Bronchioid neighbourhoods, positioned with respect to the UMAP embedding of each index cell, are coloured by the maximum correlation value across primary lung tissue neighbourhoods from four distinct anatomical locations based on the common coordinate framework (CCF) established in the HLCA [31]. e) Box plot depicting neighbourhood similarities of bronchioid neighbourhoods with primary lung tissue neighbourhoods in the respective CCF category over major cell types identified in the bronchioid scRNA-seq data. The centre line indicates the median, the box limits indicate the upper and lower quantiles, and the whiskers indicate 1.5× interquartile range. In d and e, the bronchioid neighbourhoods were constructed from cells integrated across the two samples.

Functional analysis of ciliary beating in non-COPD- and COPD-derived bronchioids

Using high-speed video microscopy [35], ciliary beating was observed in the bronchioid model (supplementary video E2). For bronchioids generated from four different patients, we determined the median CBF across different fields of acquisition (figure 5a, b). Overall, the median CBF values ranged from 6 Hz to 27 Hz (figure 5c), which is within the range of in vivo reference CBF values [36] and is also very similar to those found in bronchial epithelial strips [37], ALI cultures [38] and airway organoids [39]. Bronchioids derived from non-COPD patients were characterised by a physiological, biphasic spectrum of beating frequencies (figure 5c), as described for control ALI culture [40]. In contrast, COPD-derived bronchioids exhibited a loss of the biphasic spectrum and a shift towards lower frequencies (figure 5c).

FIGURE 5

Fast Fourier transform (FFT) analysis of ciliary movement using high-speed video microscopy analysis. a) Schematic workflow for measuring ciliary beat frequency (CBF). Representative brightfield image corresponding to a part of a bronchioid at day 20. The coloured areas correspond to the regions of interest (ROIs) determined by the plugin Stardist 2D of the Fiji software on the image generated by a maximum intensity projection over time. Scale bar: 5 µm. b) Representative measurements of the mean intensity in ROIs over time. The main frequency was determined using a FFT in each ROI. c) Violin plots showing the CBF of bronchioids derived from four different patients (non-COPD Patients 3 and 4, COPD Patients 5 and 6). For each patient, individual values represent median CBFs, determined in 9–10 fields.

Rhinovirus infection and epithelial response in the bronchioid model

We also used the model to test for rhinovirus infection, a major trigger of asthma and COPD exacerbation [41, 42]. To track viral infection, we injected a recombinant RV16 that accommodates GFP expression into the lumen of the bronchioid (figure 6a) [23]. Considering the calculated volume of the injection, estimated between 1 µL and 3 µL for a 5 cm-long bronchioid containing between 12 000 and 40 000 cells, with measured lumen diameters of ∼200 µm, the MOI was ∼0.5 and 1.5. 48 h post-infection, 95% and 81% of the cells were GFP+ in bronchioids generated from Patients 7 and 8, respectively (figure 6b, c), suggesting infection and viral RNA replication in the vast majority of cells, in agreement with the calculated MOI. The infection was less efficient (45% of GFP+ cells) in bronchioids generated from Patient 9 (figure 6b, c). Next, we examined the effects of RV16 on cytokine secretion, including type I (IFN-β) and type III (IFN-λ1/3) interferons and the chemokine CXCL8. RV16-GFP infection increased the levels of secreted IFN-β, IFN-λ1/3 and CXCL8 in comparison with the noninfected bronchioid condition for Patients 7 and 8 (figure 6d). No induction of interferon production or a modest increase in the CXCL8 level was observed in the bronchioids generated from Patient 9, which is in agreement with the decreased infection efficiency (figure 6d).

FIGURE 6

Rhinovirus infection in the bronchioid model. a) Experimental design for testing rhinovirus type 16 (RV16)-green fluorescent protein (GFP) infection. b) Histograms sh