Genetic and non-genetic factors affecting the expression of COVID-19-relevant genes in the large airway epithelium

Abstract

Background: The large airway epithelial barrier provides one of the first lines of defense against respiratory viruses, including SARS-CoV-2 that causes COVID-19. Substantial inter-individual variability in individual disease courses is hypothesized to be partially mediated by the differential regulation of the genes that interact with the SARS-CoV-2 virus or are involved in the subsequent host response. Here, we comprehensively investigated non-genetic and genetic factors influencing COVID-19-relevant bronchial epithelial gene expression.

Methods: We analyzed RNA-sequencing data from bronchial epithelial brushings obtained from uninfected individuals. We related ACE2 gene expression to host and environmental factors in the SPIROMICS cohort of smokers with and without chronic obstructive pulmonary disease (COPD) and replicated these associations in two asthma cohorts, SARP and MAST. To identify airway biology beyond ACE2 binding that may contribute to increased susceptibility, we used gene set enrichment analyses to determine if gene expression changes indicative of a suppressed airway immune response observed early in SARS-CoV-2 infection are also observed in association with host factors. To identify host genetic variants affecting COVID-19 susceptibility in SPIROMICS, we performed expression quantitative trait (eQTL) mapping and investigated the phenotypic associations of the eQTL variants.

Results: We found that ACE2 expression was higher in relation to active smoking, obesity, and hypertension that are known risk factors of COVID-19 severity, while an association with interferon-related inflammation was driven by the truncated, non-binding ACE2 isoform. We discovered that expression patterns of a suppressed airway immune response to early SARS-CoV-2 infection, compared to other viruses, are similar to patterns associated with obesity, hypertension, and cardiovascular disease, which may thus contribute to a COVID-19-susceptible airway environment. eQTL mapping identified regulatory variants for genes implicated in COVID-19, some of which had pheWAS evidence for their potential role in respiratory infections.

Conclusions: These data provide evidence that clinically relevant variation in the expression of COVID-19-related genes is associated with host factors, environmental exposures, and likely host genetic variation.

Keywords: ACE2; Bronchial epithelium; COVID-19; SARS-CoV-2; eQTL.

Fig. 1

Study design. Graphical illustration of analyses (gray boxes) carried out to study non-genetic and genetic factors affecting the expression of COVID-19-related genes in bronchial epithelium. Input data sets for these analyses are denoted with a green box (WGS and RNA-seq) and external data sets or data resources used in these analyses are denoted with a blue box

Fig. 2

ACE2 gene expression associations in SPIROMICS. a–d Box plots showing that ACE2 log₂ gene expression (x-axis) was increased in association with current but not former smoking as compared to never smokers (a), obesity (b, validated in the MAST and SARP cohorts, Additional file 3: Figure S2a-b), hypertension (c, adjustments include anti-hypertensive treatment, validated in SARP, Additional file 3: Figure S3a, data not collected in MAST), and female sex (d, not replicated in either MAST or SARP, Additional file 2: Table S1A). e Scatterplots showing that ACE2 gene expression was increased in association with higher levels of our previously validated gene signatures of the airway epithelial response to interferon (left panel, replicated in SARP) and to IL-17 inflammation (right panel, replicated in MAST and SARP) after adjusting for smoking status (Additional file 2: Table S1B). f Box plots showing that ACE2 Exon 1c, which contributes to the truncated ACE2 transcript was differentially increased in association with our interferon signature while Exons 1a and 1b that contribute to the full length ACE2 transcript were not. P values indicated by: **** < 0.0001, *** < 0.001, ** < 0.01, * < 0.05, ns = not significant in linear models adjusted for covariates. In a–d and f, the boxes denote the interquartile range, the center line denotes the median, and whiskers denote the interquartile range × 1.5

Fig. 3

COVID-19-related gene set enrichment analyses in association with comorbidities. a–f Barcode plots in which the vertical lines represent the 100 genes most upregulated (red) or downregulated (blue) in nasal/oropharyngeal swab samples obtained from COVID-19 patients as compared to other viruses at the time of diagnosis of an acute upper respiratory infection. These gene sets are plotted against log fold gene expression changes arranged from most downregulated to most upregulated with that comorbidity (horizontal gray bar). Lines above (red) and below (blue) the bar represent the running sum statistic with a significant finding indicated when the line crosses the dashed line at either end of the plot. Genes downregulated by SARS-CoV-2 infection compared to other viruses were significantly enriched amongst genes downregulated in association with cardiovascular conditions overall (a), hypertension (b), and obesity (c), while in current (d) and former smoking (f) and in COPD (e), these downregulated genes in COVID-19 were enriched amongst upregulated genes in association with comorbidity. ** indicates FDR < 0.05. g COVID-19-related pathway gene sets were generated from an IPA analysis of the genes downregulated by SARS-CoV-2 infection compared to other viruses. Gene set enrichment scores for gene sets enriched at FDR < 0.05 (columns) are shown in the heatmap plotted against comorbidities (rows) with gene sets enriched amongst downregulated and upregulated genes indicated in blue and yellow, respectively. All pathways not enriched at FDR < 0.05 were shrunk to zero (white). Euclidean distance with average linkage was used for clustering

Fig. 4

Cis-eQTLs in bronchial epithelium. a Effect size measured as allelic fold change (aFC, log₂) of the significant cis-eQTLs for COVID-19 candidate genes. Error bars denote 95% bootstrap confidence intervals. b Comparison of the regulatory effects and the effect of SARS-CoV-2 infection on the transcription of COVID-19 candidate genes in normal bronchial epithelial cells from Blanco-Melo et al. [30]. The graph shows regulatory effects as aFC as in a and fold change (log₂) of differential expression comparing the infected with mock-treated cells with error bars denoting the 95% confidence interval. Genes with adjusted P value < 0.05 in the differential expression analysis are colored in black, genes with non-significant effect are colored in gray. Highlighted genes have eQTL effect size greater than 50% of the differential expression effect size on the absolute scale. DE—differential expression. c Replication of cis-eQTLs from bronchial epithelium in GTEx v8 using the concordance rate (proportion of gene-variant pairs with the same direction of the effect, left panel) and proportion of true positives (π₁, right panel). Upper panel shows the effect of sample size on the replication and concordance measures quantified as Spearman correlation coefficient (ρ). Lower panel shows the replication and concordance measures as the function of epithelial cell enrichment of the tissues measured as median epithelial cell enrichment score from xCell. Gray dashed line denotes median enrichment score > 0.1, which classifies tissues as enriched for epithelial cells. Wilcoxon rank sum test was used to estimate the difference in replication estimates between tissues enriched or not enriched for epithelial cells. The 16 tissues enriched for epithelial cells are outlined in the figure legend, for the full legend see Additional file 3: Figure S9a

Fig. 5

Colocalization analysis of the regulatory variants for COVID-19-related genes. a Illustration of the concept of how regulatory variants for COVID-19-related genes in bronchial epithelium can be possible candidates for genetic factors that affect infection or progression of the disease. Dotted lines denote the hypothesis we are able to create by searching for the phenotypic associations of the cis-eQTLs for COVID-19-related genes. b Heatmap of the colocalization analysis results for 20 COVID-19-related genes with eQTLs that have at least one phenotypic association belonging to the experimental factor ontology (EFO) parent categories relevant to COVID-19 (respiratory disease, hematological or pulmonary function measurement). Genes highlighted in bold indicate the loci involving COVID-19-relevant EFO categories with posterior probability for colocalization (PP4) > 0.5, suggesting evidence for shared genetic causality between eQTL and GWAS trait. In the TLE locus, the nearest genome-wide significant variant for forced expiratory volume in 1 s (FEV1) from Shrine et al. [57] is more than 1 Mb away, indicating that the association between the variant and FEV1 might be confounded by incomplete adjustment for height. c–e Regional association plot for the GWAS signal on the upper panel and cis-eQTL signal on the lower panel for IFITM3 (c), ERMP1 (d), and MEPCE (e) locus, where the eQTL for the corresponding gene colocalizes with the GWAS trait relevant to COVID-19. Genomic position of the variants is shown on the x-axis and -log10(P value) of the GWAS or eQTL association on the y-axis. The lead GWAS and eQTL variants are highlighted