Single-cell multiomic profiling of human lungs reveals cell-type-specific and age-dynamic control of SARS-CoV2 host genes

Abstract

Respiratory failure associated with COVID-19 has placed focus on the lungs. Here, we present single-nucleus accessible chromatin profiles of 90,980 nuclei and matched single-nucleus transcriptomes of 46,500 nuclei in non-diseased lungs from donors of ~30 weeks gestation,~3 years and ~30 years. We mapped candidate cis-regulatory elements (cCREs) and linked them to putative target genes. We identified distal cCREs with age-increased activity linked to SARS-CoV-2 host entry gene TMPRSS2 in alveolar type 2 cells, which had immune regulatory signatures and harbored variants associated with respiratory traits. At the 3p21.31 COVID-19 risk locus, a candidate variant overlapped a distal cCRE linked to SLC6A20, a gene expressed in alveolar cells and with known functional association with the SARS-CoV-2 receptor ACE2. Our findings provide insight into regulatory logic underlying genes implicated in COVID-19 in individual lung cell types across age. More broadly, these datasets will facilitate interpretation of risk loci for lung diseases.

Keywords: COVID-19; cis-regulatory elements; developmental biology; gene regulation; genetics; genomics; human; human sequence variants; lung; single cell RNA/ATAC-seq.

Figure 1.. Single-nucleus atlas of chromatin accessibility and transcriptomes in the human lungs.

(A) UMAP (Uniform Manifold Approximation and Projection) embedding (McInnes et al., 2018) and clustering results of snATAC-seq data from 90,980 single-nucleus chromatin profiles from ten donors: premature born (30 week^GA for gestational age, n = 3), 4-month-old (n = 1), three yo (n = 3) and 30 yo (n = 3). For library quality control see Figure 1—figure supplement 1A–D. (B) Dot plot of marker genes from snRNA-seq used for cluster annotation. For additional genes see Figure 1—figure supplement 2A. (C) UMAP embedding (McInnes et al., 2018) and clustering result of 46,500 snRNA-seq data from nine donors: premature born (30 week^GA), three yo, 30 yo, n = 3 per time point, identifies 31 clusters. Each dot represents a nucleus. Spread-out gray dots correspond to nuclei of unclassified cells. For library quality control see Figure 1—figure supplement 1E–H. (D) Dot plot of marker genes from snRNA-seq used for cluster annotation. For additional genes see Figure 1—figure supplement 2B.

Figure 1—figure supplement 1.. Quality control of snATAC-seq and snRNA-seq datasets.

(A) Nuclei with less than 1000 uniquely mapped reads were filtered from snATAC-seq datasets. (B) Number of nuclei selected for snATAC-seq libraries after quality control and removing potential doublets. (C) Average number of reads per nucleus. (D) Fraction of reads in peak regions per dataset. (E) Representative UMI barcode distribution output from CellRanger pipeline (10x Genomics) for snRNA-seq libraries from human lungs. (F) Sequencing saturation of snRNA-seq libraries. (G) Number of nuclei selected for snRNA-seq libraries after quality control and removing potential doublets. (H) Genes detected per nucleus. All data are represented as mean ± SD.

Figure 1—figure supplement 2.. Expression and chromatin accessibility at marker gene loci used for annotation.

(A) Dot plot of marker genes from snATAC-seq used for cluster annotation. Extension of dot plot in Figure 1B (B) Dot plot of marker genes from snRNA-seq used for cluster annotation. Extension of dot plot in Figure 1D.

Figure 2.. snATAC-seq analysis of human lungs reveals candidate cis-regulatory elements for ACE2 and TMPRSS2.

(A) Dot plot illustrating cluster-specific gene expression of candidate SARS-CoV-2 cell entry genes. For violin plots illustrating cluster-specific gene expression please see Figure 2—figure supplement 1A–E. (B) Dot plot illustrating cluster-specific gene body chromatin accessibility of candidate SARS-CoV-2 cell entry genes. (C) Union set of peaks (vertical lines) identified in all clusters surrounding ACE2 and 15 peaks that showed co-accessibility with the ACE2 promoter (red lines, co-accessibility score >0.05) via Cicero (Cusanovich et al., 2018). (D) Zoom into ACE2 locus and genome browser tracks of snATAC-seq signal (Robinson et al., 2011). ACE2 promoter region highlighted by red box. (E) Union set of peaks (vertical lines) identified in all clusters surrounding TMPRSS2 and 73 peaks that showed co-accessibility with the TMPRSS2 promoter (red lines, co-accessibility score >0.05) via Cicero (Cusanovich et al., 2018). (F) Zoom into TMPRSS2 locus and genome browser tracks of snATAC-seq signal (Robinson et al., 2011). TMPRSS2 promoter region highlighted by red box. For genome browser tracks of BSG, FURIN, CTSL please see Figure 2—figure supplement 1F.

Figure 2—figure supplement 1.. Gene expression and chromatin accessibility for SARS-COV-2 cell entry genes.

(A-E) Cluster-specific violin plots of gene expression of (A) ACE2, (B) TMPRSS2, (C) CTSL, (D) FURIN and (E) BSG. (F) Genome browser (Robinson et al., 2011) representation of snATAC-seq signal at SARS-COV-2 cell entry gene loci CTSL, BSG, FURIN.

Figure 3.. Age-increasing gene expression and accessible chromatin in AT2 cells exhibits signatures of immune regulation and harbors TMPRSS2-linked sites of chromatin accessibility.

(A) Differential analysis was performed on AT2 cells between three ages with replicates (n = 3 per stage). (B) Fraction of AT2 cells with expression of ACE2 (left) and TMPRSS2 (right) in 30wk^GA, 3yo and 30yo human lung samples. All data are represented as mean ± SD. p values derived from unpaired, two-tailed t-tests. For expression data of BSG, CTSL, FURIN please see Figure 3—figure supplement 1A. (C) Log normalized expression of TMPRSS2 in AT2 cells. Displayed are median expression values for AT2 cells in individual samples with at least 1 UMI (unique molecular identifier). (D) Differential analysis was performed on AT2 cells using pairwise comparisons between three ages with replicates (n = 3 per stage). (E) K-means cluster analysis (K = 5) of relative accessibility scores (see Materials and methods) for 22,845 age-dynamic peaks (FDR < 0.05, EdgeR) (Robinson et al., 2010) in AT2 cells. Clusters III and IV show increasing accessibility with age and contain nine TMPRSS2-co-accessible sites. (F) GREAT (McLean et al., 2010) analysis of elements in group cIII (left panel) and cIV (right panel) shows enrichment of immune-related gene ontology terms. (G) Transcription factor motif enrichment analysis of elements in cIII and cIV. (H) K-means cluster analysis (K = 6) of TMPRSS2-co-accessible sites based on the relative percentage of AT2 cells with at least one fragment overlapping each peak. Red bars indicate dynamic peaks identified from pair-wise differential analysis (FDR < 0.05, EdgeR) (Robinson et al., 2010). (I) Locus restricted differential analysis of TMPRSS2-linked peaks with increased accessibility in AT2 with age (top panel in 3H). Data are represented as mean ± SD. Black asterisk, p<0.05 (independent t-test); Red asterisk, FDR < 0.05 (EdgeR) (Robinson et al., 2010) from dynamic peak analysis. For additional sites and promoter accessibility of TMPRSS2 please see Figure 3—figure supplement 1B,C. (J) Genome browser representation of four TMPRSS2-linked peaks across age groups (Robinson et al., 2011).

Figure 3—figure supplement 1.. Gene expression of additional SARS-COV-2 cell entry genes and chrmatin accessibility of peak linked to TMPRSS2 during aging.

(A) Fraction of AT2 cells with expression of BSG, CTSL and FURIN in 30wk^GA,3yo and 30yo human lung samples. All data are represented as mean ± SD. (B) Promoter accessibility of TMPRSS2 in AT2 cells in 30wk^GA, 3yo, and 30yo human lung samples. (C) Fraction of cells with accessibility at peaks co-accessible with TMPRSS2, which did not show increase with age. Data are represented as mean ± SD. Red asterisks denote FDR < 0.05 (EdgeR) and black asterisks denote p<0.05 via independent t-test.

Figure 4.. Genetic variants predicted to affect age-increasing AT2 accessible chromatin are associated with respiratory phenotypes and TMPRSS2 expression.

(A) Top: genome browser view of AT2 sites linked to TMPRSS2 activity and those with age-dependent increase in accessibility. Middle: Number of non-singleton genetic variants in gnomAD v3 mapping (Karczewski et al., 2019) in each age-dependent site predicted to disrupt binding of AT2-enriched TF motifs. Bottom: Common variants (minor allele frequency >0.05 in at least one population) predicted to bind AT2-enriched TF motifs, color-coded by TF family. Motif-bound variants that also have predicted effects (FDR < 0.10) on AT2 accessible chromatin in deltaSVM models (Lee et al., 2015) highlighted in red. (B) Association of common variants with predicted AT2 effects (motif-disrupting+deltaSVM) with human phenotypes in the UK Biobank (Lab, 2020). The majority of tested variants show at least nominal evidence (p<0.005) for association with phenotypes related to respiratory disease, infection and/or medication. (C) Fine-mapping probabilities for an TMPRSS2 expression QTL in human lung samples from the GTEx project release v8 (GTEx Consortium, 2020). The variant rs35074065 has the highest casual probability (PPA = 0.42) for the eQTL, maps in an age-dynamic AT2 site and is predicted to disrupt binding of IRF and STAT TFs. Variants are colored based on r² with rs35074065 in 1000 Genomes Project data using all populations (Auton et al., 2015). (D) Estimated cell type proportions for 515 human lung samples from GTEx derived using cell-type-specific expression profiles for cell types with more than 500 cells from snRNA-seq data generated in this study. (E) Association p-values between rs35074065 genotype and TMPRSS2 lung expression after including an interaction term between genotype and estimated cell-type proportions for each sample. We observed stronger eQTL association when including an interaction with AT2 cell proportion as well as macrophage proportion.

Figure 5.. Fine-mapped risk variants at the 3p21.31 locus associated with respiratory failure in SARS-CoV-2 overlap lung cell-type chromatin sites.

(A) Genome browser view (Robinson et al., 2011) showing posterior probability (PPA) of variants in the fine-mapping credible set at the 3p21.31 association signal and lung cell-type-specific accessible chromatin profiles. Credible set variants that directly overlap lung cell-type chromatin sites are highlighted. Read depth values represent counts per million (CPM). (B) Variant rs17713054 overlaps a site active in AT2 and other epithelial cells and bound by CEBPB among other TFs, and the minor allele A is predicted to bind a CEBP motif. For in-depth analysis of rs76374459 see Figure 5—figure supplement 1. Read depth values represent counts per million (CPM). (C) Co-accessible links between the site harboring rs17713054 (blue box) and other chromatin sites, including the promoter regions of four genes SLC6A20, LIMD1, SACM1L and CCRL2 (gray box). The height of each link represents the strength of co-accessibility (Cusanovich et al., 2018). (D) Expression of SLC6A20 across lung cell types, where each dot represents a nucleus. The highest expression observed was in AT2 cells. (E) Association p-values between rs17713054 genotype and SLC6A20 lung expression after including an interaction term between genotype and estimated cell-type proportions for each sample. We observed strongest eQTL association when including an interaction with AT1/AT2-like and AT2 cell proportion.

Figure 5—figure supplement 1.. Molecular characterization of variant rs76374456.

Variant rs76374456 overlaps a cCRE accessible in erythrocytes which is bound by SPI1 among other TFs, and the minor allele C is predicted to disrupt a SPI1 motif. Read depth values represent counts per million (CPM).