Integrated Bioinformatics Analysis Reveals Key Candidate Genes and Cytokine Pathways Involved in COVID-19 After Rhinovirus Infection in Asthma Patients

Background

Respiratory virus infections are closely associated with the exacerbation of asthma, in which rhinovirus (RV) is the most common pathogen [1]. Coronavirus disease 2019 (COVID-19), caused by severe acute respiratory syndrome coronavirus 2 (SARS-COV-2), may be more serious in patients with chronic lung diseases, including asthma [2]. In asthma patients, biological factors may affect the susceptibility to SARS-COV-2 or the severity of COVID-19. Kim et al. found that RV infection was also the most common additional respiratory pathogen in 116 SARS-COV-2-positive specimens [3]. It is well known that there is a great diversity in RV, in which RV-B infection usually causes mild or asymptomatic symptoms, and RV-A or RV-C are usually associated with more severe illness or exacerbation of asthma [4]. RV infection induces the expression of interferon (IFN)-stimulated genes and subsequent production of cytokines [5].

Angiotensin-converting enzyme 2 (ACE2), a transmembrane protein with an extracellular carboxypeptidase domain located on the cell membrane of a variety of epithelial cells, was originally described as playing a role in the renin-angiotensin system [6]. The various roles of ACE2 have gradually been discovered [7]. It participates in the regulation of cardiovascular physiology, dietary amino acid homeostasis, innate immunity, and intestinal microbial ecology, and plays an important role in heart failure, hypertension, myocardial infarction, and cardiovascular complications of diabetes [8]. It is also a cellular receptor for SARS-COVs, including SARS-COV-2 which caused the COVID-19 epidemic [9]. The spike protein of virus binds to the ACE2 of the host cell, mediating viral invasion of cell membranes and transmission to other cells [10]. Therefore, we speculated that there are differences in ACE2 gene expression in RV-infected asthma patients, and it may be associated with severe COVID19 infections.

Gene expression microarray, as an efficient and large-scale technique to obtain genetic data, is widely used to study the gene expression profile of diseases and plays an important role in screening differentially expressed genes (DEGs) [11]. With the wider application of this technology, a large quantity of microarray and sequencing data has been collated and stored in the public database platforms, which are available for free download. Further analysis and integration of the data through relevant statistical software or online analysis tools can lead to more in-depth study of molecular mechanisms of disease, providing a new method for the study of disease-related genes. However, the results of different studies on the identification of DEGs are not completely consistent, and there are still some limitations in single-cohort studies. Many studies have identified the DEGs of different diseases, including ovarian cancer [12], leukemia [13], and gastric cancer [14]. The identification of DEGs has been widely used in the analysis of various types of diseases and can provide effective therapeutic targets for the disease.

In this study, the original microarray dataset GSE149273 [15] was downloaded from the NCBI-Gene Expression Omnibus (GEO), which contained a total of 90 samples. DEGs were identified in RV-A16 infections (RVA) or RV-C15 infections (RVC) compared to non-infection (control) in patients with asthma. Additionally, gene ontology (GO), Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis, and Disease ontology (DO) analysis were performed to identify the intersecting DEGs of RVA and RVC. The gene expression of ACE2 and the proportion of 22 kinds of immune cells in each sample were further analyzed to determine the effects of FVA or FVC infection on cytokines. Overall, the findings of the present study highlighted key candidate genes and cytokine pathways in RV infection in people with asthma. These may provide a reliable molecular marker to better understand RV infection and analyze its relationship with COVID-19.

Material and Methods

DATASET SOURCES:

The GSE149273 profile dataset was downloaded from the NCBI-GEO database (https://www.ncbi.nlm.nih.gov/geo/). The platform for GSE149273 is GPL21290 Illumina HiSeq 3000 (Homo sapiens), which includes 90 samples from 30 asthma patients with non-infection, RVA, and RVC. Platform and series matrix file(s) were downloaded as TXT files.

EXPRESSION ANALYSIS OF DEGS:

The original data were converted into an expression matrix, and then the mRNA expressions of 90 samples were normalized using the R software package with the robust multi-array average algorithm. The t test method in the limma R package was subsequently used to identify DEGs between RVA and RVC infected samples and non-infected samples, respectively. The log fold change >2 and P<0.05 were considered as the screening criteria for DEGs identification. The intersection of the analysis results of DEGs of RVA and RVC was demonstrated using VennDiagram in the R package, and the heatmap.2 of R package was used to draw the heatmaps of RVC, RVA, and control groups.

FUNCTIONAL ANALYSES:

To further explore the biological processes and signal pathways in which the intersecting DEGs of RVA and RVC might be involved, functional analyses of the DEGs at the intersection of RVA and RVC were performed. GO analysis was performed to demonstrate the dominant function of the intersecting DEGs from molecular function and biological process using the clusterProfiler R package. KEGG pathway analysis was performed to find the correlation between the intersecting DEGs and signal pathways using the clusterprofiler R package. DO analysis was performed to explore the relationship between the intersecting DEGs and diseases using DisGeNET data (https://www.disgenet.org/). All the results were visualized, and P<0.05 and gene counts >5 were considered to show a statistically significant difference in this functional enrichment analysis.

ANALYSIS OF ACE2 AND IMMUNE CELLS:

The t test was used to analyze the differences in gene expression of ACEs between RVA infection or RVC infection and control. P<0.05 was statistically significant. Moreover, the proportion of 22 immune cells in each sample was estimated to analyze the infiltration of immune cells after infection using cibersortx (https://cibersortx.stanford.edu/index.php), and a heatmap was drawn. The vioplot in the R package was used to visualize the change in the proportion of 22 immune cells relative to the control after RVA or RVC infection. P<0.05 was considered statistically significant.

Results

DATA INFORMATION AND IDENTIFICATION OF DEGS:

The GSE149273 profile dataset consisted of 30 patients with a median age of 35 years and a balanced sex ratio, 60% of whom were non-Hispanic white individuals. Thirty individuals were treated with non-infection, RVA, and RVC infection respectively, and 90 samples were obtained. The mRNA expression levels of the 90 samples in the microarray dataset GSE149273 were standardized, and the results are shown in Figure 1. The GSE149273 dataset was screened using the limma R package (logFC >2 and P<0.05), and 555 DEGs were obtained between RVA and control and 421 DEGs between RVC and control. Among them, 529 upregulated genes and 26 downregulated genes were identified in RVA, while 408 upregulated genes and 13 downregulated genes were classified in RVC. In addition, there were a total of 415 DEGs both in RVA and RVC, among which 406 were upregulated and 9 were downregulated (Table 1). The differential expression of multiple genes from FVA or FVC and control in the microarray is shown in Figure 2A and 2B, and the intersection of RVA and RVC DEGs is shown in Figure 2C. The cluster heatmaps of FVA and FVC are shown in Figure 3.

FUNCTIONAL ENRICHMENT ANALYSIS:

To identify which biological functions and pathways were most significantly associated with the DEGs, a total of 415 genes in the intersection of RVA and RVC were analyzed for GO and KEGG pathways. GO analysis of DEGs in this study was divided into 2 functional groups: a biological processes group and a molecular function group (Figure 4A). In the biological processes group, the DEGs were mainly enriched in “response to virus”, “defense response to virus”, and “response to interferon-gamma”. In the molecular function group, the DEGs were mainly enriched in “receptor ligand activity”, “cytokine activity”, and “cytokine receptor binding”. The top 10 GO terms for each category are listed in Table 2. KEGG pathway analysis demonstrated that DEGs were significantly enriched in “Cytokine-cytokine receptor interaction”, “Influenza A”, and “NOD-like receptor signaling pathway” (Figure 4B). All KEGG pathway terms are listed in Table 3. These results indicate that most DEGs were significantly enriched in virus and cytokine. To further clarify the relationship between DEGs and various diseases, 415 DEGs were subjected to DO analysis. The specific conditions of diseases enriched by each gene are shown in Figure 5. After enrichment analysis, it was found that the DEGs were mainly enriched in “Influenza”, “Hepatitis C, Chronic”, and “Respiratory Syncytial Virus Infections”. The enrichment results are shown in Figure 6A, and the interaction network between diseases and DEGs is shown in Figure 6B.

EXPRESSION OF ACE2 GENE AND IMMUNE CELLS:

The ACE2 gene is one of the 415 DGEs, and this study further confirmed that the expression of ACE was significantly increased after RVA or RVC infection (P<0.01), along with upregulated DGEs (Figure 7). Functional enrichment analysis showed that most DEGs were significantly related to cytokines, so we analyzed the content and proportion of 22 kinds of immune cells in 90 samples. The content of immune cells is shown in Figure 8A and the heatmap of cell proportions is shown in Figure 8B. The results showed that the proportion of “NK cells activated”, “Macrophages M1”, and “Dendritic cells activated” was increased after infection with RVA or RVC, while that of “T cells CD8” and “Mast cells resting” was decreased. Based on the correlation analysis of the immune cells, we found that “T cells follicular helper”, “Dendritic cells activated”, “Neutrophils”, and “NK cells activated” had a positive correlation with “Macrophages M1”, and “T cells CD8” was also positively correlated with “Dendritic cells resting” (Figure 8C).

To identify changes in the expression of individual immune cells after infection, we visualized the proportion of 22 kinds of immune cells relative to control after RVA or RVC infection. The results showed that the proportions of “B cells naïve”, “T cells CD8”, “NK cells resting”, “Macrophages M0”, “Macrophages M2”, “Dendritic cells resting”, and “Mast cells resting” were significantly decreased, while that of “T cells follicular helper”, “NK cells activated”, “Macrophages M1”, “Dendritic cells activated” and “Neutrophils” increased significantly after RVA infection (Figure 9A). And the proportion of “B cells naïve”, “T cells CD8”, “NK cells resting”, “Macrophages M0”, “Dendritic cells resting”, and “Mast cells resting” was significantly decreased, while “T cells gamma delta”, “NK cells activated”, “Macrophages M1”, “Dendritic cells activated”, and “Neutrophils” was significantly increased after RVC infection (Figure 9B). In general, RV infection significantly activated the expression of NK cells, dendritic cells, and neutrophils, and reduced the expression of CD8 T cells.

Discussion

SARS-COV-2, causing COVID-19, is a novel virus first detected in December 2019 in Wuhan, China and the disease is thought to be more severe in patients with chronic lung disease [16]. RV airway infection is the most common virus causing the common cold and is a major contributor to the development and exacerbation of asthma [17]. RV infection induces the expression of IFN-stimulated genes and subsequent cytokines [18], and ACE2, the receptor of SARS-COV-2, has recently been identified as an IFN-stimulated gene [19]. Therefore, RV infection may be associated with severe COVID-19 infection through key candidate genes and cytokine pathways. Studies have shown that there are large clinical variations due to differences in autoantibodies among individuals in the course of SARS-COV-2 infection, and type I IFNs in protective immunity against SARS-Cov-2 are crucial [20]. There are many immunological tests for SARS-COV-2 infection, but the protective immune genes for SARS-CoV-2 have not yet been fully clarified, and the study of genes and cytokines can provide a unique opportunity to control the infection and disease [21]. Microarrays and high-throughput sequencing have been widely used to detect the expressions of various genes in the human genome and to predict potential targets for related diseases, which can help identify target genes for diagnosis or treatment of diseases [22].

In the present study, the genomic dataset was obtained from GEO, including 90 samples with non-infection, RVA infection, and RVC infection. Bioinformatics analysis was conducted, identifying 555 DEGs between RVA and control and 421 DEGs between RVC and control. A total of 415 DEGs were involved both in RVA and RVC, among which 406 were upregulated and 9 were downregulated. Functional enrichment analysis showed that DEGs were mainly enriched in “response to virus”, “response to interferon-gamma”, “receptor ligand activity”, “cytokine activity”, “cytokines-cytokine receptor interaction”, “NOD-like receptor signaling pathway”, “Influenza A”, and “respiratory syncytial virus infection”. These results indicate that most DEGs were significantly enriched in virus and cytokine, which is consistent with normal expression after virus infection.

SARS-COV-2 infection is mediated by a Spike virus protein and ACE2, and fuses with target cells through a type I mechanism [23]. ACE2 is expressed in the alveoli, heart, esophagus, kidney, bladder, and ileum, suggesting that all of these organs may be affected by SARS-COV-2 [24,25]. The expression of ACE2 mediates SARS-COV-2 infection in host lung cells [10], so it can be reasonably concluded that increased expression of ACE2 in lung cells will increase susceptibility to SARS-COV-2 infection or lead to more severe COVID19 disease. Generally, after cellular detection of viral entry into a host cell, IFN-induced IFN-stimulated genes are essential for host antiviral defense, and studies have shown that ACE2 is an INF-stimulated gene in human airway epithelial cells [19]. Inducers of the IFN pathway, including influenza infection, can lead to increased expression of ACE2 in cells, thereby increasing susceptibility to COVID-19. Our study found that ACE2 was one of the 415 DEGs in RV and control, and there was a significant increase in the expression of ACE2 in asthma patients infected with RV, indicating that RV infection can increase susceptibility to SARS-COV-2 or lead to more severe COVID19.

Immune cells are the cells involved in or related to the immune response, and they play an important role in the human body. Studies have shown that NK cells are activated in COVID-19 [26], and can also be activated in humans after influenza virus infection, while dendritic cells are necessary for NK cell activation [27]. Dendritic cells play an important role in defense against lung infections [28]. In human dendritic cells derived from monocytes infected by SARS-COV, the chemokines are upregulated [29]. Furthermore, lymphocytosis may be due to the direct effects of cytokines on T cell populations and/or the indirect effects on other cell types such as dendritic cells [30] and neutrophils [31,32] in severe disease. Our research found that NK cells, dendritic cells, and neutrophils were activated after RV infection in asthma, which was consistent with the previous description. In addition, we found a significant decrease in CD8 T cells of people with asthma infected with RV, which was similar to earlier observations about SARS-CoV-1 infection. Some studies reported that lymphopenia was caused by a sharp decrease in the number of T cells CD4 and CD8 in COVID-19 patients [33]. T cells play an important role in viral infections, in which CD4 T cells help to produce antibodies and coordinate the response of other immune cells, while CD8 kills infected cells to reduce the viral burden, and lower numbers of CD8 cells have been shown to be correlated with severity of COVID-19 [34]. The hyperactivation of monocytes, macrophages, and dendritic cells is the driving factor of human immunopathology. Our study found a significant correlation between CD8 and dendritic cells after RV infection, indicating that asthma patients infected with RV would be more susceptible to COVID-19, which was also consistent with the higher incidence of COVID-19 among young asthmatics, accounting for 27% of the hospitalized 18–49 age group in the United States [35]. All results showed that viral infections associated with asthma exacerbations and SARS-COV-2 infections showed synergistic biomolecular interactions. RV infection can aggravate COVID-19 infection by increasing key gene and cytokine pathways associated with COVID-19.

Conclusions

Using integrated analysis of microarray datasets downloaded from the NCBI-GEO database, the present study identified 415 DEGs that may be involved in the progression of RV infection in people with asthma. These DEGs were mainly related to cytokines. Our results contribute to a better understanding of the molecular mechanism of RV infection. We found that there was significant activation of NK cells, dendritic cell, and neutrophils, and a significant decrease in CD8 T cells after RV infection in people with asthma. Moreover, the expression of ACE2, which is closely related to COVID-19, was significantly increased after RV infection. The evidence suggests that RV infection activates the cytokine pathways and increases the expression of ACE2, thus exacerbating COVID-19. These results enhance our understanding of RV infection and its correlation with COVID-19. However, further molecular biological experiments are required to confirm the results of this study.