The Transcriptome Landscape of the In Vitro Human Airway Epithelium Response to SARS-CoV-2

Airway–liquid interface cultures of primary epithelial cells and of induced pluripotent stem-cell-derived airway epithelial cells (ALI and iALI, respectively) are physiologically relevant models for respiratory virus infection studies because they can mimic the in vivo human bronchial epithelium. Here, we investigated gene expression profiles in human airway cultures (ALI and iALI models), infected or not with severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2), using our own and publicly available bulk and single-cell transcriptome datasets. SARS-CoV-2 infection significantly increased the expression of interferon-stimulated genes (IFI44, IFIT1, IFIT3, IFI35, IRF9, MX1, OAS1, OAS3 and ISG15) and inflammatory genes (NFKBIA, CSF1, FOSL1, IL32 and CXCL10) by day 4 post-infection, indicating activation of the interferon and immune responses to the virus. Extracellular matrix genes (ITGB6, ITGB1 and GJA1) were also altered in infected cells. Single-cell RNA sequencing data revealed that SARS-CoV-2 infection damaged the respiratory epithelium, particularly mature ciliated cells. The expression of genes encoding intercellular communication and adhesion proteins was also deregulated, suggesting a mechanism to promote shedding of infected epithelial cells. These data demonstrate that ALI/iALI models help to explain the airway epithelium response to SARS-CoV-2 infection and are a key tool for developing COVID-19 treatments.


Introduction
The rapid spread of severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) in humans has posed a serious global health threat. Coronaviruses are part of a large family of viruses that cause illnesses ranging from common colds to severe respiratory diseases, including COVID-19 caused by SARS-CoV-2. SARS-CoV-2 is a single positive-stranded RNA enveloped virus that is capable of replicating in epithelial cells [1]. Due to the longand short-term effects of COVID-19 on human health and the need to limit the emergence of novel virus variants, significant efforts have been dedicated to help understand the viral infection mechanisms and to develop antiviral drugs using physiologically relevant 2 of 19 in vitro culture models that mimic in vivo phenotypes. Human airway epithelial cells in culture are traditionally used for modeling respiratory diseases [2]. These cells can be obtained from lung tissue biopsies and are cultured as primary airway epithelial cells in air-liquid interface (ALI) systems that support epithelial cell differentiation and mimic key aspects of the mucosal epithelium [3]. They can also be derived by differentiation of induced pluripotent stem cells (iPSCs) in ALI systems (i.e., iALI) [4][5][6][7]. In both ALI and iALI systems, epithelial cells are cultured on a permeable membrane with the medium in the basal chamber and the epithelium exposed to air at the apical side of the membrane. In this system, cells are in contact with air and can be induced to differentiate into a functional pseudostratified epithelium. Several studies confirmed that ALI culture transcriptomic profiles are very similar to those of the in vivo airway epithelium obtained via bronchial brushing or biopsy [8,9]. ALI and iALI models can be infected by viruses and have been used to model various mechanisms of viral pathogenesis [10,11]. SARS-CoV-2 can replicate in both models [1,12], and interferon can inhibit virus infection. This demonstrated interferon therapeutic potential for COVID-19 treatments and the usefulness of these models as a high-throughput screening tool [13,14]. Additionally, ALI and iALI cultures from patients with respiratory diseases (e.g., chronic obstructive pulmonary disease) recapitulate some in vivo disease characteristics and are used to assess the impact of smoke exposure on viral infection [10,15,16]. Therefore, ALI and iALI cultures could help us to understand SARS-CoV-2 effects in the bronchial epithelium by analyzing the transcriptomic changes upon infection. ALI models are very helpful in identifying the key initiating steps of viral injury and innate epithelial cell defense that may or may not lead to cell infection and replication. Increasing the number of models and donors will improve the robustness of the identified pathways by reducing the inter-individual heterogeneity in viral susceptibility.
Various omic-based studies, including in vivo (human samples) and in vitro (model systems) transcriptome profiling studies (bulk RNA sequencing), have highlighted the molecular changes induced by SARS-CoV-2 infection [17][18][19]. The recent advent of singlecell RNA sequencing (scRNA-seq) provides a precious approach to carefully analyze gene expression and cell composition. For instance, scRNA-seq has been used to identify the cells in the human respiratory system with the highest expression of transmembrane receptors for SARS-CoV-2 [20] and to show that in ALI cultures of nasal epithelial cells, ciliated and goblet/secretory cells express progressive SARS-CoV-2 entry factors. Moreover, in infected human bronchial epithelial-derived ALI cultures, scRNA-seq [21] revealed that ciliated cells are a major SARS-CoV-2 target.
In the present study, we analyzed bulk and single-cell RNA-seq datasets to provide a detailed picture of the gene expression changes in ALI and iALI models following SARS-CoV-2 infection. This analysis highlighted the molecular mechanisms involved in the induction of the hyper-inflammatory state and innate immune response, including interferon signaling, chemokines and extracellular matrix (ECM) and identified potential regulatory microRNAs (miRNAs) for therapeutic interventions.

Cellular Landscapes of Non-Infected Human Lung Epithelium Models
Non-infected HBEC in ALI culture and iPSC-derived epithelial-like cells in iALI culture were analyzed via single-cell RNA sequencing (scRNA-seq) to generate a comprehensive portrait of the cell populations present in these models. Single-cell transcriptome data were also obtained from the lung epithelium collected through bronchial biopsy or brushing ( Figure 1A). Analysis of these scRNA-seq data from these three models (biopsy/brushing bronchial cells, ALI and iALI cultures) ( Figure 1B) showed that most epithelial cell markers, such as SCGB1A1/CC10 (secretoglobin family 1A member 1; secretory cells), TP63 (tumor protein p63; basal cells) and FOXJ1 (forkhead box J1; multi-ciliated cells), were present in the three sample types. This suggested that our non-infected ALI and iALI models recapitulate the human airway and can be used for infection studies. Comparison of the single-cell landscapes indicated that unlike the biopsy/brushing and ALI models, the iALI model expressed ASCL1, a marker of neuroendocrine cells, but not FOXI1 (ionocyte marker), and all models express POU2F3 (tuft cell marker). In addition, immune cells were an important cell fraction in biopsy/brushing samples ( Figure 1C) but were absent in the ALI and iALI models. This major difference may help to (i) understand the epithelial-intrinsic innate inflammatory response to SARS-CoV-2 infection; (ii) identify the key epithelial target cells; and (iii) compare the molecular mechanisms triggered in the absence (iALI and ALI) and presence (biopsy-brushing samples) of immune cells.

Cell-Specific Expression of Single-Stranded RNA Virus Receptors in Non-Infected Human Lung Epithelium Models
Expression of virus receptors is a key factor for virus infection and transmission of all viruses. Analysis of scRNA-seq data from the three non-infected lung epithelium models showed that CDHR3 (encoding cadherin-related family member 3, the entry receptor for rhinovirus C) [22] was specifically expressed in ciliated cells and CD55 (Coxsackievirus B1 receptor) [22] in secretory cells ( Figure 1D). ACE2 (SARS-CoV-2 receptor) was expressed by ciliated, secretory and basal cells. Moreover, in biopsy/brushing samples, ACE2 + cells were overrepresented among all epithelial cell types compared with immune cells, including macrophage, endothelial and dendritic cells ( Figure 1D). This is in line with a recent study reporting the absence of viral transcripts in bronchoalveolar fluid and peripheral blood mononuclear cell samples collected from patients with COVID-19 [23]. Collectively, scRNAseq data analysis of the three cell culture models showed expression of CDHR3, CD55 and ACE2 receptors in ciliated, secretory and basal cells, suggesting the possibility of modeling ssRNA virus infection, including SARS-CoV-2, in the airway epithelium in vitro. (A) Schematic representation of the scRNA-seq experimental workflow: airway epithelium sources (biopsy/brushing-derived cells, ALI culture of primary bronchial epithelial cells and iALI culture of iPSC-derived epithelial-like cells), generation of scRNA-seq libraries and sequencing, computational analysis to identify cell types. The brushing/biopsy epithelium scRNA-seq data were from [24]. The ALI and iALI scRNA-seq data were generated in our laboratory (see Methods); (B) contribution of each cell type in the three models. UMAP and tSNE were used to show the contribution of SCGB1A1+, TP63+, FOXJ1+, FOXI1, POU2F3+ and ASCL1+ cells (dark gray). (C) UMAP projection of the different cell types. The feature plots display the subset of airway epithelial cells obtained from biopsy/brushing of healthy donors [24]. (D) Maps of the expression of ssRNA virus receptors (rhinovirus C receptor (CDHR3), Coxsackievirus B1 receptor (CD55) and SARS-CoV-2 receptor (ACE2)) in single cells from the three models. UMAP and tSNE were used to show the cell types that express the ACE2, CDHR3 and CD55 receptors.

Cell-Specific Expression of Single-Stranded RNA Virus Receptors in Non-Infected Human Lung Epithelium Models
Expression of virus receptors is a key factor for virus infection and transmission of all viruses. Analysis of scRNA-seq data from the three non-infected lung epithelium models showed that CDHR3 (encoding cadherin-related family member 3, the entry receptor for rhinovirus C) [22] was specifically expressed in ciliated cells and CD55 (Coxsackievirus B1 receptor) [22] in secretory cells ( Figure 1D). ACE2 (SARS-CoV-2 receptor) was Figure 1. Distribution of the different airway epithelial cell types in the three lung tissue models. (A) Schematic representation of the scRNA-seq experimental workflow: airway epithelium sources (biopsy/brushing-derived cells, ALI culture of primary bronchial epithelial cells and iALI culture of iPSC-derived epithelial-like cells), generation of scRNA-seq libraries and sequencing, computational analysis to identify cell types. The brushing/biopsy epithelium scRNA-seq data were from [24]. The ALI and iALI scRNA-seq data were generated in our laboratory (see Methods); (B) contribution of each cell type in the three models. UMAP and tSNE were used to show the contribution of SCGB1A1+, TP63+, FOXJ1+, FOXI1, POU2F3+ and ASCL1+ cells (dark gray). (C) UMAP projection of the different cell types. The feature plots display the subset of airway epithelial cells obtained from biopsy/brushing of healthy donors [24]. (D) Maps of the expression of ssRNA virus receptors (rhinovirus C receptor (CDHR3), Coxsackievirus B1 receptor (CD55) and SARS-CoV-2 receptor (ACE2)) in single cells from the three models. UMAP and tSNE were used to show the cell types that express the ACE2, CDHR3 and CD55 receptors.

Conserved Expression of the Epithelial Cell Response to SARS-CoV-2 Infection in the iALI Model
Then, the bulk RNA-sequencing datasets of iPSC-derived alveolar epithelial type 2like cells (iAT2) infected by SARS-CoV-2 and mock-infected [14] at 1 and 4 dpi were analyzed to determine whether the SARS-CoV-2 bulk RNA-seq signature observed in the ALI model was present also in the iALI model. Analysis of the temporal distribution of RNAseq reads allows for identifying four major gene groups ( Figure 3A): (a) genes that were upregulated early (at 1 dpi) and the expression of which gradually increased from 1 to 4 dpi (NFKBIA, CSF1, IL32 and FOSL1); (b) genes that were upregulated early and the expression of which was not changed at 4 dpi (IL23A, CXCL10, CXCL20 and PLAUR); (c) genes that were upregulated at 1 dpi and were then downregulated at 4 dpi (CXCL3,

Conserved Expression of the Epithelial Cell Response to SARS-CoV-2 Infection in the iALI Model
Then, the bulk RNA-sequencing datasets of iPSC-derived alveolar epithelial type 2-like cells (iAT2) infected by SARS-CoV-2 and mock-infected [14] at 1 and 4 dpi were analyzed to determine whether the SARS-CoV-2 bulk RNA-seq signature observed in the ALI model was present also in the iALI model. Analysis of the temporal distribution of RNA-seq reads allows for identifying four major gene groups ( Figure 3A): (a) genes that were upregulated early (at 1 dpi) and the expression of which gradually increased from 1 to 4 dpi (NFKBIA, CSF1, IL32 and FOSL1); (b) genes that were upregulated early and the expression of which was not changed at 4 dpi (IL23A, CXCL10, CXCL20 and PLAUR); (c) genes that were upregulated at 1 dpi and were then downregulated at 4 dpi (CXCL3, CXCL5, CXCL1 and LOX); and (d) genes that were upregulated only at 4 dpi (MMP9 and MMP13). In addition, and like in the ALI model, in the COVID-19 iALI model, the interferon response was activated, as indicated by the upregulation of genes encoding interferoninduced proteins and interferon regulatory factors ( Figure 3B). This suggests that the ALI and iALI models share common interferon response regulation features. Conversely, the expression kinetics of ECM-related genes showed that some genes (ITGB6, ITGB1, GJA1, VIM and the ECM regulator PLOD2, which encode proteins of the host cell cytoskeleton structure) were progressively downregulated during SARS-CoV-2 infection in the iALI model ( Figure 3C). This indicated that the ECM and adhesion pathways are affected during cell-to-cell SARS-CoV-2 transmission. Moreover, the expression of EpCAM (epithelial cell adhesion molecule), an epithelial cell surface marker, progressively decreased, and that of the stromal cell marker COL1A1 progressively increased ( Figure 3D), suggesting epithelial mesenchymal transition. CXCL5, CXCL1 and LOX); and (d) genes that were upregulated only at 4 dpi (MMP9 and MMP13). In addition, and like in the ALI model, in the COVID-19 iALI model, the interferon response was activated, as indicated by the upregulation of genes encoding interferon-induced proteins and interferon regulatory factors ( Figure 3B). This suggests that the ALI and iALI models share common interferon response regulation features. Conversely, the expression kinetics of ECM-related genes showed that some genes (ITGB6, ITGB1, GJA1, VIM and the ECM regulator PLOD2, which encode proteins of the host cell cytoskeleton structure) were progressively downregulated during SARS-CoV-2 infection in the iALI model ( Figure 3C). This indicated that the ECM and adhesion pathways are affected during cell-to-cell SARS-CoV-2 transmission. Moreover, the expression of Ep-CAM (epithelial cell adhesion molecule), an epithelial cell surface marker, progressively decreased, and that of the stromal cell marker COL1A1 progressively increased ( Figure  3D), suggesting epithelial mesenchymal transition. Expression levels of (A) inflammatory cytokines/chemokines, including (a) genes that were upregulated early and the expression of which gradually increased upon SARS-CoV-2 infection, from 1 to 4 dpi; (b) genes that were upregulated early, the expression of which remained constant between 1 and 4 dpi; (c) genes that were upregulated specifically at 1 dpi and were then downregulated at 4 dpi; and (d) genes that became upregulated at 4 dpi. Expression of genes implicated in the (B) interferon response, (C) extracellular matrix and (D) epithelial/stromal gene upon infection. Normalized expression levels (counts per million reads) were quantified via RNA-seq using data from [14] (human iPSC-derived AT2 cells infected or not with SARS-CoV-2).

Experimental Validation of the Epithelial Cell Response to SARS-CoV-2 Infection in the iALI Model
To gain insight into the molecular basis of SARS-CoV-2 infection in the iALI model, the iALI bronchial epithelium was incubated with SARS-CoV-2 Delta (21A-Delta-B.1.617.2) at low MOI (0.05) to let the infection take hold for 4 days, and then the expression of interferon-induced genes and regulatory factors, viral entry genes and ECM genes was investigated. Viral RNA quantification confirmed infection of iALI cultures. The virus Expression levels of (A) inflammatory cytokines/chemokines, including (a) genes that were upregulated early and the expression of which gradually increased upon SARS-CoV-2 infection, from 1 to 4 dpi; (b) genes that were upregulated early, the expression of which remained constant between 1 and 4 dpi; (c) genes that were upregulated specifically at 1 dpi and were then downregulated at 4 dpi; and (d) genes that became upregulated at 4 dpi. Expression of genes implicated in the (B) interferon response, (C) extracellular matrix and (D) epithelial/stromal gene upon infection. Normalized expression levels (counts per million reads) were quantified via RNA-seq using data from [14] (human iPSC-derived AT2 cells infected or not with SARS-CoV-2).

Experimental Validation of the Epithelial Cell Response to SARS-CoV-2 Infection in the iALI Model
To gain insight into the molecular basis of SARS-CoV-2 infection in the iALI model, the iALI bronchial epithelium was incubated with SARS-CoV-2 Delta (21A-Delta-B.1.617.2) at low MOI (0.05) to let the infection take hold for 4 days, and then the expression of interferon-induced genes and regulatory factors, viral entry genes and ECM genes was investigated. Viral RNA quantification confirmed infection of iALI cultures. The virus was detected at the apical side at 1 dpi and also at 4 dpi ( Figure 4A). Tissue integrity was controlled through transepithelial electrical resistance (TEER) measurements ( Figure 4B). No difference between infected and non-infected (control) conditions was observed. Im-munofluorescence analysis of infected iALI cultures showed that SARS-CoV-2 membrane M protein colocalized with tubulin β 4a (TubIV). This confirmed the culture infection and suggested a preferential infection of ciliated cells ( Figure 4C). In non-infected cells, iALI ciliated cells were more homogeneous and intact. Moreover, RT-qPCR analysis of infected iALI cultures showed activation of the interferon signaling pathways at 1 dpi and 4 dpi. In the non-infected control, gene expression levels at day 0 and day 4 were not different ( Figure S3). Interferon-induced proteins (IFI44, IFIT1, IFIT3, IFI35), interferon regulatory factors (IRF9, MX1, ISG15), oligoadenylate synthetases (OAS1, OAS3) and chemokine ligand (CXCL10) were upregulated at 4 dpi ( Figure 4D). Conversely, viral infection did not alter the expression of genes encoding adhesion molecules, such as GJA1 and ITGB1, unlike what was observed upon infection at higher MOI (0.5) ( Figure 3C). In agreement with this observation, TEER of our own iALI bronchial model was comparable in non-infected and infected (at low MOI = 0.05) cultures from 1 to 4 dpi, indicating that at this infection level, epithelium integrity was not disrupted ( Figure 4B). Altogether, these data show that the iALI model is a trustworthy and sensitive model for SARS-CoV-2 infection.
was detected at the apical side at 1 dpi and also at 4 dpi ( Figure 4A). Tissue integrity was controlled through transepithelial electrical resistance (TEER) measurements ( Figure 4B). No difference between infected and non-infected (control) conditions was observed. Immunofluorescence analysis of infected iALI cultures showed that SARS-CoV-2 membrane M protein colocalized with tubulin β 4a (TubIV). This confirmed the culture infection and suggested a preferential infection of ciliated cells ( Figure 4C). In non-infected cells, iALI ciliated cells were more homogeneous and intact. Moreover, RT-qPCR analysis of infected iALI cultures showed activation of the interferon signaling pathways at 1 dpi and 4 dpi. In the non-infected control, gene expression levels at day 0 and day 4 were not different ( Figure S3). Interferon-induced proteins (IFI44, IFIT1, IFIT3, IFI35), interferon regulatory factors (IRF9, MX1, ISG15), oligoadenylate synthetases (OAS1, OAS3) and chemokine ligand (CXCL10) were upregulated at 4 dpi ( Figure 4D). Conversely, viral infection did not alter the expression of genes encoding adhesion molecules, such as GJA1 and ITGB1, unlike what was observed upon infection at higher MOI (0.5) ( Figure 3C). In agreement with this observation, TEER of our own iALI bronchial model was comparable in non-infected and infected (at low MOI = 0.05) cultures from 1 to 4 dpi, indicating that at this infection level, epithelium integrity was not disrupted ( Figure 4B). Altogether, these data show that the iALI model is a trustworthy and sensitive model for SARS-CoV-2 infection.  Data are the mean ± SEM of three independent experiments with three technical replicates/each (9 samples); * p < 0.05, ** p < 0.01, *** p < 0.001, **** p < 0.0001, ns: not significant (Student's t test). MOI: multiplicity of cellular infection.

Epithelial Cell Communication Networks in Response to SARS-CoV-2 Infection
ALI cultures include different cell types connected by tight and adherens junctions. Communication between epithelial cells occurs through the release of a variety of small molecules, including cytokines and chemokines. To investigate the impact of SARS-CoV-2 infection on ligand/receptor interaction between the different airway cell types, scRNAseq data from HBEC ALI cultures, infected or not with SARS-CoV-2 (MOI 0.01), were analyzed using the SingleCellSignalR method. Comparison of the summary chord diagram indicated a decrease in the number of paracrine interactions between the main epithelial cell types two dpi ( Figure 5A), which strongly suggests an impact of SARS-CoV2 infection on intercellular communications. For instance, expression of receptor tyrosine kinase (RET) and its ligand artemin (ARTN), one of the most pre-eminent interacting pairs between neuroendocrine cells and other epithelial cells, was strongly decreased at 3 dpi ( Figure 5B). This suggests that neuroendocrine cells are particularly sensitive to environmental stimuli (e.g., viral infection) and act as a rheostat to orchestrate ALI culture responses. replicates/each (9 samples); * p < 0.05, ** p < 0.01, *** p < 0.001, **** p < 0.0001, ns: not significant (Student's t test). MOI: multiplicity of cellular infection.

Epithelial Cell Communication Networks in Response to SARS-CoV-2 Infection
ALI cultures include different cell types connected by tight and adherens junctions. Communication between epithelial cells occurs through the release of a variety of small molecules, including cytokines and chemokines. To investigate the impact of SARS-CoV-2 infection on ligand/receptor interaction between the different airway cell types, scRNAseq data from HBEC ALI cultures, infected or not with SARS-CoV-2 (MOI 0.01), were analyzed using the SingleCellSignalR method. Comparison of the summary chord diagram indicated a decrease in the number of paracrine interactions between the main epithelial cell types two dpi ( Figure 5A), which strongly suggests an impact of SARS-CoV2 infection on intercellular communications. For instance, expression of receptor tyrosine kinase (RET) and its ligand artemin (ARTN), one of the most pre-eminent interacting pairs between neuroendocrine cells and other epithelial cells, was strongly decreased at 3 dpi (Figure 5B). This suggests that neuroendocrine cells are particularly sensitive to environmental stimuli (e.g., viral infection) and act as a rheostat to orchestrate ALI culture responses.

RNA Velocity Reveals Discrepancies between Non-Infected and Infected Epithelium in the ALI Model
To further characterize the ALI model transcriptional dynamics in response to virus infection, the single-cell RNA velocity was measured in SARS-CoV-2-infected and -noninfected epithelial ALI cultures using a dynamic model of the transcriptional state based on unspliced and spliced transcript counts [28]. Cell trajectory analysis revealed that non-infected ALI epithelial cells displayed two distinct bifurcation points through two different epithelial transition states. The first bifurcation point included basal cells (TP63-/KRT5+) that directly differentiate into ciliated cells, and the second bifurcation point mainly included basal cells (TP63+/KRT5+) that preferentially differentiate into secretory cells and then into ciliated cells. This placed secretory cells in an intermediate position between basal cells and mature ciliated cells (DNAH9+) ( Figure 5C). Velocity ordering analysis in infected cells revealed a change in the ciliated cell differentiation trajectory. Indeed, basal cells were all oriented to differentiate into ciliated cells through a secretory cell state. Moreover, DNAH9+-infected cell mapping within ALI epithelial cells revealed that ciliated cells were more susceptible to infection compared with other cell types, like DNAH9+ mature ciliated cells that are preferentially eliminated during SARS-CoV-2 infection ( Figure 5D).

Potential miRNA Regulators of the Epithelial Cell Intrinsic Response to SARS-CoV-2 Infection
GenGo Metacore software was then used to identify miRNAs that regulate the epithelial cell response genes after SARS-CoV-2 infection in the ALI and iALI models (i.e., potential antiviral targets). This analysis identified 167 miRNAs that regulate the key epithelial cell intrinsic genes deregulated upon viral infection ( Figure 3A,B). Among these miRNAs, 44 were regulators of interferon response genes and 123 were regulated inflammatory response genes (complete list of miRNAs in Supplementary Table S6). More than 54% of the interferon gene targets were regulated by more than one miRNA. For instance, SOCS3 and MX2 were targeted by 24 and 5 miRNAs, respectively ( Figure 6A). The inflammatory genes regulated by the highest number of miRNAs were MMP9 (n = 29 miRNAs), NFKBIA (n = 19), MMP13 (n = 19), CSF1 (n = 12), FOSL1 (n = 11) and LOX (n = 11) ( Figure 6B). MIR34a was identified as a potential regulator of MMP9, NFKBIA, FOSL1 and CXCL10, whereas MIR203 and MIR19 were common regulators of NFKBIA, MMP13 and SOCS3. MIR138-5p and MIR-326-3p are validated regulators of ISG15 and ISG20. MIR650, MIR541-3p and MIR302d-3p regulated several interferon-stimulated genes, including MX1, MX2, IRF7 and IRF9. CSF1 and LOX were targets of MIR130a and MIR29, respectively. These miRNAs can attenuate the inflammatory response and inhibit key coagulation cascade factors, thus preventing inflammatory epithelium damage. Moreover, 45 miRNAs were differentially expressed between SARS-CoV-2-infected and control lung-derived epithelial cells from an independent study (GSE148729) [29]. Among them, 17 miRNAs were significantly upregulated after 24 h of infection, and 28 were downregulated. These microRNAs differentially expressed in lung epithelial cells following SARS-CoV-2 infection were compared to the list of 129 microRNAs identified with GenGo Metacore using Venn diagrams. Among the 17 upregulated miRNAs (published datasets), 5 modulated the epithelial response to SARS-CoV-2 infection in both ALI and iALI models (MIR155-5p, MIR483-3p, MIR483-5p, MIR125b-5p and MIR4284). Conversely, downregulated miRNAs were not present on our list.

Discussion
The human large airways are lined by an epithelium with three abundant specialized cell types (basal, secretory and multi-ciliated cells) and some rare cell types (neuroendocrine cells, tuft cells and ionocytes) [30][31][32], on a stroma composed of mesenchymal cells, smooth muscle cells and immune cells [33]. Primary HBECs cultured in ALI conditions can be used to study airways in vitro. More recently, our group and others have generated a functional airway epithelium from human iPSCs (iALI models) with a marked similarity to the airway epithelium in vivo [4,6,[34][35][36]. As expected, scRNA-seq analysis of ALI and iALI models showed that immune cells, which are found in biopsy-/brushing-derived primary cultures, were absent in these models. In contrast, the stromal compartment is uniquely present in the iALI model, contributing to the faithful replication of the complex interplay between the bronchial epithelium and submucosa, which could prove particularly valuable in the context of viral infection. Similarly, neuroendocrine cells were mostly identified in the iALI model, compared with the ALI model. Our analysis also characterized rare cell types, such as ionocytes, that strongly express FOXI1 and ASCL3 [31,37].
SARS-CoV-2 infection in ALI and iALI models induces virus-induced epithelial disruption, loss of mature ciliated cells and triggers intrinsic immune responses. Despite the differences in terms of cell composition, the ALI and iALI models exhibit relevant proportions of airway cell types, express SARS-CoV-2, rhinovirus C and Coxsackievirus B1 receptors (ACE2, CDHR3, CD55) and have been used to model various mechanisms of SARS-CoV-2 pathogenesis [38]. These models provide a suitable and reliable platform for researchers to study SARS-CoV-2 pathogenesis during infection.
In this study, we analyzed the bulk and single-cell transcriptomes of ALI and iALI cultures infected with SARS-CoV-2. This highlighted the emergence of pro-inflammatory and interferon signatures, in which epithelial cells upregulated the expression of several cytokines, chemokines, interferon and downregulated ECM-related genes in response to SARS-CoV-2 at different times after infection. Most of these genes play essential roles in virus control and also in disease development. For instance, many chemokines (CCL2, Figure 6. Identification of miRNA targets as potential COVID-19 therapeutics. Schematic representation of potential miRNAs that target some genes implicated in the interferon and inflammatory responses upon SARS-CoV-2 infection of ALI cultures. Only the miRNAs experimentally validated are represented. A line with a bar indicates "microRNA-targeting".

Discussion
The human large airways are lined by an epithelium with three abundant specialized cell types (basal, secretory and multi-ciliated cells) and some rare cell types (neuroendocrine cells, tuft cells and ionocytes) [30][31][32], on a stroma composed of mesenchymal cells, smooth muscle cells and immune cells [33]. Primary HBECs cultured in ALI conditions can be used to study airways in vitro. More recently, our group and others have generated a functional airway epithelium from human iPSCs (iALI models) with a marked similarity to the airway epithelium in vivo [4,6,[34][35][36]. As expected, scRNA-seq analysis of ALI and iALI models showed that immune cells, which are found in biopsy-/brushing-derived primary cultures, were absent in these models. In contrast, the stromal compartment is uniquely present in the iALI model, contributing to the faithful replication of the complex interplay between the bronchial epithelium and submucosa, which could prove particularly valuable in the context of viral infection. Similarly, neuroendocrine cells were mostly identified in the iALI model, compared with the ALI model. Our analysis also characterized rare cell types, such as ionocytes, that strongly express FOXI1 and ASCL3 [31,37].
SARS-CoV-2 infection in ALI and iALI models induces virus-induced epithelial disruption, loss of mature ciliated cells and triggers intrinsic immune responses. Despite the differences in terms of cell composition, the ALI and iALI models exhibit relevant proportions of airway cell types, express SARS-CoV-2, rhinovirus C and Coxsackievirus B1 receptors (ACE2, CDHR3, CD55) and have been used to model various mechanisms of SARS-CoV-2 pathogenesis [38]. These models provide a suitable and reliable platform for researchers to study SARS-CoV-2 pathogenesis during infection.
In this study, we analyzed the bulk and single-cell transcriptomes of ALI and iALI cultures infected with SARS-CoV-2. This highlighted the emergence of pro-inflammatory and interferon signatures, in which epithelial cells upregulated the expression of several cytokines, chemokines, interferon and downregulated ECM-related genes in response to SARS-CoV-2 at different times after infection. Most of these genes play essential roles in virus control and also in disease development. For instance, many chemokines (CCL2, CCL3, CCL20, CXCL1, CXCL3, CXCL10, IL-8) associated with inflammatory responses were expressed in the ALI and iALI models. CCL2, CXCL10 and IL-8 are associated with airway inflammation, and high serum levels of these chemokines were found in patients with severe SARS [39]. CCL3 is also involved in viral infections [40]. Our detailed analysis of the transcriptional response to SARS-CoV-2 infection showed that ALI and iALI cultures produced an unbalanced cytokine response, with preferential upregulation of genes encoding cytokines (for instance IL-6, IL-1β and IL-33) that are mainly implicated in the defense against extracellular aggression. Several studies showed that IL-6 serum levels are increased in patients with COVID-19 [41]. IL-6 is involved in acute inflammation due to its implication in controlling the acute phase response [42]. IL-6 production is induced by TNF-α and IL-1β [43]. In animal models of SARS-CoV infection, TNF activity neutralization provides protection against SARS-CoV morbidity and mortality [44]. A large number of data showed the role of interferons in SARS-CoV-2 infection. Interferons exercise their biological functions by regulating the expression of interferon-stimulated genes (ISGs). ISG upregulation has been described in various cells from patients with severe COVID-19 [45,46]. Here, we found that infection with SARS-CoV-2 induced a strong interferon response in ALI and iALI cultures, marked by a high expression of ISG15 (key factor in the innate immune response to SARS-CoV-2 infection), ISG20 (with antiviral activity against RNA viruses), IRF-7 (the master regulator of interferon responses) [47] and of several ISGs (IFITM1, IFITM3, IRF9, IFI27, OAS2, MX1, MX2, SOCS3) involved in the regulation of the host defense responses to the virus [48]. This confirms the known ALI/iALI cultures' intrinsic response to SARS-CoV-2 infection focused on the activation of the interferon pathways. This is in line with a study showing strong expression of many ISGs in the respiratory tract of patients with COVID-19, supporting the idea that interferon-mediated immune response plays a key role in SARS-CoV-2 infection control [46].
Although SARS-CoV-2 pathogenesis has not been fully understood, it seems that excessive immune responses play a key role. Evidence suggests that immune response deregulation causes lung damage [49]. Here, we found that in infected ALI and iALI cultures, chemokines and cytokines, including IL-6, IL-1β, interferon and TNF, were upregulated to coordinate all aspects of the immunogenic response to SARS-CoV-2 infection. Additionally, SARS-CoV-2-related pneumonia with severe respiratory failure is characterized by enhanced ECM [50]. Similarly, in infected ALI cultures, cells expressing MMP9, an enzyme that participates in ECM remodeling, were markedly increased. The molecular pathways involved in MMP-9 regulation during SARS-CoV-2 infection are not known. Ueland et al. suggested that MMP-9 may be an early indicator of respiratory failure in patients with COVID-19 [51], and Hsu et al. reported an important increase in MMP-9 concentration in the plasma of patients who developed acute respiratory distress syndrome [52]. Conversely, in infected ALI cultures, cells expressing ECM-related genes (e.g., VIM, ITGB1, ITGB6, GJA1 and PLOD2) were decreased compared with control cultures. Vimentin and integrin are critical targets for SARS-CoV-2 host cell invasion [53,54]. Identifying the molecular mechanisms that lead to their regulation will be pivotal to help us understand their role in epithelium damage and reparation during SARS-CoV-2 infection.
It is important to acknowledge that different MOIs may generate different transcriptomic profiles. Dose dependence is highly likely but extremely complex to address given the heterogeneity of the currently available data. Because of the internal consistency of our analysis, we defend the hypothesis that the findings we reported here represent at least the core of SARS-CoV-2-induced genes, irrespective of the used MOI, and that variations associated with MOI changes might be outside this core.
Another important question we addressed concerned the miRNA role in the regulation of genes deregulated in infected ALI/iALI cultures. Recent investigations revealed that miRNAs are implicated in viral pathogenesis by altering the miRNA-modulated host gene regulation or the host immune system [55]. The variation in miRNA levels during virus infection and their role in modulating SARS-CoV-2 infection in human cells have been well described [15]. Thus, miRNA antagonists or mimics could be used to develop new therapeutic strategies for the treatment of patients with COVID-19 [56]. The anti-viral response by miRNAs may implicate the regulation of their mRNA targets that participate in the cellular response to viral infection. Interferon signaling was one of the primary effectors against SARS-CoV-2 infections in both ALI and iALI models. Moreover, genes encoding members of this pathway are direct targets of many miRNAs. For instance, MIR138 regulates ISG15 expression by direct binding to its 3 untranslated region (UTR) [57]. Similarly, MIR326-3p reduces the activity of an ISG20 3 UTR luciferase reporter [58]. MIR650, MIR29a-3p and MIR130a-3p regulate other ISGs, such as MX1, MX2, IFITM3 and IFITM1 [59][60][61]. Moreover, many miRNAs target the 3 UTR site of matrix metalloproteaseencoding genes. For instance, MIR34a and MIR19 target NFKBIA [62] (critical for SARS-CoV-2 entry in the cells) and MMP9 [63] (macrophage-derived biomarker associated with inflammation), whereas MIR130a targets the 3 UTR of pro-inflammatory metabolite genes (e.g., LOX [64] that interacts with SARS-CoV-2) [65]. Altogether, different culture models are necessary to validate the biological effect of these candidate miRNAs to SARS-CoV-2 infection and to understand the interactions between host ISGs, cytokines and miRNAs. More investigations are needed to determine the expression profiles of key miRNAs during viral infection, to better predict disease severity and to develop new therapeutic options based on these miRNAs for COVID-19 treatment and/or prevention.

ALI Culture of Primary Airway Epithelial Cells and iPSC-Derived Airway Epithelium for scRNA-Seq Analysis
Primary Human Bronchial Epithelial Cells (HBECs) were expanded and differentiated in ALI culture following the protocol given by StemCell Technologies. Briefly, cells were dissociated mechanically from bronchial biopsy specimens obtained via fiberoptic bronchoscopy (approval number: 2013 11 05; NCT02354677) and cultured in PneumaCult-Ex Plus expansion medium (cat #05041, StemCell Technologies, Saint-Egrève, France) for 15 days. After the expansion phase, differentiation was initiated by seeding 1.1 × 10 5 cells/insert on Transwell R polyester membranes (cat #3460, Corning, Kennebunk, ME, USA). Once epithelial cells reached confluence, the apical growth medium was removed, and the basal medium was replaced with PneumaCult™-ALI maintenance medium (cat #05002, StemCell Technologies, France) (i.e., day 0 of ALI culture). Epithelial cells were allowed to differentiate at 37 • C with 5% CO 2 for 28 days. The iPSC-derived airway epithelium on ALI (iALI model) was generated as previously described [4]. Briefly, the major stages of embryonic lung development were recapitulated as follows: stage 1, definitive endoderm (day 0-3) using activin A (100 ng/mL); stage 2, anterior foregut endoderm (day 4-8); stage 3, lung progenitor specification (day 9); and stage 4, epithelial layer (day 14). After 40 days, the airway epithelium on iALI displayed morphologic and functional similarities with primary human airway epithelial cells and included different airway cell types (basal, secretory and multi-ciliated cells).

ScRNA-Seq and Data Analysis
Non-infected ALI and iALI cultures were dissociated with trypsin into single-cell suspensions. Cell viability (>83%) and aggregation (minimizing the presence of cell aggregates by using a cell strainer, with pore sizes of 30 to 40 microns) were tested before starting the single-cell library preparation. The concentration of freshly dissociated cells was adjusted to 1000 cells/µL in HBSS/0.05% BSA, and then the 10x Chromium Controller and the Chromium Single Cell 3 Reagent kit V3.1 were used to perform the scRNAseq experiments. Library preparation was performed according to the manufacturer's instructions using the Chromium Chip B Single Cell kit and Chromium Multiplex Kit (10X Genomics, Pleasanton, CA, USA). Sequencing was performed in paired-end mode with an S1 flow cell (28/8/87 cycles) and a NovaSeq 6000 sequencer (Illumina, San Diego, CA, USA) at the MGX core facility of Montpellier, France. First, the cell ranger mkfastq and cellranger count pipelines were used for the initial quality control, sample demultiplexing, mapping and raw data quantification. Fastq files were run with the Count application using default parameters and were aligned to the human genome reference sequence GRCh38, filtered and counted. C-loop software (version 6.2.0) was used to visualize clusters and sub-clusters of transcriptionally related cells and to identify candidate genes, the expression of which was enriched in specific clusters. Clustering results were visualized with the t-distributed Stochastic Neighbor-Embedding (tSNE) technique. Bronchial cell (biopsy/brushing samples) datasets included in [24] were analyzed through a COVID-19 Cell Atlas website portal and were visualized using Uniform Manifold Approximation and Projection (UMAP).

Functional Enrichment Analysis of Bulk RNA-Seq Datasets
The unique bulk RNA-seq signature was obtained from a publicly available list of differential expressed genes between the primary human lung epithelium infected with SARS-CoV-2 (multiplicity of infection, MOI, = 2 for 24 h) and mock-treated controls [18]. The GO functional enrichment and pathway enrichment analyses were performed with ShinyGO [26] and Gene Set Enrichment Analysis (GSEA) (http://www.broadinstitute.org/ gsea/, accessed on October 2022). GO annotations were divided into three categories: biological process, molecular function and cellular component. In the enrichment analysis, the Fisher' exact test was used to test whether genes were enriched in a term, and an adjusted p value < 0.05 was set as the screening condition. GenGo Metacore software (Version 4.7) was used to identify miRNA targets. Heatmaps and gene networks were generated with Ingenuity Pathway Analysis (IPA, QIAGEN, Redwood City, CA, USA).

Analysis of Publicly Available Single-Cell RNA-Seq Datasets
A public scRNA-seq dataset (GEO accession number: GSE166766) of HBEC ALI culture samples infected with SARS-CoV-2 (MOI~0.01) was also analyzed [66]. To further investigate the dynamic nature of ALI epithelial cells before and after infection, RNA velocity analysis was measured to estimate the future state of individual cells. The Velocyto ® package was used to investigate the gene expression dynamics in the scRNA-seq data before virus infection and at day 3 post-infection (dpi). RNA velocity quantifies the change in the state of a cell over time by distinguishing unspliced and spliced mRNAs reads. To obtain counts of spliced and unspliced mRNAs reads, Velocyto used the outputs of Cell-Ranger (version 3.0.2) alignments with the command line 'run10x' and the transcriptome GRCh38.p12 (accession NCBI:GCA_000001405.27). At this step, loom files were created for the input data of the scVelo tool [28]. This python package allowed for normalizing with the scv.pp.normalize_per_cell() function and transform scv.pp.log1p(). Genes were filtered by keeping only the top 2000 highly variable genes (n_top_genes = 2000 parameters for scv.pp.filter_and_normalize function). Then, moments were calculated for each cell across its nearest neighbors, with n_neighbors set to 30 and the first 10 PCs using the scv.pp.moments() function. Then, cell velocities were estimated using scvelo.tl.velocity() based on a stochastic model of transcriptional dynamics. To visualize the velocity graph, data were projected on the UMAP space and clusters defined by the CellRanger pipeline were colored. Seurat V4 and R we used to identify potential paracrine interactions (ligand receptors) between the different cell categories in ALI cultures infected with SARS-CoV-2 at 1, 2 and 3 dpi. After cluster identification on Seurat V4 [67], the function cluster_analysis from the SingleCellSignalR package [68] was used to compute paracrine interactions between cell clusters and to predict ligand-target links between interacting cells by combining their expression level with prior knowledge on gene regulatory networks and signaling pathways. To compare ALI and iALI samples, genes identified in our samples were combined with other publicly available datasets that used infected iPSC-derived AT2 cells (iAT2) [14], specifically for the genes related to inflammatory and interferon responses and extracellular ECM.

SARS-CoV-2 Virus Stock and Titration
The hCoV-19/France/HDF-IPP11602i/2021 (21A-Delta-B.1.617.2) strain was supplied by the National Reference Centre for Respiratory Viruses hosted by Institut Pasteur (Paris, France). The human sample from which this strain was isolated was provided by Dr Guiheneuf Raphaël, CH Simone Veil, Beauvais, France. The strain was propagated in VeroE6 cells with DMEM containing 25 mM HEPES at 37 • C and 5% CO 2 , and viruses were harvested 72 h post-inoculation. Virus stocks were stored at −80 • C. Viruses from infected cell culture supernatants were titrated with the plaque assays on a monolayer of VeroE6 cells and 100 µL of virus serial dilutions. The plaque-forming unit (PFU) values were determined through crystal violet staining and then scoring the wells displaying cytopathic effects. The virus titer was determined as the number of PFU/mL, and MOI was the PFU/cell ratio.

iALI Infection by SARS-CoV-2
iALI cultures were washed 24 h before infection by adding culture medium to the apical side at 37 • C for 15 min. For SARS-CoV-2 infection, iALI apical side was covered with 350 µL of the appropriate viral dilution (2.5 × 10 5 PFU per sample) for 1 h 30 min at 37 • C. Then, the inoculum was removed, and the culture was quickly washed with culture medium. The first point of the kinetic analysis (Day 0) was collected by adding culture medium to the apical side at 37 • C for 15 min. iALI's apical side were then washed every day with culture medium to collect kinetics samples. SARS-CoV-2 infection experiments were conducted inside the BSL3 facility at the CEMIPAI institute.

Reverse Transcription-Quantitative Polymerase Chain Reaction (RT-qPCR)
RNA was extracted from cells using the QIAshredder kit (QIAGEN, Redwood city, CA, USA) and the RNeasy mini kit (Qiagen, Redwood city, CA, USA) according to the manufacturer's instructions. Viral RNA was quantified by RT-qPCR in triplicate, as described [69], using the Luna Universal One-Step RT-qPCR Kit (New England Biolabs, Ipswich, MA, USA) and a BIORAD CFX Opus 384 system. iALI cultures were washed by adding culture medium at the apical side at 37 • C for 15 min; then, viral RNA was extracted and quantified, in triplicate, via RT-qPCR (amplification of the viral envelope E gene). The relative gene expression was calculated for each triplicate by normalizing to GAPDH gene level (control) and using the ∆∆Ct method, as previously described in [70]. Primers are listed in Table S1. For four genes (GJA1, IRF9, OAS1 and CXCL10), primers were designed using the website Primer3+. Their efficiency was confirmed using 2-fold serial dilutions of a cDNA pool. As IRF9, OAS1 and CXCL10 displayed a slope of −3.3 ± 0.66 with an R 2 > 0.98 (only GJA1 R 2 = 0.96), an efficiency of 2 was used for the relative expression calculation (Supplementary Figure S2).

Measurements of Transepithelial Electrical Resistance (TEER)
Transepithelial resistance was measured daily with an EVOM2 (WPI, Friedberg, Germany), while the apical side was submerged in culture medium. Measurements were performed after 10 min of incubation at 37 • C.

Statistical Analysis
Data are presented as the mean ± SEM, unless otherwise specified. Statistical analyses were performed using GraphPad Prism 5 software, version 5.0 (Student's t-test; GraphPad, San Diego, CA, USA). The shown data were from representative experiments, with similar results in at least three independent biological replicates, unless otherwise specified. Differences were evaluated using the Student's t-test. A p value ≤ 0.05 was considered significant.

Conclusions
Here, we characterized two ALI and iALI airway models to understand SARS-CoV-2 infection pathogenesis. We identified the inflammatory and interferon profiles induced in these models in response to SARS-CoV-2 infection. We provided a more detailed molecular portrait of the genes and cell types associated with SARS-CoV-2 infection in the ALI and iALI airway models. Further investigations are needed to determine whether similar responses take place in other SARS-CoV-2 target tissues models.