Single-cell eQTL analysis of activated T cell subsets reveals activation and cell type-dependent effects of disease-risk variants

Abstract

The impact of genetic variants on cells challenged in biologically relevant contexts has not been fully explored. Here, we activated CD4⁺ T cells from 89 healthy donors and performed a single-cell RNA sequencing assay with >1 million cells to examine cell type-specific and activation-dependent effects of genetic variants. Single-cell expression quantitative trait loci (sc-eQTL) analysis of 19 distinct CD4⁺ T cell subsets showed that the expression of over 4000 genes is significantly associated with common genetic polymorphisms and that most of these genes show their most prominent effects in specific cell types. These genes included many that encode for molecules important for activation, differentiation, and effector functions of T cells. We also found new gene associations for disease-risk variants identified from genome-wide association studies and highlighted the cell types in which their effects are most prominent. We found that biological sex has a major influence on activation-dependent gene expression in CD4⁺ T cell subsets. Sex-biased transcripts were significantly enriched in several pathways that are essential for the initiation and execution of effector functions by CD4⁺ T cells like TCR signaling, cytokines, cytokine receptors, costimulatory, apoptosis, and cell-cell adhesion pathways. Overall, this DICE (Database of Immune Cell Expression, eQTLs, and Epigenomics) subproject highlights the power of sc-eQTL studies for simultaneously exploring the activation and cell type-dependent effects of common genetic variants on gene expression (https://dice-database.org).

Fig. 1.. Single-cell transcriptomes of activated CD4⁺ T cell subsets.

(A) Study overview. (B) Single-cell transcriptomes of FACS-sorted activated CD4⁺ T cells are displayed by uniform manifold approximation and projection (UMAP) based on Seurat-based clustering. The pie chart shows the fraction of cells in the indicated 19 CD4⁺ T cell subsets. (C) Plot shows average expression (color scale) and percent of expressing cells (dot size scale) for selected marker transcripts in each of the indicated cell clusters.

Fig. 2.. Single-cell eQTL analysis of activated CD4⁺ T cell subsets.

(A) Total number of eGenes identified in each cell cluster by analysis of single-cell RNA-seq data from 735,147 cells (details in Table S1); number of cells analyzed in each cluster is also shown. (B) Fraction of eGenes overlapping with significant SNPs (P_GWAS < 5 × 10⁻⁸) that emerged from the catalogue of GWAS studies (GWAS SNPs) comprising 846 unique human diseases and traits (see Material and Methods). (C) eGenes identified in resting CD4⁺ T cell types (n=8) identified previously (3) or in activated CD4⁺ T cell subsets (n=19) from this study. Venn diagram indicates overlap of eGenes from resting CD4⁺ T cell types with eGenes from activated CD4⁺ T cell subsets (left panel). The XY plot shows the most significant adj. association P value for a given eGene across all resting cell types or activated cell cluster, color coded by eGene category (right panel). (D) eGenes identified in resting CD4⁺ T cell types (n=8) identified previously (3) or in activated CD4⁺ T cell subsets (n=19) from this study. The heatmaps show the most significant adj. association P value (left panels) and absolute effect sizes (right panels) for selected eGenes across all resting or activated CD4⁺ T cell subsets. (E) Expression levels of selected eGenes in the indicated cell cluster from subjects categorized based on the genotype at the indicated peak cis-eQTL; each symbol represents the aggregated mean expression of the indicated gene from an individual subject; adj. association P value *P < 0.05, **P < 0.001, and ***P < 0.00001. (F) Volcano plot of single-cell differential gene expression analysis (DGEA) shows statistical significance (-log₁₀ adjusted P value) and log2 fold change (LFC) comparing the gene expression of cells after knockdown with control siRNA or IRF7 siRNA pools (see Material and Methods). The color scale indicates the average expression (log₂(CPM+1)) and dot size is the fraction of cells expressing a given gene (CPM > 0), both metrics derived from the condition with increased gene expression. The gray dotted lines represent the threshold values for LFC and adjusted P value. (G) UCSC Genome Browser tracks for the IRF7 locus, adj. association P value for cis-eQTLs (showing the most significant peak cis-eQTL across all cell clusters) linked to the expression of IRF7, recombination rate track (96), H3K27ac ChIP-seq track, and H3K27ac HiChIP chromatin interaction maps in naïve CD4⁺ T cells (2). (H) UCSC Genome Browser tracks for the locus harboring the genes IFITM1, IFITM2 and IFITM3, adj. association P value for cis-eQTLs (showing the most significant peak cis-eQTL across all cell clusters) linked to the expression of the indicated transcripts, recombination rate track (96), H3K27ac ChIP-seq track, and H3K27ac HiChIP chromatin interaction maps in naïve CD4⁺ T cells (2). (I) Expression levels of the eGenes P2RX4 and PDCD1 in the indicated cell cluster from subjects categorized based on the genotype at the indicated peak cis-eQTL; each symbol represents the aggregated mean expression of the indicated gene from an individual subject; adj. association P value *P < 0.05, **P < 0.001, and ***P < 0.00001.

Fig. 3.. Activation-dependent effects of genetic variants are highly cell-type restricted.

(A) Cluster-specific eGene analysis, showing row-wise z-scores of the adj. association P values (left panel) and expression levels (right panel) for identified eGenes (one per row) in the indicated cell cluster. (B) Overlap of eGenes identified in the indicated CD4⁺ T cell subsets (n=19). The pie chart shows the fraction of eGenes identified in varying numbers of activated CD4⁺ T cell subsets. (C) eGenes identified in resting CD4⁺ T cell types (n=8) identified previously (3) or in activated CD4⁺ T cell subsets (n=19) from this study. The heatmaps show the most significant adj. P value (left panels) and absolute effect sizes (right panels) for selected eGenes across all resting or activated CD4⁺ T cell subsets. (D-F) Expression levels of selected eGenes in the indicated cell cluster from subjects categorized based on the genotype at the indicated peak cis-eQTL; each symbol represents the aggregated single-cell data (mean expression) from an individual subject; adj. association P value *P < 0.05, **P < 0.001, and ***P < 0.00001. (G) Volcano plot of single-cell differential gene expression analysis (DGEA) shows statistical significance (-log₁₀ adjusted P value) and log2 fold change (LFC) comparing the gene expression of cells after knockdown with control siRNA and LSP1 siRNA pools (see Material and Methods). The color scale indicates the average expression (log₂(CPM+1)) and dot size is the fraction of cells expressing a given gene (CPM > 0), both metrics derived from the condition with increased gene expression. The gray dotted lines represent the threshold values for LFC and adjusted P value. (H) Expression levels of selected eGenes in the indicated cell cluster from subjects categorized based on the genotype at the indicated peak cis-eQTL; each symbol represents the aggregated single-cell data (mean expression) from an individual subject; adj. association P value *P < 0.05, **P < 0.001, and ***P < 0.00001.

Fig. 4.. Disease-risk variants linked to activation-induced eGenes.

(A) Venn diagram indicates overlap of GWAS eGenes from resting CD4⁺ T cell types (n=8) identified previously (3) with GWAS eGenes from activated CD4⁺ T cell subsets (n=19) from this study. (B,C) Expression levels of selected eGenes in the indicated cell cluster from subjects categorized based on the genotype at the indicated peak cis-eQTL; each symbol represents the aggregated single-cell data (mean expression) from an individual subject; adj. association P value *P < 0.05, **P < 0.001, and ***P < 0.00001. (D) Single-tissue TWAS (S-PrediXcan (89)): Gene-level Manhattan plots showing the association P value (-log₁₀) for gene expression with allergic diseases (hay fever, allergic rhinitis or eczema); results across 19 activated CD4⁺ T cell subsets is shown (see Material and Methods). The red horizontal line shows gene-level Bonferroni corrected genome-wide significant P value threshold (P < 2.47 × 10⁻⁶). (E) Single-tissue TWAS (S-PrediXcan (89)): Z scores showing the direction of effect for the genotype-inferred expression of transcripts that encode protein-coding genes in activated CD4⁺ T cell subsets. Red circles indicate genes with Bonferroni corrected genome-wide significant P value threshold (P < 2.47 × 10⁻⁶). (F) Expression levels of PHF5A in the indicated cell cluster from subjects categorized based on the genotype at the indicated peak cis-eQTL; each symbol represents the aggregated single-cell data (mean expression) from an individual subject; adj. association P value *P < 0.05, **P < 0.001, and ***P < 0.00001. (G) Volcano plots of single-cell differential gene expression analysis (DGEA) show statistical significance (-log₁₀ adjusted P value) and log2 fold change (LFC) comparing the gene expression of cells after knockdown with control siRNA and PHF5A siRNA pools (see Material and Methods). The color scale indicates the average expression (log₂(CPM+1)) and dot size is the fraction of cells expressing a given gene (CPM > 0), both metrics derived from the condition with increased gene expression. The gray dotted lines represent the threshold values for LFC and adjusted P value. (H,I) Expression levels of selected eGenes in the indicated cell subsets (left panel, bulk gene expression data of resting cells as previously published (3); expression levels are shown in TPM, transcripts per million) or activated CD4⁺ T cell subsets (right panel) from subjects categorized based on the genotype at the indicated peak cis-eQTL; each symbol represents the aggregated single-cell data (mean expression) from an individual subject; adj. association P value *P < 0.05, **P < 0.001, and ***P < 0.00001. (J) UCSC Genome Browser tracks for the TNFRSF14 locus, adj. association P value for cis-eQTLs (showing the most significant peak cis-eQTL across all cell cluster) linked to the expression of TNFRSF14, recombination rate track (96), H3K27ac ChIP-seq track, and H3K27ac HiChIP chromatin interaction maps in naïve CD4⁺ T cells (2). (F) Expression levels of selected eGenes in the indicated cell subsets (left panel, bulk gene expression data of resting cells as previously published (3); expression levels are shown in TPM, transcripts per million) or activated CD4⁺ T cell subsets (right panel) from subjects categorized based on the genotype at the indicated peak cis-eQTL; each symbol represents the aggregated single-cell data (mean expression) from an individual subject; adj. association P value *P < 0.05, **P < 0.001, and ***P < 0.00001.

Fig. 5.. Sex has a major influence on activation-induced gene expression patterns.

(A) Analysis of differentially expressed transcripts (sex-biased transcripts, each dot) in all activated CD4⁺ T cell subsets from females (n=35) versus male (n=54) subjects (see Material and Methods). Right panel, fractions of sex-biased transcripts in autosomes and sex chromosomes. (B) Overlap of sex-biased transcripts. The pie chart shows the relative cell specificity of the transcripts in the indicated CD4⁺ T cell subsets (n=19). (C) Venn diagram indicates overlap of sex-biased transcripts from resting CD4⁺ T cell types (n=8) reported previously (3) with sex-biased transcripts from activated CD4⁺ T cell subsets (n=19) in this study. (D) ConsensusPathDB (CPBD) gene term enrichment analysis showing the biological pathways enriched for sex-biased transcripts. (E) DAVID pathway analysis showing the biological pathways enriched for sex-biased transcripts in activated CD4⁺ T cell subsets.