The m7G Modification Level and Immune Infiltration Characteristics in Patients with COVID-19

Introduction

Severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) has caused widespread health concerns worldwide. The clinical symptoms of individuals with coronavirus 2019 (COVID-19) can range from asymptomatic to pneumonia, severe hypoxia, or even organ failure.¹ COVID-19 patients who are elderly or exhibit obesity, hypertension, or diabetes may have higher mortality and hospitalization rates.² However, the specific mechanism underlying disease risk disparities and individual outcomes is still unclear. It is imperative to elucidate the underlying pathophysiology and identify therapeutic approaches to treat COVID-19.

Over 170 types of RNA modifications play crucial roles in regulating mRNA translation, stability, transport, and RNA processing. In eukaryotic cells, the N6-methyladenosine (m6A) modification is the most prevalent among the messenger RNAs (mRNAs).^3–5 N7-methylguanosine (m7G) is another ubiquitous RNA modification that was first reported in the 1970s located at the 5’ caps of eukaryotic mRNA or internally within rRNA and tRNA of all species,^5–9 and it affects RNA processing, metabolism, and function.¹⁰ Furthermore, in addition to its presence at the cap position in mRNAs, m7G has been found in internal mRNA regions.^9,11,12 Increasing evidence suggests that m7G plays an important role in the development of human diseases.¹³

Modifications to viral RNA also play a critical role in the viral life cycle. m6A, 5-methylcytosine (m5C), and N4-acetylcytidine (ac4C) have been reported to be involved in viral infection regulation.^14–17 Owing to the multiple functions related to the m7G cap structure in cellular mRNAs, several viruses modify the 5’ end of their viral RNA in the same way. During the infectious cycle, most viruses use the cellular nuclear capping machinery to cap RNA. However, other viruses fetch caps from other pathways and cellular mRNAs or are encoded by viral capping enzymes.^18–20 Furthermore, viral infection also regulates the expression of host gens for m7G RNA modification.²¹

Many studies have demonstrated the impact of SARS-CoV-2 on host cell m6A methylome.^22–24 However, few studies have revealed alterations in m7G in SARS-CoV-2-infected cells. In the present study, we assessed the expression of 29 major regulators of m7G RNA modification in SARS-CoV-2-infected patients. Additionally, according to the m7G-related gene expression, the individuals were divided into two clusters by consensus clustering with different prognoses. Furthermore, a nomogram model based on the six m7G genes was developed for individual risk prediction.

Material and Methods

The RNA-Seq of COVID-19 Patients and m7G-Related Gene

The GSE157103 dataset²⁵ containing transcript profiles of COVID-19 patients was downloaded from GEO.²⁶ The GSE157103 dataset included 100 COVID-19 patients and 26 non-COVID-19 controls. Clinical information included gender, age, Charlson comorbidity index score, ventilator-free days, mechanical ventilation, and intensive care unit (ICU) status. The m7G-related genes were identified by a literature search and GSEA (http://www.gsea-msigdb.org/gsea/login.jsp) database, and a total of 29 m7G-related genes were obtained. R software was used for all data processing and analyses in this study. Statistical significance was defined as p<0.05.

Classification of Disease Subtypes

Current research mainly involves three types of classification for disease subtypes, m7G clusters, gene clusters, and m7G-score clusters. Both m7G clusters and gene clusters select the optimal consensus matrix to determine the optimal number of clusters by partitioning around the medoids algorithm. The R package “limma” was used to perform differential expression analysis. The R package “ConsensusClusterPlus” was used to group the patients with COVID-19. The screening condition for differentially expressed genes (DEGs) was set to an adjusted p-value of < 0.05. First, 15 m7G-related differentially expressed genes (m7G-DEGs) were found between COVID-19 patients and individuals not infected with COVID-19. The classification of COVID-19 patients into m7G-clusters is based on the expression level of these 15 m7G-related genes, and patients were assigned to “m7G cluster A” and “m7G cluster B”. Next, we analyzed all differentially expressed genes (DEGs) between m7G cluster A and m7G cluster B, including but not limited to m7G-related genes. The classification of gene clusters is based on the expression level of these DEGs which were detected between m7G-clusters, and then the COVID-19 patients are assigned to “gene cluster A” and “gene cluster B”. In addition, we constructed an m7G score. This score index is based on the principal component analysis (PCA) algorithm to analyze the 15 m7G DEGs. The m7G score combines the PC1 and PC2 scores to form a signature score. According to the median, the patients were divided into “m7G-score high” and “m7G-score low” groups.

The Abundance Estimates for Immune Cells

ssGSEA analysis was used to evaluate the abundance estimates for the 23 immune cell types. This was done to infer the individual immune cell abundance of the current sample using the previously published immune cell signature gene panel. Furthermore, a study was conducted to determine the correlation between m7G genes and immune cells. In order to analyze and visualize data, we used “ggplot2”, “limma”, and “GSEABase” packages.

Functional Enrichment Analysis

The present research is mainly related to the Kyoto Encyclopedia of Genes and Genomes (KEGG) and Gene Ontology (GO) enrichment analyses. Through assessing the signaling pathway and biological function in DEGs enrichment, we can further infer the possible pathogenesis. The “enrichplot”, “clusterProfiler”, and “org.Hs.eg.db” packages were used to perform the analysis.

The Nomogram Model

The machine learning of the RF and SVM classifiers was compared to determine the optimal algorithm to screen disease characteristic genes. The indicators compared included residual, reverse cumulative distribution of residual, and area under the curve (AUC) value. Next, the screened feature genes were obtained by machine learning, and a predictive nomogram was generated to assess the disease risk. In addition, clinical impact curve analysis (CICA) was used to assess clinical usability. Calibration curves and decision curve analysis (DCA) were used to evaluate the accuracy and measure the net benefit. We use the “lrmModel” function to construct the nomogram model. The packages “rms” and “rmda” were used to analyze and visualize data. One thousand bootstrap resamplings were used to validate the nomogram model.

Results

Patterns of m7G-Related Genes in COVID-19

A total of 29 m7G-related genes were retrieved from the GSEA databases and the literature, and the names of these 29 genes and their chromosomal locations are shown in Figure 1. In this study, 126 individuals, including 100 COVID-19 patients and 26 non-COVID-19 patients, were enrolled. We analyzed the difference in the expression of m7G-related genes between COVID-19 patients and non-COVID-19 patients. A total of 15 m7G-DEGs were found. The expression of NSUN2, NUDT4, NUDT4B, AGO2, EIF4E, EIF4E3, GEMIN5, LARP1, NCBP1, EIF4G3, IFIT5, and LSM1 was significantly up-regulated in COVID-19 patients (p<0.05). However, the expression levels of CYFIP1, EIF4E2, and EIF3D in COVID-19 patients were significantly down-regulated (p<0.05) (Figure 1B). In addition, a specific transcription profile was generated between COVID-19 and non-COVID-19 patients based on the expression levels of the m7G-DEGs (Figure 1C).

Figure 1 The expression and patterns of m7G genes. (A) The chromosome location and names of these 29 m7G genes. (B) Differential expression of m7G-related genes in COVID-19 and non-COVID-19 patients. (C) Heat map of m7G-related DEGs between COVID-19 and non-COVID-19. (D) The best consensus matrix (k=2). (E) m7G genes expression between m7G clusters A and B. (F) Heat map of m7G genes expression between m7G cluster A and B.

By analyzing the expression levels of m7G-DEGs to determine the subtypes of the disease, it is helpful to distinguish the disease phenotype, improve the detection rate, and improve the therapeutic effect. The PAM algorithm was used to determine the best consensus matrix. The best consensus matrix (k=2) was obtained (Figure 1D). Based on the m7G DEGs expression, 100 patients were divided into 2 categories, m7G-cluster A and B. m7G cluster A revealed high expression of NUDT4, NUDT4B, AGO2, EIF4E, EIF4E3, LARP1, EIF4G3, and IFIT5, while m7G cluster B was characterized by high expression levels of CYFIP1, EIF4E2, EIF3D, and LSM1 (Figure 1E). Furthermore, a specific transcription profile was generated between m7G cluster A and m7G cluster B (Figure 1F).

Immune Infiltration and Functional Enrichment Analysis of m7G Patterns

Principal component analysis (PCA) showed that the clusters could be distinguished based on m7G DEGs (Figure 2A). As the number of COVID-19 patients increases, it is particularly important to study the immune response mechanisms of the disease. We analyzed immune infiltration between m7G clusters A and B (Figure 2B). The results showed that m7G cluster B had high immune cell abundance, suggesting that m7G cluster B had significant immune cell infiltration. These results illustrate the differences in immune penetration among the m7G-based disease subtypes of COVID-19. To investigate the relationship between different m7G genes and immune cell infiltration, we further analyzed the correlation between 15 m7G DEGs and immune cells (Figure 2C). Surprisingly, most genes were positively correlated with immune cells, and the most significant positive correlation was found between EIF3D and immune cells. Furthermore, only six genes were negatively correlated with immune cells: NUDT4, NUDT4B, AGO2, EIF4E, EIF4E3, and EIF4G3.

Figure 2 Immune Infiltration and functional enrichment analysis of m7G patterns. (A) Principal component analysis indicated that clustering was effective. (B) The 23 kinds of immune cell infiltration between m7G cluster A and B. (C) The correlation between m7G-related genes and immune cells. (D) The KEGG enrichment analysis of DEGs between m7G cluster A and B. (E) The GO functional analysis based on DEGs between m7G cluster A and B.

In addition, to further explore the biological function between the two subtypes of disease, we first analyzed the DEGs between m7G clusters A and B and finally obtained 700 DEGs. These 700 DEGs were then analyzed using KEGG (Supplementary Table 1) and GO enrichment (Supplementary Table 2) analyses. Interestingly, the KEGG analysis results showed that enrichment was predominant in viral infection and immune-related aspects, including COVID-19, Th17 cell differentiation, Th1 and Th2 cell differentiation, and human T-cell leukemia virus 1 infection (Figure 2D). Furthermore, among the three aspects of GO enrichment analysis, the number of genes enriched in biological processes (BP) was the largest. BP was mainly enriched in immune cell regulation, including T cell activation (GO:0042110), positive regulation of cell activation (GO:0050867), positive regulation of leukocyte activation (GO: 0002696), leukocyte cell-cell adhesion (GO: 0007159), and regulation of T cell activation (GO: 0050863) (Figure 2E). Molecular functional enrichment analysis was mainly enriched in the structural constituents of ribosomes (GO:0003735), NAD+ nucleosidase activity (GO:0003953), MHC protein complex binding (GO:0023023), hydrolase activity, hydrolyzing N-glycosyl compounds (GO:0016799), MHC class II protein complex binding (GO:0023026), and NAD(P)+ nucleosidase activity (GO:0050135). Functional enrichment analysis of cellular component revealed that they were mainly enriched in ribosomal subunits (GO:0044391), cytosolic ribosomes (GO:0022626), ribosomes (GO:0005840), large ribosomal subunits (GO:0015934), cytosolic large ribosomal subunit (GO:0022625), and cytosolic small ribosomal subunit (GO:0022627). These results suggest that there are significant immune characteristics and related signals in the m7G clusters.

Clinical Feature and m7G Clusters

We sought to find out the relationship between m7G clusters and clinical features. The results showed significant differences in ventilator-free days (p<0.001), Hospital-free days post 45-day follow-up (p<0.001), ICU status (p=0.004), and mechanical ventilation status (p<0.001) between m7G clusters A and B. However, no significant differences in age, gender, and Charlson comorbidity index were observed between the two groups (Table 1).

Table 1 Clinical Features in m7G Clusters

Construction of m7G Signatures

To verify whether the classification of disease subtype based on m7G is reasonable, consensus clustering analysis was conducted again based on 700 DEGs, obtained from m7G clusters. Patients were divided into different genomic subtypes. As shown in Figure 3A, k=2 was selected as the optimal consensus matrix. Based on the expression of 700 DEGs, COVID-19 patients were grouped into two subtypes: gene cluster A and gene cluster B.

Figure 3 (A) The best consensus matrix. (B) Heat map of gene expression between gene clusters A and B. (C) The expression of m7G-related DEGs between gene clusters A and B (D) The 23 kinds of immune cell infiltration between gene cluster A and B. The m7G-score in gene clusters (E) and m7G clusters (F). (G) Sankey diagram to further visualize the relationship among the disease subtypes classification, m7G clusters, gene clusters, and m7G-score. *p-value<0.05, **p value<0.01, ***p value < 0.001.

A specific transcription profile was generated between gene clusters A and B (Figure 3B). Next, we analyzed the expression of m7G DEGs between gene clusters A and B. The results showed that m7G DEGs also showed significant differences between the gene clusters (Figure 3C). Moreover, we analyzed the abundance of immune cells between gene clusters A and B. Twelve types of immune cells were highly infiltrated in gene cluster A, whereas only four types of immune cells were highly infiltrated in gene cluster B. Gene cluster A showed a higher abundance of immune cell infiltration (Figure 3D).

In addition, we constructed a scoring system, the m7G-score, which is based on a composite score of m7G gene expression in COVID-19 patients. Each patient received a corresponding m7G-score. We compared the scores of the different disease subtypes. The results revealed a higher m7G score in gene cluster B than in cluster A (p<0.001) (Figure 3E). However, the m7G score was higher in m7G cluster A than in m7G cluster B (p<0.001) (Figure 3F). Interestingly, gene cluster B and m7G cluster A, both with higher scores, were in the group with low immune infiltration. These results suggest that the two subtypes, m7G clusters and gene clusters, have favorable immune characteristics and genetic differences.

According to the median value of the m7G-score, patients can be divided into m7G-score high and low groups. To better distinguish the subtypes of the disease in these three categories, we compared the patient populations of the three groups. As shown in Figure 3G, the patient groups obtained using these three methods showed a high degree of coincidence.

The Nomogram Model

A assessment of the risk of infection with COVID-19 helps control the infection and provides guidance for targeted treatment. We attempted to establish a prediction model based on m7G genes to evaluate disease risk. Two different machine-learning methods, RF and SVM, were performed to select the characteristic genes of the disease. By comparing the residual distribution (Figure 4A) and reverse residual cumulative distribution (Figure 4B) of the two algorithms, we evaluated which algorithm was better. The results found that RF algorithm has a higher accuracy and is suitable for screening COVID-19 feature genes. Moreover, the AUC for RF (AUC=1) was higher than that for SVM (AUC=0.981) (Figure 4C). Therefore, we screened the genes using RF. The optimal random forest model was constructed using 71 trees (Figure 4D). The importance scores for the m7G genes were then calculated. A total of six genes with importance scores higher than 2.5 were included in our next prediction model, including NUDT4B, IFIT5, LARP1, EIF4E, LSM1, and NUDT4 (Figure 4E).

Figure 4 Establishment and validation of the prediction model. Residual (A) and reverse cumulative distribution of residual (B) were used to compare the accuracy of RF and SVM algorithms. (C) RF and SVM were evaluated by ROC curve. (D) The optimal nTree was obtained to find the disease feature genes. (E) The importance score of the relevant gene was calculated. (F) The nomogram prediction model. (G) The calibration curve revealed a good degree of agreement. (H) Decision curve analysis showed a high net benefit. (I) Nomogram models were validated using clinical impact curve analysis.

Next, we combined the effects of these 6 genes on the risk of COVID-19 to build a predictive nomogram for risk evaluation. The corresponding score was assigned according to the gene expression level. The combined total score of the six scores could be applied to evaluate the risk (Figure 4F). The calibration curve for the nomogram model is shown in Figure 4G. The results revealed that the prediction model fits well with the ideal curve. DCA demonstrated that the prediction nomogram based on m7G genes could obtain more net benefits (Figure 4H). The CICA results suggested that patients with a high risk predicted by the model were highly similar to the patients who were actually positive (Figure 4I).

Discussion

In the present study, we found that the highest expression of m7G genes in COVID-19 patients was significantly up-regulated compared to that of non-COVID-19 patients, suggesting that m7G modification may play a key role in the pathogenesis of viral infection. Furthermore, two subtypes were obtained by clustering based on m7G gene expression. Analysis of immune cell infiltration and clinical characteristics showed that compared with m7G-cluster A, m7G-cluster B showed higher immune cell infiltration. The number of patients in cluster B requiring ICU and mechanical ventilation was significantly lower than that in cluster A. This suggests that clusters B had mild symptoms. These results further suggest that cluster B, with rich immune infiltration, may have a better prognosis. This hypothesis was consistent with the results reported by Wang et al.²⁷ In addition, in the functional enrichment analysis of DEGs between m7G-cluster A and B, it was found that the main signaling pathways were involved in COVID-19 and immune cell differentiation, and biological functions were also enriched in the activation and regulation of immune cells. These results suggest that m7G plays an important regulatory role in the occurrence and development of COVID-19.

After SARS-CoV-2 infects human cells, the antiviral response of the host and pathogenicity of the virus affect the final clinical outcome. The mutant strain of SARS-CoV-2 RNA has appeared in just three years, and the transmission ability has strengthened, which makes it particularly important to study the pathogenesis of viral infection and the host immune response. The inflammatory response induced by alveolar macrophages leads to acute respiratory distress syndrome, which involves a variety of inflammatory factors, such as TNF, IL12, and IL6.^28,29 The number of CD8+T cells in patients with COVID-19 is decreased, which can be used as a prognostic marker of disease severity.²⁷ However, a previous study has shown that the response of CD4+T cells to SARS-CoV-2 is more common than that of CD8+T cells, and CD4+ T cells are related to the titers of anti-SARS-CoV-2 IgG and IgA.³⁰ The spike protein, as the main antigen to recognize SARS-CoV-2, is also the main target of specific CD4+T cells.³¹ Among virus-specific CD4+ T cells, follicular helper T cells account for a large proportion of SARS-CoV-2-specific CD4 T cells in patients with acute and convalescent COVID-19.³² Although the response of CD4+ T cells to SARS-CoV-2 infection is interesting, the study unexpectedly found that peak reactive CD4+ T cells were detected in 35% of healthy people,³³ which means that the possibility of cross-reaction cannot be ruled out. The complex immune response in the body and various symptoms of clinical patients indicate that the key to the success of antiviral therapy can never rely on single immune cells alone. COVID-19 patients who received neutralizing monoclonal antibodies alone did not benefit greatly.³⁴

Some viruses, such as autologous RNA, can successfully escape the recognition of retinoic acid-Inducible gene-I through post-transcriptional m7G modification, thus achieving immune escape.³⁵ The structure of the 5’ m7G cap involves methyltransferases, RNA 5’ triphosphatase (RTPase), and guanylytransferase (GTase). However, capping of SARS-CoV-2 RNA has rarely been studied. A SARS-CoV nsp13 helicase has been shown to have RTPase activity, which may play a role in the formation of the 5′ cap structure.³⁶ In addition, the SARS-CoV-2 nsp12 has nucleotidylation activity. There is evidence that the cap core structure, GpppA, is catalyzed by SARS-CoV-2 nps12 through its GTase activity.³⁷ Replication and transcription complexes assembled by a group of nps, including nps13 and nps12, facilitate the transcription of SARS-CoV-2 mRNA.

We constructed a prediction model based on six m7G genes (NUDT4B, IFIT5, LARP1, EIF4E, LSM1, and NUDT4) to predict the risk of COVID-19. Interestingly, some of these genes have been shown to be closely related to viral infection. eIF4E has been confirmed to be involved in the translation of some viral mRNA and proliferation of infected cells.³⁸ SARS-CoV-2 has been shown to replicate effectively in target cells via the ERK/MNK1/eIF4E signaling pathway. Emetine inhibits coronavirus replication by inhibiting the interaction between viral mRNA and eIF4E.³⁹ LARP1, as a host protein that can interact with nucleocapsid proteins, can also bind to SARS-CoV-2 RNA and inhibit viral replication in vivo.⁴⁰ Although the specific mechanism of other genes in coronaviruses have not been confirmed, they are related to other viral infections. LSm1 participates in the replication of hepatitis C and dengue virus.^41,42 However, studies have confirmed that dengue virus RNA is related to 63% of the RNA interactome of SARS-CoV-2.⁴⁰ When IFIT5 is stimulated by viral infection, it promotes the antiviral response in host cells.⁴³ The prediction model based on these genes helped us further evaluate the risk of disease.

This study had some limitations. First, the analysis of immune cell infiltration was estimated, which may not fully reflect the actual infiltration of immune cells in COVID-19 patients. However, this needs to be verified using real samples. Second, limited clinical information and a lack of other basic clinical features, such as lung imaging findings, basic disease conditions, treatment plans, and vaccines, may lead to different experimental results.

Conclusion

In conclusion, the expression of m7G-related genes in COVID-19 patients was higher than that in non-COVID-19 patients. There was a significant correlation between m7G genes and immune cell infiltration. The m7G-cluster B was related to higher immune infiltration and milder symptoms. The nomogram based on the six m7G genes accurately predicted the risk of COVID-19.