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.

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.

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).

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.

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.

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.

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.

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).

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).

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.

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).

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).

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.

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.

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.