An airway organoid-based screen identifies a role for the HIF1α-glycolysis axis in SARS-CoV-2 infection.

It is urgent to develop disease models to dissect mechanisms regulating severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) infection. Here, we derive airway organoids from human pluripotent stem cells (hPSC-AOs). The hPSC-AOs, particularly ciliated-like cells, are permissive to SARS-CoV-2 infection. Using this platform, we perform a high content screen and identify GW6471, which blocks SARS-CoV-2 infection. GW6471 can also block infection of the B.1.351 SARS-CoV-2 variant. RNA sequencing (RNA-seq) analysis suggests that GW6471 blocks SARS-CoV-2 infection at least in part by inhibiting hypoxia inducible factor 1 subunit alpha (HIF1α), which is further validated by chemical inhibitor and genetic perturbation targeting HIF1α. Metabolic profiling identifies decreased rates of glycolysis upon GW6471 treatment, consistent with transcriptome profiling. Finally, xanthohumol, 5-(tetradecyloxy)-2-furoic acid, and ND-646, three compounds that suppress fatty acid biosynthesis, also block SARS-CoV-2 infection. Together, a high content screen coupled with transcriptome and metabolic profiling reveals a key role of the HIF1α-glycolysis axis in mediating SARS-CoV-2 infection of human airway epithelium.

Introduction

Severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) is the cause of the coronavirus disease 2019 (COVID-19) pandemic, among the worst global health care crises of our time. There is a great and immediate need to develop new and relevant models to study virus biology and virus-host interactions and perform drug screens. Compared with typically used Vero cells (a kidney epithelial cell line derived from the African green monkey) or transformed carcinoma cell lines, human organoids may more faithfully recapitulate key aspects of viral infection and virus-host interactions. Several organoid models have been developed to study SARS-CoV-2 infection, including lung organoids (Han et al., 2021; Katsura et al., 2020; Lamers et al., 2021; Pei et al., 2021; Salahudeen et al., 2020; Samuel et al., 2020; Tiwari et al., 2021), bronchial organoids (Fang et al., 2021), intestinal organoids (Lamers et al., 2020; Zhou et al., 2020), colonic organoids (COs; Han et al., 2021), blood vessel and kidney organoids (Monteil et al., 2020), brain organoids (Jacob et al., 2020; Ramani et al., 2020; Zhang et al., 2020), and liver organoids (Yang et al., 2020; Zhao et al., 2020).

A recent study suggested that angiotensin converting enzyme 2 (ACE2), the putative receptor of SARS-CoV-2, is highly expressed in airway ciliated cells (Ziegler et al., 2020). Over the past several years, protocols have been reported to direct human pluripotent stem cell (hPSC) differentiation into various airway lineages (Chen et al., 2017, 2019b; Dye et al., 2015; Hawkins et al., 2017, 2021; Huang et al., 2014, 2015; Hurley et al., 2020; McCauley et al., 2017; Montoro et al., 2018; Mou et al., 2012). Here, we developed an hPSC-derived airway organoid (hPSC-AO) platform for monitoring infection of SARS-CoV-2 and to screen for inhibitors. Comparing with adult organoids, hPSC-AOs facilitate scaling, which is needed for high-throughput screening of small molecules and metabolic profiling. Using hPSC-AOs, we identified GW6471, which blocks SARS-CoV-2 infection in the hPSC-AOs and hPSC-derived COs (hPSC-COs). Transcriptional and metabolic profiling reveals a key role for the HIF1α-glycolysis axis in mediating SARS-CoV-2 infection, which provides a target for anti-viral drug development.

Results

Single-cell RNA-seq characterization of hPSC-AOs

We differentiated hPSCs to hPSC-AOs incorporating minor modifications into previously reported stepwise strategies (Hawkins et al., 2021; Huang et al., 2014; McCauley et al., 2017). Immunostaining confirmed that the hPSC-AOs contain MUC5AC+ goblet cells and P63+ basal cells (Figures S1A and S1B). To fully characterize the cellular complexity of the hPSC-AOs, single-cell transcriptome profiles identified FOXJ1+ ciliated-like cells, MUC5B goblet-like cells, KRT5TP63KRT17 basal cells, and MKI67CDK1TOP2A proliferating cells (Figures 1A and 1B and S1C). The hPSC-derived ciliated-like cell population was enriched for adult human ciliated and proximal ciliated cell markers (Figure 1C). Correlation analysis of signature genes further validated identity of the hPSC-derived ciliated-like cell population in hPSC-AOs showing high similarity to adult human ciliated and proximal ciliated cells, but not alveolar epithelial type 1 or type 2 cells (Travaglini et al., 2020) (Figure 1D). Consistent with these data, the hPSC-derived ciliated-like population does not express markers for the alveolar epithelial type 2 cells, such as SFTPD, SFTPC, and NAPSA (Figure S1D). Both high-resolution phase contrast imaging (Figure 1E) and video (Video S1) validate the presence of multiple cilia in ciliated-like cells. The SARS-CoV-2 receptor, ACE2 (Hoffmann et al., 2020), is highly expressed in the ciliated-like cell cluster (Figure 1F). TMPRSS2 (Hoffmann et al., 2020) and CTSL (Ou et al., 2020), two key transmembrane proteases used for SARS-CoV-2 entry; FURIN, a pro-protein convertase that pre-activates SARS-CoV-2 (Shang et al., 2020); and NRP1, a host factor of SARS-CoV-2 infection (Cantuti-Castelvetri et al., 2020; Daly et al., 2020), are expressed more widely, but also enriched in the ciliated-like cell cluster (Figure 1F). Immunostaining analysis further validated that ACE2 is expressed in acetyl-α-tubulin (AATB)+FOXJ1+ ciliated-like cells (Figures 1G and 1H).

Figure 1

Single-cell RNA-seq analysis of hPSC-AOs

(A) Uniform manifold approximation and projection (UMAP) plot illustrating five cell clusters in the hPSC-AOs. n = 1 biological replicate.

(B) UMAP and violin plots showing the expression of genes FOXJ1, MUC5B, KRT5, and MKI67.

(C) Enrichment analysis of hPSC-AOs using genes highly expressed in adult human ciliated or proximal ciliated cells.

(D) Correlation analysis of genes with cell fates in hPSC-AOs and adult human lung cells.

(E) Phase contrast image of a representative ciliated-like cell. Scale bar, 30 μm.

(F) UMAP and violin plots showing the expression of SARS-CoV-2 entry factors, including ACE2, CTSL, FURIN, TMPRSS2, and NRP1.

(G) Representative confocal images of hPSC-AOs co-stained with antibodies recognizing ACE2 and cilia marker acetyl-α-tubulin (AATB). 4′,6-diamidino-2-phenylindole (DAPI) stains nuclei. Scale bar, 20 μm.

(H) Representative confocal images of hPSC-AOs co-stained with antibodies recognizing ACE2 and ciliated cell marker FOXJ1. DAPI stains nuclei. Scale bar, 100 μm.

See also Figure S1 and Video S1.

hPSC-AOs are permissive to SARS-CoV-2 infection

To determine their permissiveness, hPSC-AOs were infected with SARS-CoV-2 (USA-WA1/2020, MOI = 0.2). At 48 h post-infection (hpi), qRT-PCR analysis using primers targeting subgenomic N transcripts confirmed that significant viral replication could be detected at the RNA level in infected hPSC-AOs (Figure 2A). Immunostaining for SARS-CoV-2 nucleocapsid (SARS-N) protein (Figure 2B) confirmed robust SARS-CoV-2 infection of hPSC-AOs. The 3D confocal imaging further validated the presence of SARS-N in ciliated-like cells (Figure 2C; Video S2). In addition, the transcriptome profiling (Figure 2D) confirmed robust SARS-CoV-2 infection of hPSC-AOs. Moreover, plotting these RNA-sequencing (RNA-seq) datasets by principle-component analysis (PCA; Figure 2E) and clustering analysis (Figure 2F) demonstrated that the infected hPSC-AOs occupy a distinct transcriptional space compared with mock-infected hPSC-AOs. A volcano plot of SARS-CoV-2-infected hPSC-AOs compared with mock-treated organoids revealed the upregulation of SARS-CoV-2-associated transcripts and the robust induction of chemokines, including CXCL1, CXCL2, CXCL3, CXCL10, and CXCL11 (Figure 2G), as has previously been observed in SARS-CoV-2-infected lung tissues (Blanco-Melo et al., 2020; Han et al., 2021; Yang et al., 2020). Interestingly, interferon-pathway-associated genes, such as IFIT1, IFIT2, IFIT3, IFI44L, IFNB1, OAS2, and EGR1, are also upregulated in SARS-CoV-2-infected hPSC-AOs (Figure 2G).

Figure 2

hPSC-AOs are permissive to SARS-CoV-2 infection

(A) Relative SARS-CoV-2 RNA expression levels in hPSC-AOs at 48 hpi (MOI = 0.2). Total viral RNA from infected hPSC-AOs was analyzed by qRT-PCR for the presence of N transcripts relative to ACTB.

(B) Representative confocal images of hPSC-AOs at 48 hpi (MOI = 0.2) co-stained with antibodies recognizing SARS-CoV-2 nucleocapsid (SARS-N) protein and ciliated cell marker FOXJ1. Scale bar, 100 μm.

(C) Representative 3D confocal images of hPSC-AOs at 48 hpi (MOI = 0.2) co-stained with antibodies recognizing SARS-N and AATB. Scale bar, 10 μm.

(D) Read coverage on viral transcriptome in the mock- and SARS-CoV-2-infected hPSC-AOs (MOI = 0.2). Schematic shows the SARS-CoV-2 genome. Coverage is normalized per million reads.

(E and F) PCA (E) and sample clustering (F) on the mock- and SARS-CoV-2-infected hPSC-AOs.

(G) Volcano plot showing the gene expression changes between mock- and SARS-CoV-2-infected hPSC-AOs.

Data in (A) are presented as mean ± SEM (n = 3 biological replicates). The p values were calculated by unpaired two-tailed Student’s t test. ∗∗∗∗p < 0.0001.

See also Video S2.

A high content chemical screen using hPSC-AOs identifies GW6471 blocking SARS-CoV-2 infection

To identify small molecules capable of blocking SARS-CoV-2 virus infection, hPSC-AOs were deposited onto 384-well plates followed by addition of chemicals from a collection of agonists and antagonists of known signaling pathways and US Food and Drug Administration (FDA)-approved drugs. The compound information is shown in Table S1. Two hours post-treatment, organoids were infected with SARS-CoV-2 virus at MOI = 0.2. At 48 hpi, the organoids were fixed and analyzed for the percentage of SARS-N+ cells. The wells which the Z score was less than −2 were chosen as primary hit compounds (Figure 3A). Hit compounds were evaluated for efficacy and cytotoxicity at different concentrations (Figures 3B and S2A). One compound, GW6471, was confirmed to decrease the percentage of SARS-N+ cells through a dose-dependent manner, independent of cytotoxicity (half maximal effective concentration [EC50] = 2.1 μM; Figures 3B and 3C and S2B). qRT-PCR analysis showed a significant reduction of replicating virus in GW6471-treated hPSC-AOs (Figure 3D). Immunostaining confirmed a significant reduction of SARS-N+ cells detected among FOXJ1+ ciliated-like cells in hPSC-AOs treated with 10 μM GW6471 (Figures 3E and 3F). The compound was also evaluated for anti-SARS-CoV-2 activity post-infection. hPSC-AOs were infected with SARS-CoV-2 (MOI = 0.2). After 24 hpi, these hPSC-AOs were treated with 10 μM GW6471. At 48 hpi, both viral RNA (Figure 3G) and SARS-N+ cells (Figures 3H and 3I) were significantly decreased in the compound-treated hPSC-AOs. Finally, GW6471 was evaluated for its capacity to block infection of the B.1.351 SARS-CoV-2 variant. Both viral RNA (Figure 3J) and SARS-N+ cells (Figures 3K and 3L) were significantly decreased in the GW6471-treated hPSC-AOs.

Figure 3

An hPSC-AO-based high-throughput chemical screen identifies GW6471 that blocks SARS-CoV-2 infection

(A) Primary screening results. The x axis is the compound number. The y axis is the Z score. Red line indicates Z score = −2.

(B) Efficacy curve of GW6471. Data are presented as mean ± SD. n = 3 biological replicates.

(C) Representative confocal images of control or GW6471-treated hPSC-AOs at 48 hpi (MOI = 0.2). Scale bar, 100 μm.

(D) qRT-PCR analysis for viral N subgenomic RNA in hPSC-AOs, which were pretreated with control or 10 μM GW6471 at 48 hpi (MOI = 0.2). Data are presented as mean ± SEM. n = 3 biological replicates.

(E and F) Representative confocal images (E) and quantification (F) of SARS-N in hPSC-AOs, which were pretreated with control or 10 μM GW6471 at 48 hpi (MOI = 0.2). Scale bar, 100 μm. Data are presented as mean ± SEM. n = 6 biological replicates.

(G) Relative SARS-CoV-2 viral RNA expression levels at 48 hpi in hPSC-AOs infected with SARS-CoV-2 virus (MOI = 0.2) and 24 h later exposed to control or 10 μM GW6471 treatment. Data are presented as mean ± SEM. n = 3 biological replicates.

(H and I) Representative confocal images (H) and quantification (I) of SARS-N at 48 hpi of hPSC-AOs infected with SARS-CoV-2 virus (MOI = 0.2) and 24 h later exposed to control or 10 μM GW6471 treatment. Scale bar, 100 μm. Data are presented as mean ± SEM. n = 6 biological replicates.

(J) qRT-PCR analysis for viral N single-guide RNA (sgRNA) at 48 hpi of hPSC-AOs, which were pretreated with control or 10 μM GW6471 (variant B.1.351, MOI = 0.2). Data are presented as mean ± SEM. n = 3 biological replicates.

(K and L) Representative confocal images (K) and quantification (L) of SARS-N at 48 hpi of hPSC-AOs, which were pretreated with control or 10 μM GW6471 (variant B.1.351, MOI = 0.2). Scale bar, 100 μm. Data are presented as mean ± SEM. n = 6 biological replicates.

(M) qRT-PCR analysis for viral N sgRNA at 48 hpi of hPSC-COs, which were pretreated with control or 10 μM GW6471 (MOI = 3). Data are presented as mean ± SEM. n = 3 biological replicates.

(N and O) Representative confocal images (N) and quantification (O) at 48 hpi of SARS-N of hPSC-COs, which were pretreated with control or 10 μM GW6471 (MOI = 3). Scale bar, 100 μm. Data are presented as mean ± SEM. n = 6 biological replicates.

(P) Relative SARS-CoV-2 viral RNA expression levels at 48 hpi of hPSC-AOs infected with SARS-CoV-2 virus (MOI = 3) and 24 h later followed by control or 10 μM GW6471 treatment. Data are presented as mean ± SEM. n = 3 biological replicates.

(Q and R) Representative confocal images (Q) and quantification (R) of SARS-N at 48 hpi of hPSC-AOs infected with SARS-CoV-2 virus (MOI = 3) and 24 h later followed by control or 10 μM GW6471 treatment. Scale bar, 100 μm. Data are presented as mean ± SEM. n = 6 biological replicates.

The p values were calculated by unpaired two-tailed Student’s t test. ∗p < 0.05 and ∗∗∗p < 0.0001.

See also Figures S2 and S3 and Table S1.

Previously, we showed that hPSC-COs are permissive to SARS-CoV-2 infection (Han et al., 2021). We further validated GW6471 as an anti-SARS-CoV-2 compound using hPSC-COs. The hPSC-COs were pretreated with GW6471 prior to infection with SARS-CoV-2. A reduction of subgenomic viral RNA (Figure 3M) and a decrease of SARS-N+ cells were found in GW6471-treated hPSC-COs, which was similar to our observations with hPSC-AOs. (Figures 3N and 3O). To determine whether GW6471 has anti-viral activity post-infection, hPSC-COs were infected with SARS-CoV-2 (MOI = 3). At 24 hpi, hPSC-COs were treated with 10 μM GW6471. At 48 hpi, both viral RNA (Figure 3P) and the number of SARS-N+ cells (Figures 3Q and 3R) were significantly decreased in the GW6471-treated hPSC-COs.

To determine whether GW6471 affects viral particle stability, a concentrated SARS-CoV-2 virus preparation was pre-incubated with 10 μM GW6471 for 3 h. At 48 hpi, no significant difference was detected for the percentage of SARS-N+ cells comparing hPSC-AOs infected with SARS-CoV-2 that was pre-incubated with GW6471 or control (Figures S3A and S3B), suggesting that GW6471 does not affect the stability of viral particles. To determine the impact of GW6471 on viral entry, we used a luciferase-based pseudovirus as reported in our previous publication (Han et al., 2021). GW6471 does not block viral entry (Figure S3C). Consistent with this observation, ACE2 expression levels are not changed after GW6471 treatment (Figure S3D).

The HIF1α pathway and SARS-CoV-2 infection

To examine the mechanism of action for GW6471 in the context of our organoid-based SARS-CoV-2 infection model, RNA-seq analysis was carried out to compare the transcriptome profiles of control and GW6471-treated hPSC-AOs or hPSC-COs. PCA plots (Figure 4A) and clustering analysis (Figure 4B) showed that profiles of control and GW6471-treated hPSC-AOs clustered separately. Volcano plot (Figure 4C) and Ingenuity pathway analysis (IPA; Figure 4D; Table S2) highlighted downregulation of the HIF1α signaling pathway in treated cells. Likewise, RNA-seq data from control and GW6471-treated hPSC-COs had profiles that clustered separately, as indicated by PCA plots (Figure 4E) and clustering analysis (Figure 4F). Volcano plot (Figure 4G) and IPA (Figure 4H; Table S2) also highlighted downregulation of the HIF1α signaling pathway. Although several cell-cycle-associated pathways were identified in GW6471-treated hPSC-AOs, GW6471 treatment does not affect cell proliferation (Figures S3E–S3H). In the absence of SARS-CoV-2 infection, GW6471 does not affect chemokine expression (Figure S3I) or the interferon pathway (Figure S3J). In the presence of SARS-CoV-2 infection, GW6471 treatment showed a trend to decrease the expression of chemokines (Figure S3K), but no consistent trend for impacting genes involved in interferon signaling (Figure S3L).

Figure 4

GW6471 blocks SARS-CoV-2 infection by inhibiting the HIF1α pathway

(A and B) PCA (A) and sample clustering (B) on the control- and 10 μM GW6471-treated hPSC-AOs at 48 hpi (MOI = 0.2).

(C) Volcano plot showing the expression changes between control- and 10 μM GW6471-treated hPSC-AOs at 48 hpi (MOI = 0.2).

(D) Enriched pathways in GW6471- versus control-treated hPSC-AOs at 48 hpi.

(E and F) PCA (E) and sample clustering (F) on the control- and 10 μM GW6471-treated hPSC-COs at 48 hpi (MOI = 1).

(G) Volcano plot showing the expression changes between control- and 10 μM GW6471-treated hPSC-COs at 48 hpi (MOI = 1).

(H) Enriched pathways in GW6471- versus control-treated hPSC-COs at 48 hpi.

(I) qRT-PCR analysis for viral N sgRNA of hPSC-AOs treated with control or 1 μM chetomin at 48 hpi (MOI = 0.2). Data are presented as mean ± SEM. n = 3 biological replicates.

(J and K) Representative confocal images (J) and quantification (K) of SARS-N+ cells of hPSC-AOs, which were treated with control or 1 μM chetomin at 48 hpi (MOI = 0.2). Scale bar, 100 μm. Data are presented as mean ± SEM. n = 6 biological replicates.

(L) qRT-PCR analysis for viral N sgRNA of hPSC-AOs expressing shHIF1α or scrambled shRNA at 48 hpi (MOI = 0.2). Data was presented as mean ± SEM. n = 3 biological replicates.

(M and N) Representative confocal images (M) and quantification (N) of SARS-N+ cells of hPSC-AOs expressing shHIF1α or scrambled shRNA at 48 hpi (MOI = 0.2). Scale bar, 100 μm. Data are presented as mean ± SEM. n = 6 biological replicates.

The p values were calculated by unpaired two-tailed Student’s t test. ∗∗∗p < 0.001 and ∗∗∗∗p < 0.001.

See also Figures S3 and S4 and Tables S2 and S3.

To test whether decreased HIF1α signaling is functionally relevant to controlling infection, the compound chetomin, a known small molecule suppressor for transcriptional activation of the HIF1α pathway, was tested (Kung et al., 2004). Both viral RNA replication (Figure 4I) and the number of SARS-N+ cells (Figures 4J and 4K) were significantly decreased following SARS-CoV-2 infection of hPSC-AOs, when treated with 1 μM chetomin. Similar to GW6471, chetomin does not affect chemokine expression (Figure S4A) or the interferon pathway (Figure S4B) in the absence of SARS-CoV-2. In the presence of SARS-CoV-2 infection, chetomin treatment showed a trend to decrease chemokine expression levels (Figure S4C), but no consistent trend for impacting genes involved in interferon signaling (Figure S4D). We used short hairpin RNAs (shRNAs) to knockdown expression levels of HIF1α (Table S3). qRT-PCR analysis confirmed the knockdown efficiency and specificity (Figure S4C). hPSC-AOs carrying scrambled shRNA or shHIF1α were infected with SARS-CoV-2 (MOI = 0.2). At 48 hpi, both qRT-PCR and immunostaining data showed decreased viral infection in hPSC-AOs expressing shHIF1α (Figures 4L–4N).

Metabolic profiling identifies the role of glycolysis in SARS-CoV-2 infection

Metabolic profiling was performed to compare mock- versus SARS-CoV-2-infected hPSC-AOs and (following SARS-CoV-2 infection) control versus GW6471-treated hPSC-AOs. Profiles from SARS-CoV-2 and mock-infected hPSC-AOs clustered separately (Figure 5A). Fatty acids, such as oleic acid, palmitic acid, palmitoleic acid, and stearic acid, and amino acids, such as glycine, isoleucine, threonine, and serine, were all increased in SARS-CoV-2-infected hPSC-AOs (Figure 5B). In addition, citric acid and D-glucose-6-phosphate levels also increased in SARS-CoV-2-infected hPSC-AOs (Figure 5B). Profiles from hPSC-AOs that were GW6471 or control-treated following SARS-CoV-2 infection clustered separately (Figure 5C). Compared with controls, GW6471-treated hPSC-AOs showed decreased levels of fatty acids, including oleic acid, palmitic acid, palmitoleic acid, and amino acids, as well as D-glucose-6-phosphate and citric acid (Figure 5D). Transcriptome profiling confirmed that genes encoding enzymes in the glycolysis pathway are suppressed in GW6471-treated conditions (Figure 5E). This is consistent with downregulation of the glycolysis pathway in GW6471-treated hPSC-COs (Figure 4H).

Figure 5

Metabolic profiling identifies a key role for glycolysis in SARS-CoV-2 infection

(A and B) Hierarchical clustering analysis (A) and heatmap (B) of metabolic profiles for mock- or SARS-CoV-2 (MOI = 0.2)-infected hPSC-AOs at 48 hpi. n = 4 biological replicates for mock and N = 3 biological replicates for SARS-CoV-2.

(C and D) Hierarchical clustering analysis (C) and heatmap (D) of metabolic profiles for control- or 10 μM GW6471-treated hPSC-AOs at 48 hpi (MOI = 0.2). n = 3 biological replicates.

(E) Heatmap of genes encoding enzymes involved in the glycolytic pathway of control- or 10 μM GW6471-treated hPSC-AOs at 48 hpi (MOI = 0.2). n = 3 biological replicates.

HIF1α signaling controls glycolysis upon SARS-CoV-2 infection

To determine whether HIF1α signaling controls glycolysis in the hPSC-AOs, we compared transcript profiles of the control and chetomin-treated samples. PCA plots (Figure 6A) and clustering analysis (Figure 6B) showed that profiles of control and chetomin-treated hPSC-AOs clustered separately. Consistent with the results of GW6471-treated conditions (Figure 5E), genes encoding enzymes in the glycolysis pathway were suppressed in chetomin-treated conditions (Figure 6C). RNA-seq was also applied to analyze the hPSC-AOs expressing shHIF1α upon SARS-CoV-2 infection. Profiles from hPSC-AOs expressing scrambled RNA or shHIF1α clustered separately (Figures 6D and 6E). RNA-seq analysis confirmed downregulation of the HIF1α pathway in hPSC-AOs expressing shHIF1α (Figure 6F). Interestingly, the coronavirus replication pathway was also found to be downregulated in hPSC-AOs expressing shHIF1α upon SARS-CoV-2 infection (Figure 6F). Genes involved in the glycolysis pathway are also downregulated in hPSC-AOs expressing shHIF1α (Figure 6G; Table S4). The hPSC-AOs expressing shHIF1α show slightly decreased levels of ACE2 expression, but this difference seems unlikely to be sufficient to explain the significant decrease of SARS-CoV-2 infection in hPSC-AOs expressing shHIF1α (Figure S4F). Metabolic profiling also confirmed the decreased levels of fatty acids, including palmitic acid, oleic acid, and palmitoleic acid, as well as cholesterol, in hPSC-AOs expressing shHIF1α (Figures 6H and 6I). SARS-CoV-2-infected hPSC-AOs were next treated with xanthohumol, a prenylated flavonoid that suppresses fatty acid and cholesterol biosynthesis (Miyata et al., 2015) and the glycolysis pathway (Yuan et al., 2020). Both viral RNA levels (Figure 6J) and the number of SARS-N+ cells (Figures 6K and 6L) were significantly decreased in the 10 μM xanthohumol-treated hPSC-AOs. Finally, 5-(tetradecyloxy)-2-furoic acid (TOFA) (Parker et al., 1977) and ND-646 (Svensson et al., 2016), two inhibitors of acetyl-coenzyme A (CoA) carboxylase 1, the key enzyme involved in fatty acid de novo synthesis, were tested to further validate the role of fatty acid metabolism in SARS-CoV-2 infection. Both qRT-PCR (Figure 6M) and immunostaining (Figures 6N and 6O) confirmed that TOFA or ND-646 blocks SARS-CoV-2 infection in hPSC-AOs.

Figure 6

HIF1α regulates glycolysis in SARS-CoV-2-infected hPSC-AOs

(A and B) PCA (A) and sample clustering (B) analysis of control- and 1 μM chetomin-treated hPSC-AOs at 48 hpi (MOI = 0.2). n = 3 biological replicates.

(C) Heatmap of genes encoding enzymes involved in the glycolytic pathway of control- or 1 μM chetomin-treated hPSC-AOs at 48 hpi (MOI = 0.2).

(D–F) PCA (D), sample clustering (E), and IPA (F) analysis of the hPSC-AOs expressing shHIF1α or scrambled shRNA at 48 hpi (MOI = 0.2). n = 3 biological replicates.

(G) Heatmap of genes encoding enzymes involved in the glycolytic pathway in the hPSC-AOs expressing shHIF1α or scrambled shRNA at 48 hpi (MOI = 0.2).

(H and I) Hierarchical clustering analysis (H) and heatmap (I) of metabolic profiles of the hPSC-AOs expressing shHIF1α or scrambled shRNA at 48 hpi (MOI = 0.2). n = 3 biological replicates.

(J) qRT-PCR analysis for viral N sgRNA at 48 hpi of hPSC-AOs treated with control or 10 μM xanthohumol (MOI = 0.2). Data are presented as mean ± SEM. n = 3 biological replicates.

(K and L) Representative confocal images (K) and quantification (L) of SARS-N+ cells at 48 hpi of hPSC-AOs, which were treated with control or 10 μM xanthohumol (MOI = 0.2). Scale bar, 100 μm. Data are presented as mean ± SEM. n = 6 biological replicates.

(M) qRT-PCR analysis for viral N sgRNA at 48 hpi of hPSC-AOs treated with control or 3 μM TOFA or 3 μM ND-646 (MOI = 0.2). Data are presented as mean ± SEM. n = 3 biological replicates.

(N and O) Representative confocal images (N) and quantification (O) of SARS-N+ cells at 48 hpi of hPSC-AOs, which were treated with control or 3 μM TOFA or 3 μM ND-646 (MOI = 0.2). Scale bar, 100 μm. Data are presented as mean ± SEM. n = 6 biological replicates.

Data are presented as mean ± SEM. The p values were calculated by unpaired two-tailed Student’s t test. ∗p < 0.05 and ∗∗∗p < 0.001.

See also Figure S4 and Table S4.

Discussion

Morbidity and mortality from COVID-19 are largely attributed to the acute viral pneumonitis that evolves to acute respiratory distress syndrome. ACE2 is highly expressed in multiciliated cells of the human airway (Ziegler et al., 2020). We successfully developed an hPSC-AO platform, containing AATB+FOXJ1+ ciliated-like cells that express ACE2. RNA-seq analysis of infected organoids revealed upregulation of cytokine/chemokine signaling, which mimics the inflammatory changes observed in primary human COVID-19 pulmonary infections (Blanco-Melo et al., 2020). We adapted the hPSC-AOs to a high content platform to screen for chemicals that block SARS-CoV-2 infection. From a pilot screen, we identified GW6471 that blocks SARS-CoV-2 infection not only in hPSC-AOs but also in hPSC-COs, indicating that the anti-SARS-CoV-2 activity of GW6471 functions across different organ systems. Transcriptome profiles suggested that GW6471 blocks SARS-CoV-2 infection at least in part by inhibiting HIF1α, which was further validated using a known HIF1α inhibitor, chetomin. At a molecular level, chetomin disrupts the structure of the cysteine/histidine-rich 1 (CH1) domain of p300 and precludes its interaction with HIF, thereby attenuating hypoxia-inducible transcription (Kung et al., 2004). We further validated the role of HIF1α in SARS-CoV-2 infection using shRNAs. HIF1α has been shown to modulate virus replication in primary airway epithelial cells infected with respiratory syncytial virus (Morris et al., 2020). In addition, studies using monocytes showed that SARS-CoV-2 infection triggers mitochondrial reactive oxygen species production, which induced HIF1α stabilization and glycolysis. In our study, we found that GW6471 downregulates the HIF1 pathway in both hPSC-AOs (Figures 4A–4D) and hPSC-COs upon SARS-CoV-2 infection (Figures 4E–4H) and that knockdown of HIF1α (Figures 4I–4K) or use of a specific HIF1α inhibitor (Figures 4L and 4M) blocks SARS-CoV-2 infection. Bojkova et al. (2020) reported changes in the HIF1 pathway and carbon metabolism in SARS-CoV-2-infected Caco-2 cells. They further showed that blocking glycolysis with 2-deoxy-D-glucose, an inhibitor of hexokinase (the rate-limiting enzyme in glycolysis), prevented SARS-CoV-2 replication in Caco-2 cells (Bojkova et al., 2020). A recent proteo-transcriptomics analysis in SARS-CoV-2-infected Huh7 cells found a dose-dependent increase of serine/threonine kinase-mechanistic target of rapamycin kinase (Akt-mTOR) signaling, but surprisingly, a decrease in activity of HIF1α. However, while 68% of differentially expressed HIF1 pathway genes were downregulated, the other 32% were upregulated (Appelberg et al., 2020). In addition, hypoxemia is the typical characteristic in patients with severe COVID-19, which is consistent with our findings (Berenguer et al., 2020; Petrilli et al., 2020). Together, the studies suggest that the HIF1 pathway might play a cell type and context dependent role to support SARS-CoV-2 infection.

Metabolic profiling suggested decreased levels of fatty acids upon GW6471 treatment (Figure 5D) or knockdown of HIF1α (Figure 6I). GW6471 was previously reported as a PPARα antagonist (Xu et al., 2002). Saturated and unsaturated fatty acids are endogenous ligands for peroxisome proliferator activated receptor alpha (PPARα). Polyunsaturated fatty acids, ω6-polyunsaturated fatty acids, and saturated fatty acids have been shown to directly bind to PPARα. In addition, PPARα has been shown to play a critical role in hepatic fatty acid oxidation and ketogenesis based on analysis of PPARα mice (Kersten, 2014). Here, we found that GW6471-treated hPSC-AOs showed decreased levels of fatty acids compared with control-treated hPSC-AOs (Figure 5D), suggesting that there might be a cell-type-dependent feedback loop that controls fatty acid levels in hPSC-AOs. In hPSC-AOs, GW6471 inhibits the HIF1α-glycolysis axis. The pyruvate produced by glycolysis is an important intermediary in the conversion of carbohydrates into fatty acids and cholesterol, which might partially explain the decrease of fatty acids in GW6471-treated hPSC-AOs. Many viruses, including cowpea mosaic virus (Carette et al., 2000), poliovirus (Cherry et al., 2006; Guinea and Carrasco, 1990), tombusvirus (Sharma et al., 2010, 2011), and Semliki Forest virus (Perez et al., 1991), require fatty acids for their replication and assembly. Modulation of fatty acid metabolism has been shown to impact the replication of several viruses, including hepatitis C virus (Nasheri et al., 2013), and Old World alphaviruses (Bakhache et al., 2019). In addition, a preprint shows that orlistat, a drug that inhibits lipases and also fatty acid synthase, and triacsin C, which inhibits long chain acyl-CoA synthetases, block SARS-CoV-2 replication in Vero cells and Calu-3 cells (Silvas et al., 2020). Consistent with these observations, RNA-seq analysis identified downregulation of the coronavirus replication pathway in hPSC-AOs expressing shRNA against HIF1α (Figure 6F). Since fatty acids are also involved in viral particle formation (Herker et al., 2010) and other steps of viral cycle as discussed in several reviews (Chukkapalli et al., 2012; Heaton and Randall, 2011; Lorizate and Kräusslich, 2011), GW6471 might also affect multiple steps of the viral life cycle.

Finally, we performed metabolic profiling using SARS-CoV-2-infected hPSC-AOs and found increased levels of D-glucose-6 phosphate, citric acid, fatty acids, and amino acids, suggesting that infection stimulates glycolysis as part of a metabolic rewiring. Indeed, GW6471 treatment and knockdown of HIF1α, which block SARS-CoV-2 infection, lead to decreased levels of D-glucose-6 phosphate, citric acid, fatty acids, and amino acids. The data are consistent with the metabolic state being a key regulator of SARS-CoV-2 infection. The metabolic analysis of the sera of patients with COVID-19 identified increased circulating levels of free fatty acids (Appelberg et al., 2020), which is consistent with our data using hPSC-AOs. HIF1α is known to regulate glycolysis in immune cells, including monocytes (Appelberg et al., 2020), macrophages (Chen et al., 2019a), and dendritic cells (Liu et al., 2019). Consistent with these observations, glycolysis-associated genes are upregulated in SARS-CoV-2-infected hPSC-lung organoids, pancreatic endocrine cells, and lung autopsy samples from patients with COVID-19 (Figures S4G–S4I). Our data are consistent with a recent study showing that elevated glucose levels promote SARS-CoV-2 infection in human monocytes by regulating HIF1α-dependent glycolysis (Codo et al., 2020). Taken together, we demonstrated that inhibiting the HIF1α-glycolysis axis blocks SARS-CoV-2 infection, which highlights the potential to target this metabolic pathway for development of anti-viral drugs.

STAR★Methods

Key resources table

Resource availability

Lead contact

Further information and requests for resources and reagents should be directed to and will be fulfilled by the Lead Contact, Shuibing Chen (shc2034@med.cornell.edu)

Materials availability

This study did not generate new unique reagents.

Experimental model and subject details

Cell lines

The H1 and RUES2 hESCs were purchased from WiCell Research Institute Inc. The hESC H1 line was grown and maintained on 1% matrigel (Corning) coated six-well plates in homemade mTesR1 medium at 37°C under the 5% CO2 culture conditions. RUES2 hESCs were grown and maintained on 1% Matrigel-coated six-well plates in StemFlex medium (GIBCO) at 37°C with 5% CO2 culture condition. The medium was changed daily. When hESCs reached ∼90% confluence, the cells were passaged at 1:6-1:10 with ReLeSR (Stem Cell Technology).

293 T cells and Vero E6 cells were purchased from ATCC and maintained in DMEM (GIBCO) +10% FBS (GIBCO) at 37°C with 5% CO2 culture condition.

Viruses

SARS-CoV-2 variants, including USA-WA1/2020 and B.1.351, were obtained from the World Reference Center for Emerging Viruses and Arboviruses located at the University of Texas, Medical Branch via the CDC.

Method details

Airway organoid differentiation

hPSC-AOs were derived using a protocol slightly modified from previous studies (Huang et al., 2014; McCauley et al., 2017). For definitive endoderm (DE) differentiation, cells were passaged with Accutase (Innovative Cell Technology) at 1:2-1:3 for the first day. When achieving 80%–90% confluency, hESCs were treated with 3 μM CHIR99021 (CHIR, Sigma) and 100 ng/ml Activin A (R&D systems) in basal medium RPMI1640 (Cellgro) supplemented with 1X Pen-Strep (GIBCO) for 1 day, and changed to the basal medium containing 100 ng/ml Activin A and 2% FBS for the next 48 hours. To induce anterior foregut endoderm, the endoderm cells were cultured in complete serum free differentiation (cSFD) medium supplemented with 2 μM dorsomorphin dihydrochloride (R&D Systems) and 10 μM SB431542 (R&D Systems) for 24 hours, and then switched to 10 μM SB431542 and 1 μM IWP2 (R&D Systems) treatment for 24 hours. For induction of early stage lung progenitor cells (day 15), the anterior foregut endoderm was treated with 3 μM CHIR99021, 10 ng/ml human BMP4 (peprotech) and 50 nM all-trans retinoic acid (ATRA) in cSFD medium. At day 15, the lung progenitor cells were detached from 2D culture by trypsinization, mixed with matrigel, and replated to form 70 μL domes on 24-well plates. The 3D culture was maintained in cSFD containing 100 ng/ml human FGF10, 250 ng/ml human FGF2, 50 nM dexamethasone, 0.1 mM 8-bromo-cAMP (Sigma Aldrich) and 0.1 mM IBMX (Sigma Aldrich). To promote the formation of ciliated cells, the air-liquid culture system was adopted (Fulcher and Randell, 2013; Hawkins et al., 2021). D30 hPSC-AOs were switched to culturing in PneumaCult Ex Plus medium (StemCell Technologies, Inc.) supplemented with 1 μM A83-01(Tocris) and 1 μM DMH1(Tocris) for at least 10 days. At day 40, the hPSC-AOs above were dissociated into single cells and seeded into 0.4 μm transwell permeable support inserts coated with 1% Matrigel at 30,000 cells per insert to expand for another 10 days. The culture medium in both upper and lower chambers was changed to PneumaCult ALI medium (StemCell Technologies, Inc.) for inducing the formation of ciliated cells at day 49. The medium on the top was removed the next day. The cells were exposed to air for at least 14 days for further differentiation.

Colonic organoids differentiation

For definitive endoderm (DE) differentiation, cells were passaged with Accutase at 1:2-1:3 for the first day. When achieving 80%–90% confluency, hESCs were treated with 3 μM CHIR99021 and 100 ng/ml Activin A in basal medium RPMI1640 supplemented with 1X Pen-Strep for 1 day, and changed to the basal medium containing 100 ng/ml Activin A and 2% FBS for the next 48 hours. To induce CDX2+ hindgut endoderm, DE was treated with 3 μM CHIR99021 and 500 ng/ml FGF4 (Peprotech) in RPMI1640 supplemented with 1X B27 supplement (GIBCO) and 1X Pen-Strep (GIBCO) for 4 days with daily changing of fresh media. Spheroids began to bud out from the 2D culture during the hindgut differentiation process. The hindgut endoderm was then subjected to colonic lineage induction by treatment with 100 ng/ml 3 μM CHIR99021, 100 ng/ml BMP2 (Peprotech), and 100 ng/ml hEGF (Peprotech) in Advance DMEM F12 medium supplemented with 1X B27 supplement (GIBCO), 1X GlutaMax (GIBCO), 10 mM HEPES (GIBCO) and 1X Pen-Strep for 3 days with daily changing of fresh medium. After colonic fate induction, the colon progenitor spheroids were collected from the initial 2D cultures and embedded in a Matrigel dome in a 24-well plate. Differentiation to mature colonic cell types was achieved by culturing these colon progenitor spheroids in differentiation medium containing 600 nM LDN193189 (Axon), 3 μM CHIR99021 and 100 ng/ml hEGF in Advance DMEM F12 medium supplemented with 1X B27 supplement, 1X GlutaMax, 10 mM HEPES and 1X Pen-Strep. The differentiation medium was changed every 3 days for at least 40 days to achieve full colonic differentiation. The hPSC-COs were passaged and expanded every 10 – 14 days at 1:6 density. To passage hPSC-COs, the Matrigel domes were scrapped off the plate and resuspended in cold splitting media (Advance DMEM F12 medium supplemented with 1X GlutaMax, 10 mM HEPES and 1X Pen-Strep). The hPSC-COs were mechanically dislodged from the Matrigel dome and fragmented by pipetting in cold splitting media. The old Matrigel and splitting media were removed after pelleting cells and the organoids were resuspended in 100% Matrigel. 50 μL Matrigel containing fragmentized hPSC-COs were plated in one well of a pre-warmed 24-well plates.

SARS-CoV-2 infections

SARS-CoV-2 variants, including USA-WA1/2020 and B.1.351, were obtained from the World Reference Center for Emerging Viruses and Arboviruses located at the University of Texas, Medical Branch via the CDC. SARS-CoV-2 was propagated in Vero E6 cells (ATCC) in EMEM supplemented with 10% FCS, 1 mM Sodium Pyruvate and 10 mM HEPES as described previously (Liu et al., 2020).

hPSC-AOs or hPSC-COs were fragmented into small cell clusters and plated on 10% matrigel-coated plates. The infection was performed in the culture media at the indicated MOIs at 37°C. For pre-infection treatment experiments, hPSC-AOs or hPSC-COs were pretreated with DMSO (control), 10 μM GW6471, 10 μM xanthohumol, 1 μM chetomin, 3 μM ND-646 or 3 μM TOFA for 4 hours prior to infection. For post-infection treatment experiments, hPSC-AOs or hPSC-COs were treated with DMSO (control) or 10 μM GW6471 at 24 hpi. At the indicated hpi, cells were washed three times with PBS. For RNA analysis cells were lysed in TRIzol (Invitrogen). For immunofluorescence staining cells were fixed in 4% formaldehyde for 60 minutes at room temperature.

All work involving live SARS-CoV-2 was performed in the CDC/USDA-approved BSL-3 facility at the Aaron Diamond AIDS Research Center located at Columbia University.

Immunohistochemistry

For FOXJ1 and ACE2 staining, the hPSC-AOs were fixed in 4% paraformaldehyde and transferred to 30% sucrose, followed by embedding and freezing in O.C.T (Fisher Scientific, Pittsburgh, PA). The hPSC-AOs or hPSC-COs were fixed in 4% formaldehyde, followed with permeabilization in PBS containing 0.1% Triton X-100. For immunofluorescence, cells or tissue sections were immuno-stained with primary antibodies at 4°C overnight and secondary antibodies at room temperature for 1 hour. The information for primary antibodies and secondary antibodies are provided in Table S5. Nuclei were counterstained by DAPI.

qRT-PCR

Total RNA samples were prepared from organoids using TRIzol and the Direct-zol RNA Miniprep Plus kit (Zymo Research) according to the manufacturer’s instructions. To quantify viral replication, measured by the accumulation of subgenomic N transcripts, one-step quantitative real-time PCR was performed using SuperScript III Platinum SYBR Green One-Step qRT-PCR Kit (Invitrogen) with primers specific for the TRS-L and TRS-B sites for the N gene as well as ACTB as an internal reference as described previously (Yang et al., 2020). Quantitative real-time PCR reactions were performed on an Applied Biosystems QuantStudio 6 Flex Real-Time PCR Instrument. Delta-delta-cycle threshold (ΔΔCT) was determined relative to the ACTB and mock infected /treated samples. Error bars indicate the standard deviation of the mean from three biological replicates. The sequences of primers/probes are provided in Table S6.

High throughput chemical screening

hPSC-AOs were dissociated using TrypLE for 10 min in a 37°C water bath and replated into 10% Matrigel-coated 384-well plates at 10,000 cells/50 μl medium/well. After overnight plating, compounds from a library containing signaling pathway regulators and FDA-approved drugs were added at 10 μM. The detailed chemical information is listed as Table S1. DMSO treatment was used as a negative control. hPSC-AOs were then infected with SARS-CoV-2 (MOI = 0.2). After 48 hpi, hPSC-AOs were fixed for Immunofluorescence assays using the ImageXpressMICRO Automated High-Content Analysis System to determine the percentage of SARS-N+ cells (infection rate).To calculate EC50, the infection rate was normalized to the average of the infection rate under DMSO-treated conditions. The efficacy curves were calculated using Prism GrapdPad. The relative survival rate was calculated by normalizing the cell number to the average of cell numbers at control conditions.

SARS-CoV-2 pseudovirus

Recombinant Indiana VSV (rVSV) expressing SARS-CoV-2 spikes were generated as previously described (Zhao et al., 2017). HEK293T cells were grown to 80% confluency before transfection with pCMV3-SARS-CoV-2-spike (kindly provided by Dr. Peihui Wang, Shandong University, China) using FuGENE 6 (Promega). Cells were cultured overnight at 37°C with 5% CO2. The next day, medium was removed and VSV-G pseudo-typed ΔG-luciferase (G∗ΔG-luciferase, Kerafast) was used to infect the cells in DMEM at an MOI of 3 for 1 hour before washing the cells with 1 × DPBS three times. DMEM supplemented with anti-VSV-G antibody (I1, mouse hybridoma supernatant from CRL-2700; ATCC) was added to the infected cells and they were cultured overnight as described previously (Liu et al., 2020). The next day, the supernatant was harvested and clarified by centrifugation at 300 g for 10 minutes and aliquots stored at −80°C.

hPSC-AOs were seeded in 384-well plates and treated with chemicals, SARS-CoV-2- pseudovirus was added at the indicated MOI = 1 and the plate centrifuged at 1200 g for 1 hour. Then, the organoids were cultured at 37°C with 5% CO2. At 24 hpi, organoids were harvested for luciferase assays following the Luciferase Assay System protocol (E1501, Promega)

Viral particle stability assay

To determine whether chemical affects viral particle stability, a concentrated SARS-CoV-2 virus preparation was pre-incubated with chemical for 3 hours. After a 1000-fold dilution, the pre-treated SARS-CoV-2 was added to hPSC-AOs at a final MOI = 0.2 along with chemical or control. Under these conditions, we can separate the impact of chemical on viral particles rather than on host cells.

Sequencing and gene expression UMI counts matrix generation

The 10X Libraries were sequenced on the Illumina NovaSeq6000 sequencer with pair-end reads (28 bp for read 1 and 91 bp for read 2). The sequencing data were primarily analyzed by the 10X cellranger pipeline (v3.0.2) in two steps. In the first step, cellranger mkfastq demultiplexed samples and generated fastq files; and in the second step, cellranger count aligned fastq files to the 10X pre-built human reference (GRCh38 v3.0.0) and extracted gene expression UMI counts matrix.

Single cell RNA-seq data analysis

We filtered cells with less than 500 or more than 6000 genes detected, cells with less than 1000 or more than 30000 UMIs detected, as well as cells with mitochondrial gene content greater than 15%, and used the remaining 7170 cells for downstream analysis.

We normalized the gene expression UMI counts using a deconvolution strategy implemented by the R scran package (v.1.14.1). In particular, we pre-clustered cells using the quickCluster function; we computed size factor per cell within each cluster and rescaled the size factors by normalization between clusters using the computeSumFactors function; we normalized the UMI counts per cell by the size factors and took a logarithm transform using the normalize function.

We identified highly variable genes using the FindVariableFeatures function in the R Seurat package (v3.1.0) (Stuart et al., 2019), and selected the top 3000 variable genes after excluding mitochondrial genes, ribosomal genes and dissociation-related genes. The list of dissociation-related genes was originally built on mouse data (van den Brink et al., 2017); we converted them to human ortholog genes using Ensembl BioMart. We scaled the data using the ScaleData function and performed PCA on the top 3,000 variable genes using the RunPCA function. We ranked the PCs based on the percentage of variance explained by each component using the ElbowPlot function and selected the top 21 PCs for downstream visualization and clustering analysis.

We ran Uniform Manifold Approximation and Projection (UMAP) dimensional reduction using the RunUMAP function in the R Seurat package with the number of neighboring points setting to 35 and training epochs setting to 1000. We clustered cells into nine clusters by constructing a shared nearest neighbor graph and then grouping cells of similar transcriptome profiles using the FindNeighbors function and FindClusters function (resolution set to 0.5) in the R Seurat package. After reviewing the clusters, we merged them into five clusters representing ciliated-like cells, goblet-like cells, basal cells and two proliferating cell clusters, for further analysis. We identified marker genes for the merged five clusters using the FindMarkers function in the R Seurat package, and selected top 10 positive marker genes per cluster for heatmap plot using the R pheatmap package. We generated UMAP and violin plots highlighting expressions of selected genes using the R ggplot2 package.

To evaluate cell fates of hPSC-derived cells, we collected unique marker genes for fifteen reference cell types in human adult lung based on a recent study (Travaglini et al., 2020). We compared the marker genes identified in our five clusters and those in the fifteen reference cell types by calculating the fraction of overlapping marker genes. We performed z-score transformation on raw fraction values across reference cell types and generated the cell-type similarity heatmap plot using the R pheatmap package. We examined the expression of adult human ciliated and proximal ciliated cell marker genes in the hPSC-derived cells. In particular, we assigned each cell a score by comparing the average expression of these marker genes to that of a random set of background genes using the AddModuleScore function in the R Seurat package.

RNA-Seq on hPSC-AOs

Total RNA was extracted in TRIzol (Invitrogen) and DNase I treated using the Directzol RNA Miniprep kit (Zymo Research) according to the manufacturer’s instructions. RNaseq libraries of polyadenylated RNA were prepared using the TruSeq RNA Library Prep Kit v2 (Illumina) or TruSeq Stranded mRNA Library Prep Kit (Illumina) according to the manufacturer’s instructions. cDNA libraries were sequenced with pair-end 51 bps using an Illumina NovaSeq6000 platform. The resulting reads were checked for quality using FastQC v0.10.1(https://www.bioinformatics.babraham.ac.uk/projects/fastqc) and were trimmed for adaptor sequences and low quality bases using cutadapt v1.18(Kechin et al., 2017). In order to measure viral gene expression, the trimmed reads were aligned to the human reference genome (GRCh37) combined with SARS-CoV-2 genome (NC_045512.2) using STAR aligner v.2.5.2b (Dobin et al., 2013). Raw gene counts were quantified using HTSeq-count v0.11.2(Anders et al., 2015).

We observed large fractions of viral reads in several SARS-CoV-2 infected hPSC-AOs, and as a result, the human gene counts are significantly lower than those in the mock and treated samples. To make a fair comparison, we down-sampled the raw gene counts in the mock and treated samples, including Mock hPSC-AOs, GW6471 treated hPSC-AOs at 48 hpi (MOI = 0.2) and 1 μM chetomin treated hPSC-AOs at 48 hpi (MOI = 0.2), such that the human gene counts in these samples are comparable to the median human gene counts in the SARS-CoV-2 infected hPSC-AOs.

We performed differential expression analysis on gene counts using R DESeq2 package v1.26.0 (Love et al., 2014), and selected differentially expressed genes with adjusted p value less than 0.05 and absolute log2 fold change greater than 1 for pathway analysis using QIAGEN Ingenuity Pathways Analysis (IPA). We applied regularized logarithm transformation to the counts data and performed principle component analysis (PCA) using the plotPCA function in the DESeq2 package. We measured the sample-to-sample distances by applying the R dist function to the transpose of the transformed counts data, and performed a hierarchical clustering based on the distance using the R hc function. We generated heatmap plots based on the transformed counts data using Heatmapper. We generated volcano plots using the R ggplot2 package.

We calculated read coverage on the viral transcriptome by extracting the number of reads at each genomic position using samtools (Li et al., 2009) coverage. We normalized the coverage by the total number of aligned reads in each sample and plotted the coverage per million reads using R ggplot package.

RNA-Seq on hPSC-COs

cDNA libraries were sequenced with pair-end 51 bps using an Illumina NovaSeq6000 platform. The resulting reads were checked for quality using FastQC v0.10.1 and were trimmed for adaptor sequences and low quality bases using cutadapt v1.18 (Dobin et al., 2013). In order to measure viral gene expression, the trimmed reads were aligned to the human reference genome (GRCh37) combined with SARS-CoV-2 genome (NC_045512.2) using STAR aligner v.2.5.2b. Raw gene counts were quantified using HTSeq-count v0.11.2(Anders et al., 2015).

We performed differential expression analysis on the gene counts using R DESeq2 (Love et al., 2014) package v1.26.0, and selected differentially expressed genes with adjusted p value less than 0.05 and absolute log2 fold change greater than 1 for pathway analysis using QIAGEN Ingenuity Pathways Analysis (IPA). We applied regularized logarithm transformation to the counts data and performed principle component analysis (PCA) using the plotPCA function in the DESeq2 package. We measured the sample-to-sample distances by applying the R dist function to the transpose of the transformed counts data, and performed a hierarchical clustering based on the distance using the R hc function. We generated volcano plots using the R ggplot2 package.

Metabolic profiling

hPSC-AOs were collected into 150 mM ammonium acetate. After centrifugation, cell pellets were collected for metabolite profiling. Metabolite extraction was performed with 200 μL of 80% methanol (20% water) with internal standards (U13C succinate, and U13C citrate, and heptadecanoate). The extracts were dried under gentle nitrogen flow. The dried samples were subjected to a two-step derivatization procedure (methoximation and silylation), and analyzed with GC-TOF/MS (GCT Premier, Waters Corporation) as described in our previous publication (Qiu et al., 2016).

In brief, the samples were methoximized with 50 μl of methoxyamine hydrochloride (MOA, 15 mg/mL in pyridine) at 30°C for 90 minutes. The silylation step was done with 50 μl of N,O-Bis(trimethylsilyl)trifluoroacetamide (BSTFA, containing 1% TMCS, Supelco, PA, USA) at 70°C for 60 minutes. After the samples cooled to room temperature, they were analyzed by GC-TOF/MS. Separation was performed with a 30-m DB-5MS column coupled with 10-m guard column. Helium was used as carrier gas at a flow of 1 ml/min. The oven program is as the following: started at 60°C for 1 minute, 10°C /min to 320 minutes and kept for 3 minutes.

Raw data files (.raw) generated from GC-TOF/MS were analyzed in Refiner MS (Genedata Expressionist, Basel, Switzerland) software. The output from Refiner MS contains retention time, quantification mass, and peak volume for each metabolite. Metabolite annotation was performed by searching against the mass spectrometry using our in-house library, Fiehn library and NIST 14 libraries in Refiner MS. To minimize the variations from cell number difference among different samples, the data was normalized to the sum of all the annotated metabolites. The normalized data was used for the statistical analysis.

HIF1a knockdown in hPSC-AOs

Two shRNAs against HIF1α and one scrambled negative control non-effective shRNA in lentiviral GFP vectors were purchased from OriGene company. The shRNA sequences are shown in Table S3. To generate lentivirus particles, shRNA vectors and lentivirus packaging vectors were co-transfected into 293T cells. The day 15 lung progenitor cells were infected with the collected lentivirus (MOI = 0.5) with 8 μg/ml polybrene. To increase the infection efficiency, the cells were centrifuged at 2300 rpm for 1 hour at 30°C. At 24 hpi, the cells were selected with 1 μg/ml puromycin for an additional 72 hours. After selection, the cells were dissociated into single cells and seeded into 24 well plates at 300-400 cells/μl for 3D organoid differentiation.

Quantification and statistical analysis

N = 3 independent biological replicates were used for all experiments unless otherwise indicated. n.s. indicates a non-significant difference. P-values were calculated by unpaired two-tailed Student’s t test unless otherwise indicated. ∗p < 0.05, ∗∗p < 0.01, ∗∗∗p < 0.001 and ∗∗∗∗p < 0.001.

REAGENT or RESOURCE	SOURCE	IDENTIFIER
Antibodies
Human ACE-2 Antibody	R & D Systems	#AF933 RRID:AB_355722
CDX2	Biogenex	#MU392A-UC RRID:AB_2650531
Cytokeratin-20	Santa Cruz	#sc-56522 RRID:AB_831485
FOXJ1	Thermo Fisher Scientific	#14-9965-82 RRID:AB_1548835
Acetyl-α-Tubulin	Cell Signaling	#5335 RRID:AB_2798556
acetyl-alpha tubulin	Sigma-Aldrich	#MABT868 RRID:AB_2819178
MUC5AC	Thermo Fisher Scientific	#MA5-12178 RRID:AB_10978001
P63	Biocare	#CM 163 A RRID:AB_10582730
SARS-CoV/SARS-CoV-2 Nucleocapsid Antibody	Sino Biological	#40143-R001 RRID:AB_2827974
Ki67 Antibody (SP6)	THERMO FISHER	#RM-9106-S1 RRID:AB_2314700
Donkey anti-Mouse IgG (H+L) Highly Cross-Adsorbed Secondary Antibody, Alexa Fluor 488	Thermo Fisher Scientific	#A-21202 RRID:AB_141607
Donkey anti-Rabbit IgG (H+L) Secondary Antibody, Alexa Fluor 594 conjugate	Thermo Fisher Scientific	#A-21207 RRID:AB_141637
Donkey anti-Rabbit IgG (H+L) Secondary Antibody, Alexa Fluor 647 conjugate	Thermo Fisher Scientific	#A-31573 RRID:AB_2536183
Donkey anti-Goat IgG (H+L) Secondary Antibody, Alexa Fluor 647	Thermo Fisher Scientific	#A-31571 RRID:AB_162542
Donkey anti-Goat IgG (H+L) Cross-Adsorbed Secondary Antibody, Alexa Fluor 647	Thermo Fisher Scientific	#A-21447 RRID:AB_141844
DAPI	Santa Cruz	#sc-3598
Chemicals, peptides, and recombinant proteins
Y-27632	MedchemExpress	#HY-10583
CHIR99021	Cayman Chemical	#13122
Retinoic acid	Sigma Aldrich	#R2625-500MG
LDN 193189 hydrochloride - DM 3189 hydrochloride	Axon Medchem	#Axon 1509
L-Ascorbic acid	Sigma Aldrich	#A4544-100G
SB431542	R&D Systems	#1614/50
Dorsomorphin dihydrochloride	R&D Systems	#3093/50
IWP2	R&D Systems	#3533/50
Dexamethasone	Sigma-Aldrich	#D4902
8-Bromo-cAMP	Sigma-Aldrich	#B5386
IBMX	Sigma-Aldrich	#I5879
A83-01	Tocris Bioscience	#2939
DMH1	Tocris Bioscience	#4126
GW6471	Selleck	#S2798
Xanthohumol	MedchemExpress	#HY-N1067
Chetomin	Sigma-Aldrich	#C9623
ND-646	MedchemExpress	#HY-101842
TOFA	MedchemExpress	#HY-101068
Activin A	R&D Systems	#338-AC-500/CF
Recombinant Human FGF-10 Protein	Peprotech	#100-26-500UG
Recombinant Human bFGF Protein	Peprotech	#100-18B-500UG
Recombinant Human BMP-4 Protein	R & D Systems	#314-BP
Recombinant Human EGF	Peprotech	#AF-100-15-500UG
Recombinant Human BMP-2	Peprotech	#AF-120-02
Recombinant Human FGF-4	Peprotech	#AF-100-31
Deposited data
scRNA-seq	This paper	GEO: GSE160231
RNA-seq	This paper	GEO: GSE160231
Metabolism Profiling	This paper	NMDR: ST001921
Experimental models: Cell lines
hESC line H1	WiCell	0043
RUES2	WiCell	0013
293T	ATCC	#CRL-11268
Vero E6	ATCC	#CRL-1586
Software and algorithms
Code for data analysis	This paper	zenodo: https://zenodo.org/record/5533014
Cell Ranger	10X Genomics	https://support.10xgenomics.com/single-cell-gene-expression/software/overview/welcome
Scran	Lun et al., 2016	https://bioconductor.org/packages/release/bioc/html/scran.html
Rstudio	Rstudio	https://www.rstudio.com/
Seurat R package v3.1.0	Butler et al., 2018	https://satijalab.org/seurat/
DAVID6.8	LHRI	https://david.ncifcrf.gov/home.jsp
Adobe illustrator CC2017	Adobe	https://www.adobe.com/products/illustrator.html
Graphpad Prism 6	Graphpad software	https://www.graphpad.com
ToppCell Atlas	Toppgene	https://toppgene.cchmc.org/
Other
StemFlex	GIBCO Thermo Fisher	#A3349401
mTeSR1 Complete Kit	Stem Cell Technologies	#85850
F12	GIBCO Thermo Fisher	#31765035
Penicillin-Streptomycin (5,000 U/mL)	GIBCO Thermo Fisher	#15070063
IMDM	GIBCO Thermo Fisher	#21056023
GlutaMAX Supplement	Thermo Fisher Scientific	#35050079
Accutase	StemCell Technologies, Inc.	# 07920
ReleSR	StemCell Technologies, Inc.	# 05872
B27	Thermo Fisher Scientific	# A3582801
Matrigel	Corning	#354234
N2 supplement-B	Stem Cell Technologies	#7156
DMEM/F12	Thermo Fisher Scientific	# 10565018
Advanced DMEM/F-12	GIBCO Thermo Fisher	# 12-634-028
RPMI 1640	GIBCO Thermo Fisher	#10-040-CV
Monothioglycerol	Sigma Aldrich	#M6145
PneumaCult-Ex Plus Medium	StemCell Technologies, Inc.	#05040
PneumaCult-ALI Medium	StemCell Technologies, Inc.	#05001