ACAT1 and Metabolism-Related Pathways Are Essential for the Progression of Clear Cell Renal Cell Carcinoma (ccRCC), as Determined by Co-expression Network Analysis

Kidney cancer ranks as one of the top 10 causes of cancer death; this cancer is difficult to detect, difficult to treat, and poorly understood. The most common subtype of kidney cancer is clear cell renal cell carcinoma (ccRCC) and its progression is influenced by complex gene interactions. Few clinical studies have investigated the molecular markers associated with the progression of ccRCC. In this study, we collected microarray profiles of 72 ccRCCs and matched normal samples to identify differentially expressed genes (DEGs). Then a weighted gene co-expression network analysis (WGCNA) was conducted to identify co-expressed gene modules. By relating all co-expressed modules to clinical features, we found that the brown module and Fuhrman grade had the highest correlation (r = −0.8, p = 1e-09). Thus, the brown module was regarded as a clinically significant module and subsequently analyzed. Functional annotation showed that the brown module focused on metabolism-related biological processes and pathways, such as fatty acid oxidation and amino acid metabolism. We then performed a protein-protein interaction (PPI) network to identify the hub nodes in the brown module. It is worth noting that only one candidate, acetyl-CoA acetyltransferase (ACAT1), was considered to be the final target most relevant to the Fuhrman grade of ccRCC, by applying the intersection of hub genes in the co-expressed network and the PPI network. ACAT1 was subsequently validated using another two external microarray datasets and the TCGA dataset. Intriguingly, validation results indicated that ACAT1 was negatively correlated with four grades of ccRCC, which was also consistent with our results from qRT-PCR analysis and immunohistochemistry staining of clinical samples. Overexpression of ACAT1 inhibited the proliferation and migration of human ccRCC cells in vitro. In addition, the Kaplan-Meier survival curve showed that patients with a lower expression of ACAT1 showed a significantly lower overall survival rate and disease-free survival rate, indicating that ACAT1 could act as a prognostic and recurrence/progression biomarker of ccRCC. In summary, we found and confirmed that ACAT1 might help to identify the progression of ccRCC, which might have important clinical implications for enhancing risk stratification, therapeutic decision, and prognosis prediction in ccRCC patients.

Kidney cancer ranks as one of the top 10 causes of cancer death; this cancer is difficult to detect, difficult to treat, and poorly understood. The most common subtype of kidney cancer is clear cell renal cell carcinoma (ccRCC) and its progression is influenced by complex gene interactions. Few clinical studies have investigated the molecular markers associated with the progression of ccRCC. In this study, we collected microarray profiles of 72 ccRCCs and matched normal samples to identify differentially expressed genes (DEGs). Then a weighted gene co-expression network analysis (WGCNA) was conducted to identify co-expressed gene modules. By relating all co-expressed modules to clinical features, we found that the brown module and Fuhrman grade had the highest correlation (r = −0.8, p = 1e-09). Thus, the brown module was regarded as a clinically significant module and subsequently analyzed. Functional annotation showed that the brown module focused on metabolism-related biological processes and pathways, such as fatty acid oxidation and amino acid metabolism. We then performed a protein-protein interaction (PPI) network to identify the hub nodes in the brown module. It is worth noting that only one candidate, acetyl-CoA acetyltransferase (ACAT1), was considered to be the final target most relevant to the Fuhrman grade of ccRCC, by applying the intersection of hub genes in the co-expressed network and the PPI network. ACAT1 was subsequently validated using another two external microarray datasets and the TCGA dataset. Intriguingly, validation results indicated that ACAT1 was negatively correlated with four grades of ccRCC, which was also consistent with our results from qRT-PCR analysis and immunohistochemistry staining of clinical samples. Overexpression of ACAT1 inhibited the proliferation and migration of human ccRCC cells in vitro. In addition, the Kaplan-Meier survival curve showed that patients with a lower expression of ACAT1 showed a significantly lower overall survival rate and disease-free survival rate, indicating that ACAT1 could act as a prognostic and recurrence/progression biomarker of ccRCC. In summary, we found and confirmed that ACAT1 might help to identify the progression of ccRCC, which might have important clinical implications for enhancing risk stratification, therapeutic decision, and prognosis prediction in ccRCC patients.

INTRODUCTION
Kidney cancer is one of the most common malignancies of the urinary system (1), ∼90% of which is renal cell carcinoma (RCC). Clear cell RCC (ccRCC) accounts for between 70 and 85% of RCC, and has the highest rate of mortality (2). In the past decades, kidney cancer patients had few treatment options other than surgery, and 5-year survival was <20% once a metastatic disease developed (3,4). A growing body of evidence shows a strong link between cancer and alerted metabolism. It is clear that many key oncogenic signaling pathways converge to accommodate tumor cell metabolism to support their growth and survival. In addition, some of these metabolic changes appear to be necessary for malignant transformation. In light of these basic findings, many researchers suggest that changes in cellular metabolism should be considered as an important marker of cancer (5). Previous studies have shown that seven known kidney cancer genes, VHL, MET, FLCN, TSC1, TSC2, FH, and SDH, are involved in pathways that respond to metabolic stress or nutrient stimulation, suggesting that kidney cancer is a disease of dysregulated cellular metabolism (6). There is, therefore, great significance in determining the effective metabolismrelated biomarkers responsible for the genesis and development of ccRCC.
Currently, microarray and high-throughput sequencing technology have been widely applied to screen biomarkers of cancer (7,8). The latest studies indicate that molecular biomarkers can improve the predictive accuracy of bladder cancer progression (9). Therefore, it is possible that we can identify such biomarkers that can predict the progression of renal cancer. Most of the studies focus on screening differentially expressed genes, but they have considerable limitations in ignoring the high correlations between genes. This finding may be functionally related between genes with similar expression patterns (10). The weighted gene co-expression network analysis (WGCNA) concentrates on the associations between genes, which is a systems biology method used to identify gene clusters associated with certain biological features (11). At present, similarly, this analysis is also used to identify tumor biomarkers, by correlating gene clusters with clinical features that can indicate tumor progression, such as tumor stage, grade, and metastasis (12,13). Thus, we intend to use the WGCNA method to find gene clusters related to ccRCC progression, which are indicated by clinical features, such as grade, stage, and metastasis. Some key metabolism-related genes and pathways can be identified from the gene clusters, which may illustrate the metabolic alteration during the progression of ccRCC (14,15).

Data Preparation
Raw data were downloaded from public Gene Expression Omnibus (GEO, http://www.ncbi.nlm.nih.gov/geo/). Differentially expressed genes (DEGs) were identified from 72 ccRCC and 72 normal kidney samples using dataset GSE53757 performed on Affymetrix HG U133 Plus 2.0 (16). Thirty-nine samples from dataset GSE29609 performed on Agilent Whole Human Genome Oligo Microarray G4112A (Agilent-012391) were used to perform weighted gene co-expression networks. Two additional external data GSE40355 and GSE73731 (17) were used for validation. To obtain reliable results, we also applied the Oncomine database (http://www.oncomine.org/) and GEPIA (18) (Gene Expression Profiling Interactive Analysis) database based on TCGA (The Cancer Genome Atlas) for validation.

Differentially Expressed Genes (DEGs)
Probes should be annotated first. We used "limma" (19) to screen the DEGs between ccRCC and normal samples. FDR (false discovery rate) <0.01 and |log 2 (FC)| ≥1 were regarded as the cut-off criteria.

Weighted Gene Co-expression Network
Detailed steps were described in our previous studies (20,21). Briefly, the "WGCNA" (11) R package was used to conduct a co-expression network (GSE29609). To ensure the reliability of the constructed network, outlier samples were excluded after plotting a clustering dendrogram using the "flushClust" R package. After excluding outlier samples, a weighted adjacency matrix was generated by the formula amn = |cmn| β (cmn represents Pearson's correlation between genes, amn represents adjacency between genes, the β parameter is based on the standard scale-free network and was used to magnify the correlation between genes). Then, a topological overlap matrix (TOM) was generated (22) and modules were identified by hierarchically clustering genes (23).

Clinically Significant Modules and Module Functional Annotation
After the co-expressed network was generated, the correlations between modules and external clinical information were calculated using the Pearson correlation method. To further clarify the underlying mechanism of module genes in corresponding clinical features, genes of the interest module were enriched using the "clusterProfiler" (24) R package for functional and pathway enrichment analysis. A false discovery rate (FDR) <0.01 was considered to be statistically significant.

Identifying Hub Genes
Based on the theory of WGCNA, hub genes had the highest degree of connectivity in a module, by which the biological significance of the module was determined. Connectivity and clinical trait relationship were measured by module connectivity and clinical trait relationship, as determined by the absolute value of the Pearson's correlation (cor.geneModuleMembership >0.8 and cor.geneTraitSignificance >0.2) (12). In addition, a proteinprotein interaction (PPI) network of the clinically significant module was constructed using STRING (25) (Search Tool for the Retrieval of Interacting Genes, https:// www.string-db.org/). The criterion for selecting hub nodes is a combined score of ≥0.8 and a connectivity degree of ≥20. To ensure the reliability of the identified hub genes, one possible strategy was to apply the intersection as the final target. Therefore, the "real" hub genes were taken as the intersection of the hub genes in the co-expression network and the hub nodes in the PPI network.

Hub Gene Validation
To assess the correlation of hub gene expression in four distinct Fuhrman grades, we conducted a linear regression analysis using two additional independent validation datasets GSE40435 and GSE73731. In addition, RNA-seq data were obtained from the GEPIA (Gene Expression Profiling Interactive Analysis, http:// gepia.cancer-pku.cn/) database to verify the association of hub gene expression with ccRCC progression. Survival plots for hub genes were also generated using the GEPIA database. Oncomine (https://www.oncomine.org) and the Human Protein Atlas database (26) (http://www.proteinatlas.org/) were used to verify mRNA and protein expression between tumor and normal samples.

Cell Culture and Transfection
The human ccRCC cell lines ACHN and Caki1 were purchased from the Chinese Academy of Sciences in Shanghai. ACHN and Caki1 were cultured in Minimum Essential Medium (Gibco, China) and McCoy's 5A Medium (Gibco, China), respectively, both containing 10% fetal bovine serum (FBS) (Gibco, Australia) in a humidified atmosphere with 5% CO 2 at 37 • C. Human ccRCC cells were transfected with lentiviral vector (Control) and lentiviral ACAT1 (ACAT1) using Lipofectamine 3000 (Invitrogen, USA). The lentiviral vector (Catalog# EX-EGFP-Lv105) and lentiviral ACAT1 (Catalog# EX-C1085-Lv105) were purchased from GeneCopoeia, USA. Cells were transfected with 1 µg of the plasmid and incubated for 48 h before harvesting.

Total RNA Isolation and qRT-PCR From Kidney Tissues
Total RNA from ccRCC and normal kidney tissues was isolated using a HiPure Total RNA Mini Kit (Magen, Shanghai, China), and quantified by NanoDrop. First-strand cDNA was synthesized using 1 µg of total RNA with a ReverTra Ace qPCR RT Kit (Toyobo, China). Each reaction was performed using 1 µg of the cDNAs with iQTM SYBR R Green Supermix (Bio-Rad, Shanghai, China) in a final volume of 20 µl. ACAT1 primer: The annealing temperature of the two primers was 60 • C.

MTT Assay and Cloning Formation Assay
For the MTT assay, 3,000-5,000/200 µl medium transfected cells were seeded into 96-well dishes. Then, 20 µl of 5 mg/ml MTT reagent (Sigma-Aldrich) was added to each well, followed by incubation at 37 • C for 4 h. After incubation, 100 µl of DMSO was added to each well to dissolve the precipitate. Absorbance at 570 nm was measured by a microplate reader (SpectraMax M2; Molecular Devices, USA). For the cloning formation assay, 1,000-1,500 transfected cells were seeded into six-well dishes to grow for 2 weeks. The cells were fixed with 4% PFA for 30 min, followed by staining with 0.1% crystal violet for 15 min.

Transwell Assay
A 24-well plate Transwell chamber system was used to perform the Transwell assay. Next, 5-10 × 10 4 transfected ccRCC cells in serum-free medium were seeded into the upper chamber (Corning, Inc. USA), while 10% FBS medium was added to the lower chamber. After incubation for 24 h at 37 • C, the cells in the upper chamber were removed. Cells on the lower side of the chamber were fixed with 4% PFA for 30 min, followed by staining with 0.1% crystal violet, and the cells were photographed by microscopy.

Western Blot Analysis
After transfection, cells were harvested and washed with cold PBS. Total protein was isolated using RIPA buffer with a protease inhibitor. A total of 20 µg of total protein was separated using 10% SDS-PAGE and transferred to PVDF membranes (Millipore, USA). After blocking with 5% bovine serum albumin (BSA), the membrane was incubated with primary antibodies overnight at 4 • C. The membrane was then incubated with the second antibodies for 2 h at room temperature. We visualized the band using an enhanced chemiluminescence (ECL) kit (Bio-Rad) with a Bio-Rad ChemiDoc MP Imaging System (Bio-Rad, USA). GAPDH was used as a loading control. The primary antibodies were as follows: ant-ACAT1, 1:1,000 (Proteintech, catalog# 16215-1AP), anti-GAPDH, 1:1,000 (Santa Cruz, catalog# sc-365062).

Immunofluorescence Staining and Evaluation for ccRCC Cells
Cells were seeded on coverslips after transfection. After washing with cold PBS, the cells were fixed with 4% PFA for 30 min and then treated with 0.1% Triton X-100 for 15 min. After blocking with 5% BSA for 30 min, the cells were incubated with Ki67 antibody (Novus, catalog# NBP2-19012) for 2 h at room temperature. After washing with PBS, the cells were incubated with Cy3-labeled secondary antibody for 1 h at room temperature. Immunofluorescence staining was visualized using a fluorescence microscope (Olympus, Japan) after nuclei were labeled with DAPI.

Immunohistochemistry (IHC) Staining for ccRCC
The tissue microarray (TMA) of ccRCC was collaborated with Shanghai OutdoBiotech (Shanghai, China). The tissue microarray (catalog# HKidE180Su02) contained 150 ccRCC specimens and 30 adjacent normal tissues. The survival and clinical correlation analyses were based on the detailed clinical data of these 150 cases. Briefly, paraffin sections were hydrated and embedded, followed by incubation with 3% H 2 O 2 for 15 min. Tissues were then incubated with citrate buffer for antigen retrieval. Tissues were incubated with ACAT1 antibody (Proteintech, catalog# 16215-1AP) after blocking with 5% BSA followed by incubation with biotinylated secondary antibody. Tissues were incubated with HRP substrate solution for 30 min, followed by incubation with DAB substrate chromogen solution. Tissue slides were counterstained with hematoxylin, dehydrated, and mounted. The ACAT1 staining signal was calculated based on the scores of the staining intensity and staining positive rate. The staining intensity is scored as follows: 0 points (negative), 1 point (weak), 2 points (moderate), and 3 points (strong). The staining positive rate is scored based on the positive cells as follows: 0 points (negative), 1 point (1-25%), 2 points (26-50%), 3 points (51-75%), and 4 points (76-100%). The ACAT1 total staining score is calculated by the formula: total score = staining intensity score × staining positive rate score. Total scores of ≥6 were regarded as a high expression and scores of <6 were regarded as a low expression.

Gene Set Enrichment Analysis (GSEA)
To explore the functional role of ACAT1 in the progression of ccRCC, we performed gene set enrichment analysis (27) (GSEA, http://software.broadinstitute.org/gsea/index.jsp) using GSE73731. Based on the median expression value of ACAT1, 265 ccRCC tissue samples were divided into two groups. The reference sets were selected as c2.cp.kegg.v5.2.symbols.gmt to annotate all gene sets. In addition, the cut-off criteria for significantly enriched KEGG pathways were a gene size of ≥30, FDR of <0.05, and of an enrichment score (ES) |> 0.65.

Differentially Expressed Genes in ccRCC Tissue Samples
Under the threshold of a false discovery rate (FDR) <0.01 and |log 2 (FC)| ≥1, a total of 2,572 DEGs were identified between 72 ccRCC samples and 72 normal kidney samples. The volcano plot for all DEGs is shown in Figure 1A. All DEGs were selected for subsequent co-expression network construction. The flow chart of the study was shown in Figure S1.

Weighted Co-expression Network
A co-expression network analysis was performed using 39 ccRCC samples in GSE29609 (Figure 1B). The 2,572 DEGs were included by adopting the "WGCNA" R package. β = 4 (scale free R 2 = 0.85) was selected as the soft-thresholding power  (Figures 2A-D). Ten modules were identified with a minimum size (gene group) of 30 for the gene dendrogram and a cut line of 0.25 for the module dendrogram ( Figure 3A).

Identification of Clinically Significant Module
Identifying the module most significantly associated with clinical features has considerable biological implications. The brown module has the highest correlation with the Fuhrman rank (r = −0.8, p = 1e-09, Figure 3B). In addition, the brown module also showed the highest gene significance associated with the Fuhrman grade ( Figure 3C). Therefore, we chose the brown module as the module of interest and analyzed it.
A total of 364 genes in the brown module were enriched for Gene Ontology (GO) and pathway analysis. Biological processes of the brown module were focused on fatty acid beta-oxidation, oxidation-reduction process and lipoprotein metabolic process (FDR < 0.01). KEGG pathways of the brown module were significantly enriched in protein and carbon metabolism-related pathways (FDR < 0.01, Figure 4A).

Identification of Hub Genes
Thirty genes with high connectivity (cor.geneModuleMembership >0.8 and cor.geneTraitSignificance >0.2) in the brown module were selected as hub genes. Table 1 lists the hub genes that significantly correlated with Fuhrman grade. To identify the "real" hub genes, a protein-protein interaction (PPI) network was also constructed. The PPI network was visualized by Cytoscape (28) (Figure 4B). Genes with more than 20 nodes were regarded as hub nodes. Intriguingly, only one common hub gene, ACAT1, in both the co-expression network and PPI network was identified as a "real" hub gene and was further validated ( Table 1).

Hub Gene Validation
Another three independent datasets were used to verify the expression of ACAT1 in different Fuhrman grades. In test set GSE40435 and GSE73731, linear regression analyses showed that ACAT1 expression had a strongly negative correlation with a Fuhrman grade of ccRCC (Figures 5A,B). The results from the clinical samples also suggested that ACAT1 mRNA and protein were both decreased in ccRCC (Figures 5D-F, 6A,B, S2A-D). In addition, the results from the qPCR and tissue microarray of ccRCC showed that ACAT1 expression at both the mRNA and protein level is lower in high stage ccRCC than in low stage ccRCC (Figures 5C,F, 6C,D,G, Table 2). Furthermore, ccRCC patients with lower ACAT1 expression had significantly shorter overall survival (OS) and disease-free survival (DFS) times ( Figure 6E and Figures S2E,F).

Overexpression of ACAT1 Inhibited Proliferation and Migration of ccRCC Cells
To investigate the effect of ACAT1 on the viability and proliferation of ccRCC cells, ACHN and Caki1 cells were treated with lentiviral control and lentiviral ACAT1 for 48 h and determined by MTT assay and cloning formation assay, suggesting that overexpression of ACAT1 drastically restrained ccRCC cell proliferation (Figures 7A-E). Immunofluorescence staining revealed that the ACAT1-overexpression group showed fewer Ki-67-positive cells than the control group ( Figure 7D).

Functional Annotation for the Hub Gene
To determine the functional role of ACAT1 in ccRCC progression, we performed a GSEA analysis using 265 ccRCC samples. Under the cut-off criteria of gene size ≥30, FDR <0.05, and |enrichment score (ES) |>0.75, six pathways were significantly enriched. Interestingly, the pathways were all involved in lipid and fatty acid metabolism processes, including "fatty acid metabolism, " "valine leucine and isoleucine degradation, " "PPAR signaling pathway, " "lysine degradation, " "butanoate metabolism, " and "citrate cycle TCA cycle" (Figure 8).

DISCUSSION
ccRCC is the most common subtype of kidney cancer and its progression is affected by complex gene interactions . Many studies suggest that altered cellular metabolism is an important marker of cancer (5). Previous studies have shown that kidney cancer is a disease of dysregulated cellular metabolism (5). Determining effective metabolismrelated biomarkers responsible for the genesis and development of ccRCC is strongly warranted. Molecular biomarkers associated with ccRCC progression have been predicted in many studies. By comparing ccRCC at different histological levels, eight genes were identified to distinguish between different levels of ccRCC (30). This analytical method might, however, result in large false positive results because it did not use a global level system biological analysis method. By directly comparing gene expression, other studies found that EphA1 (31), EphA2 (32), and VEGFR-1 (33) were indicators in different stages of ccRCC but lacked large sample support. Unlike these researches, we used a systems biology method combined with a large number of samples to screen specific biomarkers of ccRCC. Most importantly, molecular biomarkers should be well-differentiated between tumor and normal tissues. On this basis, the co-expression network was performed by the dynamic tree cutting method, and 10 co-expression modules were identified. Correlation analysis showed that the brown module had the highest correlation with the Fuhrman grade among the 10 modules. Hub genes with the highest connectivity were identified from the brown module, which determined the characteristics of the module to a large extent. In addition, a PPI network was constructed based on the genes in the brown module, and hub nodes were also identified. Finally, only one common hub gene, ACAT1, was selected as the real hub gene in the co-expression network and PPI network for further validation. Further verification also confirmed that ACAT1 was negatively correlated with the Fuhrman grade of ccRCCs, and its expression was also related to overall survival and the disease-free survival of ccRCC patients.
The ACAT1 enzyme carries out the last step in ketone breakdown during fasting or starving. The enzyme could catalyze the reversible formation of acetoacetyl-CoA from two molecules of acetyl-CoA (34). Ketone bodies are an important source of energy during fasting in normal cells. ACAT1 gene mutation induced a deficiency in mitochondrial acetoacetyl-CoA thiolase, which was also called ketothiolase deficiency (35). The most important fuel that supports tumor growth is carbohydrates. It therefore seems reasonable that a low-carbohydrate diet could reduce the progression of cancer. Many mouse studies have shown that dietary restriction reduced tumor size and growth rate (36), prolonged survival (37), and increased sensitivity to radiotherapy (36). Therefore, we can assume that the decreased ACAT1 expression in high-grade ccRCC may be caused by metabolic changes, as the invasive tumors cannot obtain enough energy from ketolysis and fatty acid oxidation to support their growth.
Compared with benign prostate tissues, ACAT1 expression was significantly increased in prostate cancer tissues (38,39). Upregulation of ACAT1 was correlated with more aggressive pancreatic cancer. Interestingly, Zhao et al. (40) and White et al. (41) both performed quantitative proteomic analysis and identified that the ACAT1 protein was decreased in ccRCC compared with adjacent tissues, which was consistent with our findings.
Previous studies reported that altered pathways, such as metabolic pathways (42), glycolysis, and fatty acid oxidation (43), were confirmed by our results from the brown module, revealing that metabolism changes were important for ccRCC progression. Lipid droplets were often found in ccRCC cytoplasm (4); therefore, ACAT1 may participate in cholesterol metabolism, similar to cytosolic acetyl-CoA acetyltransferase 2 (ACAT2), which indicated that lipid and fatty acid metabolism were different between kidney cancer and prostate cancer, as well as pancreatic cancer. In our study, ACAT1 at the transcriptional and translational levels were significantly decreased in ccRCC tissues (Figures 5D, 6). Moreover, ACAT1 had a strong negative correlation with the four grades of ccRCC, indicating that ACAT1 was closely associated with the progression of ccRCC. Notably, when ccRCC patients had lower expressions of ACAT1, they exhibited a significantly shorter OS and DFS rate. Furthermore, an in vitro study indicated that overexpressed ACAT1 inhibited the proliferation and migration of renal cell ACHN and Caki1 cells, suggesting that ACAT1 might be a favorable prognostic marker in ccRCC (Figures 5E,F).
We should also consider some of the limitations of this study. Our findings should be validated using a larger number of clinical samples. In addition, the mechanisms governing the impact of ACAT1 on the progression of ccRCC should be elucidated by molecular biology experiments, which is our next research plan. In conclusion, we built a co-expression network and identified the ACAT1 related to the progression of ccRCC, which might have important clinical significance in improving risk stratification, treatment decision-making and prognosis prediction in ccRCC patients.

ETHICS STATEMENT
Written informed consent was obtained from all participants of this study. All participants were over 16 years old. The study using human kidney tissue samples for RNA isolation was approved by the Ethics Committee of Zhongnan Hospital of Wuhan University (approval number: 2015029) (44).

AUTHOR CONTRIBUTIONS
LC, TP, YL, YX, and XW conceived and designed the study. LC, TP, YL, and YX performed the analysis procedures. LC, FZ, GW, and YX analyzed the results. YL, KQ, and YX contributed analysis tools. LC and TP performed the experiments. LC, TP, YL, and YX contributed to the writing of the manuscript. All authors reviewed the manuscript.

FUNDING
This study was supported in part by grants from Hubei Province Health and Family Planning Scientific Research Project (Grant number WJ2017H0002) and National Natural Science Foundation of China (Grant number 81772730). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

ACKNOWLEDGMENTS
We acknowledge the Gene Expression Omnibus (GEO), Oncomine, The Human Protein Atlas, and GEPIA databases for free use.