Recapitulation of SARS-CoV-2 infection and cholangiocyte damage with human liver ductal organoids

The emerging pandemic of coronavirus, SARS-CoV-2 (previously named 2019-nCoV), has posed significant threats to global public health (Wu et al., 2020). The dominant symptoms of coronavirus disease 2019 (COVID-19) are fever and cough (Chen et al., 2020; Huang et al., 2020). However, a proportion of patients showed multi-organ damage and dysfunction (Chen et al., 2020; Huang et al., 2020; Zhu et al., 2020). Of note, liver damage is emerging as a co-existed symptom reported in patients with COVID-19. A recent epidemiologic study in Shanghai (China) reported that 75 out of 148 (50.7%) COVID-19 patients had liver function abnormality, indicated by key liver function parameters above the normal range, including alanine aminotransferase (ALT), aspartate aminotransferase (AST), alkaline phosphatase (ALP) or total bilirubin (TBIL) (Fan et al., 2020). A nationwide clinical study collecting 1,099 COVID-19 patients revealed that around 20% of patients had elevated ALT and AST and around 10% of patients had elevated TBIL. Especially, the percentage of patients with liver damage is much higher among severe COVID-19 patients than patients with mild symptoms (Huang et al., 2020). Although clinical correlation has been implicated, it remains unclear whether the liver damage is caused by direct virus infection in the liver or by systematic dysfunction such as cytokine storm. Viruses bind to host receptors to initiate infection. Recent studies have demonstrated that SARS-CoV-2 uses the SARS-CoV receptor angiotensin I converting enzyme 2 (ACE2) for host cell entering and transmembrane serine protease 2 (TMPRSS2) for viral spike (S) protein priming (Kuhn et al., 2004; Hoffmann et al., 2020; Wan et al., 2020; Zhu et al., 2020). It has been shown that ACE2 expression is widely distributed across human tissues, including lung, liver, kidney and multiple digestive tract organs (Qi et al., 2020; Zhang et al., 2020; Zhao et al., 2020). Significant enrichment of ACE2+ population in cholangiocytes compared to hepatocytes in the human healthy liver was reported recently (Chai et al., 2020), implying that SARS-CoV-2 might directly target ACE2+ cholangiocytes in patients. However, whether the virus indeed infects human cholangiocytes thus causes local damage has not been addressed yet. At present, due to the lack of suitable research models, the mode of virus transmission and tissue tropism is not well established yet. Studies on mechanisms of SARS-CoV-2 pathogenesis mainly depend on bioinformatics analysis, clinical characteristics and rare autopsy reports (Xu et al., 2020). Here we report the use of human organoids as a tool to investigate the SARS-CoV-2 infection and virus-induced tissue damage ex vivo at the cellular and molecular levels. In a three-dimensional (3D) culture system with defined culture medium, liver bile duct-derived progenitor cells embedded in Matrigel can self-assemble into long-term expandable 3D structures that termed “liver ductal organoids”, which retain their tissue-of-origin commitment and genetic stability during self-renewing (Huch et al., 2015). To establish the SARS-CoV-2 infection model with human liver ductal organoids, we first determined whether the long-term organoid culture could preserve the cholangiocytes expressing ACE2 and TMPRSS2 ex vivo. We processed single-cell RNA sequencing (scRNA-seq) to interrogate the transcriptome of cholangiocytes in human liver ductal organoids. A total number of 7,978 cells were analyzed and cell populations were visualized by t-distributed stochastic neighbor embedding (t-SNE), partitioning the cells into 7 clusters (Fig. 1A). The common cholangiocyte markers epithelial cell adhesion molecule (EPCAM) and keratin 19 (KRT19) were uniformly highly expressed in all the 7 clusters, indicating the heterogeneity of cholangiocytes in these organoids was relatively low (Fig. 1B). Notably, we identified the SARS-CoV-2 receptor gene ACE2 expressed sparsely among cluster 0–5 in unbiased preferences and was detectable in 2.92% cells (233 out of 7,978) (Fig. 1C and 1D). Anti-ACE2 immunostaining further verified the presence of ACE2+ cholangiocytes in human liver ductal organoids (Fig. 1E). Besides, TMPRSS2 expressed uniformly across all the clusters and accounted for 51.45% of the total cells (4,105 out of 7,978), it is worth mentioning that 68.24% of the ACE2+ cells were co-expressing TMPRSS2 (159 out of 233) (Fig. 1C and 1D), making this cell population potentially highly vulnerable to SARS-CoV-2 infection. Interestingly, we

collected by pipetting and washed once with PBS, then repeated the centrifugation and removed supernatant. The organoids were embedded with Matrigel, followed by seeding on a 24-well plate. After polymerization, culture medium was added.
Quantitative RT-PCR. Total RNA was isolated from organoids by RNeasy Mini kit (QIAGEN,74106) and reverse-transcribed into cDNA with M-MLV Reverse Transcriptase (Invitrogen, 28025013) in biosafety level 2 facility strictly following the regulations. Quantitative real-time PCR was performed on CFX384 Touch System (Bio Rad). Primers used were listed in Supplementary Table 2. The SARS-CoV-2 primer and probe sets were   provided by Integrated DNA Technologies (IDT, 10006606).
Single-cell RNA seq and data analysis. Single-cell RNA sequencing was performed using the 10x Genomics Chromium System. Human liver ductal organoids were derived from a patient who underwent resection, cultured for 3 passages as described above. Mouse primary liver ductal organoids were cultured from biliary ducts isolated from an 8-weekold C57BL/6 mouse. Briefly, organoids were dissociated with 1× TrypLE Select Enzyme (Gibco, 12563011) to obtain single cell suspension. A total of around 8,000 cells per sample were captured on a 10×Chromium device and library preparation was carried out using Single Cell 3' Reagent Kits v2 according to the manufacturer's instructions (10× Genomics). Libraries were sequenced on an Illumina NovaSeq 6000 platform.
Cell Ranger (version 3.1) with default parameters was used to process sequencing data to generate feature-barcode matrices. The human dataset was analyzed using the standard workflow on the Seurat R Package (version 3.1.3). For the feature-barcode matrix of 8,094 cells from the human dataset, we removed cells with less than 500 genes and more than 6,000 genes as well as cells with high fraction of mitochondrial UMIs (> 20%). 7,978 high quality cells and 17,447 expressed genes were remained for downstream analysis. The cell populations were clustered using the 'FindClusters' function and After quality control, clean reads were aligned to human reference genome (GRCh38) using HISAT2 (version 2.1.0). The alignments were then passed to StringTie (version 1.3.5) to assemble and quantify the transcripts in each sample. Differentially expressed genes (DEGs) was identified by the R package edgeR (version 3.28.1). Genes were defined as DEGs if it possesses the following characteristics: 1) gene expression (FPKM) > 1 in any sample, 2) absolute log2 (fold change) > 2 and 3) p-value < 0.01. Visualization and hierarchical clustering of log2-transformed FPKM was generated by pheatmap (version 1.0.12). GO analysis was performed using metascape (http://metascape.org). Gene set enrichment analysis was performed with GSEA v3.0 software (available from the Broad Institute).

Statistical analysis.
We employed Student's t-test or ANOVA test to analyze the parametric experimental results. Significant differences were noted with asterisks.