Single-cell meta-analysis of SARS-CoV-2 entry genes across tissues and demographics

Abstract

Angiotensin-converting enzyme 2 (ACE2) and accessory proteases (TMPRSS2 and CTSL) are needed for severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) cellular entry, and their expression may shed light on viral tropism and impact across the body. We assessed the cell-type-specific expression of ACE2, TMPRSS2 and CTSL across 107 single-cell RNA-sequencing studies from different tissues. ACE2, TMPRSS2 and CTSL are coexpressed in specific subsets of respiratory epithelial cells in the nasal passages, airways and alveoli, and in cells from other organs associated with coronavirus disease 2019 (COVID-19) transmission or pathology. We performed a meta-analysis of 31 lung single-cell RNA-sequencing studies with 1,320,896 cells from 377 nasal, airway and lung parenchyma samples from 228 individuals. This revealed cell-type-specific associations of age, sex and smoking with expression levels of ACE2, TMPRSS2 and CTSL. Expression of entry factors increased with age and in males, including in airway secretory cells and alveolar type 2 cells. Expression programs shared by ACE2⁺TMPRSS2⁺ cells in nasal, lung and gut tissues included genes that may mediate viral entry, key immune functions and epithelial–macrophage cross-talk, such as genes involved in the interleukin-6, interleukin-1, tumor necrosis factor and complement pathways. Cell-type-specific expression patterns may contribute to the pathogenesis of COVID-19, and our work highlights putative molecular pathways for therapeutic intervention.

Main

Coronavirus disease 2019 (COVID-19), caused by severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) infection, can manifest with pathologies in multiple systems, including the lungs and airways, gastrointestinal tract, kidney, liver and heart, and multi-organ failure^1,2,3. SARS-CoV-2 RNA has been found in nasal and throat secretions, saliva and stool specimens⁴.

Virion infection of host cells is initiated by the viral spike (S) protein binding to ACE2. ACE2 expression has been correlated with increased viral load in human cell lines^5,6 and in mice⁷. Viral infection further requires proteolytic cleavage of the S protein, and TMPRSS2 or cathepsin L, encoded by the CTSL gene, can provide this role for cellular entry⁸.

There is substantial variation in the clinical consequences of infection across individuals, from asymptomatic illness to death. Disease severity and mortality rise with age^9,10, with a slightly higher incidence and mortality in men². Children are significantly less likely to develop severe acute disease¹¹. Smoking may be associated with more severe disease¹². Finally, adults with preexisting cardiovascular disease may have higher rates of disease acuity and death².

Identifying specific cell types that can be infected by SARS-CoV-2 and relating SARS-CoV-2 entry factors to key covariates like age or sex could inform our understanding of COVID-19 tropism and heterogeneity in disease outcomes. The Human Cell Atlas (HCA) community has generated single-cell atlases of diverse tissues in healthy individuals, which can now be leveraged to enable such studies. Early analyses of HCA data revealed that some of the cells of the nasal passages, airways, lung parenchyma and gut express ACE2 and TMPRSS2 (refs. ^13,14), most notably nasal goblet cells and multiciliated cells¹³ in the airways and AT2 cells in the distal lung^13,15,16, and identified ACE2 and TMPRSS2 expression in colonic enterocytes^13,17.

Here, we chart the cell-type-specific expression patterns of ACE2 and accessory proteases by integrated analysis of 116 single-cell and single-nucleus RNA-sequencing (scRNA-seq and snRNA-seq) studies, including 31 studies of the lung and airways, and 85 studies of other diverse tissues. With the lung and airway studies, to our knowledge, we performed the first single-cell meta-analysis of atlas datasets associating cell-type-specific changes in expression level with age, sex and smoking status. We identify cross-tissue and tissue-specific gene programs enriched in immune-associated genes in ACE2⁺TMPRSS2⁺ cells and highlight other proteases that are significantly coexpressed with ACE2 and could play a role in infection.

Results

Double-positive ACE2 ⁺ TMPRSS2 ⁺ cells across the lung, airways and other organs associated with COVID-19

We enumerated the proportion of double-positive ACE2⁺TMPRSS2⁺ cells and ACE2⁺CTSL⁺ cells across 92 human scRNA-seq or snRNA-seq studies, including 7 of the lung and airways (Fig. 1, Methods and Supplementary Table 1 and 2). We surveyed published datasets, assigning cells to five broad categories (Fig. 1a,b, Extended Data Figs. 1 and 2 and Supplementary Table 1), and analyzed more finely annotated published and unpublished datasets (Methods, Fig. 1c,d and Supplementary Tables 1 and 3).

Fig. 1: A cross-tissue survey of ACE2⁺TMPRSS2⁺ cells shows enrichment in cells at reported sites of disease transmission or pathogenesis.

a,b, Double-positive cells were more prevalent in epithelial organs and cells. a, Proportion of ACE2⁺TMPRSS2⁺ cells per dataset (dots) from 21 tissues and organs (rows). b, Proportion of ACE2⁺TMPRSS2⁺ cells within cell clusters (dots) annotated by broad cell-type categories (rows) within each of the top seven enriched datasets. SMCs, smooth muscle cells. c,d, Significant coexpression of ACE2⁺TMPRSS2⁺ or ACE2⁺CTSL⁺ highlights cells from tissues implicated in transmission or pathogenesis. Significance of coexpression (dot size; −log₁₀ adjusted (adj.) P value), by two-sided Wald test (Methods); red border: false discovery rate (FDR) < 0.1) of ACE2⁺TMPRSS2⁺ (c) or ACE2⁺CTSL⁺ (d) and effect size (dot color, color bar) for finely annotated cell classes (columns) from diverse tissues (rows). Only tissues and cells in at least one significant coexpression relationship are shown (Methods). PDAC, pancreatic ductal adenocarcinoma; CD–PC, collecting duct principal cell; PEC, parietal epithelial cell; PCT, proximal convoluted tubule; TA, transit amplifying. e–h, In situ validation of double-positive cells in the lung, airways and SMGs (n = 3 donors per experiment, images of three randomly chosen areas per donor). Proximity ligation in situ hybridization (PLISH) and immunostaining (e and g) and quantification (error bars: standard errors; f and h) in human adult lung alveoli for ACE2 (white), TMPRSS2 (green) and CTSL (e) (red; total of 1,487 DAPI-positive cells examined for quantification (f)) and ACE2 (white), TMPRSS2 (green) and HTII-280 (g) (red; total of 482 HTII-280-positive cells examined for quantification (h)).

ACE2⁺TMPRSS2⁺ epithelial cells were most prevalent (in order) within the ileum, liver, lung, nasal mucosa, bladder, testis, prostate and kidney (Fig. 1a). Consistent with previous reports¹⁸, double-positive ACE2⁺TMPRSS2⁺ cells in the nose and airways were largely secretory goblet and multiciliated cells, and double-positive cells in the distal lung were largely alveolar type 2 (AT2) cells (Fig. 1c and Extended Data Fig. 3a). ACE2 and TMPRSS2 expression in secretory and AT2 cells is also supported by single-cell sequencing assay for transposase-accessible chromatin (scATAC-seq) from the primary carina and subpleural parenchyma of one adult individual, respectively, as well as secretory and multiciliated cells, and to a lesser extent some basal and tuft cells (Supplementary Fig. 1a–d; n = 3 samples per location, n = 1 patient; Methods). In a larger aggregation of lung and nasal datasets (Methods), we observed ACE2⁺TMPRSS2⁺ cells in various lung epithelial cells in pediatric samples (Extended Data Fig. 3b,c), also supported by single-cell chromatin accessibility by transposome hypersensitive sites sequencing (scTHS-seq)¹⁹ (Extended Data Fig. 4 and Methods). Significant double-positive ACE2⁺TMPRSS2⁺ cells in other tissues included enterocytes, pancreatic ductal cells, prostate luminal epithelial cells, brain oligodendrocytes, kidney proximal tubular cells and principal cells of the collecting duct, inhibitory enteric neurons, heart fibroblasts/pericytes, and fibroblasts and pericytes in multiple tissues (Fig. 1a–c). Notably, some of the cell types in which there were double-positive cells (including brain oligodendrocytes, multiciliated cells of the upper respiratory tract and sustentacular cells in olfactory epithelium) are cell types that also express MYRF (albeit not always significant triple expressors; Supplementary Fig. 2). MYRF is a transcription factor that induces expression of myelin basic protein and myelin oligodendrocyte glycoprotein²⁰. Autoimmune reactions against these proteins are known to potentially induce neurological symptoms (Discussion).

ACE2⁺CTSL⁺ coexpressing cells were enriched among AT1 and AT2 cells, enterocytes, ventricular cardiomyocytes and heart macrophages, as well as fibroblasts and pericytes in multiple tissues, including the placenta, heart, lung, kidney and enteric nervous system (ENS; Fig. 1d). We did not observe substantial ACE2 mRNA expression in scRNA-seq profiles in the bone marrow or cord blood (Fig. 1a,b), although there was ACE2 expression in alveolar and heart macrophages (Extended Data Fig. 5). Notably, in human placenta^21,22,23, ACE2 was expressed (1.4%) in maternal decidual/stromal cells, maternal pericytes and fetal extravillous trophoblasts, cytotrophoblasts and syncytiotrophoblasts in both first-trimester and term placenta (Fig. 1d). While there was little expression of TMPRSS2 (0.2%), CTSL was expressed in most cells (56%), and there were ACE2⁺CTSL⁺ double-positive cells (1.3%).

Cell-type-specific expression of additional proteases that may be relevant to infection

SARS-CoV-2 infects cells in the absence of TMPRSS2 (ref. ⁸), so additional proteases likely play roles in proteolytic cleavage of viral proteins for entry and egress. To predict such proteases, we tested the coexpression of ACE2 with each of 625 annotated human protease genes²⁴ in a declined donor transplant dataset (‘regev/rajagopal’; Supplementary Table 1). TMPRSS2 was significantly coexpressed in multiple lung epithelial cell types (Fig. 2a and Supplementary Tables 4 and 5), as were multiple members of the proprotein convertase subtilisin kexin (PCSK) family (Fig. 2a,b), including FURIN, PCSK2, PCSK5, PCSK6 and PCSK7 in AT2 cells. Proprotein convertases have known roles in coronavirus S-protein priming. We obtained similar results in an independent dataset from 40 samples (Extended Data Fig. 6a,b, Supplementary Table 1 and datasets ‘kropski’, ‘lafyatis/rojas’, ‘misharin_new’, ‘nawijn/teichmann’, ‘northwestern_misharin_ 2018reyfman’ and ‘sanger_meyer_2019madissoon’; collectively referred to as ‘aggregated lung’). As previously reported²⁵, the SARS-CoV-2 S protein has a polybasic motif in the S1/S2 region (Extended Data Fig. 6c) that corresponds to cleavage motifs of PCSK family proteases (Extended Data Fig. 6d)²⁵ and an additional site at the S2′ position (Extended Data Fig. 6e)²⁶.

Fig. 2: ACE2–protease coexpression and SARS-CoV-2 S-protein cleavage sites suggest a possible role for additional proteases in infection.

a, Multiple proteases were coexpressed with ACE2 in human lung scRNA-seq data. Scatterplot of significance (−log₁₀ adj. P value), by two-sided Wald test (Methods) and effect size of coexpression of each protease gene (dot) with ACE within each indicated epithelial cell type (color). Dashed line: significance threshold. TMPRSS2 and PCSK proteases that were significantly coexpressed with ACE2 are marked. b, ACE2–protease coexpression with PCSKs, TMPRSS2 and CTSL across lung cell types. Significance (dot size; −log₁₀ (adj. P value), by two-sided Wald test. (Methods)) and effect size (color) for coexpression of ACE2 with selected proteases (columns) across cell types (rows). NK, natural killer. c,d, Multiple proteases were expressed across lung cell types. c, Distribution of non-zero expression for AC2, PCSK and TMPRSS2 across lung cell types. White dot: median non-zero expression. d, Proportion of cells expressing ACE2, PCSK or TMPRSS2 across lung cell types, ordered by compartment. e, ACE2⁺PCSK⁺ double-positive cells across lung cell types. Fraction of different ACE2⁺PCSK⁺ or ACE2⁺TMPRSS2⁺ double-positive cells across lung cell types. Dots: different samples; line: median of non-zero fractions. f, ACE2–protease coexpression analysis for the 20 most significant human proteases in AT2 cells. Significance (dot size; −log₁₀ (adj. P value), by two-sided Wald test (Methods)) and effect size (color) for coexpression of ACE2 with different proteases (columns) across cell types (rows). g, Additional protease expression in ACE2⁺TMPRSS2⁺ double-positive cells. Significance (−log₁₀ adj. P value, by two-sided Wald test (Methods)) and fold change of differential expression for each human protease between ACE2⁺TMPRSS2⁺ double-positive versus double-negative cells within each indicated epithelial cell type (color). Significantly differentially expressed proteases within AT2 cells and PCSK across all epithelial cell types are highlighted.

FURIN, PCSK5 and PCSK7 were coexpressed with ACE2 across multiple lung cell types (Fig. 2c and Extended Data Fig. 6f). PCSK1 and PCSK2 were mostly restricted to neuroendocrine cells²⁷; PCSK2 was also detected in some AT2 cells (Fig. 2d and Extended Data Fig. 6g). In AT2 cells, proximal multiciliated cells and basal cells, dual expression of PCSK proteases with ACE2 was at fractions comparable to or higher than that of ACE2⁺TMPRSS2⁺ cells (Fig. 2e and Extended Data Fig. 6h). Coexpression was significant across other tissues (Extended Data Fig. 6i,j), including liver, ileum, kidney and nasal airways.

Because different host proteases may contribute to different stages of the viral life cycle²⁶, we examined the prevalence of ACE2⁺TMPRSS2⁺PCSK⁺ triple-positive cells in the lung. ACE2⁺TMPRSS2⁺PCSK7⁺ were the main triple-positive cells in multiciliated (0.75%) and secretory (0.72%) cells of proximal airways, and ACE2⁺TMPRSS2⁺FURIN⁺ triple-positive cells were the most common within AT2 cells (0.36%; Extended Data Fig. 6k). Among all known human proteases (Fig. 2f and Supplementary Fig. 3), cathepsins (CTSB, CTSC, CTSD, CTSL and CTSS), proteasome subunits (PSMB2, PSMB4 and PSMB5) and complement proteases (C1R, C2 and CFI) were the most commonly coexpressed with ACE2 in lung epithelial cell types.

Orthogonal validation of ACE2, TMPRSS2 and CTSL expression in the lungs

As ACE2 expression was quite low, we next validated some of these patterns by fluorescence in situ hybridization and immunofluorescence in tissue sections of airways and alveoli from three healthy donor lungs that were rejected for lung transplantation. ACE2, CTSL and TMPRSS2 were coexpressed by fluorescence in situ hybridization in alveolar cells, albeit at low levels (Fig. 1e,f). Co-staining with cell-type-specific markers showed ACE2 expression and TMPRSS2 expression in some HTII-280⁺ AT2 cells (Fig. 1g,h); we confirmed the latter by TMPRSS2 protein immunostaining (Extended Data Fig. 7d). TMPRSS2 protein was expressed at low levels in some AT1 cells (identified by AGER; Extended Data Fig. 7d). Some non-epithelial cells also expressed these three genes. We further validated ACE2 expression by bulk mRNA-seq of sorted AT2 cells (Extended Data Fig. 7e). Immunohistochemistry with antibodies used previously to block cellular viral entry specifically labeled adult pro-SFTPC⁺ AT2 cells (Extended Data Fig. 7c, Supplementary Table 6 and Methods).

Previous studies revealed that ACE2 is highly enriched in nasal and intestinal mucous cells^13,14. While mucous cells are relatively rare in healthy surface airway epithelium, they are abundant in submucosal glands (SMGs). Analysis by scRNA-seq of microdissected SMGs from healthy donors showed enrichment of ACE2, TMPRSS2 and CTSL in mucous cells (Extended Data Fig. 7f). In situ analysis confirmed the presence of ACE2 transcripts in acinar epithelial cells of the SMGs (Extended Data Fig. 7g) and cells expressing ACE2 in the large airway epithelium (Extended Data Fig. 7).

Association of ACE2, TMPRSS2 and CTSL expression in lung and airway cells with age, sex and smoking

We next asked how the expression of ACE2, TMPRSS2 and CTSL in specific cell subsets relates to three key covariates associated with more severe disease: age (older individuals), sex (males) and smoking²⁸. As no single dataset to date was sufficiently large, we aggregated samples across 31 scRNA-seq and snRNA-seq studies (Supplementary Table 2; 14 published^{16,18,29,30,31,32,33,34,35,36,37,38}; 17 not yet published^39,40 at the time of writing). This analysis spanned 1,320,896 cells from 228 individuals without known lung disease or from histologically normal-appearing lung adjacent to the site of disease, across 377 nasal, lung and airway samples from brushes, scrapings, biopsies, bronchoalveolar lavages, resections or entire lungs that could not be used for transplant or postmortem examinations (Fig. 3a). From unpublished data, we only obtained single-cell expression counts for the three genes (preprocessed by each data generator), total unique molecular identifier (UMI) counts per cell, cell identity annotations (which we harmonized to three resolution levels across studies; Fig. 3a,b, Supplementary Table 2, Extended Data Fig. 8 and Methods), and age, sex and smoking status (when ascertained). We modeled the association between the expression counts of each gene and age, sex and smoking status using a generalized linear model, accounting for technical variation arising from dataset-related factors and covariate interactions (Methods). We fitted this model within each cell type to non-fetal lung data of donors for whom smoking history was known (985,420 cells, 286 samples, 164 donors, 21 datasets) and fitted a model without smoking status covariates to the full non-fetal lung data (1,096,604 cells, 309 samples, 185 donors, 24 datasets).

Fig. 3: ACE2, TMPRSS2 and CTSL expression increases with age and in men, and shows cell-type-specific associations with smoking.

a, Samples in the aggregated lung and airway dataset partitioned to several classes by their cell composition. Percentage of cells by level 2 cell annotations (annotations with a preceding ‘1’ indicate coarse annotations of cells that had no annotation at level 2) across samples. The 377 samples were ordered by sample composition clusters (Methods). b, Schematic of key lung and airway epithelial cell types. c, Distribution of normalized ACE2 and TMPRSS2 expression across level 3 lung cell types in 1,031,254 cells from 228 donors. Red shading indicates the main cell types that expressed both ACE2 and TMPRSS2. d, Age, sex and smoking status associations with expression of ACE2 (blue), TMPRSS2 (orange) and CTSL (green) in level 3 epithelial cells. The effect size of the association is given as a log fold change (sex and smoking status) or the slope of log expression per year with age. As the age effect size is given per year, it is not directly comparable to the sex and smoking status effect sizes. Positive effect sizes indicate increases with age, in males, and in smokers. Colored bars: associations with an FDR-corrected P value < 0.05 (one-sided Wald test on regression model coefficients), consistent effect direction in pseudo-bulk analysis, and consistent results using the model with interaction terms (Methods). White bars: associations that did not pass all of the three above-mentioned evaluation criteria. Error bars: standard errors around coefficient estimates. Error bars are only shown for colored bars (indications or robust trends). Number of cells and donors, respectively, for each cell type: basal: 155,877 and 105; multiciliated lineage: 37,530 and 157; secretory: 22,306 and 140; rare: 2,676 and 71; submucosal secretory: 33,661 and 45; AT1: 29,973 and 101; AT2: 155,512 and 104. EC: endothelial cell; MDC: monocyte-derived cell.

For simplicity, we treated each cell as an independent observation. This implicitly combines variability in both donors and cells, and, because cells from the same donor are not truly independent observations, can result in inflated P values, especially when there are few donors for a particular cell type. To address this, account for covariate interactions and ensure robustness, we (1) used a simple noise model (Poisson) to reduce overfitting of donor variability; (2) confirmed that effect directions of significant associations were consistent in a pseudo-bulk analysis (modeling only donor variation; Methods and Supplementary Data 1–4); (3) confirmed summarized age, sex and smoking associations with a model including interaction terms (Methods and Supplementary Data 1–4); and (4) separated significant associations that passed all above confirmations into ‘robust trends’ and ‘indications’ depending on their robustness to holding out individual datasets (Methods and Supplementary Data 1–4). We focused on trends or indications in cell types where ACE2 and TMPRSS2 were coexpressed (Fig. 3c): airway epithelial cells (basal, multiciliated and secretory cells), AT1 and AT2 cells and SMG secretory cells.

We found robust trends of ACE2 expression with age, sex and smoking status in these cell types (Fig. 3d, Extended Data Fig. 9 and Supplementary Figs. 4–6; nonsmoking model results in Supplementary Figs. 7–10): ACE2 expression increased with age in AT2 cells, and was elevated in males in airway secretory cells and AT1 and AT2 cells. ACE2 levels were higher in past or current smokers in basal and submucosal secretory cells, and lower in AT2 cells (Fig. 3d). Analysis of bulk RNA-seq data from bronchial brushings⁴¹ indicated an upregulation of both ACE2 and TMPRSS2 in current smokers compared with former smokers (Extended Data Fig. 10). Furthermore, we found indications of increased ACE2 expression with age and in males in multiciliated cells, but those relied on inclusion of the dataset with the most cells and samples (‘regev/rajagopal’; Extended Data Fig. 9 and Methods). All above trends and indications for sex and age were validated in a simplified model without smoking status on the full non-fetal lung dataset (Supplementary Fig. 7, Supplementary Data 5–8 and Methods).

Examining joint trends of ACE2 and the protease genes within the same cell type, we found robust trends of ACE2 and TMPRSS2 coexpression increasing with age in AT2 cells, in males in AT1 cells, and an indication of the two genes being elevated in males in multiciliated cells (ACE2 indication dependent on the ‘regev/rajagopal’ dataset; Fig. 3d and Extended Data Fig. 9). ACE2 and CTSL showed robust trends of joint upregulation in males in AT2 cells, and in smokers in submucosal secretory cells. Indications of joint upregulation of these genes were found in males in AT1 cells, and in smokers in basal cells (Fig. 3d, Extended Data Fig. 9 and Methods). All joint trends for age and sex covariates were confirmed on the full non-fetal lung data using the simple model without smoking covariates (Supplementary Fig. 7).

An immune gene program in ACE2 ⁺ TMPRSS2 ⁺ cells in airway, lung and gut

Our previous analyses revealed immune signaling genes that covary with ACE2 and TMPRSS2 in airway and lung cells^13,14. To explore these in a broader context, we identified tissue and cell programs related to double-positive ACE2⁺TMPRSS2⁺ cells in the nasal epithelium, lung and gut (Supplementary Tables 7–10). Tissue programs are shared across double-positive cells from different cell types in one tissue; cell programs distinguish double-positive cells from the rest of the cells of the same type (Methods).

Tissue programs were enriched in pathways related to viral infection and immune response, including phagosome structure, antigen processing and presentation, and apoptosis (Fig. 4a,b, Supplementary Fig. 11a,b (for selected genes) and Supplementary Tables 7–10). These included CEACAM5 (lung, nasal and gut programs) and CEACAM6 (ref. ⁴²; lung), surface attachment factors for coronavirus S protein; SLPI (lung and nasal)⁴³; PIGR (lung and gut; may promote antibody-dependent enhancement via IgA⁴⁴); and CXCL17 (lung and nasal)⁴⁵. Tissue programs also had genes associated with cholesterol and lipid metabolic pathways and endocytosis (DHCR24, LCN2 and FASN), major histocompatibility complex I and II pathways⁴⁶, preparation against cellular injury (interferons; extracellular RNase: PLAC8 and TXNIP), complement (C3 and C4BPA), immune modulation (BTG1) and tight junctions (DST, CLDN3 and CLDN4).

Fig. 4: Tissue- and cell-type-specific gene modules in ACE2⁺TMPRSS2⁺ cells highlight immune and inflammatory features.

a,b, Tissue programs of ACE2⁺TMPRSS2⁺ cells in lung, gut and nasal samples. a, Selected tissue program genes. Node: gene; edge: program membership. Genes were selected heuristically for visualization (Methods). b, Enrichment was tested using a hypergeometric test exactly as performed by gprofiler in scanpy.queries.enrich (−log₁₀ adj. P value) of Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway gene sets in the full tissue programs. c–e, Cell programs of ACE2⁺TMPRSS2⁺ cells. c,d, Top 12 genes from each cell program recovered for different lung (c) or nasal (d) epithelial cell types (nodes; colors). Colored concentric circles: overlap with a gene in the top 250 significant genes in other cell types. ACE2 and TMPRSS2 were included even if not among the top 12 genes. e, Enrichment (−log₁₀ adj. P value) of KEGG disease and non-disease pathway gene sets in either highly significant genes across all tissues (top) or in specific tissues (lung and nose; bottom). f, Motif activity in immune transcription factors in ACE2⁺ cells. Significance (−log₁₀ adj. P value) of the top ten differential ‘motif activity scores’ (Methods) between epithelial ACE2⁺ cells or ACE2⁻ cells. Epithelial cells are: AT1, AT2, secretory, ciliated, ionocytes and neuroendocrine cells (highlighted in the gray shaded area in Supplementary Fig. 1a). n = 2 locations: primary carina and lung lobes; n = 3 samples per location; n = 1 patient. Motifs were extracted from the JASPAR2020 database, and the motif code is shown in each row. Dashed line: threshold for significance (adj. P value of 0.05). P values were calculated by logistic regression and likelihood ratio test, adjusted through Bonferroni correction (Methods).

Cell programs (Fig. 4c,d, Supplementary Fig. 12a–c and Supplementary Tables 7–10) were enriched in many of the same genes and pathways (for example, CEACAM5, CXCL17 and SLPI), and further captured unique functions, including tumor necrosis factor (TNF) signaling in lung secretory cells (for example, RIPK3; ref. ⁴⁷), lysosomal functions in lung secretory and multiciliated cells⁴⁸, the immunoproteasome (AT1 cells; Fig. 4c), cytokines, chemokines and their receptors (nasal goblet cells: CSF3, CXCL1, CXCL3, IL19 and CCL20; AT1 cells: IL1R1) and genes that encode surfactant proteins (AT2 cells: SFTPA and SFTPA2). Cell programs from multiple tissues (Fig. 4c,d) included genes related to TNF signaling, raising the possibility that anti-TNF therapy may impact the expression of ACE2 and/or TMPRSS2. Some of the genes encode proteins that are targets of known drugs⁴⁹ (for example, in lung secretory cells: C3, HDAC9, IL23A, PIK3CA, RAMP1 and SLC7A11), and other gene products have been shown to interact with SARS-CoV-2 proteins⁵⁰, for example, GDF15 (ref. ⁵¹), a central regulator of inflammation⁵², and yet others may be related to COVID-19 pathological features, including MUC1 (ref. ⁵³; in tissue and specific cell programs), IL6ST (lung tissue and gut enterocyte programs) and IL6 (AT2 program; Supplementary Fig. 12d). Other cell types, such as heart pericytes, were enriched for cells coexpressing ACE2 with IL6R or IL6ST (Supplementary Fig. 13). The immune-like programs of ACE2⁺ epithelial cells were also reflected in the regulatory features of the ACE2 locus by scATAC-seq (Fig. 4f). Cell–cell interaction analysis⁵⁴ (Methods) predicted interactions (Supplementary Table 11) between AT2 cells (overall or ACE2⁺TMPRSS2⁺) and myeloid cells through oncostatin, complement, interleukin (IL)-1 receptor and colony-stimulating factor signaling.

Conserved expression patterns in mouse models

Preclinical studies of SARS-CoV-2 infection and treatment require model systems that approximate human physiology. Transgenic mouse models that express human ACE2 (hereafter, hACE2) have been identified as a valuable resource to evaluate diverse therapeutics for COVID-19 (ref. ⁵⁵). We thus asked whether expression patterns of SARS-CoV-2 entry factors were similar in human and mouse model cell types of interest.

Ace2⁺Tmprss2⁺ and Ace2⁺Ctsl⁺ double-positive cells were present primarily in club and multiciliated cells in the airway epithelia of healthy mice⁵⁶ (Fig. 5a), consistent with human airways (Extended Data Fig. 3a), and increased from 2 to 4 months of age (Fig. 5a,b). Moreover, the expression patterns observed in scRNA-seq data of whole lungs from mice exposed daily to cigarette smoke for 2 months (Fig. 5c–k and Methods) are consistent with our observations in human airway epithelial cells (Fig. 3d and Extended Data Fig. 9a). Upon smoke exposure, there was a significant increase in the number Ace2⁺ cells and Ace2 expression in airway secretory cell numbers, but not AT2 cells (Fig. 5f–i). There was also agreement in expression patterns between the human placenta and mouse placenta development (Figs. 1c,d and 5l and Supplementary Fig. 14).

Fig. 5: Ace2, Tmprss2 and Ctsl expression in mouse in similar cell types, and follows similar patterns with age and smoking.

a, Gradual increase in Ace2 expression by airway epithelial cell type with age. Mean expression of Ace2 in different airway epithelial cells of mice of three consecutive ages. Shown are replicate mice (dots; n = 3 for each age), mean (bar) and error bars (s.e.m.). The effect of mouse age was tested using a two-sided Wald test (P values). TPM, transcripts per million. b, Increase in proportion of Ace2⁺Ctsl⁺ goblet and club cells with age. Percentage of Ace2⁺Ctsl⁺ cells in different airway epithelial cell types of mice of three consecutive ages. The effect of mouse age was tested using a Wald test (P values). c–k, Increase in Ace2 expression in secretory cells with smoking. Mice were exposed daily to cigarette smoke (CS) or filtered air (FA) as control for 2 months after which cells from whole-lung suspensions were analyzed by scRNA-seq (Drop-seq). AM, alveolar macrophages; IM, interstitial macrophages; DC, dendritic cells; LEC, lymphatic endothelial cells; CEC, capillary endothelial cells; EC, endothelial cells; Mono, monocytes. c,d, Uniform manifold approximation and projection analysis of scRNA-seq profiles (dots) colored by experimental group (c) or by Ace2⁺ cells and indicated double-positive cells (d). AT1 and AT2 cells and airway epithelial secretory and ciliated cells are marked by the red dashed line. Macro, macrophage; mono, monocyte. e, Marker genes of AT1, AT2, multicilitated and secretory cell clusters. f, The relative frequency of Ace2⁺ cells is increased by smoking in airway secretory cells but not AT2 cells. Relative proportion of Ace2⁺ and Ace2⁻ cells in smoke-exposed and control mice of different cell types (FA: n = 9 mice; CS: n = 5 mice; error bars represent 95% confidence intervals). g,h, Expression of Ace2 was increased in airway secretory cells (FA: 187 cells; CS: 62 cells), but not in AT2 cells (FA: 3,808; CS: 1,882). Distribution of Ace2 expression in secretory (g) and AT2 (h) cells from control and smoke-exposed mice (P value derived from a Wilcoxon rank-sum test; NS, not significant). i–k, Reanalysis of published bulk mRNA-seq⁶⁹ of lungs exposed to different daily doses of CS show increased expression of Ace2 (i), Tmprss2 (j) and Ctsl (k) after 5 months of chronic exposure; n = 8 mice per condition. Bars show the mean, and error bars show the standard error (**P = 0.0046, ***P = 0.0002 and ****P < 0.0001; one-way ANOVA with Dunnett’s multiple comparisons test, compared to FA group.) l, Expression in placenta. Mean expression (color) and proportion of expressing cells (dot size) of Ace2, Tmprss2 and Ctsl along with marker genes (Supplementary Fig. 14) in single- and double-positive cells from embryonic day (E) 9.5 to E18 of mouse placenta development.

Discussion

To the best of our knowledge, this study represents the first single-cell meta-analysis. Our meta-analysis provided the required power to uncover age, sex and smoking associations at single-cell resolution. The contrasting smoking associations of ACE2 across epithelial cell types show the importance of single-cell resolution, as downregulation in AT2 cells would have been otherwise masked by increases in airway epithelial signal in bulk RNA-seq⁵⁷. Although we have aggregated over 200 donors in our dataset, effects such as race, ethnicity, genetic ancestry, cumulative smoking or healthy tissue with a distal disease site may still confound the associations we have obtained.

Our models included tested covariates, technical covariates and interaction terms, allowing us to uncover complex associations (for example, sex and smoking associations are typically stronger for younger individuals; Supplementary Fig. 5). Modeling the smoking status of a donor was important to reduce background variation and account for the unbalanced distribution of covariates. Fitting this model required aggregating many datasets, harmonized by a consistent cell-type annotation. However, the annotation remains coarse in some cases, where cell labels still aggregate over considerable diversity, and can be further refined in the future. As the HCA grows and further datasets become available, our model could be extended to allow nonlinear associations with the tested covariates. Such associations may uncover, for example, distinct effects in the particularly affected geriatric population. While there is a trend of an increased proportion of ACE2⁺TMPRSS2⁺ cells with age (Extended Data Fig. 3b,c), this cannot be modeled reliably given the compositional diversity (Fig. 3a and Supplementary Fig. 15), potential confounders and limited sample sizes. Further metadata can help to address this.

Our findings in human and mouse models are consistent with respect to smoking and age associations. In line with our human data, we find an increase in Ace2 expression in maturing mice (2–4 months). Others have reported lower expression of entry factors in aged mice (24 months), showing potential limitations of mice as a model system⁵⁸.

Our comprehensive cross-tissue analysis expands on our^13,14,16,59 and others^60,61,62 earlier efforts, identifying cell subsets across tissues that may be implicated in transmission or pathogenesis. For example, double-positive cells in the SMGs may be a reservoir for viruses that escape from expulsion associated with severe cough in the airway luminal surface. Another intriguing hypothesis is that neurological symptoms^63,64,65 and Guillain–Barré syndrome⁶⁶ may arise as an autoimmune response to myelin antigens expressed by infected ACE2⁺TMPRSS2⁺ and ACE2⁺ cells that express myelin-producing genes (Supplementary Fig. 2 and Supplementary Table 7).

ACE2 and TMPRSS2 expression in lung, nasal and gut epithelial cells is associated with programs involving key immunological genes and genes related to viral infection. Expression of IL6, IL6R and IL6ST in lung epithelial cells raises the hypothesis that infection may trigger uncontrolled cytokine expression, as IL-6 levels were reported to increase with COVID-19 severity⁶⁷. The prediction of TNF, complement and IL-1 pathways may suggest a benefit for therapies that target these axes. The accessibility of binding sites for the transcription factors STAT and IRF in scATAC-seq data is consistent with interferon regulation of ACE2 expression in epithelial cells¹⁴ and with high activity of STAT1, STAT2, IRF1, IRF2, IRF5, IRF7, IRF8 and IRF9 in macrophage states, which increased in patients with severe COVID-19 (ref. ⁶⁸). Future lines of inquiry could include investigating the impact of lysosomal genes in lung secretory and multiciliated cells on viral infection and of RIPK3 expression in airway cells on necroptosis.

Finally, the expression of other potential accessory proteases may help pursue therapeutic hypotheses related to disruption of viral processing via protease inhibition. FURIN, PCSK5 and PCSK7 are more broadly expressed than TMPRSS2 across lung cell types (Fig. 2d) and across tissues (Extended Data Fig. 6i). Viral proteins may physically interact with PCSK6 (ref. ⁵⁰), which is significantly coexpressed with ACE2 in AT2 cells (Fig. 2b and Extended Data Fig. 6b). Because PCSK proteases are localized in different membrane compartments²⁷, they might process SARS-CoV-2 S proteins at different viral stages. Altogether, this could provide SARS-CoV-2 with immense flexibility in entry and egress.

Our meta-analysis provides a detailed molecular and cellular map to aid in our understanding of SARS-CoV-2 transmission, pathogenesis and clinical associations. We herein demonstrated how this can be done despite restrictions on data sharing. As the HCA progresses, we envision such meta-analyses in the context of other diseases, for example, by combining large healthy reference atlases with both epidemiological and genetic risk factors. In parallel, as new atlases are generated from COVID-19 tissues and models, their integration will further advance our understanding of this disease.

Methods

Patient samples

Sample collection underwent Institutional Review Board (IRB) review and approval at the institutions where the samples were originally collected. ‘Adipose_healthy_manton_unpublished’ was collected under IRB no. 2007P002165/1 (ORSP-3877). Tissue samples from breast, esophagus muscularis, esophagus mucosa, heart, lung, prostate, skeletal muscle and skin, referred to as ‘tissue_healthy_regev_snRNA-seq_unpublished’, were collected under ORSP-3635. Samples referred to as ‘eye_sanes_unpublished’ were collected under Dana-Farber/Harvard Cancer Center protocol no. 13-416 and Massachusetts Eye and Ear protocol no. 18-034H. Samples referred to as ‘kidney_healthy_greka_unpublished’ were collected under Massachusetts General Hospital IRB no. 2011P002692. Samples referred to as ‘liver_healthy_manton_unpublished’ were collected under IRB no. 02-240/ORSP-1702, as well as ORSP-2630 under ORSP-2169. Lung samples from smokers and non-smokers (41 samples from ten patients, 2–6 locations each) with the suffix ‘regev/rajagopal_unpublished’ were collected under Massachusetts General Hospital IRB no. 2012P001079/ORSP-3900 under ORSP-3490. Healthy and fibrotic lung samples with the suffix ‘xavier_snRNA-seq_unpublished’ were collected under Massachusetts General Hospital IRB no. 2003P000555 (CG-5242 under ORSP-3490), and Medoff no. 2015P000319 (CG-5145 under ORSP-3490). Pancreatic ductal adenocarcinoma samples were collected under C. Fernandez-del Castillo, 2003P001289 (CG-4692) under ORSP-3490 at Massachusetts General Hospital. Samples in the dataset ‘barbry’ were derived from a study that was approved by the Comité de Protection des Personnes Sud Est IV (approval no. 17/081), and informed written consent was obtained from all participants involved. All experiments were performed during 8 months, in accordance with relevant guidelines and French and European regulations. No deviations were made from our approved protocol named 3Asc (An Atlas of Airways at a single-cell level: NCT03437122). Lungs with chronic obstructive pulmonary disease and idiopathic pulmonary fibrosis in the ‘kaminski’ dataset were obtained from patients undergoing transplant, while healthy lungs were obtained from rejected donor lung organs that underwent lung transplantation at the Brigham and Women’s Hospital or donor organs provided by the National Disease Research Interchange. Patient tissues relating to the dataset ‘krasnow’ were obtained under a protocol approved by Stanford University’s Human Subjects Research Compliance Office (IRB no. 15166), and informed consent was obtained from each patient before surgery. The study protocol was approved by the Partners Healthcare IRB (protocol no. 2011P002419). Samples in the dataset ‘kropski_banovich’ were collected under Vanderbilt IRB nos. 060165 and 171657, and Western IRB no. 20181836 (ethics approval no. 2018/769-31). ‘Meyer_b’ samples were collected by the Cambridge Biorepository for Translational Medicine, under research ethics committee (REC) approval no. 15/EE/0152. Samples in the dataset ‘linnarsson’ are covered by 2018/769-31, approved by the Swedish Ethical Review Authority. Samples in the ‘misharin’ dataset were collected under STU00056197, STU00201137 and STU00202458 and approved by the Northwestern University IRB. Samples in the ‘rawlins’ dataset were obtained from terminations of pregnancy from Cambridge University Hospitals NHS Foundation Trust under permission from NHS Research Ethical Committee (96/085) and the Joint MRC/Wellcome Trust Human Developmental Biology Resource (grant no. R/R006237/1; www.hdbr.org; REC approval nos. 18/LO/0822 and 18/NE/0290). The studies relating to datasets ‘schultze’ and ‘schultze_falk’ were approved by the ethics committees of the University of Bonn and University Hospital Bonn (local ethics vote no. 076/16) and the Medizinische Hochschule Hannover (local ethics vote no. 7414/2017). Fifteen human tracheal airway epithelia in the ‘schultze’ dataset were isolated from de-identified donors whose lungs were not suitable for transplantation. Lung specimens were obtained from the International Institute for the Advancement of Medicine and the Donor Alliance of Colorado. The National Jewish Health IRB approved the research under protocol nos. HS-3209 and HS-2240. Samples in the ‘xu/whitsett’ dataset were provided through the federal United Network of Organ Sharing via the National Disease Research Interchange and International Institute for Advancement of Medicine and entered into the National Heart, Lung and Blood Institute (NHLBI) LungMAP Biorepository for Investigations of Diseases of the Lung at the University of Rochester Medical Center, overseen by the IRB as RSRB00047606 (Supplementary Tables 1 and 2).

Integrated analysis of published datasets

Publicly available (Supplementary Table 1) single-cell RNA-seq datasets were downloaded from the Gene Expression Omnibus (GEO). We searched the GEO for datasets that met the following criteria: (1) provided unnormalized count data; (2) were generated using the 10x Genomics Chromium platform; and (3) profiled human samples. These samples spanned a wide range of tissues, including primary tissues, cultured cell lines and chemically or genetically perturbed samples. Applying these filters increases standardization of sample as the vast majority were prepared using the same 10x Chromium instrument and Cell Ranger pipelines.

Datasets comprise one or more samples (individual gene expression matrices), which often correspond to individual experiments or patient samples. In total, this yielded 2,333,199 cells from 469 samples from 64 distinct datasets (Supplementary Table 1). To allow comparison across samples and datasets, we mapped genes using a common dictionary of gene symbols and excluded unrecognized symbols. If a gene from an aggregated master list was not found in a sample, the expression was considered to be zero for every cell in that sample.

After all datasets were collected, we quantified the percentage of cells with >0 UMIs for both ACE2 and TMPRSS2 or ACE2 and CTSL. For further analyses with broad cell classes, we only used datasets with more than 15 double-positive cells yielding 252,871 cells from 40 samples.

For integration across datasets, we used two levels of annotations. When possible, every sample was annotated with its tissue of origin based on the available metadata from the GEO. We excluded any sample for which tissue was not specified. For the smaller subset of 252,871 cells, we manually annotated cell clusters with broad cell-type classes using marker genes. These clusters were generated using the harmony-pytorch Python implementation (v0.1.1; https://github.com/lilab-bcb/harmony-pytorch/) of the Harmony scRNA-seq integration method⁷⁰ for batch correction and leiden clustering from the Scanpy package (v1.4.5). Clusters without clear markers distinguishing types were excluded from further analysis.

Data were processed using Scanpy. Individual datasets were log normalized (UMIs/10,000 + 1) by column sum and the log1p function (ln(10,000 × g_ij+ 1), where a gene’s expression profile, g, is the result of the UMI count for each gene, i, for cell j, normalized by the sum of all UMI counts for cell j. This data normalization step was only used for generating the clusters and cell-type annotations.

All other statistical tests for the integrated analysis were performed on the cell’s binary classification as double positive or not. For example, for a cell to be considered ACE2⁺, it has >0 ACE2 transcripts. Double-positive cells have >0 transcripts for both genes of interest. We used Fisher’s exact test to determine the statistical dependence between the expression of ACE2 and TMPRSS2 or CTSL and corrected for multiple testing using the Benjamin–Hochberg method over all tests for each gene pair.

Bronchial brushings from current and former smokers

Bronchial brushings were obtained from high-risk individuals undergoing lung cancer screening at ~1-year intervals by white light and autofluorescence bronchoscopy and computed tomography (n = 137 brushings from n = 50 patients; GSE109743) and profiled via RNA-seq as described previously⁴¹. Differential expression analysis of entry factors in former and current smokers was performed via voom-limma⁵¹ using the model:

$$Y_i \sim {\mathrm{smoking}} + {\mathrm{batch}} + {\mathrm{TIN}} + \left( {1\left| {\mathrm{patient}} \right.} \right),$$

where smoking denotes the encoded smoking status (‘current’ or ‘former’), batch refers to the experimental batch effect derived from the sequencing run, TIN represents the RNA integrity score, and (1|patient) is a random intercept per patient. Multiple-testing correction was performed via Benjamini–Hochberg to obtain an FDR-corrected P value.

Integrated coexpression analysis of high-resolution cell annotations across tissues

We compiled a compendium of published and unpublished datasets consisting of 2,433,890 cells from 21 tissues and/or organs including adipose, bone marrow, brain, breast, colon, cord blood, ENS, esophagus mucosa, esophagus muscularis, anterior eye, heart, kidney, liver, lung, nasal, olfactory epithelium, pancreas, placenta, prostate, skeletal muscle and skin. After the harmonization of cell-type annotations, ACE2-TMPRSS2 and ACE2-CTSL expression were assessed using a logistic mixed-effect model:

$$Y_i \sim ACE{\it{2}} + \left( {{1|{\rm{sample}}\_{\rm{id}}}} \right)$$

(1)

where Y_i was the binarized expression level of either TMPRSS2 or CTSL, and covariates were binarized ACE2 expression in cell i and a sample-level random intercept.

Models were fit separately for each cell type in each dataset. To avoid spurious associations in cell types with very few ACE2⁺ cells and due to very low expression of ACE2, we subsampled ACE2⁻ cells to the number of ACE2⁺ cells within each cell type and discarded cell types containing fewer than five cells expressing either ACE2 or the other gene being tested after the subsampling procedure. The significance of the association between ACE2 and TMPRSS2/CTSL was controlled for 10% FDR using the statsmodels Python package (v0.11.1)⁷¹. Data processing was performed using Scanpy (v1.4.6)⁷², and logistic models were fit using lme4 R package (v1.1.21)⁷³.

Single-cell ATAC-sequencing analysis

Library generation and sequencing

We performed single-cell ATAC-seq from primary carina and subpleural parenchyma of one individual (n = 3 samples per location). Libraries were generated using the 10x Chromium Controller and the Chromium Single Cell ATAC Library & Gel Bead Kit (1000111) according to the manufacturer’s instructions (CG000169-Rev C; CG000168-Rev B) with unpublished modifications relating to cell handling and processing. Briefly, human lung-derived primary cells were processed in 1.5 ml DNA LoBind tubes (Eppendorf), washed in PBS via centrifugation at 400g for 5 min at 4 °C and lysed for 3 min on ice before washing via centrifugation at 500g for 5 min at 4 °C. The supernatant was discarded and lysed cells were diluted in 1× diluted nuclei buffer (10x Genomics) before counting using trypan blue and a Countess II FL Automated Cell Counter to validate lysis. If large cell clumps were observed, a 40-µm Flowmi cell strainer was used before the tagmentation reaction, followed by Gel Bead-In-Emulsion generation and linear PCR as described in the protocol. After breaking the emulsion, the barcoded tagmented DNA was purified and further amplified to enable sample indexing and enrichment of scATAC-seq libraries. The final libraries were quantified using a Qubit dsDNA HS Assay kit (Invitrogen) and a High Sensitivity DNA chip run on a Bioanalyzer 2100 system (Agilent).

All libraries were sequenced using NextSeq High Output Cartridge kits and a NextSeq 500 sequencer (Illumina), and 10x scATAC-seq libraries were characterized by paired-end sequencing (2 × 72 cycles).

Initial data processing and quality control

Fastq files were demultiplexed using 10x Genomics Cell Ranger ATAC mkfastq (v1.1.0). We obtained peak-barcode matrices by aligning reads to GRCh38 (CR v1.2.0 pre-built reference) using Cell Ranger ATAC count. Peak-barcode matrices from six channels were normalized per sequencing depth and pooled using the Cell Ranger ATAC command ‘aggr’.

The aggregated, depth-normalized, filtered dataset was analyzed with Signac (v0.1.6; https://github.com/timoast/signac/), a Seurat⁷⁴ extension developed for the analysis of scATAC-seq data. All the analyses in Signac were run with a random number generator seed set as 1234. Cells that appeared as outliers in quality-control metrics (peak_region_fragments ≤ 750 or peak_region_fragments ≥ 20,000 or blacklist_ratio ≥ 0.025 or nucleosome_signal ≥ 10 or TSS.enrichment ≤ 2) were excluded from the analysis.

Normalization and dimensionality reduction

The aggregated dataset was processed with latent semantic indexing⁷⁵, that is, datasets were normalized using term frequency–inverse document frequency. Next, singular value decomposition, ran on all binary features, was used to embed cells in low-dimensional space. UMAP⁶⁹ was then applied for visualization, using the first 30 dimensions of the singular value decomposition space.

Gene activity matrix and differential motif activity analysis

A gene activity matrix was calculated as the chromatin accessibility associated with each gene locus (extended to include 2 kb upstream of the transcription start site, as described in the vignette ‘analyzing PBMC scATAC-seq’ (March 13, 2020; https://satijalab.org/signac/articles/pbmc_vignette.html), using as gene annotation the genes.gtf file provided together with Cell Ranger’s ATAC GRCh38-1.2.0 reference genome. For the motif analysis, we note that because epithelial cells with an accessible ACE2 locus tend to have a higher number of fragments in peaks than cells with inaccessible ACE2 (Supplementary Fig. 1e), consistent also with higher UMIs in scRNA-seq, some of the cells with inaccessible ACE2 could be false negatives, thus reducing our power.

Clusters were annotated using label transfer from matching scRNA samples or by literature/expert search of marker ‘active’ (that is, accessible) genes. Differential motif activity analysis was performed using Signac’s implementation of ChromVAR⁷⁶, with motif position frequency matrices from JASPAR2020 (ref. ⁷⁷; http://jaspar.genereg.net/) selecting transcription factor motifs from human (species = 9606), broadly following the vignette ‘motif analysis with Signac’ (https://satijalab.org/signac/articles/motif_vignette.html). Cells were identified as positive for ACE2 and/or TMPRSS2 (that is, with the loci accessible) if at least one fragment was overlapping with the gene locus or within 2 kb upstream. Differential activity scores between epithelial cells positive for ACE2 (with the above-mentioned definition of ‘positive’) and nonexpressing ACE2 was performed with the FindMarkers function of Seurat (v3.1.1), using the test function set to ‘LR’ (that is, logistic regression) and the number of counts per peak as the latent variable. The function constructs a logistic regression model predicting group membership based on each motif score individually and compares this to a null model with a likelihood ratio test. An adjusted P value is the result of Bonferroni correction.

Immunohistochemistry and proximity ligation in situ hybridization

PLISH was performed as described previously⁷⁷. Briefly, frozen human trachea and distal lung sections were fixed with 4% paraformaldehyde for 20 min, treated with protease (20 μg ml⁻¹ proteinase K for lung or pepsin for trachea for 9 min) at 37 °C, and dehydrated with ethanol. The sections were incubated with gene-specific oligonucleotides (Supplementary Table 6) in hybridization buffer (1 M sodium trichloroacetate, 50 mM Tris (pH 7.4), 5 mM EDTA and 0.2 mg ml⁻¹ heparin) for 2 h at 37 °C. Common bridge and circle probes were added to the section and incubated for 1 h followed by T4 ligase reaction for 2 h. Rolling circle amplification was performed by using phi29 polymerase (30221, Lucigen) for 12 h at 37 °C. Fluorophore-conjugated detection probe was applied and incubated for 30 min at 37 °C. For the combination of PLISH and immunostaining, sections were incubated with primary antibody for HTII-280 (Terrace Biotech, TB-27AHT2-280), pro-SFTPC (Millipore, ab3786) or ACTA2 (Sigma, F3777) for 1 h at room temperature. Sections were incubated with goat anti-mouse IgM secondary antibody (Thermo Fisher Scientific, A21044) or donkey anti-rabbit IgG secondary antibody (Thermo Fisher Scientific, A32795) for 45 min at room temperature, and then sections were mounted in medium containing DAPI. We imaged three representative areas per patient for three patients in total for the images and quantification shown in Fig. 1 and one representative area for a single patient for Extended Data Fig. 7a,c,d,g. Images were captured using an Olympus confocal microscope FV3000 with Olympus FLUOVIEW FV31S-SW (v2.1.1.98) using a ×20 or ×60 objective.

Bulk mRNA sequencing of sorted cell populations from human lung

Human lung tissue was received from New England Donor Services under the Massachusetts General Hospital approved institutional review board protocol. Tissue was used from three individual patients with no history of lung disease or smoking. Primary bronchus and a piece of right lung lobe were manually dissected out. Single cells were dissociated in HBSS media (Sigma, 55021C) containing collagenase (225 units ml^–1), dispase (2.5 units ml^–1), elastase (2 units ml^–1), pronase (1 mg ml^–1), DNAse (1 unit ml^–1) and Y-27632 (5 μM) inhibitor. Cell suspension was treated with ACK lysis buffer (Thermo Fisher Scientific A1049201) for 2 min on ice to remove red blood cells. Large airway basal cells were isolated from primary bronchus tissue while small airway basal cells and alveolar type 2 (AT2) cells were isolated from lung lobe tissue. Basal cells from both tissue suspensions were isolated using the anti-human CD271 MicroBeads kit (Miltenyl Biotech, 130-099-023) following the manufacturer’s protocol. AT2 cells were isolated with anti-HT2-280 antibody (Terrace Biotech, TB-27AHT2-280) and anti-mouse IgM MicroBeads kit (Miltenyl Biotech, 130-047-301) following the manufacturer’s protocol. For mRNA-seq sample preparation, total RNA was isolated using Trizol reagent (Thermo Fisher Scientific, 15596026) following the manufacturer’s instructions. Quality and quantity of total RNA was assayed using Nanodrop and an Agilent Bioanalyzer. mRNA-seq libraries were prepared using the TrueSeq protocol from Illumina. For mRNA-seq analysis, the raw fastq files were mapped to the human genome using Tophat (with bowtie2) (version 2.1.1). The output files were processed through Cufflinks (version 2.2.1.3)) and Cuffdiff (version 2.2.1.6) to conduct differential gene expression analysis.

Transposome hypersensitive sites sequencing on human pediatric samples

THS-seq was performed as previously reported¹⁹ on human pediatric samples (full gestation, with no known lung disease) collected at day 1 of life, and again at 14 months, 3 years and 9 years (n = 1 at each time point).

Integrated analysis for associating ACE2, TMPRSS2 and CTSL expression with age, sex and smoking status in nasal, airway and lung cells

To assess the association of age, sex, and smoking status with the expression of ACE2, TMPRSS2 and CTSL, we aggregated 31 scRNA-seq datasets of healthy human nasal and lung cells, as well as fetal samples containing the expression counts of only the three genes. Aggregation of these datasets was enabled by harmonizing the cell-type labels of individual datasets and dataset concatenation within Scanpy⁷¹ (v1.4.5.1). We harmonized annotations manually on the basis of provided cell-type labels together with data contributors using a preliminary ontology generated on the basis of five published datasets^{31,32,33,36,38} with three levels of annotations. Level 1 has the lowest resolution and distinguishes epithelial from stromal/mesenchymal, endothelial and immune cells. Level 2 breaks up each of the level 1 categories in the coarsest available further observed annotations. Level 3 in turn splits up the observed level 2 annotations where finer annotations were available (Supplementary Table 2; consent to publish was obtained from all contributors). To compare AT2 cells and their possible fetal progenitors, we mapped progenitor cells labeled ‘AT2-like’ and ‘SpC⁺ progenitors’ to the AT2 label. We further harmonized metadata by collapsing the smoking covariate into ‘has smoked’ and ‘has never smoked’ and by taking the mean age where only age ranges were given. This resulted in a dataset of 1,320,896 cells and three genes in 377 samples from 228 donors (the cell by three-gene count matrix with annotations is available on the Single Cell Portal (https://singlecell.broadinstitute.org/single_cell/study/SCP1257)). We divided the data into fetal (136,450 cells, 41 samples and 34 donors), adult nasal (57,548 cells, 20 samples and 18 donors) and adult lung (1,126,898 cells, 316 samples and 187 donors) datasets based on the metadata provided.

To get an overview of sample diversity, we clustered the samples using the proportion of cells in level 2 cell types as features. Clustering was performed using louvain clustering (resolution of 0.3; louvain package v0.6.1) on a k-nearest-neighbor graph (k = 15) computed on Euclidean distances over the top five principal components of the cell-type proportion data within Scanpy. This produced four clusters. Sample cluster labels were assigned based on cell-type compositions and metadata for anatomical location that were obtained from the published datasets and via input from the data generators.

Within non-fetal datasets, we modeled the association of age, sex and smoking status with gene expression for ACE2, TMPRSS2 and CTSL within each cell type using a generalized linear model with the log total counts per cell as offset and Poisson noise as implemented in Statsmodels⁷¹ (v0.11.1) and using a Wald test from Diffxpy (www.github.com/theislab/diffxpy/; v0.7.3, batchglm v0.7.4). Specifically, we fit the model:

$$Y_{ij} \sim {\mathrm{age}} + {\mathrm{sex}} + {\mathrm{age:sex}} + {\mathrm{smoking}} + {\mathrm{sex:smoking}} + {\mathrm{age:smoking}} + {\mathrm{dataset}},$$

(2)

which models effects of age, sex and smoking while accounting for potential interactions between covariates and the uneven distribution of covariates across the dataset. Here, Y_ij denotes the raw count expression of gene i in cell j; age, sex and smoking denote the modeled covariates; and ‘age:sex’, ‘sex:smoking’ and ‘age:smoking’ represent the interaction terms between these covariates. The interaction terms model whether there is a difference in the smoking effect in men and women, and likewise whether the age effect is different for smokers and non-smokers. We included the ‘dataset’ term to model the technical variation (for example, sampling and processing differences) between the diverse datasets, and the log total count per cell was used as an offset. Here, the total counts were scaled to have a mean of 1 across all cells before the log was taken. Due to the inclusion of interaction terms, the complex interaction model (2) fits the overall effects of age (k_age), sex (k_sex) and smoking (k_smoking) as linear functions of the other two covariates, given by the equations:

$$k_{\mathrm{age}}\left({\mathrm{sex,smoking}} \right) = \beta _{\mathrm{age}} + {\mathrm{sex}}\ \beta _{\mathrm{age:sex}} + {\mathrm{smoking}}\ \beta _{\mathrm{age:smoking}},$$

$$k_{\mathrm{sex}}\left( {\mathrm{age,smoking}} \right) = \beta _{\mathrm{sex}} + {\mathrm{age}}\ \beta _{\mathrm{age:sex}} + {\mathrm{smoking}}\ \beta _{\mathrm{sex:smoking}},$$

$$k_{\mathrm{smoking}}\left( {\mathrm{age,sex}} \right) = \beta _{\mathrm{smoking}} + {\mathrm{age}}\ \beta _{\mathrm{age:smoking}} + {\mathrm{sex}}\ \beta _{\mathrm{sex:smoking}},$$

Here, β_age and β_age:sex represent the model coefficients for age and the interaction of age and sex in model 2, respectively, and age denotes the age covariate. Sex and smoking covariates were converted into a one-hot encoded format such that sex = 0 denoted females and smoking = 0 denoted non-smokers. As linear dependencies on covariates can be summarized by showing two values per covariate, we displayed effect sizes for the overall age, sex and smoking associations by computing k_age, k_sex and k_smoking for sex ∈ {0,1}, smoking ∈ {0,1} and age ∈ {31,62} (the first and third quartiles of the age distribution). Standard errors for these effects were computed with the variance–covariance matrix Σ using $\textrm {SE} = \sqrt {C^T{\Sigma}C}$, where SE is the standard error and C is the vector of covariate values used to compute the respective overall effect (for example, k_age). P values were obtained using a Wald test, and correction for multiple testing was performed over all tests on the same cell-type data using the Benjamini–Hochberg method. To fit this model, we pruned the data to contain only datasets that had at least two donors and for which smoking status metadata were provided. This resulted in a dataset of 985,420 cells and 286 samples from 164 donors for adult lung data. Only 15 donors remained for adult nasal data after this filtering, which we deemed too few to obtain robust results. To obtain cell-type-specific associations, the above model was fit within each cell type for all cell types with at least 1,000 cells.

While cells from different donors are not truly independent observations, model 2 treats them as such and thus models cellular and donor variation jointly. As donor variation tends to be larger than single-cell variation, when most cells come from few donors (either there are few donors, or few donors contribute most of the cells), this can lead to inflation of P values. To counteract this effect, we verified that significant associations were consistent when modeling only donor variation via pseudo-bulk analysis (Supplementary Data 1–4). Furthermore, we tested whether effects were dependent on few donors by holding out datasets.

Pseudo-bulk data were generated by computing the mean for each gene expression value and the number of UMIs (nUMI) covariate for cells in the same cell type and donor. After filtering as described above, model 2 was fit to the data (Supplementary Data 1–4). In contrast to the single-cell model, pseudo-bulk analysis underestimates certainty in modeled effects as uncertainty in the pseudo-bulk means are not taken into account when estimating background variance. Thus, we used only effect directions from pseudo-bulk analysis to validate single-cell associations. In further analysis, we regarded only those associations as confirmed by pseudo-bulk analysis, where the FDR-corrected P value in the single-cell model was below 0.05, and the sign of the estimated effect was consistent in both the single-cell and the pseudo-bulk analysis.

We further separated significant associations into robust trends and indications depending on the holdout analysis. A significant association was regarded as a robust trend if the effect direction is consistent when holding out any dataset when fitting the model (without considering the P value). In the case that holding out one dataset caused the maximum-likelihood estimate of the coefficient to be reversed, we denoted this as the effect no longer being present, which characterized the association as an indication. Two dataset holdouts led to indications in our analysis: the largest declined donor transplant dataset (‘regev-rajagopal’; Supplementary Table 2; most cells and most samples; indication in ACE2 multiciliated lineage age and sex associations, and CTSL AT1 sex association) and a declined donor tracheal epithelium dataset (‘seibold’; Supplementary Table 2; most donors in the smoking analysis; CTSL basal smoking association).

At least four values for each covariate are required to describe a single association in model 2 (for example, male nonsmoker, female nonsmoker, male smoker and female smoker for the k_age effect). To summarize these effects and present a single association per covariate, we also fit the simplified model:

$$Y_{\mathrm{ij}} \sim {\mathrm{age}} + {\mathrm{sex}} + {\mathrm{smoking}} + {\mathrm{dataset}}$$

(3)

As in model 2, the logarithmized, scaled total counts per cell were used as an offset, data were filtered as described, and multiple-testing correction was performed via Benjamini–Hochberg. To increase the robustness of our reported associations, we again performed pseudo-bulk and holdout analysis. Additionally, to still account for covariate interactions, we discarded associations where the complex model 2 and the simplified model 3 results were inconsistent. Here, consistency was defined by two criteria: at least one model 2 indication or robust trend in the same direction as the model 3 effect, and no model 2 indication or robust trend in the opposite direction to the model 3 effect.

As metadata on smoking status were only available for a subset of the data, we also fitted a reduced version of models 2 and 3 without the smoking covariate on a larger dataset to confirm sex and age associations (Supplementary Data 5–8). The nonsmoking model was fit on 1,096,604 cells in 309 samples from 185 donors of adult lung data. Again, log total count (scaled) was used as an offset, pseudo-bulk and holdout analysis was performed, and associations from the simple model were tested for consistency with the complex model.

Normalizing ACE2 ⁺ TMPRSS2 ⁺ double-positive fractions of human lung samples

Proportions of ACE2⁺TMPRSS2⁺ cells (Extended Data Fig. 3a and Supplementary Fig. 15) were normalized to account for differences in total UMI counts. Normalization was completed per donor, for each cell type by calculating $\frac{{X_{i,j}}}{{N_{i,j}}} \times 10,000$, where X_i,j is the double-positive fraction of cell type i in donor j, and N_i,j represents the median total UMI count of cells of type i in donor j.

Identification of gene programs using feature importance for a random forest trained to classify ACE2 ⁺ TMPRSS2 ⁺ versus ACE2 ⁻ TMPRSS2 ⁻ cells

To infer tissue programs, we trained a random forest classifier to discriminate between double-positive and double-negative cells (excluding ACE2 and TMPRSS2; a 75:25 class-balanced test:train split), generalizing across multiple cell types in one tissue, and ranked genes according to their importance scores in the classifier. To infer cell programs, we performed differential expression analysis between double-positive and double-negative cells within each cell subset.

Importantly, these methods do not assume that ACE2⁺TMPRSS2⁺ cells form a distinct subset within each cell type. Rather, our goal is to leverage the variation among single cells within a single type to identify gene programs that are co-regulated with ACE2 and TMPRSS2 within each expressing cell subset.

For each of the lung, nasal and gut datasets, we labeled the cells with non-zero counts for both ACE2 and TMPRSS2 as double-positive cells, and the cells with zero counts for both ACE2 and TMPRSS2 as double-negative cells. Within each tissue, we identified cell types with greater than ten double-positive cells, and for each of these cell types, we selected the genes with increased expression (log fold change > 0) in double-positive cells compared to double-negative cells (to focus on important ‘positive’ features). We trained a classifier with a 75:25 train:test split to classify the double-positive cells from double-negative cells within each of these cell types using the ‘sklearn’ (v0.21.3)⁷⁸ ‘RandomForestClassifier’ function with the following parameters: ‘n_estimators’ set to 100, the ‘criterion as gini’, and the ‘class_weight’ parameter set to ‘balanced_subsample’. We first trained individual classifiers separately for each of the cell types and pooled genes with positive feature importance values (using the ‘feature_importance’⁷⁹ field in the trained RandomForestClassifier object) to train a final double-positive cell versus double-negative classifier across each tissue. We used the top 500 genes, as ranked by their feature importance scores, to define the signature for the gene expression program of double-positive cells for the tissue. This procedure was carried out in lung, nasal and gut datasets, yielding tissue-specific signatures for gene expression programs of double-positive cells from each tissue.

For visualization purposes only, we generated network diagrams using the ‘networkx’ (v2.2) tool with the ForceAtlas2 graph layout algorithm⁸⁰. We scored genes that appeared in signatures for multiple tissues by their aggregated feature importance (using a plotting heuristic method that used the sum of importance ranks for genes in individual tissues and by assigning a large valued rank (10,000) to a gene that did not appear in a particular tissue) and selected the top ten genes that were shared by each pair of tissues or shared by all tissues along with additional genes that included the ones unique to each tissue’s signature to plot in the network visualization. The GO terms enriched in the gene expression programs shared by double-positive cells across tissues were found using g:Profiler (v1.0.0)⁸¹ using the ‘scanpy.queries.enrich’ tool.

This analysis was performed in two ways: on the original data, as well as after accounting for differences in distribution of the nUMI per cell between double-positive cells and double-negative cells. This was performed by binning the nUMI distribution in the double-positive cells for each tissue into 100 bins and then randomly sampling from the nUMI distribution for the double-negative cells in each bin to match the distribution of the double-positive cells in that bin. The nUMI distributions before and after the matching procedure are shown in Supplementary Fig. 11b.

Identification of gene programs enriched in double-negative cells versus double-negative cells using regression

In parallel, we used a regression framework to recover gene modules enriched in double-positive versus double-negative cells (Fig. 4c,d and Supplementary Fig. 12a,b) in the nasal, lung and gut datasets. We first restricted our analysis to cell subsets derived from at least two donor individuals that each contained a mixture of double-negative and double-positive cells (nawijn nasal: multiciliated; goblet; regev/rajagopal lung: AT1, AT2, basal, multiciliated, and secretory; aggregated lung: AT2, multiciliated and secretory; regev/xavier colon: BEST4⁺ enterocytes, cycling TA (transit amplifying), enterocytes, immature enterocytes 2 and TA-2). For each of these cell subsets, we then used MAST (v1.8.2)⁸² to fit the following regression model to every gene with cells as observations:

$$Y_i \sim X + \left( {1|S} \right),$$

where Y_i is the expression level of gene i in cells, measured in units of log₂(transcripts per 10,000 reads (TP10K) + 1), X is the binary coexpression state of each cell (that is, double-positive versus double-negative cells), and S is the donor that each cell was isolated from. To control for donor-specific effects (that is, batch effects), we used a mixed model with a random intercept that varies for each donor. To fit this model, we subsampled cells from double-positive and double-negative groups to ensure that both the donor distribution and the cell complexity (that is, the number of genes per cell) were evenly matched between the two groups, as follows. First, for each subset, we restricted our analysis to donors containing at least two double-negative and two double-positive cells. Using these samples, we partitioned the cells into ten equally sized bins based on cell complexity and subsampled double-negative cells from each bin to match the cell complexity distribution of the double-positive cells. Finally, we fit the mixed model (above), controlling for both donor and cell complexity.

To build gene modules for double-positive cells, we prioritized genes by requiring that they be expressed in at least 10% of double-positive cells, and to have a model coefficient greater than 0 with an FDR-adjusted P value of less than 0.05 (for the combined coefficient in the hurdle model). After this filtering step, genes were ranked by their model coefficient (that is, estimated effect size). The top 12 genes were selected for network visualization within each cell type (Fig. 4c,d and Supplementary Fig. 12a,b). In three cases (gut cycling TA, TA-2 and BEST4⁺ cells), RP11 antisense genes were flagged and excluded from visualizations. To visualize overlap across each network, we indicated whether each gene was among the top 250 genes from each of the other cell types. Putative drug targets were identified by querying the DrugBank database⁴⁹. Gene-set enrichment analysis was performed using the R package Enrichr (v1.0)⁸³, selecting the top 25 genes from each cell type for the pan-tissue analysis (‘all’ category; Fig. 4e) and the top 50 genes from each cell type for the tissue-specific analyses (‘nose’ and ‘lung’ categories; Fig. 4e). We note a few limitations that may influence our results, including nonuniform sampling across donors, variation in cell compositions across regions (for example, distal lung versus carina) and additional cellular heterogeneity that the current level of broad subset annotation may not have captured.

Cell–cell interaction analysis

CellphoneDB⁵⁴ (v.2.0.0) was run with default parameters on the ten human lung samples of the regev/rajagopal dataset (41 samples from 10 patients, 2–6 locations each), analyzing the cells from each dissected region separately. For each sample (patient/location combination) and for each cell type, we distinguished double-positive cells (ACE2 > 0 and TMPRSS2 > 0) from all others. Only interactions highlighted as significant (that is, present in the ‘significant means’ output P < 0.05) from CellphoneDB were considered. AT2 cells and myeloid cells were present in lung lobe samples from all ten patients, whereas samples from five patients contained both ACE2⁺TMPRSS2⁺ double-positive AT2 cells and myeloid cells.

Coexpression patterns of additional proteases and IL6/IL6R/IL6ST

ACE2–protease coexpression (Fig. 2 and Extended Data Fig. 5) and ACE2–IL6/IL6R/IL6ST coexpression (Supplementary Fig. 13) were tested via the logistic mixed-effects model described above (model 1).

Mouse smoke exposure experiments

For these experiments, 8- to 10-week-old pathogen-free female wild-type C57BL/6 mice were obtained from Charles River and housed in rooms maintained at constant temperature and humidity with a 12-h light cycle. Animals were allowed food and water ad libitum. All animal experiments were approved by the ethics committee for animal welfare of the local government for the administrative region of Upper Bavaria (Regierungspräsidium Oberbayern) and were conducted under strict governmental and international guidelines in accordance with EU Directive 2010/63/EU. The female C57BL/6 mice (n = 5) were whole-body exposed to 100% mainstream cigarette smoke at a particle concentration of 500 mg/m³, generated from 3R4F research cigarettes (filter removed; Tobacco Research Institute, University of Kentucky), for 50 min twice daily, 5 d per week for 2 months to mimic human smoking habits⁸⁴. Control mice (n = 3) were exposed to filtered air, but exposed to the same stress as mice exposed to cigarette smoke.

Reporting Summary

Further information on research design is available in the Nature Research Reporting Summary linked to this article.