Comparative transcriptomic analyzes of human lung epithelial cells infected with wild-type SARS-CoV-2 and its variant with a 12-bp missing in the E gene

Abstract

The severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) is a novel coronavirus that caused a global outbreak of coronavirus disease 2019 (COVID-19) pandemic. To elucidate the mechanism of SARS-CoV-2 replication and immunogenicity, we performed a comparative transcriptome profile of mRNA and long non-coding RNAs (lncRNAs) in human lung epithelial cells infected with the SARS-CoV-2 wild-type strain (8X) and the variant with a 12-bp deletion in the E gene (F8). In total, 3,966 differentially expressed genes (DEGs) and 110 differentially expressed lncRNA (DE-lncRNA) candidates were identified. Of these, 94 DEGs and 32 DE-lncRNAs were found between samples infected with F8 and 8X. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyzes revealed that pathways such as the TNF signaling pathway and viral protein interaction with cytokine and cytokine receptor were involved. Furthermore, we constructed a lncRNA-protein-coding gene co-expression interaction network. The KEGG analysis of the co-expressed genes showed that these differentially expressed lncRNAs were enriched in pathways related to the immune response, which might explain the different replication and immunogenicity properties of the 8X and F8 strains. These results provide a useful resource for studying the pathogenesis of SARS-CoV-2 variants.

Keywords: COVID-19; E gene; RNA-seq; immune response; lncRNAs.

Figure 1

RNA categories of quantified transcripts in RNA-seq. The percentages of each RNA type are shown in F8/8X infected or mock-infected Calu-3 cells.

Figure 2

Volcano plots showing the DEGs between 8X and control (A), F8 and control (B), and the 8X and F8 (C). (D) Venn diagram of the identified DEGs shared among the above three groups. (E) Heatmap of the top 10 DEGs between 8X and F8 groups.

Figure 3

Dot plots displaying enriched GO terms (A) or KEGG pathways (B) of identified DEGs. Dot size represents the count of the enriched gene within each category, and the x-axis represents the gene ratio (rich factor). The GO terms/KEGG pathways are arranged on the basis of the value of the gene ratio, not on their adjusted p value, for easier visualization purposes.

Figure 4

Comparison of the expression level (A) and exon number (B) of lncRNAs and protein-coding genes.

Figure 5

Functional enrichment analysis of identified lncRNAs between 8X and F8. Representative over-represented KEGG pathway (top) and GO terms (bottom) of gene-expression clusters. BP, biological process; MF, molecular function; CC, cellular component.

Figure 6

Co-expression analysis of lncRNAs and protein-coding genes. (A) The top 10 over-represented KEGG pathways of co-expressed genes. The network of lncRNAs with (B) viral protein interaction with cytokine and cytokine receptor, and (C) TNF signaling pathway-related genes (p < 0.01). The solid line represents a positive correlation, while the dotted line represents a negative correlation.

Figure 7

The competing endogenous RNA (ceRNA) network. mRNAs, circRNAs, host lncRNAs and miRNAs in the networks are represented as squares, circles, triangles and arrows, respectively.