Co-Expression Network Analysis Identified LTF in Association with Metastasis Risk and Prognosis in Clear Cell Renal Cell Carcinoma

Introduction

As one of the most common malignant tumors of the urinary system, renal cell carcinoma (RCC) accounts for about 5% of adult malignancies.¹ Clear cell RCC (ccRCC) has the highest mortality rate, accounting for 70% to 85% of RCC incidence rate.^2,3 Due to the lacking of early diagnostic indicators and strong resistance to chemoradiotherapy, the prognosis of ccRCC patients is still poor.⁴ Especially for the metastatic patients, 5-year survival was less 20%.^5,6 Therefore, novel biomarkers related to metastasis need to be found to provide decisions for early diagnosis and treatment.

In recent years, microarray and high-throughput sequencing technology have been widely applied.^7,8 However, most researches focus on the screening of differentially expressed genes (DEGs), ignoring the high degree of correlation between genes.⁹ Based on the relationship between genes, a co-expression network was constructed by the weighted gene co-expression network analysis (WGCNA). Then significant gene modules and hub genes associated with the clinical traits were identified. Therefore, it is possible that we can identify such biomarkers that can predict the metastasis and progression of ccRCC.

Materials and Methods

Data Collection

Gene expression matrix and ccRCC patients’ clinical data were downloaded from Gene Expression Omnibus (GEO) database. Dataset from GSE105288 was used to construct expression network and identify hub genes. The microarray platform was Illumina HumanMethylation450 BeadChip (HumanMethylation450_15017482), and included 35 samples of ccRCC and 9 samples of non-tumor tissues. Then, our results were further validated by RNA sequencing data downloaded from the Cancer Genome Atlas (TCGA) database. The flow chart for the study selection for this integrated analysis is shown in Supplementary Figure 1.

Data Preprocessing

According to the distance among different samples in Pearson correlation matrix, the quality of microarray was evaluated by sample clustering. The height cut of 90 was chosen, and 3 samples (GSM2825244, GSM2825254 and GSM2825274) were detected as outliers (Supplementary Figure 2).

Differentially Expressed Genes Screening

The “limma” package was used to screen the DEGs between the tumor and nontumor tissues of ccRCC. The false discovery rate (FDA) <0.05 and |log₂ fold change (FC)| > 1 were considered to be statistically different.

Co-Expression Network Construction

Those genes with expression variance was greater than the quartile of all variances were used to construct co-expression network. The pearson’s correlation matrices were calculated. Then, a weighted adjacency matrix was constructed. Next, the adjacency was transformed into topological overlap matrix (TOM).¹⁰ Average linkage hierarchical clustering was conducted according to the TOM-based dissimilarity measure with a minimum size of 30 for the resulted dendrogram.¹¹

Identification of Clinical Significant Modules

Two important parameters were used to identify clinical significant modules. Gene significance (GS) was defined as the log10 transformation of the P value, and module significance (MS) was defined as the average GS for all the genes in a module. Of all the modules selected, the module with the largest absolute MS was generally considered to be a module related to clinical characteristics. Finally, select modules that were highly relevant to certain clinical features for further analysis.

Identification of Hub Genes

Hub genes consist of highly interconnected stacks within modules and have proven to have important functions.¹² In this study, hub genes were defined as genes with high module membership (MM) (cor. Weighted >0.8).¹³ We defined hub genes as genes with high MM (cor. Weighted >0.8). Then, the protein–protein interaction (PPI) network was also constructed based on the STRING database (https://string-db.org/). In the PPI network, genes with Top 10 Hubba nodes ranked by Maximal Clique Centrality (MCC) also defined as hub genes. The common hub genes in both co-expression network and PPI network were regarded as “real” hub genes for further analyses.

Validation Set Database

When the hub gene is determined, TCGA expression data and clinical data were used as validation set. Based on TCGA data, GEPIA (http://gepia.cancer-pku.cn/) can visually analyze the expression level and clinical relevance. The Human Protein Atlas (HPA, https://www.proteinatlas.org/) provided the location and qualitative data of protein in ccRCC.

Cell Culture and Lentiviral Vector Construction

Normal human kidney derived HK-2, and human RCC cell lines 786-O, ACHN, 769-P and Caki-1 were purchased from the Cell Bank of Type Culture Collection (CBTCC, Chinese academy of science, Shanghai, China) and Cell Bank of Wuhan University. The above cell lines were grown in RPMI 1640 medium (Gibco, Carlsbad, CA, USA) containing 10% fetal bovine serum (FBS) (Gibco, Carlsbad, CA, USA) and 1% penicillin/streptomycin solution (Gibco, Carlsbad, CA, USA) in a humidified condition with 5% CO₂ at 37 °C. All cell lines were expanded in monolayers in dishes or flasks. At 75–80% confluence, the culture medium was exchanged to serum-free solution for 24 h before their further experiments. Lentivirus containing pLVX-Puro (Vector) and pLVX- lactotransferrin (LTF)-Puro obtained from Genechem (GeneChen, China) was established for overexpression of LTF.

Flow Cytometry Analysis

The cell apoptosis was measured by an Annexin V-fluorescein isothiocyanate (FITC)/Propidium Iodide (PI) apoptosis kit. Firstly, the cells were collected after 48 hours’ transfection. Next, the transfected cells were staining with 5 µL Annexin V-FITC and 5 µL PI in the dark for 30 minutes. And then, 400 µL binding buffer was added into cell suspension. Finally, we calculated the amounts of the early and late apoptosis rate through a flow cytometer (Beckman Counter, Inc) in one hour after the reaction.

Cell Proliferation Assay

Cell Counting Kit-8 (CCK-8, Japan) was purchased to detect the cell proliferation in different time according to the protocols. The transfected cells were plated in 96-well plates at 10⁴ cells per well with 100 µL of culture medium in humidified condition at 37 °C. Then, 10 µL CCK-8 solutions were added into each well and further incubated in humidified condition at 37 °C for 1 h. Finally, the cells were measured at the absorbance of 450 nm according to the indicated time point. Each experiment repeated at least three times.

Cell Migration and Invasion Assay

The cell migration and invasion experiments were performed to detect the cell adhesive and invasive ability. As for the cell migration experiment by wound healing assay, the cells were sowed in 6-well plates (10⁶/well) with serum-free culture medium. The scratches were made through a 200 µL plastic pipette tip. Forty-eight hours later, the migrated cells would be washed with phosphate buffer solution (PBS), fixed with 4% paraformaldehyde. The migration distance would be calculated under the light microscope.

As for the invasion assay by the 8-um chamber. Firstly, the chamber was pretreated with matrigel (BD Biosciences, USA). And then, the 6 ×10⁴ cells were sowed into the upper chamber without serum, while the culture medium with 10% FBS was added into the lower chamber. Twenty-four hours later, the migrated cells on the subface would be washed with PBS, fixed with 4% paraformaldehyde and stained with crystal violet solutions. The invasion ability would be determined by calculated the migrated cells on the subface through the light microscope.

Quantitative Real-Time PCR Assay (qRT-PCR)

Total RNA from cell lines were isolated by a trizol agent and reverse transcribed to synthesize complementary DNA (cDNA) by the RevertAid First Strand cDNA Synthesis Kit. We performed the qRT-PCR with SYBR Green Mix in the RT-PCR detection system (Bio-Rad, USA) based on the manufacture’s protocols. The relative mRNA expression of LTF was measured by 2^−∆∆Ct methods and the GAPDH was served as an internal control. The primers involved in our study are as follows:

LTF forward:5ʹ- AGTCTACGGGACCGAAAGACA −3ʹ

LTF reverse: 5ʹ- CAGACCTTGCAGTTCGTTCAG −3ʹ

GAPDH forward: 5ʹ- GGACCTGACCTGCCGTCTAG −3ʹ

GAPDH reverse: 5ʹ- GTAGCCCAGGATGCCCTTGA −3ʹ

Statistical Analysis

All statistical analysis was conducted with SPSS 21.0 and the results in our study were expressed as mean ± SD. The significance of the changes between the two group was determined by Student’s t-test, and the data were considered significant when P < 0.05.

Results

DEGs Screening

Under the threshold of FDA < 0.05 and |log₂ FC| > 1, a total of 453 DEGs (165 up-regulated and 288 down-regulated in ccRCC) were selected from 41 samples. The DEGs and clustering dendrogram of tumor samples, as well the clinical traits were shown in Figure 1.

Figure 1 The differentially expressed genes (DEGs) and clustering dendrogram of tumor samples, as well the clinical traits. (A) The volcano plot for all DEGs based on GSE105288. (B) The clustering was based on DEGs between ccRCC and normal. The red color represents positive tumor metastasis.

Co-Expression Network Construction and Key Modules Identification

Using “WGCNA” package in R, the genes those expression variance were greater than the quartile of all variances were grouped into modules. Here, the power of β = 6 (scale free R² = 0.87) was selected as the soft-thresholding to ensure a scale-free network (Figure 2). Fifteen modules are considered to be significantly correlated (Figure 3A). And only the correlation between the light-green module and metastasis risk was high (P = 0.02, R² = –0.4) (Figure 3B). Thus, the light-green module with metastasis risk was identified as the clinical significant module.

Figure 2 Determine soft-thresholding power in WGCNA. (A) The scale-free fit index for various soft-thresholding powers (β). (B) The mean connectivity for various soft-thresholding powers. (C) Histogram of connectivity distribution (β = 6). (D) Checking the scale free topology (β = 6).

Figure 3 Identifying modules associated with the clinical traits of ccRCC. (A) Dendrogram of all DEGs clustered based on a dissimilarity measure (1-TOM). (B) Heatmap of the correlation between module eigengenes and clinical traits of ccRCC.

Identification of Hub Genes for Metastasis Risk in the Light Green Module

Eighty-eight genes from the light green module were enriched for Gene Ontology (GO) and pathway analysis. The results of GO of the light green module were focused on cytoplasmic vesicle lumen and vesicle lumen in CC (Supplementary Figure 3A), organic acid binding in MF (Supplementary Figure 3B) and neutrophil degranulation, neutrophil activation involved in immune response, neutrophil activation and neutrophil mediated immunity in BP (Supplementary Figure 3C). KEGG pathways of the light green module were significantly enriched in Lysosome and Staphylococcus aureus infection related pathways (Supplementary Figure 3D).

The highly linked hub genes in the module act as a key function in biological processes.¹⁴ Therefore, in the light green module, 27 genes with high connectivity (weighted correlation coefficient > 0.8) were selected as the candidate hub genes for the metastasis risk of the module (Table 1). Furthermore, we also constructed a network of PPI (Figure 4A), and top 10 Hubba nodes were identified based on Maximal Clique Centrality (MCC) though CytoHubba plug-in (Figure 4B).

Table 1 Hub Genes in the Module Related with Metastasis Risk

Figure 4 The hub genes in the PPI network of genes in the light green module. (A) PPI network of genes in the light green module. The red color represents upregulated genes and blue color represents downregulated genes. (B) The hub genes in the PPI network.

Then, combining the results of PPI, WGCNA and DEGs, LTF gene was regarded as “real” hub genes for metastasis risk was identified. Specifically, in our research, “real” hub gene was the intersection of multiple gene subgroups. (1) We used WNCNA to obtain hub module related to metastasis. The key parameter was the soft-thresholding (β = 6, scale free R2 = 0.87). (2) Then PPI network (weighted correlation coefficient > 0.8) filtered to Top 10 Hubba nodes by MCC algorithm of CytoHubba plug-in. (3) Under the threshold of FDA < 0.05 and |log2 FC| > 1, a total of 453 DEGs (165 up-regulated and 288 down-regulated in ccRCC) were selected from 41 samples. Finally, after finding the intersection of top 10 hubba nodes and DEGs, we identified LTF gene as “real” hub gene for metastatic risk.

Hub Gene Validation

The expression of LTF in ccRCC tumor tissue was significantly lower than that in normal tissue, which was verified by TCGA RNA-sequencing data (Figure 5A) and immunohistochemistry results from HPA database (Figure 5B). Similarly, the expression of LTF was significantly decreased in ccRCC cell lines (786-O, ACHN, 769-P and Caki-1) compared with the normal human kidney cell line (HK-2) (Figure 5C). We also found a poorer survival rate in the samples with low expression levels of LTF (HR = 0.66, 95% CI: 0.49–0.89, P = 0.0067; Figure 6A). In addition, the LTF expression level of the sample was negatively correlated with tumor stage (P = 0.0385, Figure 6B), and patients with low expression of LTF gene were more likely to have distant metastasis (P = 0.038, Figure 6C).

Figure 5 Validation of hub gene. (A) Differential expression of LTF gene in ccRCC tissues and normal tissues based on GEPIA (http://gepia.cancer-pku.cn/). *P < 0.05 (B) Immunohistochemistry confirmed the differential expression of LFT in ccRCC tissues and normal tissues based on the Human Protein Atlas (HPA, https://www.proteinatlas.org/). (C) Quantitative RT-PCR analysis of LTF expression in human normal kidney cell line (HK-2) and ccRCC cell lines (786-O, ACHN, 769-P, Caki-1). Data are presented as the mean ± SD. ***P < 0.001 ccRCC cell lines vs HK-2 cell line.

Figure 6 The effect of LTF on the clinical features of ccRCC. (A) Kaplan-Meier plot of overall survival for LTF expression. (B) The correlation of LTF expression with the tumor stages of ccRCC based on the TCGA data in the GEPIA database. (C) The correlation of LTF expression with the tumor metastasis of ccRCC based on the TCGA data.

LTF Prevented Proliferation and Promoted Apoptosis in vitro

The low expression of LTF in ccRCC patients manifests as an antagonistic medical behavior, and we doubt whether LTF can play a tumor-suppressive role in tumors. Therefore, we further explored the role of LTF on cell proliferation and apoptosis. 786-O and ACHN were selected for in vitro RCC cell models. The above two cell lines were transfected with lentivirus to induce stable overexpression of LTF, and verified by qRT-PCR (Figure 7A). In addition, CCK-8 and flow cytometry were used to detect cell proliferation and apoptosis. The results of CCK-8 revealed that overexpression of LTF markedly decreased the proliferation of 786-O and ACHN cells in 72 h (Figure 7B and C). Moreover, the flow cytometry analysis demonstrated that apoptotic cells were dramatically up-regulated in LTF overexpressed group (Figure 7D-G).

Figure 7 LTF prevented proliferation and promoted apoptosis in vitro. (A) LTF overexpression efficiency was estimated by RT-PCR in 786-O and ACHN cell lines. ***P < 0.001 (B and C) CCK8 assay showed LTF up-regulation suppressed cell proliferation in 786-O and ACHN cell lines.**P < 0.01 (D-G) The apoptosis of 786-O and ACHN cells improved after LTF overexpression. *P < 0.05, **P < 0.01.

Overexpressed LTF Prevented Cell Migration and Invasion

As shown in Figure 8A-D, the migration capacity of 786-O and ACHN cells was inhibited in LTF-upregulated group compared with that in the control group by wound healing assay; In addition, the invasive capacity was repressed in LTF-upregulated group compared with that in the control group by the transwell assay (Figure 8E-H). These data indicated that LTF could suppress the capacity of migration and invasion in ccRCC cell lines.

Figure 8 LTF overexpression prevented cell migration and invasion. (A-D) Wound-healing assay presented up-regulation of LTF reduced migratory ability in 786-O and ACHN cell lines.**P < 0.01 (E-H) High-expression of LTF in 786-O and ACHN cells suppressed cell invasion discovered by transwell assay. **P < 0.01.

Discussion

As we all know, ccRCC is not sensitive to radiotherapy and chemotherapy, and the recurrence rate and 5-year survival rate after radical nephrectomy are 40% and 20%, respectively. Thus, novel research on the mechanisms of the tumorigenesis and metastasis of ccRCC to identify druggable and predictable biomarkers is urgent.

In the present study, we found that the LTF was strongly correlated with metastasis risk in ccRCC, and distinguished the localized and non-localized ccRCC. Then the expression of LTF is dramatic reduced and positively related with the overall survival of ccRCC patients. We further validated its effects in ccRCC cell lines. Firstly, we found that the LTF levels were decreased in ccRCC cell lines. Then, the lentivirus with overexpressed LTF was transfected into ccRCC cell lines to test the results of bioinformatic analysis. Functionally, LTF could inhibit the cellular capacity of proliferation, migration and invasion, and promoted apoptosis. All these results suggested that LTF might be served as a potential biomarker for the diagnosis and therapy in ccRCC.

What is LTF and what does it do? LTF, an 80-kDa basic glycoprotein, is a member of the transferrin family of iron binding protein, which widely expressed in milk, trachea, saliva, nasal secretions, tears, semen, pancreatic juice and neutrophils particles.¹⁵ It plays multiple roles in mammals, such as anti-inflammation, anti-bacterial, anti-fungus, anti-tumorigenesis and iron transportation. LTF could initiate immunological cascade to increase cytokines release (such as IL-18 and IFN- γ) and cell activity (such as natural killer cell and cytotoxic T lymphocyte). Previously, several studies^16,21 reported that LTF suppressed the growth and metastasis of human carcinoma cell, such as thyroid tumors, gastric carcinoma, breast carcinoma, prostatic carcinomas, parotid gland, colorectal carcinomas, gallbladder carcinomas, endometrial carcinomas. Besides, it has been proved that LTF significantly inhibited the VEGF-mediated response in terms of the microvessel spatial extension, overall vascularity and incidence of crossover in adult rats. However, the studies of LTF in ccRCC were very few. Kevin et al,²² suggested that LTF mediated the anemia often seen in association with RCC. Giuseppe et al,²³ found lower LTF intensity-distribution score in ccRCC than non-ccRCC tissues. As we know, LTF has a high affinity for iron, which has been considered an essential nutrient for cells that are dividing quickly such as tumor cells, participating in various metabolic cellular processed such as oxidative phosphorylation and RNA and DNA synthesis. Therefore, neoplastic RCC could take use of LTF in order to have a greater availability of iron for their turnover. Although biological functions of LTF had attracted more and more attention, how did it involved in antioncogenesis was still largely unclear.

In the current study, low expression of LTF gene increases the risk of metastasis in ccRCC. Previous studies have reported that LTF regulated the risk of metastasis of soft tissue sarcoma,¹⁶ oral squamous cell carcinoma¹⁷ and nasopharyngeal carcinoma (NPC).^19,24 LTF was considered as a candidate tumor susceptibility suppressor gene located at 3p21.31–21.2 that could negatively regulate the development and metastasis of NPC²⁵ by inducing cell cycle G1/S arrest, and modulating the mitogen-activated protein kinase (MAPK) and Akt signaling pathways.^19,25 However, there is not enough evidence and research to explain the relationship between LTF and the risk of metastasis in ccRCC, which still need further study.

In summary, we have demonstrated that LTF is an antioncogene in ccRCC, impairing the ability of cells to proliferate, migrate and invade. Loss of LTF was closely related to the malignant features of ccRCC, which further proved that LTF could be a prognostic biomarker in ccRCC.