Hub Long Noncoding RNAs with m6A Modification for Signatures and Prognostic Values in Kidney Renal Clear Cell Carcinoma

Background: N6-methyladenosine (m6A)–modified long noncoding RNAs (m6A-lncRNAs) have been proven to be involving in regulating tumorigenesis, invasion, and metastasis for a variety of tumors. The present study aimed to screen lncRNAs with m6A modification and investigate their biological signatures and prognostic values in kidney renal clear cell carcinoma (KIRC). Materials and Methods: lncRNA-seq, miRNA-seq, and mRNA-seq profiles of KIRC samples and the clinical characteristics of corresponding patients were downloaded from The Cancer Genome Atlas (TCGA). The R package “edgeR” was utilized to perform differentially expressed analysis on these profiles to gain DElncRNAs, DEmiRNAs, and DEmRNAs, respectively. The results of intersection of DElncRNAs and m6A-modified genes were analyzed by the weighted gene co-expression network analysis (WGCNA) to screen hub m6A-lncRNAs. Then, WGCNA was also used to construct an lncRNA-miRNA-mRNA (ceRNA) network. The Cox regression analysis was conducted on hub m6A-lncRNAs to construct the m6A-lncRNAs prognostic index (m6AlRsPI). Receiver operating characteristic (ROC) curve was used to assess the predictive ability of m6AlRsPI. The m6AlRsPI model was tested by internal and external cohorts. The molecular signatures and prognosis for hub m6A-lncRNAs and m6AlRsPI were analyzed. The expression level of hub m6A-lncRNAs in KIRC cell lines were quantified by qRT-PCR. Results: A total of 21 hub m6A-lncRNAs associated with tumor metastasis were identified in the light of WGCNA. The ceRNA network for 21 hub m6A-lncRNAs was developed. The Cox regression analysis was performed on the 21 hub m6A-lncRNAs, screening two m6A-lncRNAs regarded as independent prognostic risk factors. The m6AlRsPI was established based on the two m6A-lncRNAs as follows: (0.0006066 × expression level of LINC01820) + (0.0020769 × expression level of LINC02257). The cutoff of m6AlRsPI was 0.96. KM survival analysis for m6AlRsPI showed that the high m6AlRsPI group could contribute to higher mortality. The area under ROC curve for m6AlRsPI for predicting 3- and 5-year survival was 0.760 and 0.677, respectively, and the m6AlRsPI was also tested. The mutation and epithelial–mesenchymal transition (EMT) analysis for m6AlRsPI showed that the high m6AIRsPI group had more samples with gene mutation and had more likely caused EMT. Finally, gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis were performed for mRNAs interacted with the two m6A-lncRNAs, showing they were involved in the process of RNA splicing and regulation of the mRNA surveillance pathway. qRT-PCR analysis showed that the two m6A-lncRNAs were upregulated in KIRC. Conclusion: In the present study, hub m6A-lncRNAs were determined associated with metastasis in KIRC, and the ceRNA network demonstrated the potential carcinogenic regulatory pathway. Two m6A-lncRNAs associated with the overall survival were screened and m6AlRsPI was constructed and validated. Finally, the molecular signatures for m6AlRsPI and the two m6A-lncRNAs were analyzed to investigate the potential modulated processes in KIRC.


INTRODUCTION
Renal cell carcinoma (RCC) is the third most common malignancy in the urogenital system (Siegel et al., 2019), with approximately 270,000 new diagnosed cases annually and 116,000 deaths around the world (Ljungberg et al., 2011). Among all renal tumors, 75% are kidney renal clear cell carcinoma (KIRC) (Rini et al., 2009), characteristically of high morbidity and mortality, for which chemotherapy, radiotherapy, and some other nonoperative treatments are difficult to obtain ideal outcomes (Rini et al., 2009). Patients with KIRC benefit majorly from surgery. There are still about 30% patients suffering from the risk of metastasis (Hsieh et al., 2017). Due to the complexity of carcinogenic regulatory network and the diversity of regulatory molecular modification, inhibition of known key targets cannot produce a marked effect. Other key factors that restrain tumor progression have become the major means for management of patients with metastatic KIRC. Therefore, it is crucial to identify other key regulatory factors for KIRC in an effort to improve the clinical outcome.
The accumulation of epigenetic alteration has become a momentous factor to advance development for tumors formatted by gene mutations. The aberrantly expressed long noncoding RNAs (lncRNAs) in specific cancers have been recognized as significant epigenetic regulatory molecules and as novel biomarkers for early diagnosis and therapy (Bach and Lee, 2018;Camacho et al., 2018). LncRNAs, defined as noncoding protein RNAs with more than 200 nucleotides are implicated in numerous bioprocesses of tumor cells, such as supporting proliferative signaling, eluding immune destruction and surveillance, capacitating replicative immortality, inducting angiogenesis, and activating invasion and metastasis (Gutschner and Diederichs, 2012;Bach and Lee, 2018). It has been found that upregulated TRPM2-AS contributes to tumor cell proliferation through TRPM2 stabilized by TAF15 in colorectal cancer (Pan et al., 2020). Qin et al. (2019) demonstrated that lncRNA MIR155HG acted as a tumorpromoting factor through inhibition of miR-802 that is a tumor suppressor in pancreatic cancer, providing a potential diagnostic and therapeutic target. Also in ovarian cancer, lncRNA HOTTIP boosts the expression of IL-6 to lead the immune escape of tumor cells by upregulating PD-L1 expression in neutrophils (Shang et al., 2019). In addition, the occurrence of tumor resistant to chemotherapy can also be caused by dysregulated lncRNA H19 through modulation of the glutathione metabolism pathway in advanced ovarian cancer (Zheng et al., 2016). Therefore, the regulatory role of lncRNAs to tumor is not less than that of oncogene or tumor suppressor gene, elucidating the significance of the landscape of lncRNAs in tumors. N6-methyladenosine (m6A) refers to methylation at the N6 position of adenosine, which has been demonstrated as the most ubiquitous, abundant, and conserved internal dynamic posttranscriptional chemical modification within mRNA (Meyer and Jaffrey, 2014;Zhao et al., 2017). Functionally, m6A modulates mRNA maturation, splicing, translation, expression, and degradation and other major bioprocesses that have not been discovered yet (Meyer and Jaffrey, 2014;Zhao et al., 2017). It is not so surprising that accumulating evidence demonstrated that m6A was implicated in regulating hallmarks of tumors with tumorigenesis, proliferation, differentiation, invasion, and metastasis (Lin et al., 2016;Chen et al., 2019;Wang T. et al., 2020). The reversible processes of m6A modification were catalyzed by methyltransferases (termed as "writers"), demethylating enzymes (termed as "erasers"), and m6A-binding proteins (termed as "readers") (Meyer and Jaffrey, 2014). METTL3, classified as writers, inhibits SOX2 degradation for conducting the progression of colorectal carcinoma in an m6A-IGF2BP2-dependent manner  and was linked with poor prognosis of gastric cancer due to METTL3-mediated HDGF of methylation facilitating cancer progression (Wang Q. et al., 2020). ALKBH5, an m6A demethylase, is capable of attenuating WIF-1 RNA methylation and mediating the Wnt signaling pathway, to restrain the progression of pancreatic cancer, causing a favorable clinical outcome for patients (Tang et al., 2020). Readers are the factors that execute functions. YTHDF2, classified as readers, causes liver cancer metastasis by elevating the expression of OCT4 in an m6A RNA methylation manner . Therefore, m6A modification offers a novel insight into the interplay between mRNAs and cancers. m6A modification also occurs within lncRNAs, with the equally significant role in regulating bioprocesses of tumor (Yi et al., 2020). lncRNA NEAT1-1 with m6A modification is capable to regulate the Pol II ser2 phosphorylation and promotes the formation of the complex CYCLINL1/CDK19/NEAT1-1 contributing to bone metastasis of prostate cancer (Wen et al., 2020). In colorectal cancer, lncRNA LINRIS modulates the expression of MYC to influence the process of glycolysis, leading to the progression of cancer by stabilizing IGF2BP2 . However, the hub lncRNAs with m6A modification (m6A-lncRNAs) that are associated with the prognosis in KIRC remain unclear. The present study aims to identify m6A-lncRNAs in KIRC and investigate their biological signatures and prognostic values, providing promising targets for further research.

Acquisition of Data About KIRC
The lncRNA-seq profiles, miRNA-seq profiles, and mRNA-seq profiles of KIRC samples and corresponding clinical characteristics including survival time, status, age, gender, grade, and TNM staging were downloaded from The Cancer Genome Atlas (TCGA) (https://portal.gdc.cancer.gov/). External cohort (GSE40914) was gained from Gene Expression Omnibus (GEO) (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi) to validate the prognostic value of screened m6A-lncRNA. The identified genes with m6A modification were retrieved from the RMVar database (http://rmvar.renlab.org/index.html) (Luo et al., 2020). Data with missing or unknown values were deleted. Then, we randomly and equally divided the parental TCGA cohort into two subgroups. There was no discrepancy between them. A subgroup was regarded as the training dataset and another as the testing dataset.

Production of Gene Modules
In the training dataset, the R package "edgeR" was utilized to perform differentially expressed analysis on the lncRNA-seq data in R software (version 4.0.2), with the statistical significance of | log2 fold change (FC) | > 1 and the false discovery rate (FDR) < 0.05. After intersecting DElncRNAs with m6A-mediated genes, the weighted gene co-expression network analysis (WGCNA) was conducted to identify hub m6A-lncRNAs in KIRC. First, the Pearson's correlation coefficient of the pairwise m6A-lncRNAs was calculated for constructing the similarity matrix. The similarity matrix was transformed into a weighted adjacency matrix after a soft threshold of β was identified. Next, the adjacency matrix was transformed into topological overlap matrix (TOM) measuring the connectivity between genes. Then, in terms of dissimilarity measure (1-TOM), average linkage hierarchical clustering was performed to cluster the genes with similar expression profiles for producing gene modules. The dissimilarity of module eigengenes (MEs) was calculated.

Identification of Significant Modules Linked With Clinical Traits
Gene significance (GS) obtained from the log10 transformation of the p value in a linear regression between the gene expression and clinical traits was computed to assess the significance of each module. Then, the average GS within a module was called as module significance (MS), which can be used to analyze the association between a module and clinical characteristics. The module correlated to a clinical trait was the one with the largest MS among all selected modules. The correlation between the MEs and clinical characteristics including grade and TNM staging was also calculated to determine the relevant module associated with clinical information. These processes were called as module-trait relationship (MTR) analysis (Langfelder and Horvath, 2008).

Obtaining Hub m6A-lncRNAs
A network was constructed according to the edges between two m6A-lncRNAs with weight >0.5 after screening m6A-lncRNAs correlated with clinical modules (the red and turquoise modules). On the basis of the plug-in unit "cytoHubba" within the Cytoscape software (version 3.7.2), m6A-lncRNAs in the network with Maximal Clique Centrality (MCC) ≥ 2 were identified as the hub m6A-lncRNAs.
Construction of an lncRNA-miRNA-mRNA Network Based on Hub m6A-lncRNAs The R package "edgeR" was applied to conduct a differentially expressed analysis on the miRNA-seq and mRNA-seq of KIRC samples in R software. The selection criteria were the same as described previously. Then, weighted correlation network analysis (WGCNA) was performed on hub m6A-lncRNAs-DEmiRNAs profile to calculate the connectivity between molecules. The analysis process was described previously. An optimal soft threshold of β was equal to three. The edges between hub m6A-lncRNAs and DEmiRNAs with a threshold weight >0.3 were selected. Next, WGCNA was also performed on DEmiRNAs-DEmRNAs profile. The soft threshold of β was set to five. The edges between DEmiRNAs and DEmRNAs with a threshold weight >0.3 were opted. Finally, took the DEmiRNAs as a bridge, an lncRNA-miRNA-mRNA (ceRNA) network was constructed and visualized by the Cytoscape software (version 3.7.2). The Kaplan-Meier (KM) survival analysis along with the logrank test on the regulatory molecules in the ceRNA network was performed.

Development of m6A-lncRNAs Prognostic Index
The univariate Cox regression analysis was conducted on the hub m6A-lncRNAs to screen the m6A-lncRNAs associated with overall survival. Based on the results of the univariate analysis, the multivariate Cox regression analysis was utilized to identify the independent prognostic factors. The m6A-lncRNAs prognostic index was constructed with the hub m6A-lncRNAs of p < 0.05 in the multivariate Cox regression model through multiplying the expression values of them by their coefficient in the model and then adding them together. The patients were then divided into high-m6AlRsPI and low-m6AlRsPI groups based on the median m6AlRsPI value. After screening the hub m6A-lncRNAs correlated with the overall survival through the Cox regression analysis, we performed the KM survival analysis on the relevant m6A-lncRNAs with statistical significance and the m6AlRsPI to estimate their prognostic power. The performance of m6AlRsPI for predicting 3-and 5-year overall survival of KIRC patients was analyzed by receiver operating characteristic (ROC) curve.

Validation for m6AlRsPI
Internal and external cohorts were applied to validate the prognostic signature of m6AlRsPI. In the testing dataset, the analysis processes were described previously. The m6AlRsPI was constructed and the prognostic and predictive signatures for m6AlRsPI in the testing cohort were estimated. External cohort (GSE40914) was also utilized to verify prognostic values of screened m6A-lncRNAs.

Analysis of Gene Mutation in Different m6AlRsPI Subgroups
The mutation analysis for KIRC samples was investigated to display the discrepant gene mutation in different m6AlRsPI subgroups. The major mutant genes and variant classifications in different m6AlRsPI subgroups were our concern.

Association of m6AlRsPI Subgroups With Epithelial-Mesenchymal Transition
The correlation between different m6AlRsPI subgroups and epithelial-mesenchymal transition (EMT) was explored to show the prognostic power in different m6AlRsPI subgroups, by analyzing the expression of CDH1, VIM, SNAI1, SNAI2, and CHD2 in high m6AlRsPI and low m6AlRsPI groups. The expression of PDCD1 in different m6AlRsPI subgroups was also calculated.

STATISTICAL ANALYSIS
The continuous variables were compared by independent t-test. Univariate and multivariate prognostic analysis were performed in the Cox regression model. The survival analysis was performed by the Kaplan-Meier survival analysis with the log-rank test. The R package "survival ROC" was utilized to plot the time-dependent ROC curve. The threshold of statistical significance was p < 0.05.

Obtaining m6A-lncRNAs for KIRC
Differentially expressed analysis performed on the training dataset (271 tumors vs. 35 normal samples), the sum of DElncRNAs was 3,836, including 1,221 downregulated lncRNAs and 2,615 upregulated lncRNAs, as shown in Figure 1A. A total of 1,643 differentially expressed m6A-lncRNAs were identified following intersection of 3,836 DElncRNAs with 39,136 m6A-modified genes, as shown in Figure 1B. After clearing invalid data, the data of clinical characteristics are presented in Supplementary Table S1.

Identification of Hub m6A-lncRNAs Related to Clinical Traits
In the training cohort, WGCNA was conducted on 1,643 differentially expressed m6A-lncRNAs to obtain the hub m6A-lncRNAs in KIRC. The optimal soft threshold value was selected as 14 (scale free R 2 0.96) to ensure a scale-free network (Figures  2A,B). We further analyzed the reliability of the scale-free topology in the light of soft threshold equal to 14 ( Figure 2C). The logarithm log10 (k) of the node with connectivity K was negatively associated with the logarithm log10 [p(k)] of the probability of the node. According to the similar expression pattern, a total of eight modules were determined by average linkage clustering ( Figure 2D). Red and turquoise modules were found to have a correlation with tumor M stage by MTR analysis ( Figure 2E). There were 22 m6A-lncRNAs within red module and 25 m6A-lncRNAs within turquoise module, with a threshold weight >0.5. Then, m6A-lncRNAs within the two modules were put into the Cytoscape software, to construct a network ( Figure 2F) and to calculate the MCC value which is listed in Supplementary Table S2, and m6A-lncRNAs with MCC ≥2 were identified as the hub m6A-lncRNAs. Finally, the total of hub m6A-lncRNAs was equal to 21.
Construction of an lncRNA-miRNA-mRNA Network After performing differential expression analysis on the miRNA-seq and mRNA-seq for KIRC samples, a total of 173 DEmiRNAs and 722 DEmRNAs were obtained, respectively. First, WGCNA was conducted on the hub m6A-lncRNA-DEmiRNAs profile to determine 19 hub m6A-lncRNAs associated with 107 DEmiRNAs. Similarly, 77 DEmiRNAs correlated with 459 DEmRNAs were obtained by performing WGCNA on DEmiRNAs-DEmRNAs profile. Regarding the common DEmiRNAs on both profiles as the connection point, the lncRNA-miRNA-mRNA network (ceRNA network) comprising of two hub m6A-lncRNAs, eight DEmiRNAs, and 18 DEmRNAs was constructed ( Figure 3A). The potential regulatory pathways based on ceRNA network contain MIR663AHG/hsa-mir-141/SLC4A1 etc. and PRRT3-AS1/hsa-mir-141/SLC4A1 etc. Then, other lncRNA-miRNA correlation networks and miRNA-mRNA correlation networks are shown in Supplementary Figures S1A,B. The survival analysis on the molecules within the ceRNA network was conducted. However, the two hub m6A-lncRNAs (MIR663AHG and PRRT3-AS1) and miRNA (has-mir-141) had no impact on the overall survival of patients with KIRC ( Supplementary Figures S2A-C). In addition, four mRNAs (ATP6V1C2, CLCNKA, PVALB, and RHCG) were found to be associated with the prognosis for patients ( Figures 3B-E).

Establishment of m6A-lncRNAs Prognostic Index Based on Cox Regression Analysis
The Cox regression analysis was conducted on the 21 hub m6A-lncRNAs to screen the prognostic factors for patients with KIRC. The outcome of the univariate Cox regression analysis was illustrated in Figure 4A. The hub m6A-lncRNAs with statistical significance were included in the multivariate Cox regression analysis ( Figure 4B). It was found that LINC01820 (HR: 1.0006; 95% CI: 1.0003-1.0009; p < 0.001) and LINC02257 (HR: 1.0021; 95% CI: 1.0001-1.0041; p 0.0423) were closely correlated with the poor prognosis of patients with KIRC (risky hub m6A-lncRNAs). The m6A-lncRNAs prognostic index (m6AlRsPI) was established based on the two m6A-lncRNAs as follows: m6AlRsPI (0.0006066 × expression level of LINC01820) + (0.0020769 × expression level of LINC02257). The median value of m6AlRsPI was 0.96, as the cutoff to separate KIRC samples into the high and the low m6AlRsPI group. The Kaplan-Meier survival analysis for the two hub m6A-lncRNAs and m6AlRsPI in the training cohort were carried out to investigate the association with the overall survival. The results demonstrated that the high expression level of LINC02257 had an unfavorable overall survival as compared with low expression level of LINC02257 ( Figure 4D). However, there was no survival discrepancy between low expression level and high expression level of LINC01820 ( Figure 4C). The high m6AlRsPI group could contribute to higher mortality than that of the low m6AlRsPI group ( Figure 4E). The area under ROC curve for m6AlRsPI for predicting 3-and 5-years survival was 0.760 and 0.677, respectively Figure 4F, demonstrated a moderate performance for predictive prognosis.

Validation of the Reliability for the m6AlRsPI
For internal validation, the testing cohort was used to validate the m6AlRsPI model developed by the training cohort (Supplementary Table S1). The m6A-lncRNAs in the testing cohort were screened by WGCNA as described previously (Supplementary Figures S3A-D). The hub m6A-lncRNAs with p < 0.05 within the univariate Cox regression analysis were taken into multivariate Cox regression analysis ( Figures 5A,B) to identify hub m6A-lncRNAs correlated with prognosis in the testing cohort. Then, the m6AlRsPI model based on LINC02257 and LINC01820 was also constructed, and the KM survival analysis showed high m6AlRsPI group with poor prognosis ( Figure 5C). Moreover, in the testing cohort, the area under ROC curve for predicting 3-and 5-years survival was 0.676 and 0.741, respectively ( Figure 5D), with moderate predictive power, consistent with the result above. The KM survival analysis for LINC02257 within GEO cohort demonstrated a high expression level of LINC02257 which increased the mortality ( Figure 5E) and no relevant external dataset was found to validate the prognostic value of LINC01820.

Gene Mutation Analysis for m6AlRsPI Subgroups
The gene mutation was analyzed to gain a novel insight into the molecular nature in m6AlRsPI subgroups. The gene mutation frequency in the high m6AlRsPI group was higher than that of the low m6AlRsPI group (89.22 vs. 85.45%). The top 10 mutational genes in the high m6AlRsPI group contain VHL, PBRM1, SETD2, TTN, BAP1, ATM, MTOR, MUC16, ADGRV1, and DNAH2 ( Figure 6A). The mutation rate of more than 10% genes includes VHL, PBRM1, SETD2, TTN, and BAP1. The top 10 mutational genes in the low m6AlRsPI group comprise of VHL, PBRM1, TTN, SETD2, BAP1, DNAH9, HMCN1, KDM5C, MTOR, and MUC4 ( Figure 6B). VHL, PBRM1, and TTN are the genes with the mutation rate of more than 10% genes. The most common mutation type in both groups is missense-mutation, followed by frameshift deletions and nonsense mutation (Figures 6C,D).

Association of Epithelial-Mesenchymal Transition With m6AlRsPI Subgroups
As the hub m6A-lncRNAs are closely linked with the clinical trait of M stage, we made further analysis of the association of EMT with m6AlRsPI subgroups. E-cadherin is regarded as epithelial biomarker, which is coded by CDH1 gene. It was found that in the high m6AlRsPI group, the expression level of CDH1 gene decreased as compared with the low m6AlRsPI group ( Figure 7A). Mesenchymal biomarkers include N-cadherin, vimentin, SNAI1, and SNAI2. N-cadherin and vimentin are coded by CDH2 and VIM gene, respectively. The analysis results showed the expression level of VIM, SNAI1, and SNAI2 in the high m6AlRsPI group were higher than that in the low m6AlRsPI group (Figures 7B-D). There was no statistical difference for the expression of CDH2 between the two groups ( Figure 7E). In addition, the expression of PDCD1 between the two groups is statistically different ( Figure 7F).

GO and KEGG Enrichment Analysis for Proteins Interacted with the Two Hub m6A-lncRNAs
After looking for RNA-binding proteins (RBPs) in The Encyclopedia of RNA Interactomes (ENCORI), RBPDB database, and catRAPID database, a total of 153 mRNAs interacted with the two hub m6A-lncRNAs were identified for GO and KEGG enrichment analysis. The RBP-lncRNA interaction in ENCORI was supported by CLIP-seq Data, and in RBPDB and catRAPID omics was predicted by RNA sequencing. The outcome from the GO analysis showed they were involved in processes of RNA splicing and mRNA metabolism ( Figures 8A,B). KEGG analysis demonstrated they are available to modulate RNA transport and degradation and mRNA surveillance signaling pathway ( Figures 8C,D).

The Expression Levels of Hub m6A-lncRNAs in KIRC
The expression levels of hub m6A-lncRNAs in KIRC were quantified by qRT-PCR. Compared to the normal kidney cell line (HK-2), the expression level of LINC01820 was elevated in KIRC cell lines (ACHN and 769-P), with statistical significance ( Figure 9A). The expression level of LINC02257 was upregulated in 769-P and 786-O cell lines ( Figure 9B). These results were consistent with our analysis results based on TCGA KIRC cohort, showing LINC01820 and LINC02257 were significant molecules in modulating KIRC progression. Frontiers in Molecular Biosciences | www.frontiersin.org July 2021 | Volume 8 | Article 682471 9 DISCUSSION Metastatic KIRC is the major cause of poor clinical outcome and higher mortality for patients with kidney cancer. Although traditional targeted therapy and novel immunotherapy were utilized to treat patients with metastatic KIRC, no satisfactory lasting responses of drugs to tumors and 5-year overall survival rate for patients could not be obtained. For those patients who have lost the opportunity of surgery, inhibition of key regulatory targets for metastasis may be the most effective therapy (Thomas and Kabbinavar, 2015). Due to discovery of implicating in the diverse considerable bioprocesses for tumor, lncRNAs are regarded as a novel regulatory factor for investigating the mechanism of cancer metastasis (Gutschner and Diederichs, 2012). A variety of evidence has manifested the pivotal role of lncRNAs in regulating tumor metastasis. LncRNA MALAT1 has previously been regarded as a metastasis-promoting factor in various tumors, but it is a metastasis-suppressing factor in breast cancer through binding and inactivating the pro-metastatic transcription factor TEAD (Kim et al., 2018). LncRNA RP11-390F4.3 is induced by hypoxia/HIF-1α to facilitate EMT through modulation of multiple EMT-associated factors, leading to tumor metastasis . In addition, m6A has been identified as the most common chemical modification manner for various RNAs in eukaryotic cells, and it is involved in the progression of tumor through relative key regulatory factors with m6A modification (Wang Q. et al., 2020). The brain metastasis of lung cancer results from m6A-mediated matured miR-143-3p upregulating the expression of VASH1 to increase ubiquitylation of VEGFA . BATF2 with m6A modification suppresses gastric tumor cell metastasis by stabilizing p53 protein to inhibit phosphorylation of ERK, causing a favorable prognosis . Furthermore, the process of EMT correlated with tumor metastasis can be triggered by m6A-related molecules. Li et al. (2020) demonstrated the overexpression of METTL3 induced EMT of cancer cells through modulation of expression and secretion of TGFβ1. Hua et al. (2018) found METTL3 promoted EMT by upregulating the receptor tyrosine kinase AXL in ovarian cancer. The crucial regulatory role of m6A modification in tumor metastasis ought to be paid attention. Therefore, the hub m6A-lncRNAs are vital targets for diagnosis, surveillance, and therapy of metastatic tumor, including KIRC. In colorectal cancer, m6A methylation can induce the upregulation of lncRNA RP11 to modulate Siah1-Fbxo45/Zeb1 complex to promote the dissemination of tumor cells . Zheng et al. (2019) demonstrated lncRNA FAM225A is an upregulated oncogenic lncRNA in nasopharyngeal carcinoma and m6Amediated lncRNA FAM225A gives rise to tumor metastasis by sponging miR-590-3p and miR-1275 to upregulate ITGB3 and activate the FAK/PI3K/Akt signaling pathway. m6A-mediated LINC00470 facilitates gastric cancer cell distant metastasis by accelerating degradation of PTEN mRNA (Yan et al., 2020). In cervical cancer, lncRNA ZFAS1 with m6A modification conduces to unfavorable clinical outcome by suppressing miR-647 to lead to tumor metastasis . Currently, the significant biomarkers correlated with metastatic KIRC and relative regulatory mechanism have remained unknown. It is urgent to look for effective factors implicated in the modulation of metastasis to improve the overall survival of patients with KIRC. Additionally, hub m6A-lncRNAs can be regarded as extremely important factors for investigating the complex mechanism of metastasis as a result of abundant regulation functions.
Based on WGCNA, DElncRNAs with m6A modification were divided into modules with the correlation analysis. After calculating the correlation between modules and clinical characteristics, we found two modules that were linked with the clinical trait of M stage. Although they were weakly correlated, the hub m6A-lncRNAs with a potential regulatory role were of interest. Therefore, m6A-lncRNAs within the two modules were considered related with metastasis of KIRC and might be treated as potential major molecules for modulating tumor metastasis. Then, construction of a ceRNA network aims to investigate the potential regulation pathways of 21 hub m6A-lncRNAs. Cytoscape software was utilized to screen the hub m6A-lncRNAs, and the Cox regression analysis was performed to  identify hub m6A-lncRNAs correlated closely with the prognosis. However, prognosis-associated m6A-lncRNAs were not embodied in the ceRNA network. The m6A-lncRNAs within the ceRNA network were not correlated with overall survival by the KM survival analysis and Cox regression analysis. So it was shown that prognosis-associated m6A-lncRNAs may influence biological behaviors of KIRC through other manners (such as RNA-binding proteins manner) instead of lncRNA-miRNA-mRNA manner. Identified LINC02257 was shown to be associated with poor prognosis of patients with colorectal cancer based on bioinformatics analysis Huang et al., 2020), but it was not clear whether m6A modification involved in the regulation of colorectal progression. In consideration of its unfavorable role in colorectal cancer, LINC02257 is considered as a novel significant regulator in metastatic KIRC, particularly in m6A manner. There was no survival discrepancy between low expression and high expression of LINC01820 through KM survival analysis. We thought KM survival analysis on LINC01820 may be disturbed by other confounding factors, leading to LINC01820 was not related with prognosis, and multivariate Cox regression analysis could solve the problem. The result from multivariate Cox regression analysis showed that LINC01820 was an independent factor associated with prognosis. Based on the result, we preferred to believe LINC01820 was a prognosis-associated m6A-lncRNA. Currently, there were no articles reporting that LINC01820 is correlated with prognosis in other tumors. ROC curve illustrated a moderate predictive signature for the m6AlRsPI. The m6AlRsPI developed by the training cohort was tested by self-validation and the prognostic signature of LINC0257 was validated by external GEO cohort. These results supported m6AlRsPI was a robust and reasonable model for predicting prognosis. The analysis of biological signatures for m6AlRsPI based on two hub m6A-lncRNAs showed that the high m6AlRsPI group had more samples with gene mutation, more likely to cause epithelial-mesenchymal transition, and led to worse prognosis. As LINC01820 and LINC02257 did not regulate tumor behavior by lncRNA-miRNA-mRNA manner, RNA-binding proteins which connected with the two hub m6A-lncRNAs were investigated. Then, we found that modulated by the two hub m6A-lncRNAs, the mRNA surveillance signaling pathway, termed nonsensemediated mRNA decay (NMD), has been found to be associated with mutation in human cancer by NMD-triggering manner (Lindeboom et al., 2016). This signaling pathway may play an important role in tumorigenesis and development of KIRC. The over-expression of the two hub m6A-lncRNAs in KIRC cell lines validated by qRT-PCR further indicated that they may influence tumor behavior, contributing to poor prognosis. These findings demonstrated that the two hub m6A-lncRNAs and m6AlRsPI were crucial for patients with metastatic KIRC to predict the prognosis and investigate the complex metastatic regulatory mechanism. Accordingly, they were significantly promising targets for treatment of patients with metastatic KIRC.
However, some limitations in our study should be taken into consideration. First, WGCNA was only performed on differentially expressed m6A-lncRNAs, other lncRNAs were excluded, causing some m6A-lncRNAs loss. Second, the red and turquoise modules had a weak correlation with the M stage. Third, the prognostic values of LINC01820 and m6AlRsPI cannot be validated by the external cohort currently. Fourth, as the database on m6A gene is updating, the screening of m6A-lncRNAs in our study was based on the current m6A dataset. Fifth, downstream molecular mechanisms for the two hub m6A-lncRNAs are unknown. It is required to conduct relevant experiments on the two hub m6A-lncRNAs to elucidate regulatory signaling pathways for metastatic KIRC.

CONCLUSION
In the present study, hub m6A-lncRNAs in KIRC were screened based on WGCNA, the m6A-lncRNAs prognostic index (m6AlRsPI) was constructed and validated according to two hub m6A-lncRNAs linked with prognosis, and the biological signatures and prognostic values were investigated. The two hub m6A-lncRNAs can be looked upon as critical regulatory molecules in metastatic KIRC and as potential therapeutic targets.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/ Supplementary Material.