The Landscape of Coronavirus Disease 2019 (COVID-19) and Integrated Analysis SARS-CoV-2 Receptors and Potential Inhibitors in Lung Adenocarcinoma Patients

Introduction

The outbreak of coronavirus disease 2019 (COVID-19) caused by severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) began in Wuhan, China, in December 2019 (Gorbalenya et al., 2020) and then spread rapidly around the world (Berry et al., 2020; Vaid et al., 2020; Wang et al., 2020). As of Oct 2nd, 2020, the World Health Organization (WHO) has reported more than 34.40 million confirmed positive cases and over 1.02 million deaths¹. COVID-19 is a potentially fatal disease that has aroused public health concerns worldwide, and its high morbidity and mortality pose a major challenge as a public health problem for China and many other countries.

The severity of COVID-19 can range from mild to moderate with infection symptoms such as fever, cough, and fatigue to severe and fatal with symptoms characterized by respiratory dysfunction and/or multiple organ failure (Huang et al., 2020; Rothan and Byrareddy, 2020). Currently, the identification of the progression of COVID-19 mainly depends on clinical manifestations. There are no established biomarkers that can effectively predict disease progression. The treatment strategy for COVID-19 patients mainly focuses on providing supportive care, such as oxygenation, ventilation, and infusion management, and there are no vaccines or recognized specific antiviral drug regimens for the treatment of critical patients (Kupferschmidt and Cohen, 2020; Tobaiqy et al., 2020).

Coronaviruses use their spike proteins to bind and enter target cells through specific receptors. Severe acute respiratory syndrome coronavirus (SARS-CoV) was confirmed to use angiotensin converting enzyme 2 (ACE2) as a receptor to enter cells, and SARS-CoV-2 and SARS-CoV have 79.5% homologous sequences. Both bioinformatics modeling and in vitro experiments demonstrated that ACE2 was recognized as the receptor for SARS-CoV-2 entry into target cells, and the cellular protease TMPRSS2 for SARS-CoV-2 spike protein priming is also essential for target cell entry and spread (Hoffmann et al., 2020; Yan et al., 2020; Zhou et al., 2020). It is worth noting that cancer patients are more likely to be infected with COVID-19 than patients without cancer (Longbottom et al., 2016; Liang et al., 2020). ACE2 is also aberrantly expressed in many tumors, such as LUAD (Chai et al., 2020). SARS-CoV-2 enters target cells through ACE2-mediated endocytosis, and AP2-associated protein kinase 1 (AAK1) is one of the known regulators of endocytosis (Richardson et al., 2020). The destruction of AAK1 may interrupt the spread of the virus to cells and the intracellular assembly of virus particles (Lu et al., 2020). Considering the crucial role of ACE2, TMPRSS2, and AAK1 in the SARS-CoV-2 invasion of target cells and intracellular transmission, performing an analysis of the expression and distribution characteristics and related biological processes of ACE2, TMPRSS2, and AAK1 in cancer patients can help us better understand the pathogenesis of their susceptibility to COVID-19 and better explore potential novel treatment strategies.

In this study, we first performed enrichment pathway analysis and protein-protein interaction (PPI) network analysis on genes with significantly different expression characteristics between COVID-19 patients and normal subjects by analyzing the GSE147507 dataset in order to explore the genes and their functions that change significantly in the progression of COVID-19. Then, we focused on the expression characteristics and functions of SARS-CoV-2 Receptors including ACE2, TMPRSS2, and AAK1 in LUAD as well as their relationship with the level of immune infiltration and explored the potential therapeutic targets and drugs for ACE2, TMPRSS2, and AAK1, then further indicated the anti-inflammatory effect of SARS-CoV-2 Receptors inhibitor bexarotene, providing promising ideas for the treatment of COVID-19.

Materials and Methods

Acquisition of Differentially Expressed Genes (DEGs) Between Healthy Subjects and COVID-19 Patients

The mRNA sequencing data of COVID-19 patients were obtained from the GSE147507 dataset (containing 2 normal samples and 2 COVID-19 patient samples) in the Gene Expression Omnibus (GEO) database (Blanco-Melo et al., 2020). Differential expression analysis is performed through the R package “edgeR.” Before performing the differential expression analysis, the variance of gene expression will be estimated; the degree of internal gene expression difference is estimated to see whether the difference in gene expression between groups is greater than the internal difference. If it is, the gene is defined as a differentially expressed gene. Then we use the package “edgeR” to standardize and the standardization method is as follows: In edgeR, the TMM method computes normalization factors that represent sample-specific biases. These factors are multiplied by the library size to yield the effective library size, i.e., the library size that we would have gotten if those biases were not present. The effective library sizes can then be used for various normalization purposes, most frequently as offsets in generalized linear models. We identified the obtained false discovery rate (FDR) as the adjusted P-value. The DEGs between COVID-19 patients and normal subjects were determined under the condition of absolute log2-fold change (FC) > 1 and the adjusted P-value < 0.05.

Functional Enrichment Analysis

Gene Ontology (GO) analysis was performed in The Database for Annotation, Visualization and Integrated Discovery (DAVID)². Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis of the DEGs was performed using ClueGO and CluePedia in Cytoscape software (version 3.5.1) (Shannon et al., 2003). P < 0.05 was considered to represent statistical significance.

Protein-Protein Interaction Network Analysis

The BioGRID database³ was used to perform PPI network analysis on the proteins encoded by the identified DEGs (Stark et al., 2006). The PPI network was constructed in an online gene function analysis website- Cytoscape software⁴, and the key PPI network modules were selected by the Molecular Complex Detection (MCODE) plugin (Bader and Hogue, 2003).

Estimation of the Infiltration of Multiple Types of Immune Cells

Cell-type Identification By Estimating Relative Subsets Of RNA Transcripts (CIBERSORT) analysis was applied to analyze the absolute abundance of 21 immune cells in heterogeneous tissues to evaluate the correlation between the level of immune cell infiltration and expression of ACE2, TMPRSS2, and AAK1 (Gentles et al., 2015). The R package “CIBERSORT” was used to convert the mRNA data of non-tumor cells into the tumor microenvironment infiltration level. For each sample, the sum of all estimated immune cell type scores is equal to 1. And Tumor Immune Estimation Resource (TIMER)⁵ (Li et al., 2017) was also used for interactively explore the associations between immune infiltrates and gene expression.

Weighted Gene Co-expression Network Analysis (WGCNA)

The DEGs between LUAD patient samples and normal samples in the TCGA database were used to construct a gene co-expression network with the R package “WGCNA” (Langfelder and Horvath, 2008). First, the transcription data of the DEGs were entered to establish an expression matrix. The best “soft thresholding power” β was graphically determined to ensure the scale-free distribution network, and a co-expression matrix was constructed using the expression matrix and the β value. Using the DynamicTreeCut algorithm, genes with similar expression levels were classified into the same gene module, resulting in co-expression modules. The module Eigengenes function of the WGCNA R package was used to calculate the differences of the module eigengenes (ME) and the relevance of the modules to ACE2, TMPRSS, and AAK1. A heatmap was used to visualize the relevance of each module. Pearson’s coefficient was calculated to assess the correlations between the genes in the modules and ACE2, TMPRSS, and AAK1.

Determination of the Co-regulated miRNA Factor for ACE2, TMPRSS, and AAK1

miRNAs targeted ACE2, TMPRSS2, and AAK1 were predicted in miRWalk⁶, miRanda⁷, RNA22⁸, Targetscan⁹, and RNAhybrid¹⁰ with different algorithms. The databases to predict targeted miRNA with different algorithms predict a set of miRNA IDs rather than specific regulated sequences. miRNAs which are repeatedly predicted in at least 2 databases are identified as the targeted miRNAs of the genes. Then we overlapped the IDs of the three sets of miRNAs predicted by the three genes to obtain the IDs of the co-regulated miRNAs of the three genes, and subsequently we used the RNAhybrid¹⁰ database to determine the specific sequences of co-regulated genes regulating these three genes.

Exploration of the Potential Targeted Drug for ACE2, TMPRSS2, and AAK1

Potential targeted drugs for ACE2, TMPRSS2, and AAK1 were explored in the GDSC database¹¹. We selected drugs whose response sensitivity is positively correlated with gene expression levels as the targeted drugs for the genes. And the targeted drugs for ACE2, TMPRSS2, and AAK1 were gathered to determine the potential targeted drug that can co-regulate these three genes.

Results

Characterization of DEGs Between Healthy Subjects and COVID-19 Patients

The expression profile of the DEGs between healthy subjects and COVID-19 patients is shown in Figure 1A. Hierarchical clustering analysis was adopted to categorize genes with similar expression profiles into two groups. Through KEGG analysis of these DEGs, we found that the genes are mainly enriched in inflammation and immune-related signaling pathways, such as “lymphocyte activation,” “response to virus,” “antiviral mechanism by IFN-stimulated genes,” “regulation of immune effector process,” and “regulation of inflammatory response” (Figure 1B).

FIGURE 1

Figure 1. Hierarchical clustering analysis and KEGG analysis of DEGs. (A) Hierarchical clustering analysis showed the expression of DEGs. Each row represents a single gene and each column represents a sample. (B) KEGG analysis presented the main enrichment pathways of DEGs.

Protein-Protein Interaction Network Analysis of DEGs

We performed PPI network analysis on the proteins encoded by the identified DEGs in the BioGRID database and constructed the PPI network in the Metascape website (Figure 2A; Zhou et al., 2019). The MCODE plugin in Metascape was used to determine important modules and related hub genes in the constructed PPI network, and chemokines such as CCR1, CCL19, CXCL10, and CXCL16 were identified to present the most obvious changes in COVID-19 disease progression (Figure 2B).

FIGURE 2

Figure 2. PPI network and MCODE components were associated with DEGs. (A) PPI network of proteins encoded by DEGs. (B) The essential modules identified by MCODE from the PPI network. Level >6 is set as the cutoff criterion.

The Expression Profiles of ACE2, AAK1, and TMPRSS2 in LUAD Patients and COVID-19

The expression characteristics of SARS-CoV-2 receptors including ACE2, AAK1 and TMPRSS2 between LUAD patients and normal subjects were analyzed in the Oncomine database. ACE2 and TMPRSS2 exhibited significant overexpression in LUAD patients (Figures 3A,B), while there was no significant difference in the expression level of AAK1 between LUAD patients and normal subjects (Figure 3C). Then, we evaluated the expression of ACE2 and TMPRSS2 in lung cancer patients through IHC staining, and markedly upregulated expression of the two genes was presented in lung cancer patients compared with normal samples (Figures 3D,E). The finding that ACE2 and TMPRSS2 are highly expressed in patients with lung cancer may be part of the reason why LUAD patients are more susceptible to SARS-CoV-2 infection. The expression patterns of ACE2, AAK1, and TMPRSS2 were also evaluated in COVID-19 patients in GSE147507 to explore gene changes in the progression of the disease (Figures 3F–H). The results showed that compared with normal subjects, the expression levels of AAK1 and TMPRSS2 in patients who had been infected with SARS-CoV-2 were significantly down-regulated, suggesting that COVID-19 patients may consume lots of AAK1 and TMPRSS2 in the advanced stage of the disease.

FIGURE 3

Figure 3. Expression characteristics of ACE2, AAK1, and TMPRSS2 in LUAD patients and COVID-19. (A–C) Differences in expression of ACE2 (A), TMPRSS2 (B), and AAK1 (C) between LUAD patients and normal subjects. (D–E) IHC staining showed expression patterns of ACE2 (D) and TMPRSS2 (E) in LUAD patients and normal samples. (F–H) Expression levels of ACE2 (F), TMPRSS2 (G), and AAK1 (H) in COVID-19 patients compared to normal subjects.

The Expression Levels of ACE2, TMPRSS2, and AAK1 Are Closely Related to the Infiltration of Immune Cells

The correlations between the infiltration of 21 immune cells and the expression levels of ACE2, TMPRSS2, and AAK1 performed through CIBERSORT analysis are presented in Figures 4A,C,E. The immune cells and their infiltration significantly affected by the expression levels of ACE2 and TMPRSS2 were relatively consistent; the infiltration levels of naive B cells, activated dendritic cells, CD4 memory T cells and CD4 regulatory cells were up-regulated in the ACE2 high-expression group and TMPRSS2 high-expression group (Figures 4B,D). The immune cell infiltration of patients with different AAK1 expression levels is presented in Figure 4F. The immune infiltration analysis results of TIMER2 were also relatively consistent with CIBERSORT analysis (Supplementary Figure 2A). The infiltration levels of B cells and CD4 T cells significantly increased in the groups with high expression of ACE2 and high expression of TMPRSS2 (Supplementary Figure 2B,C). And the AAK1 high-expression group also showed a high level of CD4 T cell infiltration (Supplementary Figure 2D). The cBioportal dataset¹² (Cerami et al., 2012) showed ACE2 level was positively correlated with TMPRSS2 level (Figure 4G), the result maybe partially explained the consistence of immune cells infiltration. However, the correlation of other two genes was not significant difference (Figures 4H,I).

FIGURE 4

Figure 4. The correlation between expression patterns of ACE2, TMPRSS2, and AAK1 and the infiltration of immune cells in LUAD patients. (A,C,E) The expression levels of ACE2 (A), TMPRSS2 (C), and AAK1 (E) are closely related to the infiltration characteristics of immune cells. (B,D,F) Differences in the infiltration of 21 immune cells between the high and low expression groups of ACE2 (B), TMPRSS2 (D), and AAK1 (F). (G–I) Correlation between the expression characteristics of ACE2 and TMPRSS2 (G), ACE2 and AAK1 (H), and AAK1 and TMPRSS2 (I). The “high” and “low” expression of ACE2 and TMPRSS2 are defined using the first and last quartile of the expression level.

WGCNA of DEGs Between LUAD Patients and Normal Subjects

To further explore to potential mechanism and signaling pathway of SARS-CoV-2 receptors including ACE2, AAK1, and TMPRSS2 in LUAD patients. WGCNA analysis was used to identify the co-expressive modules. A total of 15415 mRNA-sequencing profiles (535 LUAD samples and 59 non-tumor samples) were obtained from The Cancer Genome Atlas (TCGA) (up to May 15, 2020). In addition, 4674 DEGs between LUAD patient samples and normal samples were selected with an absolute log2-fold change (FC) > 1 and an adjusted P-value < 0.05. WGCNA (Langfelder and Horvath, 2012) of the DEGs was performed to identify highly synergistically changing gene sets and to explore the associations between these gene sets and ACE2, TMPRSS2, and AAK1 (Supplementary Figure 1). We found that the “MEbrown” and “MEblue” modules were significantly related to both TMPRSS2 and AAK1, while none of the modules had a clear correlation with ACE2 (Figure 5A). Then, we explored the correlations of the genes in the modules with TMPRSS2 and AAK1. These genes were determined to be not only highly related to their corresponding modules but also significantly related to TMPRSS2 and AAK1 (Figures 5B,C). GO analysis was further performed on these genes in the two modules. The results indicated that the genes in the “MEbrown” module were enriched in biological functions and pathways such as “cell-substrate adherens junction assembly,” “focal adhesion assembly,” “adherens junction assembly,” “cell-substrate junction assembly” and “adherens junction,” and the genes in the “MEblue” module were enriched in “vasculogenesis,” “endothelium development,” “regulation of body fluid levels,” “endothelial cell differentiation” and “cell-substrate adhesion” (Figures 5D,E).

FIGURE 5

Figure 5. WGCNA of DEGs between LUAD patients and normal subjects for exploring potential mechanism and signaling pathway of ACE2, AAK1, and TMPRSS2. (A) Correlation of co-expression modules with ACE2, TMPRSS2 and AAK1. (B) Genes in the blue module are closely associated with TMPRSS2. (C) Genes in the brown module are highly correlated with AAK1. (D) GO analysis to explore signal pathways associated with TMPRSS2. (E) GO analysis to explore signal pathways associated with AAK1.

miR-432-5p Is a Coregulated Factor of ACE2, TMPRSS2, and AAK1

MicroRNAs (miRNAs) are a class of small non-coding RNA molecules that regulate gene expression by targeting the mRNA 3′UTR (Fabian et al., 2010). To explore the potential inhibitor of SARS-CoV-2 receptors including ACE2, TMPRSS2, and AAK1, we performed bioinformatics analysis and used more than two different types of prediction algorithms to clarify that ACE2, TMPRSS2, and AAK1 were regulated by miR-432-5p (Figure 6A) (Supplementary Table 1). The targeting sites of miR-432-5p and the three genes are presented in Figures 6B–D. Then, we explored the signaling pathways with which miR-432-5p mainly interfered with physiological and pathological processes such as inflammation, oncogenesis and negative regulation via Gene Set Enrichment Analysis (GSEA) (Hänzelmann et al., 2013). The top 5 positively regulated pathways are shown in Figure 6E, and the top 5 negatively regulated pathways are shown in Figure 6F. The result suggested that miR-432-5p was potential agent of SARS-CoV-2 receptors and playing an important role in tumor progression in the LUAD patients.

FIGURE 6

Figure 6. Bioinformatics identified miR-432-5p as a co-regulatory factor of ACE2, TMPRSS2, and AAK1. (A) Venn diagram exhibits miR-432-5p was a co-regulatory factor of ACE2, TMPRSS2, and AAK1. (B–D) Targeting sites of miR-432-5p and ACE2 (B), TMPRSS2 (C), and AAK1 (D). (E,F) GSEA analysis presented top 5 positive regulatory pathways (E) and top 5 negatively regulatory pathways of miR-432-5p (F).

Bexarotene Is a Potential Target Drug of ACE2, TMPRSS2, and AAK1

In further exploring the clinical therapeutic value of targeting ACE2, TMPRSS2, and AAK1, we identified bexarotene as a potential targeted drug of all three genes in the GDSC database¹³ (Yang et al., 2013; Figure 7A and Supplementary Table 2–4). Bexarotene is a special retinoic acid X receptor (RXR) agonist, and RXRs play a key role in the proliferation, differentiation, apoptosis and angiogenesis of tumor cells through transcriptional regulation (Shen et al., 2018). The 3D structure of bexarotene is shown in Figure 7B. We used AutoDock software (version 3.6.1) to dock bexarotene with the three genes and found that bexarotene effectively docked with specific sites of ACE2 (Figures 7C,D). The docking analysis in AutoDock mainly showed hydrogen bonding between the gene and bexarotene. Bexarotene may bind to TMPRSS2 and AAK1 through ionic bonding, van der Waals forces or hydrophobic interactions, so AutoDock did not show the docking sites of bexarotene with these two genes. We further studied the specific mechanism of bexarotene in the Cancer Treatment Response Portal (CTRP) database¹³ and identified 87 genes that can be regulated by bexarotene. The correlation of these genes is shown in Figure 7E. GO analysis was performed to explore the signaling pathways in which the 87 genes were mainly enriched. The results demonstrated that these genes are mainly enriched in “cellular response to organic cyclic compound,” “response to acid chemical,” “signaling by nuclear receptors,” “response to reactive oxygen species,” and especially “response to lipopolysaccharide” signaling pathways (Figure 7F), indicating that bexarotene may participate in regulating inflammatory response. Figure 7G shows the core gene clusters within the 87 genes that have changed significantly during the development of LUAD.

FIGURE 7

Figure 7. Bexarotene is a potential target drug of ACE2, TMPRSS2 and AAK1. (A) Venn diagram shows a Bexarotene is a potential targeted drug for all the three genes. (B) 3D structure of Bexarotene. (C,D) Possible docking site between Bexarotene and ACE2. (E) Correlation between genes regulated by Bexarotene. (F) The mainly signal pathways enriched by 87 genes regulated by bexarotene. (G) PPI network analysis of proteins encoded by the 87 genes.

Bexarotene Treatment Inhibited LPS-Induced Inflammatory Response in the Macrophage

To further explore the role of bexarotene in the inflammatory response, bexarotene was used to treated LPS-induced macrophage such as RAW264.7 and peritoneal macrophage. Then qPCR was used to measure the level of pro-inflammatory cytokines including IL-1β, IL-6 and TNF-a. The result showed the mRNA level of IL-1β, IL-6, and TNF-a was significantly inhibited in LPS-induced macrophage with bexarotene treatment both RAW264.7 macrophage and peritoneal macrophage (Figures 8A–F). In addition, DCFH-DA kit was used to detect reactive oxygen species (ROS) production in the LPS-induced macrophage. The result indicated the ROS generation was obviously suppressed with bexarotene administration in the LPS-activated macrophage including RAW264.7 macrophage and peritoneal macrophage (Figures 8G,H). The result suggested bexarotene treatment played an role in important in anti-inflammatory response, and the anti-inflammatory effect of bexarotene in the COVID-19 patients will also be further validated in our future research.

FIGURE 8

Figure 8. Bexarotene inhibits LPS-induced macrophage inflammation. (A–C) Bexarotene significantly inhibited mRNA levels of IL-1β (A), IL-6 (B), and TNF-α (C) in LPS-induced RAW264.7 macrophages. (D–F) Bexarotene showed significant inhibition on mRNA levels of IL-1β (D), IL-6 (E), and TNF-α (F) in LPS-induced peritoneal macrophage. (G,H) The ROS generation was obviously suppressed with bexarotene administration in the LPS-activated RAW264.7 macrophage (G) and peritoneal macrophage (H).

Discussion

The global outbreak of COVID-19 caused by the new coronavirus SARS-CoV-2 has become the largest health threat in almost all countries, bringing an unprecedented burden on the global medical system (Sun et al., 2020). The natural and adaptive immunity of immune cells is activated after SARS-CoV-2 infects tissues to prevent further invasion of the virus. Most patients with COVID-19 only show mild to moderate infection symptoms, but some patients have a disorder between natural and adaptive immunity, which aggravates local tissue lesions and leads to excessive inflammation (Lee et al., 2020). In particular, the “cytokine storm” caused by the release of large amounts of cytokines and chemokines can lead to acute lung injury and respiratory distress syndrome, endangering the lives of patients (Chen et al., 2020; Xu et al., 2020).

In the current study, when we explored the signaling pathways of DEGs between COVID-19 and normal samples, we found that these DEGs were mainly enriched in “response to virus,” “lymphocyte activation,” “antiviral mechanism by IFN-stimulated genes” and “regulation of inflammatory response,” which were related to the physiological processes of inflammation and immunity, demonstrating that genes with altered expression in patients infected with SARS-CoV-2 were associated with the regulation of inflammation and immune response. Recently a study used a large number of bulk RNA sequences and single-cell sequencing data to show that the immune response of COVID-19 patients is driven by IFN and causes a large number of immune cell infiltration (Lee et al., 2020), which is also consistent with our results. And the subsequent PPI analysis results in our study also confirmed that chemokines related to inflammatory responses, such as CCR1, CCL19, and CXCL10, and the immune response-related gene FPR1 presented significant expression changes in the progression of COVID-19.

Previous studies confirmed that SARS-CoV-2 invades host cells through the ACE2 receptor (Rivellese and Prediletto, 2020), and ACE2 is expressed not only in the lung but also in the gastrointestinal tract, kidney, heart, blood vessels and other organs and tissues, causing SARS-CoV-2 to damage not only the lung but also the heart and kidneys, as well as physiological systems such as the immune system and blood system, causing death in patients due to multiple organ failure and acute respiratory distress syndrome (Bourgonje et al., 2020). TMPRSS2 has also been studied and demonstrated to cleave at the spike proteins S1/S2 and S2 of the virus to initiate the conformational change of the spike protein and promote the entrance of the virus into the host cell through membrane fusion or endocytosis (Stahlmann and Lode, 2020). The expression of ACE2 and TMPRSS is the key factor affecting the first step of SARS-CoV-2 infection. It is worth noting that previous studies revealed that many infected COVID-19 patients were patients with malignant tumors, suggesting that tumor patients are susceptible people who need important attention in the prevention and control of COVID-19 (Liang et al., 2020). In this study, we analyzed the expression characteristics of ACE2, TMPRSS and the endocytosis regulator AAK1 in LUAD patients and found that the expression levels of ACE2 and TMPRSS in LUAD patients were significantly higher than those in normal subjects, indicating that the high expression of these genes in LUAD patients may be one of the reasons they are more susceptible to SARS-CoV-2 infection. However, some studies found that ACE2 and TMPRSS2 will decrease during disease progression (Kong et al., 2020), which may be related to the massive invasion of viruses and the consumption of these receptors. And our analysis results also consistently showed that the TMPRSS level of COVID-19 patients is significantly lower than that of normal subjects.

In the subsequent analysis, the infiltration characteristics of immune cells in the tumor tissue of LUAD patients were revealed to be closely related to the expression levels of ACE2, TMPRSS, and AAK1. Patients with high expression of ACE2 and patients with high expression of TMPRSS had relatively identical infiltration changes of immune cells. Compared with patients with low expression of these genes, patients with high expression exhibited higher infiltration fractions of B cells, activated dendritic cells and CD4 T cells. CD4 T cells play a critical role in attenuating or suppressing overactive innate immune responses during viral infection to maintain self-tolerance and immune homeostasis (Shaw et al., 2013). The significantly increased infiltration of CD4 T cells and B cells in LUAD patients with high expression levels of ACE2 and TMPRSS2 may cause an imbalance in immune regulation, making patients more susceptible to SARS-CoV-2 invasion, and the inflammatory response after infection is aggravated. While in the advanced stage of COVID-19, the CD4 cells exhibit functional exhaustion and the absolute number of them decreased significantly (Qin et al., 2020; Tan et al., 2020), and the expression of ACE2 decreased as well.

The gene sets identified by the WGCNA of DEGs between LUAD patients and normal subjects with high correlations with ACE2, TMPRSS2, and AAK1 and highly synergistic changes had biological functions and signaling pathways mainly concentrated in “cell-substrate junction assembly,” “cell-substrate adhesion,” “vasculogenesis” and “endothelium development,” which were related to the development of tissue cells and cell-mediated immunity (Tozeren, 1990).

In further exploring potential current treatment options for ACE2, TMPRSS2 and AAK1, we identified miR-432-5p as a potential targeting molecule with docking sites and regulatory functions for all three genes. MiR-432-5p plays a key role in tumor progression by regulating gene expression at the posttranscriptional level. MiR-432-5p has been revealed to inhibit the progression of LUAD by inhibiting the cell cycle to suppress cell proliferation (Chen et al., 2016). Studies have also found that miR-432 can negatively regulate the production of MCP-1 and the migration of monocytes to inhibit the inflammatory response (Liu et al., 2017). The specific regulatory effect of miR-432-5p on the expression of ACE2, TMPRSS2, and AAK1 is not yet clear, and we will further explore the mechanism of miR-432-5p in future studies. In addition, pharmacological analysis suggested that bexarotene is a potential targeted drug of these three genes. Bexarotene is a synthetic high-affinity retinoid X receptor agonist that is currently used clinically for the treatment of T-cell lymphoma (Talpur et al., 2002). In preclinical mouse models, bexarotene has also been confirmed to inhibit the occurrence of lung tumors (Zhang et al., 2019). Some studies have demonstrated that bexarotene exerts anti-inflammatory effects by downregulating the expression of interleukin (IL)-6, IL-8, monocyte chemoattractant protein 1 (MCP-1) and high mobility group box-1 (Li et al., 2019). In addition, we found bexarotene treatment obviously inhibited LPS-activated inflammatory response and ROS generation in the macrophage. These results suggest bexarotene maybe regard as a promising therapeutic strategy for COVID-19. Next, the anti-inflammatory mechanism of bexarotene and its regulating effects on ACE2, TMPRSS2, and AAK1 will also be further explored and validated in our follow-up studies.

Inevitably, the small dataset of COVID-19 with only 2 patients and 2 normal subjects is indeed the limitation of our study, which may affect the credibility of the results. Regrettably, as of now, the GSE147507 that we analyzed in this study is the only publicly available database containing COVID-19 patient information. In the future, we will conduct more in-depth analysis after obtaining more information of COVID-19 patients. Besides, whether bexarotene plays a role in patients with COVID-19 requires a large amount of preclinical data support, which is also the work we need to further explore.

Data Availability Statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author/s.

Author Contributions

JJ conceived and designed the experiments and edited the manuscript. BT, JZ, and YC performed the experiments. YC, WY, and YW analyzed the data. CK, WC, and YZ contributed analysis tools. BT and JZ wrote the manuscript. All authors read and approved the final manuscript.

Funding

This study was supported by the National Natural Science Foundation of China (No. 81803778), The Key Research and Development Project of Zhejiang Province (No. 2018C03024), the Public Welfare Research Program of Zhejiang Province (Nos. LQ20H160056 and LGD19h160002), the Natural Science Foundation of Zhejiang Province (No. LYQ20H280003), and the Science and Technology Project of Lishui City (2020KJGG10 to YZ).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

The authors wish to thank Y. Y. Zhang from the General Hospital of PLA Rocket Force for his help.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fcell.2020.577032/full#supplementary-material

Supplementary Figure 1 | WGCNA for identifying co-expression modules. (A) Volcano Plot exhibited DEGs between LUAD patients and normal samples. (B) Identification of the optimal soft threshold. (C–E) Construction and visualization of co-expression modules.

Supplementary Figure 2 | The infiltration of immune cells of LUAD patients with different expression characteristics of ACE2, TMPRSS2 and AAK1. (A) The correlation between expression of ACE2, TMPRSS2 and AAK1 and the infiltration of immune cells. (B–D) The infiltration of immune cells in the high and low expression groups of ACE2 (B), TMPRSS2 (C), and AAK1 (D).

Supplementary Table 1 | miRNA-target predicted by miRWalk, miRnda, RNA22, Targetscan, or RNAhybrid (predicted = 2 algorithms).

Supplementary Table 2 | Sensitive drugs associated with ACE2 expression were predicted by GDSC database.

Supplementary Table 3 | Sensitive drugs associated with TMPRSS2 expression were predicted by GDSC database.

Supplementary Table 4 | Sensitive drugs associated with AAK1 expression were predicted by GDSC database.

Abbreviations

ACE2, angiotensinconverting enzyme 2; COVID-19, Novel coronavirus disease 2019; CTRP database, Cancer Treatment Response Portal database; DEGs, differentially expressed genes; LUAD, lung adenocarcinoma; TCGA database, The Cancer Genome Atlas database; WGCNA, weighted gene coexpression network analysis; WHO, World Health Organization.

Footnotes

1. ^ https://covid19.who.int/
2. ^ https://david.ncifcrf.gov
3. ^ https://thebiogrid.org/
4. ^ https://metascape.org/gp/index.html
5. ^ http://cistrome.dfci.harvard.edu/TIMER/
6. ^ http://mirwalk.umm.uni-heidelberg.de/
7. ^ http://www.miranda.org/
8. ^ https://cm.jefferson.edu/rna22/
9. ^ http://www.targetscan.org/vert_72/
10. ^ https://bibiserv.cebitec.uni-bielefeld.de/rnahybrid/
11. ^ https://www.cancerrxgene.org/
12. ^ http://www.cbioportal.org/
13. ^ http://portals.broadinstitute.org/ctrp/

References

Bader, G. D., and Hogue, C. W. (2003). An automated method for finding molecular complexes in large protein interaction networks. BMC Bioinform. 4:2. doi: 10.1186/1471-2105-4-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Berry, I., Soucy, J. R., Tuite, A., and Fisman, D. (2020). Open access epidemiologic data and an interactive dashboard to monitor the COVID-19 outbreak in Canada. Can. Med. Assoc. J. 192:E420. doi: 10.1503/cmaj.75262

PubMed Abstract | CrossRef Full Text | Google Scholar

Blanco-Melo, D., Nilsson-Payant, B. E., Liu, W. C., Uhl, S., Hoagland, D., and Møller, R. (2020). Imbalanced host response to SARS-CoV-2 drives development of COVID-19. Cell 181, 1036–1045. doi: 10.1016/j.cell.2020.04.026

PubMed Abstract | CrossRef Full Text | Google Scholar

Bourgonje, A. R., Abdulle, A. E., Timens, W., Hillebrands, J. L., Navis, G. J., Gordijn, S. J., et al. (2020). Angiotensin-converting enzyme 2 (ACE2), SARS-CoV-2 and the pathophysiology of coronavirus disease 2019 (COVID-19). J. Pathol. 251, 228–248. doi: 10.1002/path.5471

PubMed Abstract | CrossRef Full Text | Google Scholar

Cerami, E., Gao, J., Dogrusoz, U., Gross, B. E., Sumer, S. O., Aksoy, B. A., et al. (2012). The cBio cancer genomics portal: an open platform for exploring multidimensional cancer genomics data. Cancer Discov. 2, 401–404. doi: 10.1158/2159-8290.CD-12-0095

PubMed Abstract | CrossRef Full Text | Google Scholar

Chai, P., Yu, J., Ge, S., Jia, R., and Fan, X. (2020). Genetic alteration, RNA expression, and DNA methylation profiling of coronavirus disease 2019 (COVID-19) receptor ACE2 in malignancies: a pan-cancer analysis. J. Hematol. Oncol. 13:43. doi: 10.1186/s13045-020-00883-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, L., Kong, G., Zhang, C., Dong, H., Yang, C., Song, G., et al. (2016). MicroRNA-432 functions as a tumor suppressor gene through targeting E2F3 and AXL in lung adenocarcinoma. Oncotarget 7, 20041–20053. doi: 10.18632/oncotarget.7884

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, N., Zhou, M., Dong, X., Qu, J., Gong, F., Han, Y., et al. (2020). Epidemiological and clinical characteristics of 99 cases of 2019 novel coronavirus pneumonia in Wuhan, China: a descriptive study. Lancet 395, 507–513. doi: 10.1016/S0140-6736(20)30211-7

CrossRef Full Text | Google Scholar

Fabian, M. R., Sonenberg, N., and Filipowicz, W. (2010). Regulation of mRNA translation and stability by microRNAs. Annu. Rev. Biochem. 79, 351–379. doi: 10.1146/annurev-biochem-060308-103103

PubMed Abstract | CrossRef Full Text | Google Scholar

Gentles, A. J., Newman, A. M., Liu, C. L., Bratman, S. V., Feng, W., Kim, D., et al. (2015). The prognostic landscape of genes and infiltrating immune cells across human cancers. Nat. Med. 21, 938–945. doi: 10.1038/nm.3909

PubMed Abstract | CrossRef Full Text | Google Scholar

Gorbalenya, A. E., Baker, S. C., Baric, R. S., de Groot, R. J., Drosten, C., Gulyaeva, A. A., et al. (2020). Severe acute respiratory syndrome-related coronavirus: The species and its viruses – a statement of the Coronavirus study group. bioRxiv [Preprint]. doi: 10.1101/2020.02.07.937862

CrossRef Full Text | Google Scholar

Hänzelmann, S., Castelo, R., and Guinney, J. (2013). GSVA: gene set variation analysis for microarray and RNA-seq data. BMC Bioinform. 14:7. doi: 10.1186/1471-2105-14-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Hoffmann, M., Kleine-Weber, H., Schroeder, S., Krüger, N., Herrler, T., Erichsen, S., et al. (2020). SARS-CoV-2 cell entry depends on ACE2 and TMPRSS2 and is blocked by a clinically proven protease inhibitor. Cell 181, 271–280. doi: 10.1016/j.cell.2020.02.052

PubMed Abstract | CrossRef Full Text | Google Scholar

Huang, C., Wang, Y., Li, X., Ren, L., Zhao, J., Hu, Y., et al. (2020). Clinical features of patients infected with 2019 novel coronavirus in Wuhan, China. Lancet 395, 497–506. doi: 10.1016/S0140-6736(20)30183-5

CrossRef Full Text | Google Scholar

Kong, Q., Xiang, Z., Wu, Y., Gu, Y., Guo, J., and Geng, F. (2020). Analysis of the susceptibility of lung cancer patients to SARS-CoV-2 infection. Mol. Cancer 19:80. doi: 10.1186/s12943-020-01209-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Kupferschmidt, K., and Cohen, J. (2020). Race to find COVID-19 treatments accelerates. Science 367, 1412–1413. doi: 10.1126/science.367.6485.1412

PubMed Abstract | CrossRef Full Text | Google Scholar

Langfelder, P., and Horvath, S. (2008). WGCNA: an R package for weighted correlation network analysis. BMC Bioinform. 9:559. doi: 10.1186/1471-2105-9-559

PubMed Abstract | CrossRef Full Text | Google Scholar

Langfelder, P., and Horvath, S. (2012). Fast R functions for robust correlations and hierarchical clustering. J. Statist. Softw. 46:i11. doi: 10.18637/jss.v046.i11

CrossRef Full Text | Google Scholar

Lee, J. S., Park, S., Jeong, H. W., Ahn, J. Y., Choi, S. J., Lee, H., et al. (2020). Immunophenotyping of COVID-19 and influenza highlights the role of type I interferons in development of severe COVID-19. Sci. Immunol. 5: eabd1554. doi: 10.1126/sciimmunol.abd1554

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, T., Fan, J., Wang, B., Traugh, N., Chen, Q., Liu, J. S., et al. (2017). TIMER: a web server for comprehensive analysis of tumor-infiltrating immune cells. Cancer Res. 77, e108–e110. doi: 10.1158/0008-5472.CAN-17-0307

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, Y., Xing, Q., Wei, Y., Zhao, L., Zhang, P., Han, X., et al. (2019). Activation of RXR by bexarotene inhibits inflammatory conditions in human rheumatoid arthritis fibroblast-like synoviocytes. Intern. J. Mol. Med. 44, 1963– 1970. doi: 10.3892/ijmm.2019.4336

PubMed Abstract | CrossRef Full Text | Google Scholar

Liang, W., Guan, W., Chen, R., Wang, W., Li, J., Xu, K., et al. (2020). Cancer patients in SARS-CoV-2 infection: a nationwide analysis in China. Lancet Oncol. 21, 335–337. doi: 10.1016/S1470-2045(20)30096-6

CrossRef Full Text | Google Scholar

Liu, S. C., Chiu, C. P., Tsai, C. H., Hung, C. Y., Li, T. M., Wu, Y. C., et al. (2017). Soya-cerebroside, an extract of Cordyceps militaris, suppresses monocyte migration and prevents cartilage degradation in inflammatory animal models. Sci. Rep. 7:43205. doi: 10.1038/srep43205

PubMed Abstract | CrossRef Full Text | Google Scholar

Longbottom, E. R., Torrance, H. D., Owen, H. C., Fragkou, P. C., Hinds, C. J., Pearse, R. M., et al. (2016). Features of postoperative immune suppression are reversible with interferon gamma and independent of interleukin-6 pathways. Ann. Surg. 264, 370–377. doi: 10.1097/SLA.0000000000001484

PubMed Abstract | CrossRef Full Text | Google Scholar

Lu, R., Zhao, X., Li, J., Niu, P., Yang, B., Wu, H., et al. (2020). Genomic characterisation and epidemiology of 2019 novel coronavirus: implications for virus origins and receptor binding. Lancet 395, 565–574. doi: 10.1016/S0140-6736(20)30251-8

CrossRef Full Text | Google Scholar

Qin, C., Zhou, L., Hu, Z., Zhang, S., Yang, S., Tao, Y., et al. (2020). Dysregulation of immune response in patients with COVID-19 in Wuhan, China. Clin. Infect. Dis. 71, 762–768. doi: 10.1093/cid/ciaa248

PubMed Abstract | CrossRef Full Text | Google Scholar

Richardson, P., Griffin, I., Tucker, C., Smith, D., Oechsle, O., Phelan, A., et al. (2020). Baricitinib as potential treatment for 2019-nCoV acute respiratory disease. Lancet 395, e30–e31. doi: 10.1016/S0140-6736(20)30304-4

CrossRef Full Text | Google Scholar

Rivellese, F., and Prediletto, E. (2020). ACE2 at the centre of COVID-19 from paucisymptomatic infections to severe pneumonia. Autoimmun. Rev. 19:102536. doi: 10.1016/j.autrev.2020.102536

PubMed Abstract | CrossRef Full Text | Google Scholar

Rothan, H. A., and Byrareddy, S. N. (2020). The epidemiology and pathogenesis of coronavirus disease (COVID-19) outbreak. J. Autoimmun. 109:102433. doi: 10.1016/j.jaut.2020.102433

PubMed Abstract | CrossRef Full Text | Google Scholar

Shannon, P., Markiel, A., Ozier, O., Baliga, N. S., Wang, J. T., Ramage, D., et al. (2003). Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. 13, 2498–2504. doi: 10.1101/gr.1239303

PubMed Abstract | CrossRef Full Text | Google Scholar

Shaw, A. C., Goldstein, D. R., and Montgomery, R. R. (2013). Age-dependent dysregulation of innate immunity. Nat. Rev. Immunol. 13, 875–887. doi: 10.1038/nri3547

PubMed Abstract | CrossRef Full Text | Google Scholar

Shen, D., Yu, X., Wu, Y., Chen, Y., Li, G., Cheng, F., et al. (2018). Emerging roles of bexarotene in the prevention, treatment and anti-drug resistance of cancers. Expert Rev. Anticancer Ther. 18, 487–499. doi: 10.1080/14737140.2018.1449648

PubMed Abstract | CrossRef Full Text | Google Scholar

Stahlmann, R., and Lode, H. (2020). Medication for COVID-19-an overview of approaches currently under study. Deutsches Arzteblatt Intern. 117, 213–219. doi: 10.3238/arztebl.2020.0213

PubMed Abstract | CrossRef Full Text | Google Scholar

Stark, C., Breitkreutz, B. J., Reguly, T., Boucher, L., Breitkreutz, A., and Tyers, M. (2006). BioGRID: a general repository for interaction datasets. Nucleic Acids Res. 34, D535–D539. doi: 10.1093/nar/gkj109

PubMed Abstract | CrossRef Full Text | Google Scholar

Sun, J., He, W. T., Wang, L., Lai, A., Ji, X., Zhai, X., et al. (2020). COVID-19: epidemiology, evolution, and cross-disciplinary perspectives. Trends Mol. Med. 26, 483–495. doi: 10.1016/j.molmed.2020.02.008

PubMed Abstract | CrossRef Full Text | Google Scholar

Talpur, R., Ward, S., Apisarnthanarax, N., Breuer-Mcham, J., and Duvic, M. (2002). Optimizing bexarotene therapy for cutaneous T-cell lymphoma. J. Am. Acad. Dermatol. 47, 672–684. doi: 10.1067/mjd.2002.124607

PubMed Abstract | CrossRef Full Text | Google Scholar

Tan, M., Liu, Y., Zhou, R., Deng, X., Li, F., Liang, K., et al. (2020). Immunopathological characteristics of coronavirus disease 2019 cases in Guangzhou, China. Immunology 160, 261–268. doi: 10.1111/imm.13223

PubMed Abstract | CrossRef Full Text | Google Scholar

Tobaiqy, M., Qashqary, M., Al-Dahery, S., Mujallad, A., Hershan, A. A., and Kamal, M. A. (2020). Therapeutic management of COVID-19 patients: a systematic review. Infect. Prevent. Pract. 2:100061. doi: 10.1016/j.infpip.2020.100061

CrossRef Full Text | Google Scholar

Tozeren, A. (1990). Cell-cell, cell-substrate adhesion: theoretical and experimental considerations. J. Biomechan. Eng. 112, 311–318. doi: 10.1115/1.2891189

CrossRef Full Text | Google Scholar

Vaid, S., Cakan, C., and Bhandari, M. (2020). Using machine learning to estimate unobserved COVID-19 infections in North America. J. Bone Joint Surg. 102:e70. doi: 10.2106/JBJS.20.00715

CrossRef Full Text | Google Scholar

Wang, C., Horby, P. W., Hayden, F. G., and Gao, G. F. (2020). A novel coronavirus outbreak of global health concern. Lancet 395, 470–473. doi: 10.1016/S0140-6736(20)30185-9

CrossRef Full Text | Google Scholar

Xu, Z., Shi, L., Wang, Y., Zhang, J., Huang, L., Zhang, C., et al. (2020). Pathological findings of COVID-19 associated with acute respiratory distress syndrome. Lancet Respir. Med. 8, 420–422. doi: 10.1016/S2213-2600(20)30076-X

CrossRef Full Text | Google Scholar

Yan, R., Zhang, Y., Li, Y., Xia, L., Guo, Y., and Zhou, Q. (2020). Structural basis for the recognition of SARS-CoV-2 by full-length human ACE2. Science 367, 1444–1448. doi: 10.1126/science.abb2762

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, W., Soares, J., Greninger, P., Edelman, E. J., Lightfoot, H., Forbes, S., et al. (2013). Genomics of drug sensitivity in cancer (GDSC): a resource for therapeutic biomarker discovery in cancer cells. Nucleic Acids Res. 41, D955–D961. doi: 10.1093/nar/gks1111

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, Q., Lee, S. B., Chen, X., Stevenson, M. E., Pan, J., Xiong, D., et al. (2019). Optimized bexarotene aerosol formulation inhibits major subtypes of lung cancer in mice. Nano Lett. 19, 2231–2242. doi: 10.1021/acs.nanolett.8b04309

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhou, P., Yang, X. L., Wang, X. G., Hu, B., Zhang, L., Zhang, W., et al. (2020). A pneumonia outbreak associated with a new coronavirus of probable bat origin. Nature 579, 270–273. doi: 10.1038/s41586-020-2012-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhou, Y., Zhou, B., Pache, L., Chang, M., Khodabakhshi, A. H., Tanaseichuk, O., et al. (2019). Metascape provides a biologist-oriented resource for the analysis of systems-level datasets. Nat. Commun. 10:1523. doi: 10.1038/s41467-019-09234-6

PubMed Abstract | CrossRef Full Text | Google Scholar