key: cord-255440-ls1l2mlg
authors: Tindle, Courtney; Fuller, MacKenzie; Fonseca, Ayden; Taheri, Sahar; Ibeawuchi, Stella-Rita; Beutler, Nathan; Claire, Amanraj; Castillo, Vanessa; Hernandez, Moises; Russo, Hana; Duran, Jason; Crotty Alexander, Laura E.; Tipps, Ann; Lin, Grace; Thistlethwaite, Patricia A.; Chattopadhyay, Ranajoy; Rogers, Thomas F.; Sahoo, Debashis; Ghosh, Pradipta; Das, Soumita
title: Adult Stem Cell-derived Complete Lung Organoid Models Emulate Lung Disease in COVID-19
date: 2020-10-18
journal: bioRxiv
DOI: 10.1101/2020.10.17.344002
sha: 
doc_id: 255440
cord_uid: ls1l2mlg

SARS-CoV-2, the virus responsible for COVID-19, causes widespread damage in the lungs in the setting of an overzealous immune response whose origin remains unclear. We present a scalable, propagable, personalized, cost-effective adult stem cell-derived human lung organoid model that is complete with both proximal and distal airway epithelia. Monolayers derived from adult lung organoids (ALOs), primary airway cells, or hiPSC-derived alveolar type-II (AT2) pneumocytes were infected with SARS-CoV-2 to create in vitro lung models of COVID-19. Infected ALO-monolayers best recapitulated the transcriptomic signatures in diverse cohorts of COVID-19 patient-derived respiratory samples. The airway (proximal) cells were critical for sustained viral infection, whereas distal alveolar differentiation (AT2→AT1) was critical for mounting the overzealous host immune response in fatal disease; ALO monolayers with well-mixed proximodistal airway components recapitulated both. Findings validate a human lung model of COVID-19, which can be immediately utilized to investigate COVID-19 pathogenesis and vet new therapies and vaccines. GRAPHIC ABSTRACT HIGHLIGHTS Human lung organoids with mixed proximodistal epithelia are created Proximal airway cells are critical for viral infectivity Distal alveolar cells are important for emulating host response Both are required for the overzealous response in severe COVID-19 IN BRIEF An integrated stem cell-based disease modeling and computational approach demonstrate how both proximal airway epithelium is critical for SARS-CoV-2 infectivity, but distal differentiation of alveolar pneumocytes is critical for simulating the overzealous host response in fatal COVID-19.

SARS-CoV-2, the virus responsible for COVID-19, causes widespread inflammation and injury in the lungs, giving rise to diffuse alveolar damage (DAD) (Arrossi and Farver, 2020; Borczuk et al., 2020; Damiani et al., 2020; Li et al., 2020; Roden et al., 2020) , featuring marked infection and viral burden leading to apoptosis of alveolar pneumocytes (Hussman, 2020) , along with pulmonary edema (Bratic and Larsson, 2013; Carsana et al., 2020) . DAD leads to poor gas exchange and, ultimately, respiratory failure; the latter appears to be the final common mechanism of death in most patients with severe COVID-19 infection. How the virus causes so much damage remains unclear. A particular challenge is to understand the out-of-control immune reaction to the SARS-CoV-2 infection known as a cytokine storm, which has been implicated in many of the deaths from COVID-19. Although rapidly developed pre-clinical animal models have recapitulated some of the pathognomonic aspects of infection, e.g., induction of disease, and transmission, and even viral shedding in the upper and lower respiratory tract, many failed to develop severe clinical symptoms (Lakdawala and Menachery, 2020) . Thus, the need for pre-clinical models remains both urgent and unmet.

To address this need, several groups have attempted to develop human pre-clinical COVID-19 lung models, all within the last few months (Duan et al., 2020; Mulay et al., 2020; Salahudeen et al., 2020) . While a head-to-head comparison of the key characteristics of each model can be found in Table 1 , what is particularly noteworthy is that none recapitulate the heterogeneous epithelial cellularity of both proximal and distal airways, i.e., airway epithelia, basal cells, secretory club cells and alveolar pneumocytes. Also noteworthy is that models derived from iPSCs lack propagability and/or cannot be reproducibly generated for biobanking; nor can they be scaled up in cost-effective ways for use in drug screens. Besides the approaches described so far, there are a few more approaches used for modeling COVID-19-(i) 3D organoids from bronchospheres and tracheospheres have been established before (Hild and Jaffe, 2016; Rock et al., 2009; Tadokoro et al., 2016) and are now used in apical-out cultures for infection with SARS-COV-2 (Suzuki et al., 2020); (ii) the most common model used for drug screening is the air-liquid interphase (ALI model) in which pseudo-stratified primary bronchial or small airway epithelial cells are used to recreate the multilayered mucociliary epithelium (Mou et al., 2016; Randell et al., 2011) ; (iii) several groups have also generated 3D airway models from iPSCs or tissue-resident stem cells (Dye et al., 2015; Ghaedi et al., 2013; Konishi et al., 2016; McCauley et al., 2017; Miller et al., 2019; Wong et al., 2012) ; (iv) others have generated AT2 cells from iPSCs using closely overlapping protocols of sequential differentiation starting with definitive endoderm, anterior foregut endoderm, and distal alveolar expression (Chen et al., 2017; Gotoh et al., 2014; Huang et al., 2014; Jacob et al., 2017; Jacob et al., 2019; Yamamoto et al., 2017) . (v) Finally, long term in vitro culture conditions for pseudo-stratified airway epithelium organoids, derived from healthy and diseased adult humans suitable to assess virus infectivity (Sachs et al., 2019; van der Vaart and Clevers, 2020; Zhou et al., 2018) have been pioneered; unfortunately, these airway organoids expressed virtually no lung mesenchyme or alveolar signature. What remains unclear is if any of these models accurately recapitulate the immunopathologic phenotype that is seen in the lungs in COVID-19.

We present a rigorous transdisciplinary approach that systematically assesses an adult lung organoid model that is propagable, personalized and complete with both proximal airway and distal alveolar cell types against existing models that are incomplete, and we cross-validate them all against COVID-19 patient-derived respiratory samples. Findings surprisingly show that cellular crosstalk between both proximal and distal components are necessary to emulate how SARS-CoV-2 causes diffuse alveolar pneumocyte damage; the proximal airway mounts a sustained viral infection, but it is the distal alveolar pneumocytes that mount the overzealous host response that has been implicated in a fatal disease.

To determine which cell types in the lungs might be most readily infected, we began by analyzing a human lung single-cell sequencing dataset (GSE132914) for the levels of expression of angiotensin-converting enzyme-II (ACE2) and Transmembrane Serine Protease 2 (TMPRSS2), the two receptors that have been shown to be the primary sites of entry for the SARS-CoV-2 (Hoffmann et al., 2020) . The dataset was queried with widely accepted markers of all the major cell types (see Table 2 ). Alveolar epithelial type 2 (AT2), ciliated and club cells emerged as the cells with the highest expression of both receptors (Fig 1A; Fig S1A) . These observations are consistent with published studies demonstrating that ACE2 is indeed expressed highest in AT2 and ciliated cells (Jia et al., 2005; Mulay et al., 2020; Zhao et al., 2020) . In a cohort of deceased COVID-19 patients, we observed by H&E (Fig S1B) that gas-exchanging flattened AT1 pneumocytes are virtually replaced by cuboidal cells that were subsequently confirmed to be AT2-like cells via immunofluorescent staining with the AT2-specific marker, surfactant protein-C (SFTPC; Fig 1B upper panel; Fig S1C; top). We also confirmed that club cells express ACE2 (Fig S1C; bottom) , underscoring the importance of preserving these cells in any ideal lung model of COVID-19. When we analyzed the lungs of deceased COVID-19 patients, the presence of SARS-COV-2 in alveolar pneumocytes was also confirmed, as determined by the colocalization of viral nucleocapsid protein with SFTPC (Fig 1B; lower panel; Fig S1D) . Immunohistochemistry studies further showed the presence of SARS-COV-2 virus in alveolar pneumocytes and in alveolar immune cells (Fig S1E) . These findings are consistent with the gathering consensus that alveolar pneumocytes support the interaction between the epithelial cells and inflammatory cells recruited to the lung; via mechanisms that remain unclear, they are generally believed to contribute to the development of acute lung injury and acute respiratory distress syndrome (ARDS), the severe hypoxemic respiratory failure during COVID-19 (Hou et al., 2020; Spagnolo et al., 2020) . Because prior work has demonstrated that SARS-CoV-2 infectivity in patient-derived airway cells is highest in the proximal airway epithelium compared to the distal alveolar pneumocytes (AT1 and AT2) (Hou et al., 2020) , and yet, it is the AT2 pneumocytes that harbor the virus, and the AT1 pneumocytes that are ultimately destroyed during diffuse alveolar damage, we hypothesized that both proximal airway and distal (alveolar pneumocyte) components might play distinct roles in the respiratory system to mount the so-called viral infectivity and host immune response phases of the clinical symptoms observed in COVID-19 (Chen and Li, 2020) .

Because no existing lung model provides such proximodistal cellular representation (Table 1) , and hence, may not recapitulate with accuracy the clinical phases of COVID-19, we first sought to develop a lung model that is complete with both proximal and distal airway epithelia using adult stem cells that were isolated from deep lung biopsies. Lung organoids were generated using the protocol outlined in Fig 1C and 

Methods. Organoids grown in 3D cultures were subsequently dissociated into single cells to create 2Dmonolayers (either maintained submerged in media or used in ALI model) for SARS-CoV-2 infection, followed by RNA Seq analysis. Primary airway epithelial cells and hiPSC-derived alveolar type-II (AT2) pneumocytes were used as additional models (Fig 1D; left panel) . Each of these transcriptomic datasets was subsequently used to cross-validate our ex-vivo lung models of SARS-CoV-2 infection with the human COVID-19 autopsy lung specimens (Fig 1D; right panel) to objectively vet each model for their ability to accurately recapitulate the gene expression signatures in the patient-derived lungs.

Three lung organoid lines were developed from deep lung biopsies obtained from the normal regions of lung lobes surgically resected for lung cancer; both genders, smokers and non-smokers were represented (Fig S2A; Table 3). Three different types of media were compared (Fig S2B) ; the composition of these media was inspired either by their ability to support adult-stem cell-derived mixed epithelial cellularity in other organs (like the gastrointestinal tract (Miyoshi and Stappenbeck, 2013; Sato et al., 2009; Sayed et al., 2020) , or rationalized based on published growth conditions for proximal and distal airway components (Gotoh et al., 2014; Sachs et al., 2019; van der Vaart and Clevers, 2020) . A growth condition that included conditioned media from L-WRN cells which express Wnt3, R-spondin and Noggin, supplemented with recombinant growth factors, which we named as 'lung organoid expansion media' emerged as superior compared to alveolosphere media-I and II (Jacob et al., 2019; Yamamoto et al., 2017) (details in the methods), based on its ability to consistently and reproducibly support the best morphology and growth characteristics across multiple attempts to isolate organoids from lung tissue samples. Three adult lung organoid lines (ALO1-3) were developed using the expansion media, monitored for their growth characteristics by brightfield microscopy and cultured with similar phenotypes until P10 and beyond (Fig S2C-D) . The 3D morphology of the lung organoid was also assessed by H&E staining of slices cut from formalin-fixed paraffin-embedded (FFPE) cell blocks of HistoGel-embedded ALO1-3 (Fig S2E) .

To determine if all the 6 major lung epithelial cells (illustrated in Fig 2A) are present in the organoids, we analyzed various cell-type markers by qRT-PCR (Fig 2B-H) . All three ALO lines had a comparable level of AT2 cell surfactant markers (compared against hiPSC-derived AT2 cells as positive control) and a significant amount of AT1, as determined using the marker AQP5. ALOs also contained basal cells (as determined by the marker ITGA6), ciliated cells (as determined by the marker FOXJ1), club cells (as determined by the marker SCGB1A1) and stem cells (as determined by marker p75-NGFR). As expected, the primary human bronchial epithelial cells (NHBE) had significantly higher expression of basal cell markers than the ALO lines (hence, served as a positive control), but they lacked stemness and club cells (hence, served as a negative control).

The presence of all cell types was also confirmed by assessing protein expression of various cell types within organoids grown in 3D cultures. Two different approaches were used-(i) slices cut from FFPE cell blocks of HistoGel-embedded ALO lines (Fig 2I-J) or (ii) ALO lines grown in 8-well chamber slides were fixed in matrigel (Fig 2K) , stained, and assessed by confocal microscopy. Such staining not only confirmed the presence of all cell types in each ALO line but also demonstrated the presence of more than one cell type (i.e., mixed cellularity) of proximal (basal-KRT5) and distal (AT1/AT2 markers) within the same organoid structure. For example, AT2 and basal cells, marked by SFTPB and KRT5, respectively, were found in the same 3D-structure (Fig 2J, interrupted curved line). Similarly, ciliated cells and goblet cells stained by Ac-Tub and Muc5, respectively, were found to coexist within the same structure (Fig 2J, interrupted box; Fig 2K, arrow) . The presence of heterogeneous cellularity was documented in all three ALO lines (see multiple additional examples in Fig S3) .

To model respiratory infections such as COVID-19, it is necessary for pathogens to be able to access the apical surface. Many groups have modeled such access by dissociating 3D organoids into single cells and plating them as 2D-monolayers (Duan et al., 2020; Han et al., 2020; Huang et al., 2020; Mulay et al., 2020; Sachs et al., 2019; Zhou et al., 2018) . Because the loss of dimensionality can have a major impact on cellular proportions and impact disease-modeling in unpredictable ways, we assessed the impact of the 3D-to-2D conversion on cellularity by RNA seq analyses. Two commonly encountered methods of growth in 2D-monolayers were tested: (i) monolayers polarized on trans-well inserts but submerged in growth media, and ii) monolayers were grown at the air-liquid interface (popularly known as the 'ALI model '(Dvorak et al., 2011; Prytherch et al., 2011) ) for 21-days to differentiate into the mucociliary epithelium (see Fig 3A; Fig S4A-E) . The submerged 2D-monolayers had several regions of organized vacuolated-appearing spots (Fig S4A, D; arrow), presumably due to morphogenesis and cellular organization even in 2D. The ALI-monolayers appeared to be progressively hazier with time after air-lift, likely due to the accumulation of secreted mucin (Fig S4E) , and formed an intact epithelial barrier, as determined by trans-epithelial electrical resistance (TEER) (Fig S4C) . RNA Seq datasets were analyzed using the same set of cell markers, as we used in Fig 1A ( Table 2 ). Cell-type deconvolution of our dataset using CIBERSORTx (https://cibersortx.stanford.edu/runcibersortx.php) showed that cellular proportions in the human lung tissues were also relatively well-preserved in organoids grown in 3D over several passages (Fig 3B; left) ; both showed a mixed population of simulated alveolar, basal, club, ciliated and goblet cells. When 3D organoids were dissociated and plated as 2D monolayers on transwells, the AT2 signatures were virtually abolished with a concomitant and prominent emergence of AT1 signatures, suggesting that growth in 2D-monolayers favor differentiation of AT2 cells into AT1 cells (Shami and Evans, 2015) (Fig 3B; middle) . A compensatory reduction in proportion was also observed for the club, goblet and ciliated cells. The same organoids, when grown in long-term 2D culture conditions in the ALI model, showed a strikingly opposite pattern; alveolar signatures were almost entirely replaced by a concomitant increase in ciliated and goblet cells (Fig 3B; right). These findings are consistent with the well-established notion that ALI conditions favor growth as pseudo-stratified mucociliary epithelium (Dvorak et al., 2011; Prytherch et al., 2011) . As an alternative model for use as monolayers for viral infection, we developed hiPSC-derived AT2 cells and alveolospheres (Fig 3C) , using established protocols (Huang et al., 2020) . Because they were grown in the presence of CHIR99021 (an aminopyrimidine derivate that is a selective and potent Wnt agonist) (Abdelwahab et al., 2019; Jacob et al., 2019; Yamamoto et al., 2017) , which probably inhibits the AT2→AT1 differentiation, these monolayers were enriched for AT2 and devoid of AT1 cells (Fig 3D) .

The multicellularity of lung organoid monolayers was also confirmed by immunofluorescence staining and confocal microscopy of the submerged and ALI monolayers, followed by the visualization of cell markers in either max-projected z-stacks (Fig 3E; left) or orthogonal views of the same (Fig 3E; right) . As expected, markers for the same cell type (i.e., SFTPB and SFTPC, both AT2 markers) colocalize, but markers for different cell types do not. Submerged monolayers showed the prominent presence of both AT1 (AQP5-positive) and AT2 cells.

Compared to the submerged monolayers, the ALI model showed a significant increase in the ciliated epithelium (as determined by Ac Tub; compare Ac Tub stained panels in Fig 3E with 3F) . This increase was associated with a concomitant decrease in both KRT5-stained basal and SFTPC-positive AT2 cells (Fig 3F) .

Taken together, the immunofluorescence images are in agreement with the RNA seq dataset; both demonstrate that the short-term submerged monolayer favors distal differentiation (AT2→AT1), whereas the 21day ALI model favors proximal mucociliary differentiation. It is noteworthy that these distinct differentiation phenotypes originated from the same 3D-organoids despite the seeding of cells in the same basic media composition (i.e., PneumaCult TM ) prior to switching over to an ALI-maintenance media for the prolonged growth at Air-Liquid interface; the latter is a well-described methodology that promotes differentiation into ciliated and goblet cells (Rayner et al., 2019) .

Because the lung organoids with complete proximodistal cellularity could be differentiated into either proximalpredominant monolayers in submerged short-term cultures or distal-predominant monolayers in long-term ALI cultures, this provided us with an opportunity to model the respiratory tract and assess the impact of the virus along the entire proximal-to-distal gradient. We first asked how viral infectivity varies in the two models. Because multiple groups have shown the importance of the ciliated airway cells for infectivity (i.e., viral entry, replication and apical release (Hou et al., 2020; Hui et al., 2020; Milewska et al., 2020; Zhu et al., 2020) ), as positive controls, we infected monolayers of human airway epithelia (see Fig S4F-I) . AT2 cells, which express high levels of viral entry receptors ACE2 and TMPRSS2 (Fig 1A, Fig S1A) have been shown to be proficient in viral entry, but are least amenable to sustained viral release and infectivity (Hou et al., 2020; Hui et al., 2020) . To this end, we infected monolayers of hiPSC-derived homogeneous cultures of AT2 cells as secondary controls (see Fig S4J- L). Infection was carried out using the Washington strain of SARS-CoV-2, USA-WA1/2020 [BEI Resources NR-52281 (Rogers et al., 2020) ]. As expected, the 2D-lung monolayers we generated, both the submerged or the ALI models, were readily infected with SARS-CoV-2 (Fig S5A) , as determined by the presence of the viral envelope gene (E-gene; Fig 3G) ; however, the kinetics of viral amplification differed. When expressed as levels of E gene normalized to the peak values in each model (Fig 3G) , the kinetics of the ALI-monolayer model mirrored that of the primary airway epithelial monolayers; both showed slow beginning (0 -48 hpi) followed by an exponential increase in E gene levels from 48 to 72 hpi. The submerged monolayer model showed sustained viral infection during the 48-72 hpi window (Fig S5A; left) . In the case of AT2 cells, the 48-72 hpi window was notably missing in monolayers of hiPSC-derived AT2 cells (Fig 3G; Fig S5A; right). When we specifically analyzed the kinetics of viral E gene expression during the late phase (48-72 hpi window), we found that proximal airway models [human Bronchial airway Epi (HBEpC)] were more permissive than distal models [human Small Airway Epi (HSAEpC) and AT2] to viral replication (Fig S5B) ; the ALO monolayers showed intermediate sustained infectivity (albeit with variability). All models showed extensive cell death and detachment by 96 hrs and, hence, were not analyzed. Confocal imaging of infected ALO monolayers with anti-SARS-COV-2 nucleocapsid protein antibody showed that submerged ALO monolayers did indeed show progressive changes during the 48 to 72h window after infection ( Taken together, these findings show that sustained viral infectivity is best simulated in monolayers that resemble the proximal mucociliary epithelium, i.e., 2D-monolayers of lung organoids grown as ALI models and the primary airway epithelia. Because prior studies conducted in patient-derived airway cells (Hou et al., 2020) mirror what we see in our monolayers, we conclude that proximal airway cells within our mixed-cellular model appear to be sufficient to model viral infectivity in COVID-19.

Next, we asked if the newly generated lung models accurately recapitulate the host immune response in COVID-19. To this end, we analyzed the infected ALO monolayers (both the submerged and ALI variants) as well as the airway epithelial (HSAEpC) and AT2 monolayers by RNA seq and compared them all against the transcriptome profile of lungs from deceased COVID-19 patients. A publicly available dataset (GSE151764) (Nienhold et al., 2020) , comprised of lung transcriptomes from victims deceased either due to non-infectious causes (controls) or due to COVID-19, was first analyzed for differentially expressed genes (Fig 4A-B) . This cohort was chosen as a test cohort over others because it was the largest one available at the time of this study with appropriate postmortem control samples. Differentially expressed genes showed an immunophenotype that was consistent with what is expected in viral infections (Fig 4C; Table 4; Fig S6) , and showed overrepresentation of pathways such as interferon, immune, and cytokine signaling (Fig 4D; Table 5 ; Fig S7) . Differentially expressed gene signatures and reactome pathways that were enriched in the test cohort were fairly representative of the host immune response observed in patient-derived respiratory samples in multiple other validation cohorts; the signature derived from the test cohort could consistently classify control (normal) samples from COVID-19samples (ROC AUC 0.89 to 1.00 across the board; Fig 4E) . The most notable finding is that the patient-derived signature was able to perfectly classify the EpCAM-sorted epithelial fractions from the bronchoalveolar lavage fluids of infected and healthy subjects (ROC AUC 1.00; GSE145926-Epithelium (Liao et al., 2020) , suggesting that the respiratory epithelium is a major site where the host immune response is detected in COVID-19. When compared to existing organoid models of COVID-19, we found that the patient-derived COVID-19-lung signature was able to perfectly classify infected vs. uninfected late passages (> 50) of hiPSC-derived AT1/2 monolayers (GSE155241) (Han et al., 2020) and infected vs. uninfected liver and pancreatic organoids (Fig 4F) . The COVID-19-lung signatures failed to classify commonly used respiratory models, e.g., A549 cells and bronchial organoids, as well as intestinal organoids (Fig 4F) . A similar analysis on our own lung models revealed that the COVID-19lung signature was induced in submerged monolayers with distal-predominant AT2→AT1 differentiation, but not in the proximal-predominant ALI model (ROC AUC 1.00 and 0.50, respectively; Fig 4G) . The ALI model and the small airway epithelia, both models that mimic the airway epithelia (and lack alveolar pneumocytes; see Fig 3B) , failed to mount the patient-derived immune signatures (Fig 4H; left) . These findings suggested that the presence of alveolar pneumocytes is critical for emulating host response. To our surprise, induction of the COVID-19-lung signature also failed in hiPSC-derived AT2 monolayers (Fig 4H; right), indicating that AT2 cells are unlikely to be the source of such host response. These findings indicate that both proximal airway and AT2 cells, when alone, are insufficient to induce the host immune response that is encountered in the lungs of COVID-19 patient.

Next, we analyzed the datasets from our ALO monolayers for differentially expressed genes when challenged with SARS-COV-2 (Fig 5A-B) . Genes and pathways upregulated in the infected lung organoidderived monolayer models overlapped significantly with those that were upregulated in the COVID-19 lung missing components, the model-derived DEG signature was sufficient to consistently and accurately classify diverse cohorts of patient-derived respiratory samples (ROC AUC ranging from 0.88 to 1.00; Fig 5F) ; the modelderived DEG signature was significantly induced in COVID-19 samples compared to normal controls (Fig 5G-H ). Most importantly, the model-derived DEG signature was significantly induced in the epithelial cells recovered from bronchoalveolar lavage (Fig 5I) .

Taken together, these cross-validation studies from disease to model (Fig 4) and vice versa (Fig 5) provide an objective assessment of the match between the host response in COVID-19 lungs and our submerged ALO monolayers. Such a match was not seen in the case of the other models, e.g., the proximal airway-mimic ALI model, HSAEpC monolayer, or hiPSC-derived AT2 models. Because the submerged ALO monolayers contained both proximal airway epithelia (basal cells) and promoted AT2→AT1 differentiation, findings demonstrate that mixed cellular monolayers can mimic the host response in COVID-19. A subtractive analysis revealed that the cell type that is shared between models which showed induction of host response signatures (i.e., ALO submerged monolayers and GSE155241 (Han et al., 2020); Fig 5F) but is absent in models that do not show such response (hu Bronchial organoids, small airway epi, ALI-model of ALO) is AT1. We conclude that distal differentiation from AT2→AT1, a complex process that is comprised of distinct intermediates (Choi et al., 2020), is essential for modeling the host immune response in COVID-19.

Both proximal and distal airway epithelia are required to mount the overzealous host response in We next asked which model best simulated the overzealous host immune response that has been widely implicated in fatal COVID-19. To this end, we relied upon a recently described artificial intelligence (AI)-guided definition of the nature of the overzealous response in fatal COVID-19 (Sahoo et al., 2020) . Using ACE2 as a seed gene, a 166-gene signature was identified and validated as an invariant immune response that was shared among all respiratory viral pandemics, including COVID-19 (Fig 6A) . A subset of 20 genes within the 166-gene signature was subsequently identified as a determinant of disease severity/fatality; these 20 genes represented translational arrest, senescence, and apoptosis. These two signatures referred to as ViP (166-gene) and severe ViP (20-gene) signatures, were used as a computational framework to first vet existing SARS-CoV-2 infection models that have been commonly used for therapeutic screens (Fig 6B-D) . Surprisingly, we found that each model fell short in one way or another. For example, the Vero E6, which is a commonly used cultured cell model, showed a completely opposite response; instead of being induced, both the 166-gene and 20-gene ViP signatures were suppressed in infected Vero E6 monolayers (Fig 6B) . Similarly, neither ViP signature was induced in the case of SARS-CoV-2 challenged human bronchial organoids(Suzuki et al., 2020) ( Fig 6C) . Finally, in the case of the hiPSC-derived AT1/2 organoids, which recapitulated the COVID-19-lung derived immune signatures (in Fig 4F) , the 166-gene ViP signature was induced significantly (Fig 6D; top) , but the 20-gene severity signature was not (Fig 6D; bottom) . These findings show that none of the existing models capture the overzealous host immune response that has been implicated in a fatality.

Our lung models showed that both the 166-and 20-gene ViP signatures were induced significantly in the submerged ALO-derived monolayers that had distal differentiation (Fig 6E; left), but not in the proximal-mimic ALI model (Fig 6E; right) . Neither signatures were induced in monolayers of small airway epithelial cells (Fig 6F) or hiPSC-derived AT2 cells (Fig 6G) .

Taken together with our infectivity analyses, these findings demonstrate that although the proximal airway epithelia and AT2 cells may be infected, and as described by others (Dye et al., 2015; Hou et al., 2020) , may be vital for mounting a viral response and for disease transmission, these cells alone cannot mount the overzealous host immune response that is associated with the fatal disease. Similarly, even though the alveolar pneumocytes, AT1 and AT2 cells, are sufficient to mount the host immune response, in the absence of proximal airway components, they too are insufficient to recapitulate the severe ViP signature that is characterized by cellular senescence and apoptosis. However, when both proximal and distal components are present, i.e., basal, ciliated and AT1 cells, the model mimicked the overzealous host immune response in COVID-19 (Fig 6H) .

The most important discovery we report here is the creation of adult lung organoids that are complete with both proximal airway and distal alveolar epithelia; these organoids can not only be stably propagated and expanded in 3D cultures but also used as monolayers of mixed cellularity for modeling viral and host immune responses during respiratory viral pandemics. Furthermore, an objective analysis of this model and other existing SARS-CoV-2-infected lung models against patient-lung derived transcriptomes showed that the model which most closely emulates the elements of viral infectivity, lung injury, and inflammation in COVID-19 is one that contained both proximal and distal alveolar signatures (Fig 6H) , whereas, the presence of just one or the other fell short.

There are three important impacts of this work. First, successful modeling of the human lung organoids that are complete with both proximal and distal signatures has not been accomplished before. The multicellularity of the lung has been a daunting challenge that many experts have tried to recreate in vitro; in fact, the demand for perfecting such a model has always remained high, not just in the current context of the COVID-19 pandemic but also with the potential of future pandemics. We have provided the evidence that the organoids that were created using our methodology retain proximal and distal cellularity throughout multiple passages and even within the same organoid. Although a systematic design of experiment (DoE) approach (Bukys et al., 2020) was not involved in getting to this desirable goal, a rationalized approach was taken. For example, a Wnt/Rspondin/Noggin-containing conditioned media was used as a source of the so-called 'niche factors' for any organoid growth (Sato and Clevers, 2015) . This was supplemented with recombinant FGF7/10; FGF7 is known to help cell proliferation and differentiation and is required for normal branching morphogenesis (Padela et al., 2008) , whereas FGF10 helps in cell maturation (Rabata et al., 2020) and in alveolar regeneration upon injury (Yuan et al., 2019) . Together, they are likely to have directed the differentiation toward distal lung lineages (hence, the preservation of alveolar signatures). The presence of both distal alveolar and proximal ciliated cells was critical: proximal cells were required to recreate sustained viral infectivity, and the distal alveolar pneumocytes, in particular, the ability of AT2 cells to differentiate into AT1 pneumocytes was essential to recreate the host response. It is possible that the response is mediated by a distinct AT2-lineage population, i.e., damageassociated transient progenitors (DATPs), which arise as intermediates during AT2→AT1 differentiation upon injury-induced alveolar regeneration (Choi et al., 2020) . Although somewhat unexpected, the role of AT1 pneumocytes in mounting innate immune responses has been documented before in the context of bacterial pneumonia (Wong and Johnson, 2013; Yamamoto et al., 2012) . In work (Huang et al., 2020) that was published during the preparation of this manuscript, authors used long-term ALI models of hiPSC-derived AT2 monolayers (in growth conditions that inhibit AT2→AT1 differentiation, as we did here for our AT2 model) and showed that SARS-CoV-2 induces iAT2-intrinsic cytotoxicity and inflammatory response, but failed to induce type 1 interferon pathways (IFN α/β). It is possible that prolonged culture of iAT2 pneumocytes gives rise to some DATPs but cannot robustly do so in the presence of inhibitors of AT1 differentiation. This (spatially segregated viral and host immune response) is a common theme across many lung infections (including bacterial pneumonia and other viral pandemics (Chan et al., 2013; Hou et al., 2020; Taubenberger and Morens, 2008; Weinheimer et al., 2012) and hence, this mixed cellularity model is appropriate for use in modeling diverse lung infections and respiratory pandemics to come.

Second, among all the established lung models so far, ours features 4 key properties that are desirable whenever disease models are being considered for their use in HTP modes for rapid screening of candidate therapeutics and vaccines -(i) reproducibility, propagability and scalability, (ii) cost-effectiveness, (iii) personalization, and (iv) modularity, with the potential to be used in multi-dimensional complex lung models with other cell types if/when needed. We showed that the protocol we have optimized supports isolation, expansion and propagability at least up to 12-15 passages (at the time of submission of this work), with documented retention of proximal and distal airway components up to P8 (by RNA seq). Feasibility has also been established for scaling up for use in 384-well HTP assays. We also showed that the protocols for generating lung organoids could be reproduced in both genders and regardless of the donor's smoking status, consistency in outcome and growth characteristics were observed across all isolation attempts. The ALOs are also cost-effective; the need for exclusive reliance on recombinant growth factors was replaced at least in part with conditioned media from a commonly used cell line (L-WRN/ ATCC® CRL-2647 cells). Such media has a batch to batch stable cocktail of Wnt, R-Spondin, and Noggin and has been shown to improve reproducibility in the context of GI organoids in independent laboratories (VanDussen et al., 2019) . The use of conditioned media may not only have made our model more cost-effective than others but also likely improved the rigor and reproducibility in lung modeling and research. Because the model is propagable, repeated iPSC-reprogramming (another expensive step) is also eliminated, further cutting costs compared to many other models. As for personalization, our model is derived from adult lung stem cells from deep lung biopsies; each organoid line was established from one patient. By avoiding iPSCs or EPSCs, this model not only captures genetics but also retains organ-specific epigenetic programming in the lung, and hence, potentially additional programming that may occur in disease (such as in the setting of chronic infection, injury, inflammation, somatic mutations, etc.). The ability to replicate donor phenotype and genotype in vitro allows for potential use as pre-clinical human models for Phase '0' clinical trials.

As for modularity, by showing that the 3D lung organoids could be used as polarized monolayers on transwells to allow infectious agents to access the apical surface (in this case, SARS-CoV-2), we demonstrate that the organoids have the potential to be reverse-engineered with additional components in a physiologically relevant spatially segregated manner: for example, immune and stromal cells can be placed in the lower chamber to model complex lung diseases that are yet to be modeled and have no cure (e.g., Idiopathic pulmonary fibrosis, etc).

Third, the value of the ALO models is further enhanced due to the availability of companion readouts/ biomarkers (e.g., ViP signatures in the case of respiratory viral pandemics, or monitoring the E gene, or viral shedding, etc.) that can rapidly and objectively vet treatment efficacy based on set therapeutic goals. Of these readouts, the host response, as assessed by ViP signatures, is a key vantage point because an overzealous host response is what is known to cause fatality. Recently, a systematic review of the existing pre-clinical animal models revealed that most of the animal models of COVID-19 recapitulated mild patterns of human COVID-19; no severe illness associated with mortality was observed, suggesting a wide gap between COVID-19 in humans (Spagnolo et al., 2020) and animal models (Ehaideb et al., 2020) . The model revealed here, in conjunction with the ViP signatures described earlier (Sahoo et al., 2020) , could serve as a pre-clinical model with companion diagnostics to identify drugs that target both the viral and host response in pandemics.

Milewska, A., Kula-Pacurar, A., Wadas, J., Suder, A., Szczepanski, A., Dabrowska, A., Owczarek, K., Marcello, A., Ochman, M., Stacel, T., et al. (2020) . Replication of Severe Acute Respiratory Syndrome Coronavirus 2 in Human Respiratory Epithelium. J Virol 94. Sayed, I.M., Suarez, K., Lim, E., Singh, S., Pereira, M., Ibeawuchi, S.R., Katkar, G., Dunkel, Y., Mittal, Y., Chattopadhyay, R., et al. (2020) . Host engulfment pathway controls inflammation in inflammatory bowel disease. Zhu, N., Wang, W., Liu, Z., Liang, C., Wang, W., Ye, F., Huang, B., Zhao, L., Wang, H., Zhou, W., et al. (2020) . Morphogenesis and cytopathic effect of SARS-CoV-2 infection in human airway epithelial cells. Nat Commun 11, 3910.

Whisker plots display relative levels of ACE2 expression in various cell types in the normal human lung. The cell types were annotated within a publicly available single-cell sequencing dataset (GSE132914) using genes listed in Table S1 . p values were analyzed by one-way Anova and Tukey post hoc test. B. FFPE sections of the human lung from normal and deceased COVID-19 patients were stained for SFTPC, alone or in combination with 

To generate adult healthy lung organoids, fresh biopsy bites were prospectively collected after surgical resection of the lung by the cardiothoracic surgeon. Before collection of the lung specimens, each tissue was sent to a gross anatomy room where a pathologist cataloged the area of focus, and the extra specimens were routed to the research lab in Human Transport Media (HTM, Advanced DMEM/F-12, 10 mM HEPES, 1X Glutamax, 1X

penicillin-streptomycin, 5 M Y-27632) for cell isolation. Deidentified lung tissues obtained during surgical resection, that were deemed excess by clinical pathologists, were collected using an approved human research protocol (IRB# 101590; PI: Thistlethwaite). Isolation and biobanking of organoids from these lung tissues were carried out using an approved human research protocol (IRB# 190105: PI Ghosh and Das) that covers human subject research at the UC San Diego HUMANOID Center of Research Excellence (CoRE). For all the deidentified human subjects, information including age, gender, and previous history of the disease, was collected from the chart following the rules of HIPAA and described in the Table. A portion of the same lung tissue specimen was fixed in 10% Zinc-Formalin for at least 24hrs followed by submersion in 70% EtOH until embedding in FFPE blocks.

The lung specimens from COVID-19 positive human subjects were collected through autopsy (the study was IRB Exempt). All donations to this trial were obtained after telephone consent followed by written email confirmation by the next of kin/power of attorney per California state law (no in-person visitation could be allowed into our COVID-19 ICU during the pandemic). The team member followed the CDC guidelines for COVID19 and the autopsy procedures ((CAP), 2020; (CDC), 2020)). Lung specimens were collected in 10% Zinc-formalin and stored for 72 hrs before processing for histology. Patient characteristics are listed in the Table. Autopsy #2 was a standard autopsy performed by anatomical pathology in the BSL3 autopsy suite. The patient expired and his family consented for autopsy. After 48 hrs, the lungs were removed and immersion fixed whole in 10% formalin for 48 hrs and then processed further. Lungs were only partially fixed at this time (about 50% fixed in thicker segments) and were sectioned into small 2-4 cm chunks and immersed in 10% formalin for further investigation.

Autopsy #4 and #5 were collected from rapid post-mortem lung biopsies. The procedure was performed in the Jacobs Medical Center ICU (all of the ICU rooms have a pressure-negative environment, with air exhausted through HEPA filters [Biosafety Level 3 (BSL3)] for isolation of SARS-CoV-2 virus). Biopsies were performed 2-4 hrs after patient expiration. The ventilator was shut off to reduce the aerosolization of viral particles at least 1 hr after the loss of pulse and before sample collection. Every team member had personal protective equipment in accordance with the University policies for procedures on patients with COVID-19 (N95 mask + surgical mask, hairnet, full face shield, surgical gowns, double surgical gloves, booties). Lung biopsies were obtained after L-thoracotomy in the 5th intercostal space by the cardiothoracic surgery team. Samples were taken from the left upper lobe (LUL) and left lower lobe (LLL) and then sectioned further.

A previously published protocol was modified to isolate lung organoids from 3 human subjects (Sachs et al., 2019; Zhou et al., 2018) . Briefly, normal human lung specimens were washed with PBS/4X penicillinstreptomycin and minced with surgical scissors. Tissue fragments were resuspended in 10 mL of wash buffer (Advanced DMEM/F-12, 10 mM HEPES, 1X Glutamax, 1X penicillin-streptomycin) containing 2 mg/ml Collagenase Type I (Thermo Fisher, USA) and incubated at 37°C for approximately 1 hr. During incubation, tissue pieces were sheared every 10 min with a 10 mL serological pipette and examined under a light microscope to monitor the progress of digestion. When 80-100% of single cells were released from connective tissue, the digestion buffer was neutralized with 10 mL wash buffer with added 2% Fetal Bovine Serum; the suspension was passed through a 100-µm cell strainer and centrifuged at 200 rcf. Remaining erythrocytes were lysed in 2 ml red blood cell lysis buffer (Invitrogen) at room temperature for 5 min, followed by the addition of 10 mL of wash buffer and centrifugation at 200 rcf. Cell pellets were resuspended in cold Matrigel (Corning, USA) and seeded in 25 µl droplets on a 12 well tissue culture plate. The plate was inverted and incubated at 37°C for 10 min to allow complete polymerization of the Matrigel before the addition of 1 mL Lung Expansion Medium per well. Lung expansion media was prepared by modifying the GI-organoid media (50% conditioned media, prepared from L-WRN cells with Wnt3a, R-spondin, and Noggin, ATCC-CRL-3276 TM ) (Ghosh et al., 2020; Sayed et al., 2020a; Sayed et al., 2020b; Sayed et al., 2020c) with a proprietary cocktail from the HUMANOID CoRE containing B27, TGF-β receptor inhibitor, antioxidants, p38 MAPK inhibitor, FGF 7, FGF 10 and ROCK inhibitor.

The lung expansion media was compared to alveolosphere media I (IMDM and F12 as the basal medium with B27, low concentration of KGF, Monothioglycerol, GSK3 inhibitor, Ascorbic acid, Dexamethasone, IBMX, cAMP and ROCK inhibitor) and II (F12 as the basal medium with added CaCl 2 , B27, low concentration of KGF, GSK3

inhibitor, TGF-β receptor inhibitor Dexamethasone, IBMX, cAMP and ROCK inhibitor) modified from previously published literature (Jacob et al., 2019; Yamamoto et al., 2017) . Neither alvelosphere media contain any added Wnt3a, R-spondin, and Noggin. The composition of these media was developed either by fundamentals of adultstem cell-derived mixed epithelial cellularity in other organs (like the gastrointestinal tract (Miyoshi and Stappenbeck, 2013; Sato et al., 2009; Sayed et al., 2020c) , or rationalized based on published growth conditions for proximal and distal airway components (Gotoh et al., 2014; Sachs et al., 2019; van der Vaart and Clevers, 2020) . Organoids were maintained in a humidified incubator at 37°C/5% CO2, with a complete media change performed every 3 days. After the organoids reached confluency between 7-10 days, organoids were collected in PBS/0.5 mM EDTA and centrifuged at 200 rcf for 5 min. Organoids were dissociated in 1 mL trypLE Select (Gibco, USA) per well at 37°C for 4-5 min and mechanically sheared. Wash buffer was added at a 1:5, trypLE to wash buffer ratio. The cell suspension was subsequently centrifuged, resuspended in Matrigel, and seeded at a 1:5 ratio. Lung organoids were biobanked and passage 3-8 cells were used for experiments. Subculture was performed every 7-10 days.

Lung-organoid-derived monolayers were prepared using a modified protocol of GI-organoid-derived monolayers (Ghosh et al., 2020; Sayed et al., 2020a; Sayed et al., 2020b; Sayed et al., 2020c) . Briefly, transwell inserts (6.5 mm diameter, 0.4 um pore size, Corning) were coated in Matrigel diluted in cold PBS at a 1:40 ratio and incubated for 1hr at room temperature. Confluent organoids were collected in PBS/EDTA on day 7 and dissociated into single cells in trypLE for 6-7 min at 37°C. Following enzymatic digestion, the cell suspension was mechanically sheared through vigorous pipetting with a 1000 µl pipette and neutralized with wash buffer. The suspension was Human iPSC-derived alveolar epithelial type 2 cells (iHAEpC2) were obtained from Cell Applications Inc. and cultured in growth media (i536K-05, Cell Applications Inc.) according to the manufacturer's instructions. All the primary cells were used within early passages (5-6) to avoid any gradual disintegration of the airway epithelium with columnar epithelial structure and epithelial integrity.

Lung organoid-derived monolayers or primary airway epithelial cells either in 384 well plates or in transwells were washed twice with antibiotic-free lung wash media. 1E5 PFU of SARS-CoV-2 strain USA-WA1/2020 (BEI Resources NR-52281) in complete DMEM was added to the apical side of the transwell and allowed to incubate for 24, 48, 72 and 96 hrs at 34°C and 5% CO2. After incubation, the media was removed from the basal side of the transwell. The apical side of the transwells was then washed twice with (antibiotic-free lung wash media) and

then twice with PBS. TRIzol™ Reagent (Thermo Fisher 15596026) was added to the well and incubated at 34˚C and 5% CO2 for 10 min. The TRIzol™ Reagent was removed and stored at -80 ˚C for RNA analysis.

Organoids and monolayers used for lung cell type studies were lysed using RNA lysis buffer followed by RNA

CoV2 studies were lysed in TRI-Reagent and RNA was extracted using Zymo Research Direct-zol RNA Miniprep.

Organoid and monolayer cell-type gene expression was measured by qRT-PCR using 2x SYBR Green qPCR Master Mix. cDNA was amplified with gene-specific primer/probe set for the lung cell type markers and qScript cDNA SuperMix (5x). qRT-PCR was performed with the Applied Biosystems QuantStudio 5 Real-Time PCR System. Cycling parameters were as follows: 95 °C for 20 s, followed by 40 cycles of 1 s at 95 °C and 20 s at 60 °C. All samples were assayed in triplicate and eukaryotic 18S ribosomal RNA was used as a reference.

ITGA6 ITGA6 F 'CGAAACCAAGGTTCTGAGCCC' ITGA6 R 'CTTGGATCTCCACTGAGGCAGT' Goblet Muc5AC Muc5AC F ' GGAACTGTGGGGACAGCTCTT' Muc5AC R ' GTCACATTCCTCAGCGAGGTC' Cilia FoxJ1 FoxJ1 F 'ACTCGTATGCCACGCTCATCTG' FoxJ1 R 'GAGACAGGTTGTGGCGGATTGA' Club cell SCGB1A1 SCGB1A1 F ' CAAAAGCCCAGAGAAAGCATC' SCGB1A1 R ' CAGTTGGGGATCTTCAGCTTC' Alveolar type 1 AQP5, PDPN, P2RX4 AQP5 F ' TACGGTGTGGCACCGCTCAATG' AQP5 R ' AGTCAGTGGAGGCGAAGATGCA' PDPN F ' GTGCCGAAGATGATGTGGTGAC' PDPN R ' GGACTGTGCTTTCTGAAGTTGGC' P2RX4 F ' GTGGCGGATTATGTGATACCAGC' P2RX4 R ' CACACAGTGGTCGCATCTGGAA' Alveolar type 2 ABCA3, SFTPA1, SFTPC, LAMP3 ABCA3 F ' CTTGACAGTCGCAGAGCACCTT' ABCA3 R ' CTCCGTGAGTTCCACTTGTCCT' SFTPA1 F ' CACCTGGAGAAATGCCATGTCC' SFTPA1 R ' AAGTCGTGGAGTGTGGCTTGGA' SFTPC F ' GTCCTCATCGTCGTGGTGATTG' SFTPC R ' AGAAGGTGGCAGTGGTAACCAG' LAMP3 F ' TGGGAGCCTATTTGACCGTCTC' LAMP3 R ' GCTGACAACTGGAGGCTCTGTT'

Assessment of SARS-CoV-2 infectivity test was determined by qPCR using TaqMan assays and TaqMan

Universal PCR Master Mix as done before (Corman et al., 2020; Lamers et al., 2020) . cDNA was amplified with gene-specific primer/probe set for the E gene and QPCR was performed with the Applied Biosystems rRNA. Cycling parameters were as follows: 95 °C for 20 s, followed by 40 cycles of 1 s at 95 °C and 20 s at 60 °C. All samples were assayed in triplicate and gene eukaryotic 18S ribosomal RNA was used as a reference.

Organoids and Lung organoid-Derived Monolayers were fixed by either: (1) 4% PFA at room temperature for 30 min and quenched with 30 mM glycine for 5 min, (2) ice-cold 100% Methanol at -20°C for 20 min, (3) ice-cold 100% Methanol on ice for 20 min. Subsequently, samples were permeabilized and blocked for 2 hrs using an inhouse blocking buffer (2 mg/mL BSA and 0.1% Triton X-100 in PBS); as described previously (Lopez-Sanchez et al., 2014) . Primary antibodies were diluted in blocking buffer and allowed to incubate overnight at 4°C;

Secondary antibodies were diluted in blocking buffer and allowed to incubate for 2 hrs in the dark. Antibody dilutions are listed in the Supplementary Key Resource Table. ProLong Glass was used as a mounting medium.

#1 Thick Coverslips were applied to slides and sealed. Samples were stored at 4°C until imaged. FFPE embedded Organoid and Lung Tissue sections underwent antigen retrieval as previously described in methods for Immunohistochemistry staining. After antigen retrieval and cooling in DI water; samples were permeabilized and blocked in blocking buffer and treated as mentioned above for immunofluorescence. Images were acquired at room temperature with Leica TCS SPE confocal and with DMI4000 B microscope using the Leica LAS-AF Software. Images were taken with a 40× oil-immersion objective using 405-,488-, 561-nm laser lines for excitation. Z-stack images were acquired by successive Z-slices of 1µm in the desired confocal channels. Fields of view that were representative and/or of interest were determined by randomly imaging 3 different fields. Zslices of a Z-stack were overlaid to create maximum intensity projection images; all images were processed using FIJI (Image J) software.

Organoids were seeded on a layer of Matrigel in 6-Well plates and grown for 7-8 days. Once mature, organoids were fixed in 4% PFA at room temperature for 30 min and quenched with 30 mM glycine for 5 min. Organoids were gently washed with PBS and harvested using a cell scraper. Organoids were resuspended in PBS using wide-bore 1000 µL pipette tips. Organoids were stained using Gill's hematoxylin for 5 min for easier FFPE embedding and sectioning visualization. Hematoxylin stained organoids were gently washed in PBS and centrifuged and excess hematoxylin was aspirated. Organoids were resuspended in 65°C histogel and centrifuged at 65°C for 5 min. Histogel embedded organoid pellets were allowed to cool to room temperature and stored in 70% ethanol at 4°C until ready for FFPE embedding by LJI Histology Core. Successive FFPE embedded organoid sections were cut at a 4 µm thickness and fixed on to microscope slides.

For SARS CoV-nucleoprotein (np) antigen retrieval, slides were immersed in Tris-EDTA buffer (pH 9.0) and boiled for 10 min at 100°C inside a pressure cooker. Endogenous peroxidase activity was blocked by incubation with 3% H 2 O 2 for 10 minutes. To block non-specific protein binding 2.5% goat serum was added. Tissues were then incubated with a rabbit SARS CoV-NP antibody (Sino Biological, See Supplementary Key Resource Table) for 1.5 hrs at room temperature in a humidified chamber and then rinsed with TBS or PBS 3 times, for 5 min each. Sections were incubated with horse anti-rabbit IgG secondary antibodies for 30 min at room temperature and then washed with TBS or PBS 3 times for 5 min each. Sections were incubated with DAB and counterstained with hematoxylin for 30 sec. annotations were computed using tximport and biomaRt R package. A custom python script was used to organize the data and log reduced using log2(TPM+1). The raw data and processed data are deposited in Gene

Expression Omnibus under accession no GSE157055, and GSE157057.

Publicly available COVID-19 gene expression databases were downloaded from the National Center for Biotechnology Information (NCBI) Gene Expression Omnibus website (GEO) (Barrett et al., 2005; Barrett et al., 2013; Edgar et al., 2002) . If the dataset is not normalized, RMA (Robust Multichip Average) (Irizarry et al., 2003a; Irizarry et al., 2003b) is used for microarrays and TPM (Transcripts Per Millions) (Li and Dewey, 2011; Pachter, 2011) is used for RNASeq data for normalization. We used log2(TPM+1) to compute the final log-reduced expression values for RNASeq data. Accession numbers for these crowdsourced datasets are provided in the figures and manuscript. All of the above datasets were processed using the Hegemon data analysis framework (Dalerba et al., 2011; Dalerba et al., 2016; Volkmer et al., 2012) .

DESeq2 (Love et al., 2014) was applied to uninfected and infected samples to identify Up-and Down-regulated genes. A gene signature score is computed using both the Up-and Down-regulated genes which are used to order the sample. To compute the gene signature score, first, the genes present in this list were normalized according to a modified Z-score approach centered around StepMiner threshold (formula = (expr -SThr)/3*stddev (Fabregat et al., 2018) . Reactome identifies signaling and metabolic molecules and organizes their relations into biological pathways and processes. Violin, Swarm and Bubble plots are created using python seaborn package version 0.10.1.

Single Cell RNASeq data from GSE145926 was downloaded from GEO in the HDF5 Feature Barcode Matrix

Format. The filtered barcode data matrix was processed using Seurat v3 R package (Stuart et al., 2019) and scanpy python package (Wolf et al., 2018) . Pseudo bulk analysis of GSE145926 data was performed by adding counts from the different cell subtypes and normalized using log2(CPM+1). Epithelial cells were identified using SFTPA1, SFTPB, AGER, AQP4, SFTPC, SCGB3A2, KRT5, CYP2F1, CCDC153, and TPPP3 genes using SCINA algorithm (Zhang et al., 2019) . Pseudo bulk datasets were prepared by adding counts from the selected cells and normalized using log (CPM+1).

CIBERSORTx (https://cibersortx.stanford.edu/runcibersortx.php) was used for cell-type deconvolution of our dataset (which was normalized by CPM). As reference samples, we first used the Single Cell RNASeq dataset (GSE132914) from Gene Expression Omnibus (GEO). Next, we analyzed the bulk RNA seq datasets for the identification of cell types of interest using relevant gene markers (see Table 2 (ACE2,TMPRSS2) were identified using SCINA algorithm. Then, normalized pseudo counts were obtained with the CPM normalization method. The cell-type signature matrix was derived from the Single Cell RNASeq dataset, cell-types, and gene markers of interest. It was constructed by taking an average from gene expression levels for each of the markers across each cell type.

All experiments were repeated at least three times, and results were presented either as one representative experiment or as average ± S.E.M. Statistical significance between datasets with three or more experimental groups was determined using one-way ANOVA including a Tukey's test for multiple comparisons. For all tests, a P-value of 0.05 was used as the cutoff to determine significance (*P < 0.05, **P < 0.01, ***P < 0.001, and ****P < 0.0001). All experiments were repeated a least three times, and P-values are indicated in each figure. All statistical analyses were performed using GraphPad prism 6.1. A part of the statistical tests was performed using R version 3.2.3 (2015-12-10). Standard t-tests were performed using python scipy.stats.ttest_ind package (version 0.19.0). Marker used in this work for Immunofluorescence (IF)

Marker used in this work for qPCR

Obscure markers (Not a lot of research relative to lung) BRCA2  CYBB  KRT5  C1QB  FCGR1A  IL10  IL6  CD44  CD276   XAGE1B  CCR5  CCR2  ALOX15B  CMKLR1  MX1  TNFRSF17  CCR1  CXCR3   CDK1  GBP1  HLA-G  IDO1  ISG20  LAG3  MAD2L1  CXCL9  MKI67   SNAI2  IFITM1  GZMB  CD163  CD38  BST2  BUB1  CCL20  CCNB2   CXCL11  IFI27  IFI35  TDO2  GZMA  OAS3  POU2AF1  CXCL13  GNLY   DMBT1  DDX58  TNFAIP8  LAMP3  KIAA0101  MELK   SLAMF8  IL21  FOXM1  IFIH1  IFI6  PDCD1LG2   IFIT2  IFIT1  CXCL10  IRF4  PSMB9  CCL18   TNFSF18  ISG15  CDKN3  C1QA  OAS1  OAS2   IFIT3  TOP2A  LILRB1  HERC6  TNFSF13B  IFI44L  STAT1 

Collection and Submission of Post-mortem Specimens from Deceased Persons with Known or Suspected COVID-19

NCBI GEO: mining millions of expression profiles--database and tools

NCBI GEO: archive for functional genomics data sets--update

Detection of 2019 novel coronavirus (2019-nCoV) by real-time RT-PCR

Single-cell dissection of transcriptional heterogeneity in human colon tumors

CDX2 as a Prognostic Biomarker in Stage II and Stage III Colon Cancer

Gene Expression Omnibus: NCBI gene expression and hybridization array data repository

The Reactome Pathway Knowledgebase

The stress polarity signaling (SPS) pathway serves as a marker and a target in the leaky gut barrier: implications in aging and cancer

Generation of alveolar epithelial spheroids via isolated progenitor cells from human pluripotent stem cells

Summaries of Affymetrix GeneChip probe level data

Exploration, normalization, and summaries of high density oligonucleotide array probe level data

Derivation of self-renewing lung alveolar epithelial type II cells from human pluripotent stem cells

SARS-CoV-2 productively infects human gut enterocytes

RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome

GIV/Girdin is a central hub for profibrogenic signalling networks during liver fibrosis

Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2

In vitro expansion and genetic modification of gastrointestinal stem cells in spheroid culture

Models for transcript quantification from RNA-Seq

Long-term expanding human airway organoids for disease modeling

Single Lgr5 stem cells build crypt-villus structures in vitro without a mesenchymal niche

The DNA Glycosylase NEIL2 Suppresses Fusobacterium-Infection-Induced Inflammation and DNA Damage in Colonic Epithelial Cells

Helicobacter pylori infection downregulates the DNA glycosylase NEIL2, resulting in increased genome damage and inflammation in gastric epithelial cells

Host engulfment pathway controls inflammation in inflammatory bowel disease

Comprehensive Integration of Single-Cell Data

Three differentiation states risk-stratify bladder cancer into distinct subtypes

SCANPY: large-scale single-cell gene expression data analysis

Long-term expansion of alveolar stem cells derived from human iPS cells in organoids

SCINA: A Semi-Supervised Subtyping Algorithm of Single Cells and Bulk Samples

Differentiated human airway organoids to assess infectivity of emerging influenza virus

Name pValue FDR Interferon Signaling 1.1102230246251565e-16 1.1102230246251565e-14 Interferon alpha/beta signaling 1.1102230246251565e-16 1.1102230246251565e-14 Cytokine Signaling in Immune system 1.1102230246251565e-16 1.1102230246251565e-14 Immune System 1.1102230246251565e-16 1.1102230246251565e-14 Interleukin-10 signaling 9.845457782375888e-13 7.87636622590071e-11 Interferon gamma signaling 9.256706512417168e-12 6.109426298195331e-10 Chemokine receptors bind chemokines 1.0827139185209944e-10 6.171469335569668e-9 Signaling by Interleukins 6.813459352272844e-9 3.406729676136422e-7 Insulin-like Growth Factor-2 mRNA Binding Proteins (IGF2BPs/IMPs/VICKZs) bind RNA IFI35  EPSTI1  AMIGO2  IFITM2  SLC4A11  CMPK2  WARS1  FAAP100  APOL1  OASL  IFI27  ISG15  OAS3  IFI44L  CD14  SLC35F6  IFIT3  IFI44  SAMD9L  IFIT2 PARP9 SRP9P1 

A. Whisker plots display relative levels of ACE2 expression in various cell types in the normal human lung. The cell types were annotated within a publicly available single-cell sequencing dataset (GSE132914) using genes listed in Table S1 . A. Schematic lists the various markers used here for qPCR and immunofluorescence to confirm the presence of all cell types in the 3D lung organoids here and in 2D monolayers later (in Fig 3) . B-H. Bar graphs display the relative abundance of various cell type markers (normalized to 18S) in adult lung organoids (ALO), compared to the airway (NHBE) and/or alveolar (AT2) control cells, as appropriate. p-values were analyzed by one-way Anova. Error bars denote S.E.M; n = 3-6 datasets. I. 3D organoids grown in 8-well chamber slides were fixed, immunostained and visualized by confocal microscopy, as in Fig 2K. Scale bar = 50 µm. C. Monolayers of ALO1-3 were challenged with SARS-CoV-2 for indicated time points prior to fixation and staining for KRT5 (red) and viral nucleocapsid protein (green) and DAPI (blue; nuclei) and visualized by confocal microscopy. Representative images are shown, displaying various cytopathic effects. Scale bar = 15 µm.

Publicly available RNA Seq datasets (GSE151764) of lung autopsies from patients who were deceased due to COVID-19 or noninfectious causes (normal lung control) were analyzed for differential expression of genes and displayed as a heatmap. Figure 7 . Reactome pathway analysis of differentially expressed genes in lung autopsies (normal vs. COVID-19). Reactome-pathway analysis of the differentially expressed genes shows the major pathways upregulated in COVID-19-affected lungs. Top: Visualization as flattened (left) and hierarchical (right, insets) reactome. Bottom: Visualization of the same data as tables with statistical analysis indicative of the degree of pathway enrichment. Figure 8 . Differential expression analysis of RNA Seq datasets from lung organoid monolayers, infected or not, with SARS-CoV-2. Adult lung organoid monolayers infected or not with SARS-CoV-2 were analyzed by RNA seq and differential expression analysis. Differentially expressed genes are displayed as a heatmap. Figure 9 . Reactome pathway analysis of differentially expressed genes in lung organoid monolayers infected with SARS-CoV-2. Reactome-pathway analysis of the differentially expressed genes shows the major pathways upregulated in SARS-CoV-2-infected lung organoid monolayers. Top: Visualization as flattened (left) and hierarchical (right, insets) reactfoam. Bottom: Visualization of the same data as tables with statistical analysis indicative of the degree of pathway enrichment.