A preferential ARE-mRNA upregulation in cell lines infected with SARS-CoV-2 and other respiratory viruses
Several recent studies used RNA-seq-based global transcriptomic analysis to investigate differential gene expression resulting from SARS-CoV-2 and other viral infections. In such a study, the human lung cancer cell line A549 was used as a model to assess the global effect on gene regulation of the following viruses: SARS-CoV-2, respiratory syncytial virus (RSV), human parainfluenza virus type 3 (HPIV3) or influenza A virus (IAV). In the same study, another lung cancer cell line Calu-3 was also infected with SARS-CoV-2 [18]. The authors concluded that compared to other viral infections, host response to SARS-CoV-2 is deficient in interferon response, but is excessive in inflammatory cytokine and chemokine production [18]. Here, we determined the fraction of ARE-containing mRNAs within the groups of infection-dependent upregulated or downregulated genes by crossing the list of differentially expressed genes (DEGs) with the AU-rich element (ARED-plus) database [12]. The statistical significance for differential expression was set to adjusted P value < 0.05. Interestingly, in all cases where cell lines were used, the fraction of ARE-containing mRNAs within the overexpressed genes was higher than the ARE-fraction in the downregulated mRNAs with very high statistical significance (Fig. 1A; Additional file 1: Table S1: Upregulation). For instance, 37% of the upregulated mRNAs in SARS-CoV-2-infected A549 cells contained AREs and only 16% of the downregulated. Similar results were observed in A549 expressing the SARS-CoV-2 receptor protein ACE2 with high multiplicity of SARS-CoV-2 infection (MOI) (37 vs. 24%) or low MOI (40 vs. 24%). A similar observation was also made in Calu-3 cells infected with SARS-CoV-2 (42 vs. 21%). Other viruses transfecting A549 cells like RSV or HPIV3 or IAV had similar but more pronounced effects. Most notably was RSV (43 and 13%). The median of the Log2 Fold Change between transfected versus un-transfected was always clearly and significantly higher in ARE-mRNAs compared to the non-ARE mRNAs. Due to the large sample sizes of ARE and non-ARE mRNAs, the results were always very highly significant with very low P values typically less than 10–17 and consequently with a narrow range of 95% CI (data not shown) and low margin of error.
Importantly, the difference between the median of Log2 Fold Change of ARE and that of the non-ARE was higher in A549 infected with RSV, HPIV3, IAV compared to the same cell line and Calu-3 infected with the SARS-CoV2, suggesting that SARS-CoV-2 drives a weaker ARE-dependent response than the other viruses (Fig. 1B, C).
Volcano plots of ARE-mRNAs in SARS-CoV-2-infected A549 and Calu-3 and RSV-infected A549 are displayed and show that all of the 20 most significantly differentially expressed ARE-mRNAs are overexpressed and none is under-expressed (Fig. 1D). Similar observations were made in all other infections (data not shown). Among the recurring, most significantly upregulated ARE-mRNAs are typically cytokines and chemokines like CXCL1, CXCL2, TNF and IL6 but also regulatory genes like TNFAIP3, NFKBIA, DUSPs and JUN (Fig. 1D).
Similarly, the transcriptomic effects of SARS-CoV-1 and MERS-CoV infection of the human fetal lung fibroblast cells MRC5 were analyzed (GSE56192) and [18]. A notable 51% of overexpressed mRNAs in the SARS-CoV-1-infected MRC5 cells contained AREs, whereas only 13% of the downregulated (Fig. 2A). A significant but less pronounced effect on MRC5 was observed in MERS-CoV (33 and 20%) (Fig. 2B). Volcano plots show that all or most of the 20 most significantly differentially expressed mRNAs in MERS-CoV-1 and SARS-CoV-1 are ARE-containing (Fig. 2C). Taken together, these in vitro observations clearly state that cell lines infected with respiratory viruses tend to overexpress ARE-containing mRNA more than non-ARE, and this observation is weaker with SARS-CoV-2 compared to other viruses.
ARE-mRNA levels in SARS-CoV-2-infected clinical samples
Similar to the investigations performed on cell lines above, a possible higher upregulation of the ARE-mRNA containing fraction from clinical samples was investigated. First, we analyzed transcriptomic data originating from post-mortem lung samples [18]. SARS-CoV-2 infection did not have a significant enrichment in the ARE-gene fraction within the overexpressed group of mRNAs and subsequently no higher global levels of ARE-mRNAs could be detected (Fig. 3A, B). Xiong et al. [19] investigated the differential expression of mRNAs in bronchoalveolar lavage fluid (BLAF) and PBMCs of SARS-CoV2-infected patients. Our analysis shows no significant enrichment of the ARE-mRNA fraction in the upregulated group; in fact, a tendency to the opposite, especially in BALF (16 vs. 24%), could be observed; however, this did not result in a significantly higher global upregulation of ARE-mRNAs (Fig. 3A, B). Analysis of transcriptomic studies on neutrophils, epithelial cells and macrophages from infected individuals compared to healthy donors resulted in a comparatively weak preferential enrichment and upregulation of ARE-mRNAs in neutrophils only. This, however, did not result in a global preferential upregulation of ARE-mRNAs [20]. In the same study, the authors investigated the effect of IVA infection in the same type of cells. The enrichment of the ARE-signature within the overexpressed genes was observed and is significant only in neutrophils infected with IAV (35 vs. 26%), and this resulted in a comparatively weak enrichment and upregulation of ARE-mRNAs (Fig. 3A, B). Sarma et al. [22] performed transcriptomic studies on tracheal aspirate from COVID-19 patients with acute respiratory distress syndrome (ARDS) with uninfected controls without ARDS. We found a significant enrichment of the ARE-mRNAs within the overexpressed genes (36 vs. 22%) along with a higher upregulation of ARE-mRNAs compared with those that do not contain AREs (P = 0.0004) (Fig. 4A, B).
Overall, the results of the investigation of clinical samples may suggest that the higher upregulation of ARE-mRNAs observed in SARS-CoV-2-infected cell lines is non-existent or weak in clinical samples compared to cell lines. However, it is important to note that the Blanco Mello et al. and Xiong et al. transcriptomic data were based on only two and three SARS-CoV-2 patients and healthy donors and are likely to be lacking the statistical power to identify transcripts with moderate but significant differential expression. Similarity, the Gao et al. study was based on five SARS-CoV-2 and IAV samples. The transcriptomic data from Sarma et al. however, was based on 15 COVID-19 patients and 5 controls and resulted in a significant ARE-mRNA upregulation signature.
Preferential ARE-mRNA upregulation in blood cells of patients with mild forms of COVID-19
Bibert et al. [21] conducted a large comparative transcriptomic profile study that is based on blood cells from 103 patients with different severity levels of COVID-19 along with 27 healthy and 22 influenza-infected individuals. The authors classified COVID-19 patients into four groups according to the level of respiratory failure; 23 with no oxygen support requirements (OXY0), 40 who received oxygen but no need for mechanical ventilation (OXY1), 15 who required mechanical ventilation within the first 7 days in hospital (TUBE early) while 25 were sampled after 7 days (TUBE late). Log2 Fold Change for influenza infection and each of the four severity groups was determined using DESeq2 at a transcriptomic level compared to the healthy donors. The influenza A and B virus-infected patients displayed the typical ARE-over-representation and expression signature found in cell lines (Fig. 4A). Most interesting was the OXY0 group of SARS-CoV-2-infected individuals: a strong and highly significant enrichment of ARE-mRNAs within the upregulated mRNAs of infected patients compared to the downregulated group could be observed (36 vs. 17%). This correlated with a clear and significant higher upregulation of ARE-mRNAs compared to non-ARE, median of Log2 Fold Change of ARE-mRNAs 0.5 versus 0.241 for non-ARE. (Fig. 4B, C). With the large sample sizes of ARE and non-ARE mRNAs, the result was very highly significant (P value = 9.26 × 10–30) and the 95% CIs were narrow (0.48–0.52 for ARE-mRNAs and 0.185–0.28 for non-ARE). Strikingly, this type of signature became reduced with increasing severity of the disease (Fig. 4B, C). For instance, for OXY1, the enrichment was 35% versus 19% and the medians 0.38 versus 0.27, while for TUBE early, the enrichment was comparatively weakly significant, 27 versus 24% (P = 0.0004), and the medians, 0.53 versus 0.44. For TUBE late, the enrichment was 26 versus 20%, and the difference in the Log2 Fold Change between ARE and non-ARE mRNAs was 0.4058 versus 0.3676 and no longer significant (P = 0.0622) (Fig. 4B, C). Similar to results from cells line, a volcano plot from OXY0 and OXY1 patients shows that most of the 20 most significantly differentially expressed ARE-mRNAs are overexpressed (Fig. 4D).
ARE-mRNAs that are upregulated in the mild form of the disease
Since patients with mild disease have a higher overall signature of upregulated ARE-mRNAs, a Venn diagram was plotted to find ARE-mRNAs that are unique to patients with mild disease. For comparison, another diagram was plotted for non-ARE mRNAs (Fig. 5A). The largest proportion of upregulated ARE and non-ARE mRNAs was common to OXY0, OXY1, TUBE early and TUBE late, 41.3 and 34.1%, respectively. This suggests the existence of common regulatory networks that are activated in all levels of case severity of SARS-CoV-2. Interestingly, however, is a relative large proportion of upregulated ARE-mRNAs that are unique to the mildest form of COVID-19, OXY0, 12.6% (258 mRNAs) versus 5.4% (281 mRNAs) for non-ARE (Fig. 5A). Also interestingly and in contrast, the fractions of upregulated non-ARE mRNAs that are unique to TUBE early and TUBE late are twice the fractions in ARE-mRNAs. Precisely, 6.2 + 8.3 + 10.6 = 25.1% for non-ARE and 3.3 + 3.8 + 5.8 = 12.9% for ARE. These observations show that some upregulated ARE-mRNAs tend to be unique to the mild forms of the disease while some upregulated non-ARE mRNAs tend to be unique to the severe forms of the disease.
Gene ontology analysis of ARE-mRNAs that are differentially up-regulated in mild COVID-19
The DESeq2 software was used to find mRNAs that are upregulated in the mildest form of the disease (OXY0) relative to patients in need of intubation (TUBE early and late). The median of the levels of ARE-mRNA in OXY0 relative to TUBE was higher than that of non-ARE mRNAs confirming a better ARE-response in the mild form of the disease (Fig. 5B). Next, the most significantly upregulated ARE-mRNAs in OXY0/TUBE (padj < 0.0001, n = 56) were subjected to gene ontology clustering in STRING portal (string-db.org). 27 (48%) turned out to be related to defense against viral infections and immune responses under the biological processes category (Additional file 2: Table S2). However, and similarly, 66 of 156 (42%) upregulated non-ARE mRNAs also clustered in same and similar categories that relate to defense against viral infections (Additional file 3: Table S3). The above observations indicate that the significant ARE-response is not alone responsible for the stronger interferon and immune response found in the mildest form of the disease.
Next, the 258 unique ARE-mRNAs that were upregulated in OXY0 (Fig. 5A) were analyzed using the STRING portal (string-db.org/). In the biological processes category, the highest strength of enrichment was weak, 0.77 in nucleotide excision repair (nine proteins of 111 of network, FDR 0.0159). All other enrichments in categories like molecular function, biological process, KEGG pathway and cellular component resulted in strengths of less than 0.5 indicating that the unique genes do not particularly cluster in common functional categories. Next, we focused on the 849 mRNAs that are upregulated in all four severity levels of COVID-19. Specifically, we searched for a possible cluster of ARE-mRNAs that has a higher upregulation level in the mild forms of the disease compared to the severe forms. The Log2 Fold Change values of all four forms were subjected to a K-means clustering in JMP resulting in four cluster. Cluster (3) consisted of 233 genes that are upregulated in OXY0 and OXY1 compared to TUBE early and TUBE late (Fig. 5C). The 233 genes clustered into three clusters using K-means in the string portal. Cluster 1 consisted of 81 genes and showed the highest enrichment (p value < 1.0e–16). PTEN signaling was central in this cluster and is in agreement with a previous observation [22] (Fig. 5D). KEGG pathway enrichments of the 81 genes were ranked according to strength and start with mTOR signaling but also contain other pathways including the anti-inflammatory PD-L1/PD-1 checkpoint pathway, viral infection pathways like human cytomegalovirus infection, human papillomavirus infection and chemokine signaling pathway (Additional file 4: Table S4).