key: cord-0909973-y0nmlkdp authors: Han, Lei; Wei, Xiaoyu; Liu, Chuanyu; Volpe, Giacomo; Wang, Zhifeng; Pan, Taotao; Yuan, Yue; Lei, Ying; Lai, Yiwei; Ward, Carl; Yu, Yeya; Wang, Mingyue; Shi, Quan; Wu, Tao; Wu, Liang; Liu, Ya; Wang, Chunqing; Zhang, Yuanhang; Sun, Haixi; Yu, Hao; Zhuang, Zhenkun; Tang, Tingting; Huang, Yunting; Lu, Haorong; Xu, Liqin; Xu, Jiangshan; Cheng, Mengnan; Liu, Yang; Wong, Chi Wai; Tan, Tao; Ji, Weizhi; Maxwell, Patrick H.; Yang, Huanming; Wang, Jian; Zhu, Shida; Liu, Shiping; Xu, Xun; Hou, Yong; Esteban, Miguel A.; Liu, Longqi title: Single-cell atlas of a non-human primate reveals new pathogenic mechanisms of COVID-19 date: 2020-04-23 journal: bioRxiv DOI: 10.1101/2020.04.10.022103 sha: e1e6416bb022439e5f15b3c20efc2ca444886247 doc_id: 909973 cord_uid: y0nmlkdp Stopping COVID-19 is a priority worldwide. Understanding which cell types are targeted by SARS-CoV-2 virus, whether interspecies differences exist, and how variations in cell state influence viral entry is fundamental for accelerating therapeutic and preventative approaches. In this endeavor, we profiled the transcriptome of nine tissues from a Macaca fascicularis monkey at single-cell resolution. The distribution of SARS-CoV-2 facilitators, ACE2 and TMRPSS2, in different cell subtypes showed substantial heterogeneity across lung, kidney, and liver. Through co-expression analysis, we identified immunomodulatory proteins such as IDO2 and ANPEP as potential SARS-CoV-2 targets responsible for immune cell exhaustion. Furthermore, single-cell chromatin accessibility analysis of the kidney unveiled a plausible link between IL6-mediated innate immune responses aiming to protect tissue and enhanced ACE2 expression that could promote viral entry. Our work constitutes a unique resource for understanding the physiology and pathophysiology of two phylogenetically close species, which might guide in the development of therapeutic approaches in humans. Bullet points We generated a single-cell transcriptome atlas of 9 monkey tissues to study COVID-19. ACE2+TMPRSS2+ epithelial cells of lung, kidney and liver are targets for SARS-CoV-2. ACE2 correlation analysis shows IDO2 and ANPEP as potential therapeutic opportunities. We unveil a link between IL6, STAT transcription factors and boosted SARS-CoV-2 entry. 7 cluster 13, pancreatic acinar cells in cluster 26 and parotid acinar cells in cluster 30 (Fig. 1c, 161 Extended Data Fig. 1c) . We next performed clustering and differential gene expression 162 analysis to dissect the cellular composition of each individual organ. These analyses confirmed 163 the typical patterns of cell heterogeneity for all the organs/tissues. When examining the lung 164 tissue, we defined 10 major clusters with specific molecular markers, including ciliated cells, 165 macrophages, cycling macrophages, smooth muscle cells, fibroblasts, pericytes, pulmonary 166 alveolar (pneumocytes) type 1 and type 2, endothelial and club cells (Extended Data Fig. 2a) . 167 The kidney consisted of 11 clusters, those being podocytes, thick ascending limb cells, Fig. 2h) . Finally, our data 186 showed the largest population of the pancreas to be acinar cells, while smaller clusters were 187 comprised of stromal, ductal, and islet cells (alpha and beta), together with a population that 188 could not be assigned to any known cell type (Extended Data Fig. 2i) . 189 In conclusion, we have successfully profiled the transcriptome of nine organs at single-190 cell resolution in monkey, which could assist in the study of 192 8 ACE2 and TMPRSS2 single-cell expression landscape in a non-human primate . 193 Recent studies have reported that, similarly to SARS-CoV, the capacity of SARS-CoV-2 virus to 194 infect host cells relies on viral spike (S) protein binding to ACE2 entry receptor 9, 10 , which is 195 involved in the control of blood pressure through the renin-angiotensin system 19 . This 196 phenomenon is primed by the multifunction serine protease TMPRSS2 20 . Accordingly, double 197 positive (ACE2 + /TMPRSS2 + ) cells have higher risk of infection by SARS-CoV-2. Although 198 immunohistological studies have demonstrated localization of these two proteins in the 199 respiratory tract 21 , it is unclear which cell subtypes express these genes and how homogenous 200 the expression among a specific cell subtype is. Also, comprehensive information about other 201 cell types and organs that express these two proteins and could be targeted by the virus in 202 human or monkey is lacking. 203 We inspected our data to see how widespread and homogenous ACE2 expression was 204 in the monkey tissues. As expected, ACE2 was detected in several lung clusters, mainly ciliated 205 cells, club cells and pulmonary alveolar type 2 cells (Fig. 2a, 2d, 3a upper panel) , whereas in 206 the kidney, ACE2 was primarily present in proximal tubule cells (Fig. 2a, 2d, 3b upper panel) . 207 The latter is consistent with reports describing that a significant number of COVID-19 patients 208 display altered kidney function 15, 22 . Interestingly, ACE2 expression was heterogenous among 209 these cell subtypes in both lung and kidney. In the liver, ACE2 was mostly expressed in 210 cholangiocytes, with a smaller degree of expression also found in hepatocytes (Fig. 2a, 2d, 3c 211 upper panel). Notably, the closely related SARS-CoV caused liver injury due to hepatitis in 212 some patients 23 , suggesting that the liver may also be a direct target for SARS-CoV-2. A small 213 proportion of ACE2 + was also observed in pancreatic islet cells (Fig. 2a, 2d, Extended Data Fig. 214 3a). In contrast, little or no expression was observed in thyroid, neocortex, parotid and PBMC 215 ( Fig. 2a, 2d, Extended Data Fig. 3a) . Negligible ACE2 expression in the neocortex suggests that 216 other tissues may be affected by SARS-CoV-2 that cause loss of taste and smell, regarding the 217 latter in particular the olfactory epithelium. 218 TMPRSS2 displayed more broadly expressed across cell types in multiple tissues, 219 although it was highest in kidney cells. However, in contrast to ACE2, its expression was 220 highest in the distal convoluted tubule, thin limb, intercalated and principal cell 1 and 2 kidney 221 clusters (Fig. 2b, 2d, 3b lower panel, Extended Data Fig. 3b) . Additionally, significant 222 TMPRSS2 was observed in both parotid and pancreatic acinar cells, thyroid follicular cells, 223 cholangiocytes and in several lung clusters (Fig. 2b, 2d, Extended Data Fig. 3) . We then 9 determined which cells co-expressed both genes (ACE2 + /TMPRSS2 + ). Notably, the largest 225 overlap between ACE2 and TMPRSS2 was observed in the ciliated and club cell clusters of the 226 lung and to a lesser extent the proximal tubule cells of the kidney (Fig. 2c, 2e) . A smaller 227 overlap was also observed in cholangiocytes and in pancreatic islet cells (Fig. 2c, 2e) . 228 Therefore, our data show that ACE2 and TMPRSS2 are expressed in a variety of cell 229 types, mainly epithelial cells, within the nine monkey organs/tissues (Supplementary Given the heterogeneous nature of ACE2 and TMPRSS2 expression within monkey tissues, we 237 decided to investigate similarities and differences between human and monkey. For this 238 purpose, we retrieved publicly available data from single-cell studies in human (see methods). 239 TMPRSS2 distribution was similar in cell subtypes of lung, kidney and liver between human 240 and monkey (Fig. 3d-3f) . However, strikingly, ACE2 showed distinct patterns among cell 241 subtypes in all three organs between the two species ( Fig. 3d-3f) . The biggest differences 242 were observed in ciliated cells of the lung, which had the highest expression of ACE2 in 243 monkey, and pulmonary alveolar type 2 cells, which had the highest expression of ACE2 in 244 human. The function of ciliated cells is to move substances (e.g. cell debris and toxic material) 245 across the surface of the respiratory tract and are commonly targeted by respiratory viruses, 246 whereas pulmonary alveolar type 2 cells have regenerative properties, are crucial for alveolar 247 homeostasis and produce surfactant 24, 25 . In the kidney of both monkey and human, the 248 highest ACE2 expression was in proximal tubule cells (Fig. 3e) , which are responsible for 249 electrolyte and nutrient reabsorption. However, renal endothelial cells had higher expression 250 in monkey compared to human. In liver, cholangiocytes had similarly high ACE2 expression in 251 monkey and human, but hepatocytes showed higher expression and more positive cells in the 252 human (Fig. 3f) . Considering the heterogenous expression of ACE2 within the proximal tubule 253 cells in monkey, we revisited the previously analyzed data and were able to sub cluster this 254 population of cells into two (S1 and S3) based on the expression of SLC5A2 and SLC7A13 26 255 (Extended Data Fig. 4, Supplementary Table 2b ). These two genes are sodium and glucose 256 cotransporters involved in glucose reabsorption in the kidney 27, 28 . We did not include thyroid, 257 pancreas or aorta in these analyses because of lack of high-quality available human single-cell 258 datasets. As for the neocortex and PBMC, they have little to no expression of ACE2 in human 259 (data not shown). 260 These differences in ACE2 expression across cell subtypes in the lung, kidney and liver 261 in monkey and human raise the possibility that infection with SARS-CoV-2 in the two species 262 will have different effects. 263 264 ACE2 correlation analysis across cell types reveals potential therapeutic targets. 265 To shed light on potential mechanisms that could facilitate ACE2-mediated SARS-CoV-2 266 infection, we performed an analysis of the Pearson's correlation coefficient, based on gene 267 expression in the 44 cell subtypes, to determine what genes are co-regulated with ACE2 in 268 monkey tissues. Correlated genes were considered as those displaying a coefficient higher 269 than 0.6 with an adjusted P value < 0.001. Using these criteria, we observed several genes 270 with marked correlation, including genes that belong to metabolic and developmental 271 pathways and genes involved in the cellular response to xenobiotic stimuli (Fig. 4a, b) . The 272 highest correlation was observed for transmembrane protein 27 (TMEM27, cor = 0.84), a 273 protein involved in trafficking amino acid transporters to the apical brush border of kidney 274 proximal tubules 29 . This is unsurprising considering that TMEM27 is an important paralog of 275 ACE2, and high expression was restricted to kidney cells. DnaJ heat shock protein family 276 (Hsp40) member C12 (DNAJC12, cor = 0.78), a gene with a role in immune response 277 processes 30 , had a distribution like TMEM27. Importantly, we also observed high correlation 278 with Indoleamine 2,3-dioxygenase 2 (IDO2, cor = 0.77), a gene with abundant expression in 279 kidney and liver cells that was also expressed in the lung and other organs. IDO2 functions 280 during the early phases of immune responses and promotes inflammatory autoimmunity 31, 281 32 . ANPEP, which encodes for alanyl aminopeptidase N, was also co-expressed with ACE2 in 282 kidney, liver and to a lesser extent in lung too (cor = 0.64), like IDO2 (Fig. 4c, d) . To understand whether epigenetic mechanisms underlie the heterogeneity of ACE2 294 expression in the kidney, as representative for other organs, we employed DNBelab C4 295 technology to perform high-throughput single-cell ATAC-seq (Fig. 5a) . After filtering, 6,353 296 nuclei were used for downstream analysis (Extended Data Fig. 5a , b, Supplementary Table 297 4). We integrated these data with the kidney transcriptomic data described in Fig. 1 and 298 proceeded to perform Louvain clustering to map all the different cell types within the dataset 299 ( Fig. 5b) . Consistent with the transcriptomic data, our epigenomic mapping identified thick 300 ascending limb cells and proximal tubule cells as the largest kidney clusters (Extended Data 301 ACE2 locus, with the highest signal detected in proximal tubule cells S1 and S3, which are also 305 the highest ACE2-expressing cells (Fig. 5d) . Our approach failed to detect significant signal 306 enrichment in the ACE2 locus in endothelial cells, possibly related to the low level of 307 expression (Fig. 5d) . Within the cells of the kidney we observed the highest percentage of 308 ACE2 + cells in the proximal tubule S3, with a lower percentage in the proximal tubule S1 and 309 endothelial cells (Fig. 5e) . Motif analysis within the open chromatin regions in ACE2 + cells 310 demonstrated that these regions were preferentially enriched in STAT1 and 3 and IRF1 311 binding sites (Fig. 5f) . These findings suggested that tissue protective cytokines including IL5, 312 IL6, EGF and interferons are acting on these proximal tubule cells S3 to induce ACE2. We 313 focused on IL6 because a recent clinical trial has been started that uses anti-IL6 receptor (IL6R) 314 antibodies in the treatment of COVID-19 315 (http://www.chictr.org.cn/showprojen.aspx?proj=49409). IL6 is a potent regulator of 316 immune responses and can be produced by a variety of interstitial cells including fibroblasts, 317 endothelial cells and more importantly tissue macrophages 36 . Interestingly, we also noticed 318 that distribution of IL6R correlated well with ACE2 in proximal tubule cells (Fig. 5g , Extended 319 Data Fig. 5d) . In human kidney a similar co-expression pattern was detected (Fig. 5h) . Our observations suggest a potential positive feedback loop between IL6 and ACE2 321 expression that can exacerbate COVID-19 disease progression due to increased viral entry and 322 dissemination. In this study, we have generated a single-cell transcriptomic atlas of nine organs (liver, 337 kidney, lung, pancreas, neocortex, aorta, parotic gland, thyroid and peripheral blood) from 338 cynomolgus monkey. We used this dataset not only to provide fundamental information 339 about the cellular composition of the different tissues tested but also as a platform to dissect 340 the overall expression distribution of the SARS-CoV-2 entry receptor, ACE2, and its serine 341 protease coactivator TMPRSS2 9, 10 . Interestingly, ACE2 was expressed in multiple epithelial 342 tissues besides the lung, especially the kidney and liver. Other organs of epithelial origin such 343 as the gut have also been implicated in the pathogenesis of the disease 38 . A consequence of 344 this is that the SARS2-CoV-2 virus could infect these organs too, which would explain some of 345 the reported clinical manifestations of COVID-19 2 . By comparing our dataset with publicly 346 available human single-cell RNA-seq data, we have also uncovered significant differences in 347 cell subtypes expressing ACE2 between human and monkey. We showed different expression 348 patterns for ACE2 in the lung, where the highest levels were detected in ciliated cells in 349 monkey and pulmonary alveolar type 2 cells in human. Similarly, we observed marked 350 differences in liver, in which monkey hepatocytes displayed higher ACE2 and a larger number 351 of positive cells compared to the human. We do not know whether these differences will 352 13 affect the pathogenesis of COVID-19 between these two species. Nevertheless, this is a 353 relevant finding considering that monkeys are a preferred model for studying the 354 effectiveness of drug treatments and of vaccines against the impending COVID-19 pandemic. (http://www.chictr.org.cn/showprojen.aspx?proj=49409). IL6 has been related to aging and 384 14 tissue damage 42 , and this may explain why elderly individuals and those with underlying 385 inflammatory conditions have more severe reactions to SARS-CoV-2 infection (Fig. 6) . 386 Importantly, high IL6 levels have been detected in plasma from COVID-19 patients 43 c 01 02 03 04 05 06 07 08 09 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 UCP1 APOL6 SLC1A3 SLC1A2 MFGE8 MS4A1 CD74 LEF1 CD4 GZMB GZMK NKG7 KRT19 SCGB3A1 PIGR FOXJ1 TUBB4B SCGB1A1 SLC8A1 RHCG SLC12A3 WNK1 TMEM52B FLT1 VWF CLDN5 EMCN SLC17A7 SLC17A6 TG TPO IYD SEC16B HAL APOA1 AMBP GAD1 SLC6A1 SLC32A1 GAD2 FOXI1 SLC4A1 DMRT2 SLC4A9 AIF1 CD163 C1QB PTPRC ADGRE1 DCN LAMA2 PKHD1L1 MTUS2 HES2 MYLK MOBP MOG MAG MBP OPALIN PDGFRA CPA2 PLA2G1B CLPS CPA1 GCG INS PRR27 LPO PLA2R1 PODXL PTPRO NPHS2 FXYD4 HSD11B2 NR4A2 LRP2 SLC16A9 SLC34A1 AGER CLDN18 VEGFA SFTPB SFTPC ATP6V1B1 EGF SLC4A7 SLC12A1 COCH ACTA1 MYH10 MYH11 Searching for the source of 649 Ebola: the elusive factors driving its spillover into humans during the West African 650 outbreak of 2013-2016 A Novel Coronavirus from Patients with Pneumonia in China Origin and evolution of pathogenic coronaviruses Genetic Recombination, and Pathogenesis of Coronaviruses Severe acute respiratory syndrome coronavirus-like virus in Chinese 658 horseshoe bats Middle East respiratory syndrome coronavirus neutralising serum 660 antibodies in dromedary camels: a comparative serological study A familial cluster of pneumonia associated with the 2019 novel 663 coronavirus indicating person-to-person transmission: a study of a family cluster Transmission of 2019-nCoV Infection from an Asymptomatic Contact 666 in Germany SARS-CoV-2 Cell Entry Depends on ACE2 and TMPRSS2 and Is 668 Blocked by a Clinically Proven Protease Inhibitor Structure, Function, and Antigenicity of the SARS-CoV-2 Spike 670 Glycoprotein. Cell, Epub ahead of print Comparative pathogenesis of COVID-19, MERS, and SARS in a 672 nonhuman primate model. Science,Epub ahead of print Reinfection could not occur in SARS-CoV-2 infected rhesus macaques. 674 bioRxiv Construction of a human cell landscape at single-cell level Tissue distribution of ACE2 protein, the functional receptor for SARS 677 coronavirus. A first step in understanding SARS pathogenesis The Novel Coronavirus 2019 epidemic and kidneys Induction of pro-inflammatory cytokines (IL-1 and IL-6) and lung 681 inflammation by Coronavirus-19 (COVI-19 or SARS-CoV-2): anti-inflammatory 682 strategies Coincidence of COVID-19 epidemic and olfactory dysfunction 684 outbreak. medRxiv A portable and cost-effective microfluidic system for massively parallel 686 single-cell transcriptome profiling ACE2: from vasopeptidase to SARS virus 688 receptor Phenotypic analysis of mice lacking 690 the Tmprss2-encoded protease ACE2 receptor expression and severe acute respiratory syndrome 692 coronavirus infection depend on differentiation of human airway epithelia Caution on Kidney Dysfunctions of COVID-19 Patients. medRxiv Coronavirus as a possible cause of severe acute respiratory syndrome The surfactant system of the adult lung: physiology 699 and clinical perspectives Type II alveolar cell. Defender of the alveolus Deep Sequencing in Microdissected Renal 703 Tubules Identifies Nephron Segment-Specific Transcriptomes Familial renal glucosuria and SGLT2: from a mendelian trait to 706 a therapeutic target Abnormal expression and dysfunction of novel SGLT2 mutations identified 708 in familial renal glucosuria patients Renal Collectrin Protects against Salt-Sensitive Hypertension and Is 710 Downregulated by Angiotensin II Immune response profiling identifies autoantibodies specific to 712 Moyamoya patients Indoleamine 2,3-714 dioxygenase-2; a new enzyme in the kynurenine pathway. The international journal of 715 biochemistry & cell biology 41 Dioxygenase in Hepatitis C Virus Infection Aminopeptidase N is a major receptor for the entero-pathogenic 719 coronavirus TGEV Molecular determinants of species specificity in the 721 coronavirus receptor aminopeptidase N (CD13): influence of N-linked glycosylation The origin, transmission and clinical therapies on coronavirus disease 724 2019 (COVID-19) outbreak -an update on the status Interleukin-6 and the acute phase response. The 726 Mapping the Mouse Cell Atlas by Microwell-Seq COVID-19 in gastroenterology: a clinical perspective. Gut, 730 Epub ahead of print Functional exhaustion of antiviral lymphocytes in COVID-19 patients Epub ahead of print Mesenchymal stem (stromal) cells for treatment of ARDS: a phase 1 734 clinical trial Tocilizumab for induction and maintenance of remission in giant 736 cell arteritis: a phase 2, randomised, double-blind, placebo-controlled trial Senescence promotes in vivo 739 reprogramming through p16(INK)(4a) and IL-6 Characteristics of lymphocyte subsets and cytokines in peripheral blood 741 of 123 hospitalized patients with 2019 novel coronavirus pneumonia (NCP). medRxiv Structural basis of receptor recognition by SARS-CoV-2 Single-nucleus and single-cell transcriptomes compared in matched 745 cortical cell types STAR: ultrafast universal RNA-seq aligner Sambamba: fast processing 749 of NGS alignment formats SCANPY: large-scale single-cell gene expression 751 data analysis Comprehensive Integration of Single-Cell Data Identifying ChIP-seq enrichment using 755 MACS SeqKit: A Cross-Platform and Ultrafast Toolkit for 757 FASTA/Q File Manipulation Spatiotemporal immune zonation of the human kidney scRNA-seq assessment of the human lung, spleen, and esophagus 761 tissue stability after cold preservation A human liver cell atlas reveals heterogeneity and epithelial 763 progenitors