Single cell gene expression profiling of nasal ciliated cells reveals distinctive biological processes related to epigenetic mechanisms in patients with severe COVID-19

Abstract

Objective: To explore the molecular processes associated with cellular regulatory programs in patients with COVID-19, including gene activation or repression mediated by epigenetic mechanisms. We hypothesized that a comprehensive gene expression profiling of nasopharyngeal epithelial cells might expand our understanding of the pathogenic mechanisms of severe COVID-19.

Methods: We used single-cell RNA sequencing (scRNAseq) profiling of ciliated cells (n = 12,725) from healthy controls (SARS-CoV-2 negative n = 13) and patients with mild/moderate (n = 13) and severe (n = 14) COVID-19. ScRNAseq data at the patient level were used to perform gene set and pathway enrichment analyses. We prioritized candidate miRNA-target interactions and epigenetic mechanisms.

Results: We found that mild/moderate COVID-19 compared to healthy controls had upregulation of gene expression signatures associated with mitochondrial function, misfolded proteins, and membrane permeability. In addition, we found that compared to mild/moderate disease, severe COVID-19 had downregulation of epigenetic mechanisms, including DNA and histone H3K4 methylation and chromatin remodelling regulation. Furthermore, we found 11-ranked miRNAs that may explain miRNA-dependent regulation of histone methylation, some of which share seed sequences with SARS-CoV-2 miRNAs.

Conclusion: Our results may provide novel insights into the epigenetic mechanisms mediating the clinical course of SARS-CoV-2 infection.

Keywords: COVID-19; Histones; Methylation; Mitochondria; SARS-CoV-2; Single-cell transcriptomics.

Graphical abstract

Fig. 1

Methodology flow chart Step 1: The specific type of cells is selected from each patient (ciliated cells in our case). Step 2: Normalized counts for each gene are averaged (arithmetic mean) over the selected cells population of each patient, obtaining the average expression profile of ciliated cells for each patient. Step 3: Signal2Noise metric is computed among the patients of the two conditions of interest. μA(μB) and σA(σB) are the mean and standard deviation of gene expression levels computed in the previous step for patients under condition A (B). Then, the genes are ranked by decreasing order of the metric value. Step 4: Finally, it is determined if the ranked gene set, associated with a given biological process (BP), shows statistically significant, concordant differences between the two conditions. This last step can be performed with different BPs or even with gene set associated with molecular functions or miRNAs.

Fig. 2

Gene set enrichment analysis derived from ciliated cells in mild/moderate COVID-19 patients vs. healthy controls The Figure illustrates a network of GO (Gene Ontology) terms corresponding to biological processes (BPs) significantly enriched in genes upregulated (blue nodes) and downregulated (red nodes) in mild-moderate COVID-19 (COVID-M) patients concerning to healthy control (HC) individuals (Bonferroni adjusted p-value, FDR 0.05). The node's size is proportional to the number of genes associated with that BP. The edges represent gene overlap between gene sets related to different GO terms. Connected nodes are organized in clusters of interconnected BPs obtained by the MCL algorithm, which considers similarity among gene sets to assign the edges, with a similarity score threshold of 0.5.

Fig. 3

Gene set enrichment analysis derived from ciliated cells in severe vs. mild/moderate COVID-19 patients The Figure illustrates enrichment analysis of biological processes (BPs) upregulated (blue nodes) and downregulated (red nodes) in patients with severe COVID-19 vs. patients with mild/moderate disease(Bonferroni adjusted p-value, FDR 0.05). BPs were clustered by using the MCL algorithm with a similarity score threshold of 0.5.

Fig. 4

Heat map of expression genes belonging to the pathways associated with mitochondrial membrane permeabilization, protein targeting to mitochondria, and lysosome and protein targeting the membrane. Blue squares correspond to downregulated genes and red squares to upregulated genes. We selected only the subset of genes most upregulated (the percentile 0.1 of the ranked genes)in mild/moderate COVID-19 vs. healthy patients. The colour scale represents the RNA abundance relative to the media of transcript levels in all subjects and then log2 transformed.

Fig. 5

The role of epigenetic mechanisms in COVID-19: a regulatory network of miRNAs-target genes A. Rank of genes according to the Signal2Noise metric as implemented in the GSEA platform (COVID-19 S vs. COVID-19 M) using the media of the normalized counts of the whole cell population in each individual. The top 10% of downregulated genes (rank over the dashed line) was used in the subsequent analysis shown in panels B and C. B. Venn Diagram showing the overlap between the 3076 gene targets of the 11 miRNAs associated with the target sets enriched in COVID-19 S vs. COVID-19 M and the 1167 top-ranked genes downregulated in COVID-19 S (as explained in panel A). The intersection shows the 336 genes used in panel C. C. A heat map of the 336 genes that result from the intersection mentioned in panel B. Squares in blue and red represent genes with lower and higher RNA abundance relative to the media of transcript level of all subjects and then log2 transformed, respectively. COVID-19 S: severe disease; COVID-19 M: mild/moderate disease.

Fig. 6

Network of the 336 genes downregulated in COVID-S and their associated GO terms for biological processes (BP) The network was constructed using the ToppGene/ToppCluster resource (https://toppgene.cchmc.org/) with an FDR of 0.05. Hexagons in red and squares in light blue stand for genes and BPs, respectively, which were distributed using the edge-weighted Spring embedded layout and slightly modified for readability using Cytoscape v3.4.0. In general, the upper left region shows BP associated with DNA and histone modifications, regulation of transcription, and chromatin or chromosome remodelling. In contrast, the right middle part shows BPs associated with organ morphogenesis, particularly cardiac and vascular. The lower left region shows metabolic processes, primarily catabolic.

figs1

The regulatory network of miRNAs associated with severe COVID-19: correlation analysis The comparative analysis of GSEA between cells derived from mild and severe COVID-19 shows that gene expression enriched in the group of mild patients (FDR <0.05) are the targets of 29 miRNAs. This Figure shows the fraction of target genes in common among the associated target genes of each of the 29 miRNAs. The miRNAs associated with each target gene set show high redundancy (green clusters). Therefore, there are only 11 target gene sets with low redundancy among themselves. We selected the corresponding 11 miRNAs listed to the right for the subsequent analyses

figs2

miRNA-target pairs involved in epigenetic mechanisms The heat plot shows the percentage of overlap (green scale) between the target genes associated with each miRNA, with the genes involved in each downregulated BP in patients with severe COVID-19 shown in Fig. 2 (FDR <0.05) and listed in Supplementary Table 2. We only took those genes within the 10% most enriched genes in mild patients. It is observed that most of the genes linked to histone methylation processes (highlighted in blue font) and that were downregulated in patients with severe COVID-19 are also target genes of miRNAs, which suggests that the presence of these miRNAs (or their viral homologs) deregulates histone methylation in these patients