FAM50A: A Novel Prognostic Marker That Modulates the Proliferation of Colorectal Cancer Cells

1. Introduction

Colorectal cancer (CRC) ranks as the third most common malignant tumor type globally, and the second leading cause of cancer-related death . The incidence and mortality rates of CRC have stabilized and declined in some Western countries because of nationwide screening programs and the widespread use of colonoscopy . However, the incidence of CRC remains high, resulting in a significant financial burden to society . Owing to continued economic progress and changes in dietary practices in developing countries, the number of new CRC cases is expected to increase to approximately 2.4 million globally by 2035 . Despite recent progress in medical technology and the advent of multidisciplinary and comprehensive treatments for CRC, the high mortality rate associated with this disease has not improved significantly . The 5-year overall survival (OS) rate is 90% in patients diagnosed with early-stage CRC, 70% in those with locally advanced CRC, and 15% in patients with metastatic CRC . The majority of patients with CRC are diagnosed with an advanced stage of disease. While treatments such as chemotherapy, molecular targeted therapy, and immunotherapy have increased survival rates, some patients develop initial insensitivity or resistance to these drugs over time, leading to primary and secondary drug resistance [7,8]. Thus, the development of drug resistance remains a significant clinical challenge in the treatment of CRC . The primary cause of drug resistance is typically abnormal metabolism, transport, or targeting of anti-tumor drugs . Cell death pathways, cytokine pathways, oncogenic signals, compensatory feedback loop signaling pathways, and the tumor immune microenvironment can contribute to various drug resistance mechanisms . Together with tumor recurrence and metastasis, drug resistance is therefore a major factor in the mortality rates related to CRC . A thorough understanding of the mechanisms underlying the occurrence and progression of CRC is essential for the discovery of novel diagnostic biomarkers and treatment approaches. This will enhance our understanding of CRC and facilitate the development of more efficient clinical diagnosis and treatment methods.

Tumor progression is characterized by gradual cellular changes over time . Tumor cells evolve into subtypes that have clear proliferation advantages. In vivo, they progressively lose the normal controls for cell proliferation and subsequently transition to a new phase of proliferation . As cell proliferation and DNA mutations increase simultaneously, the cells escape their normal regulatory mechanisms. The morphology and metabolism of tumor cells are also altered. In a review published in 2011, Hanahan and Weinberg further elaborated the defining features of tumors. Cancer cells possess the ability to evade apoptosis, overcome growth suppression, induce angiogenesis, and sustain proliferative signaling, while exhibiting immortalization and unlimited replication potential . Many scientists have since focused their research on tumor cell proliferation, resulting in a deeper understanding of this process [15,16].

The family of sequence similarity (FAM) refers to a group of genes that share similar sequences and are involved in various diseases . Owing to the widespread expression of FAM20C in various cancer types, it has been linked to tumorigenesis in humans . Furthermore, deletion of FAM13A following hypoxia was reported to decrease the proliferation and metastatic capabilities of tumor cells, suggesting that it may be a promising target for cancer therapy . Sequence similarity family 50 member A (FAM50A) and FAM50B are two distinct members of the sequence similarity 50 gene family. FAM50A is known as X chromosome-associated protein 5 (Xap5) . This protein has a nuclear localization sequence, and can act as a DNA-binding protein or transcription factor. FAM50A is cytogenetically positioned on human chromosome Xq28 and comprises five exons that encode 339 amino acids with a molecular weight of 40 kDa . It may serve as a splicing factor in the processing of RNA precursors, and may also be associated with the development of X-chromosome-related intellectual disability . However, the function of FAM50A in CRC is not well understood, and there is limited research on its expression and role in this tumor type.

The aim of this study was, therefore, to investigate FAM50A expression during the onset and progression of CRC. First, data from The Cancer Genome Atlas (TCGA) public database was extracted and analyzed to study the potential tumorigenic effects of FAM50A in CRC. Immunohistochemical (IHC) staining was used to quantify the expression levels of FAM50A protein in CRC samples. Additionally, the prognostic significance of FAM50A in CRC was evaluated by analyzing patient survival. Finally, given the importance of cell proliferation in tumor growth, several experimental methods (Cell Counting Kit-8 [CCK-8], 5-ethynyl-2′-deoxyuridine [EdU], and colony-formation assays) were used to investigate the impact of FAM50A on cell proliferation.

2. Materials and Methods

2.1. Extraction and Analysis of Data from Public Databases

Data from public databases, TCGA and GTEx, were used to compare the mRNA expression profile of FAM50A between pan-cancer tissues and corresponding normal tissues . Subsequently, RNA-seq data and corresponding clinical information from unpaired and paired samples of CRC in TCGA database were systematically collected and processed based on the research interests of our research . The expression level of FAM50A was statistically analyzed following log₂ transformation. The analysis and mapping of FAM50A expression in tumors was then conducted using “limma” and other R (v4.2.2) packages, with statistical significance set at P < 0.05 .

2.2. Immunohistochemical Analysis

Building upon a previous study , IHC analysis was utilized to confirm the differential expression levels of FAM50A in CRC samples. Following surgical resection, the tumor diameter was measured, signs of bleeding and necrosis were noted, and relevant tissue was collected for further analysis. The samples were fixed in 10% formalin, embedded, sectioned, stained with hematoxylin and eosin (HE), and observed under a light microscope. IHC analysis was carried out on suitable samples using a rabbit anti-FAM50A antibody (Bioss, Woburn, MA, USA; bs-8208R, 1:200) using a previously described procedure . The stained sections were photographed under a microscope at either 40× or 200× magnification. The expression level of FAM50A was quantified based on the intensity of staining and the proportion of positive cancer cells, as described in our previous study . This study was approved by the Ethics Committee of Bozhou Hospital, Anhui Medical University (No. 202223). Samples were collected from newly diagnosed patients with CRC who had not received any anti-tumor treatment. Informed consent was obtained from patients prior to specimen collection.

2.3. Clinical Correlation and Prognostic Analysis

The Wilcoxon test was utilized to investigate correlations between FAM50A expression levels and various clinicopathological characteristics, including T stage, N stage, clinicopathological stage, and age. Patients with CRC were stratified into high- and low-expression FAM50A groups according to the FAM50A expression level cut-off value, as determined by receiver operating characteristic (ROC) analysis. Data collection and analysis was carried out by two dedicated individuals. Following data entry and review, Kaplan–Meier analysis was employed to examine patient survival. Cox regression analysis (univariate and multivariate) was utilized to investigate the relationship between the FAM50A expression level and patient survival. A nomogram was utilized to examine the ability of FAM50A levels to determine CRC prognosis .

2.4. Cell Culture and Construction of Cell Lines With Stable FAM50A Knockdown and Overexpression Plasmids

CRC cell lines (SW480, HCT-8, HCT-116, and RKO) and human normal intestinal epithelial cell lines, FHC and HEK293T, were obtained from the American Type Culture Collection. All cells were cultured in Dulbecco’s modified Eagle medium with a high glucose concentration and 10% fetal bovine serum at 37°C and 5% CO₂ . The lentiviral vector plasmids pLKO.1-Scramble and pLKO.1-shFAM50A were generated from the lentiviral vector plasmid pLKO.1-Puro; moreover, pCDH-puro-control and pCDH-puro-FAM50A plasmids were designated sequentially . The target sequences for FAM50A knockdown were: sh-FAM50A-1, 5′-CCGGCGAGATCCTTCGGAAAGACTTGGATCCAAGTCTTTCCGAAGGATCTCGTTTTTG-3′ (forward) and 5′-AATTCAAAAA CGAGATCCTTCGGAAAGACTT GGATCC AAGTCTTTCCGAAGGATCTCG-3′ (reverse); sh-FAM50A-2, 5′-CCGG GAGCTGGTACGAGAAGAACAAGGATCCTTGTTCTTCTCGTACCAGCTCTTTTTG-3′ (forward) and 5′-AATTCAAAAAGAGCTGGTACGAGAAGAACAAGGATCCTTGTTCTTCTCGTACCAGCTC-3′ (reverse). To generate the FAM50A-overexpression system, a gene fragment was obtained using the following primers: 5′-GCGAATTC ATGGCTCAATACAAGGGCGC-3′ (forward) and 5′-GCGGATCC TCAGCGGATCGTGTACTTGTCC-3′ (reverse). HEK293T cells were used for virus packaging and to collect viruses .

2.5. Quantitative Reverse Transcription Polymerase Chain Reaction (qRT-PCR)

Total RNA was extracted from cells using a commercially available RNA extraction kit (MN-MS-RNA-250; Bioscience, Shanghai, China). The qRT-PCR assay was performed using a Prime Script RT and SYBR Premix Ex Taq kit (Vazyme, Nanjing, China) according to the manufacturer’s instructions. Gene expression levels were normalized to that of β-actin. The qRT-PCR results were analyzed to obtain the Ct value of gene amplification products, and the 2^-ΔΔCt method was used for further analysis. The specific primers used for qRT-PCR were as follows: FAM50A, 5′-TCCCTCACCATCACAGCTTC-3′ (forward) and 5′-TCACTGAGCAACCGCACATC-3′ (reverse) and β-actin, 5′-GGCACCCAGCACAATGAAGA-3′ (forward) and 5′-ACTCCTGCTTGCTGATCCAC-3′ (reverse).

2.6. Western Blotting

Cells were lysed using radioimmunoprecipitation assay buffer (P0013B; Beyotime Biotechnology, Shanghai, China) containing proteinase and phosphatase inhibitors. The protein concentration of the samples was measured using a bicinchoninic acid (BCA) kit (P0010; Beyotime Biotechnology, Shanghai, China) according to the manufacturer’s protocol. Equal amounts of protein sample were separated using 10% sodium dodecyl sulfate-polyacrylamide gel electrophoresis and then transferred to a nitrocellulose membrane. Subsequently, the membrane was blocked with 5% skim milk powder for 1 h at room temperature (20-25°C) and incubated with the following primary antibodies overnight at 4°C: anti-FAM50A (Proteintech, Rosemont, IL, USA; 19849-1-AP, 1:1,000) and anti-β-tubulin (Proteintech, 10068-1-AP, 1:5,000). The membrane was then incubated with the appropriate horseradish-peroxidase-conjugated secondary antibody (Beyotime Biotechnology, A0208 and A0216, 1:5,000) for 30 min at room temperature. Finally, the bands were visualized using an enhanced chemiluminescence system and quantified using Image J software (National Institutes of Health, Bethesda, MD, USA).

2.7. CCK-8 Assay

Cells were seeded (1,000 cells per well) in a 96-well plate, with three replicate wells. Following inoculation, they were incubated overnight in a controlled environment. Subsequently, 100 μL of CCK-8 solution was added to the cells and they were incubated for 1 h at five different time points (after seeding 0 [24h], 1 [24h], 2 [72h], 3 [96h] and 4 [120h ] days). The optical density at 450 nm was then determined using a microplate reader (Thermo Fisher Scientific, Waltham, MA, USA) .

2.8. EdU Assay

Cells were cultured in six-well plates and incubated with a final concentration of 10 μM EdU for 2 h. They were then immobilized prior to Click Additive Solution buffer, DNA staining, and image acquisition as described previously .

2.9. Colony-Formation Assay

For each experimental group, 1,000 cells per well were seeded into six-well plates and cultured for 10 days with careful monitoring and regular changes of the growth medium. The resulting colonies were then stained and photographed as described previously .

2.10. Statistical Analysis

Descriptive and statistical analyses were conducted using Statistical Package for the Social Sciences (SPSS) version 22.0 (IBM; Armonk, New York, USA) or R software (v4.2.2). The Wilcoxon rank-sum test was utilized to assess the differential expression levels of FAM50A between the two groups. Correlations between FAM50A expression levels and clinicopathological characteristics were determined using Spearman’s rank correlation coefficient (Spearman’s r). Kaplan–Meier analysis was used to examine patient survival, and groups were compared utilizing the log-rank test. One-way analysis of variance was used to compare multiple groups, and the Student’s t test was used for pairwise comparisons.

3. Results

3.1. FAM50A was Highly Expressed in CRC

We first extracted FAM50A expression data for 33 tumor types from TCGA database. Analysis revealed that FAM50A was significantly overexpressed in various tumors, including bladder urothelial carcinoma, colon cancer, rectal cancer, and other common tumor types (P < 0.05, Figure 1A and 1B). Since FAM50A was highly expressed in both colon and rectal cancer, and given that CRC is our long-term research focus, this gene was selected for further detailed investigation. Subsequent analysis revealed a notable upregulation of FAM50A expression levels in CRC tumor tissue than in the adjacent normal tissue. Figure 1C depicts the results for unpaired CRC samples, while Figure 1D depicts the results for paired CRC samples. Following the observation of high FAM50A expression levels in CRC from a public database, we confirmed its expression in CRC tissue using IHC staining. A total of 145 patients with CRC were studied, comprising 67 women and 78 men. The baseline characteristics of this CRC cohort are presented in Table 1. Consistent with results of the analysis of data from TCGA database, FAM50A expression levels were higher in tumor tissue than in the adjacent normal tissue (P < 0.05, Figure 1E and 1F). Based on the above results, we conducted a series of in-depth studies of FAM50A.

3.2. FAM50A Expression was Related to Clinical Pathological Features of CRC

Subsequently, we investigated correlations between FAM50A expression levels and tumor clinical characteristics. Using ROC curve analysis, 2.9 was determined as the cut-off value for distinguishing patients with CRC revealing a high (95 patients, 65.5%) or low (50 patients, 34.5%) level of FAM50A expression. Positive correlations were observed between FAM50A expression levels and several clinical characteristics, including N stage (N0 vs. N1–N2; P < 0.001; Figure 2H), M stage (M0 vs. M1; P < 0.05; Figure 2I), and TNM stage (I–II vs. III–IV; P < 0.001; Figure 2J). No statistically significant correlations were observed between other features (all P > 0.05; Figure 2A–2G). Further results are presented in Table 1. These findings suggest that FAM50A plays a crucial role in driving tumor progression and they support the hypothesis that FAM50A may act as a proto-oncogene in CRC.

3.3. The Prognostic Value of FAM50A in CRC

Patient prognosis is an important factor in clinical practice. Kaplan–Meier survival analysis was utilized to investigate the association between FAM50A expression level and prognosis in CRC. The mean OS duration after surgery was 69.50 months (95% confidence interval [CI]: 65.24–73.75; Figure 3A). Kaplan–Meier and univariate regression analyses revealed that an advanced T stage (III–IV vs. I–II; hazard ratio [HR]: 2.015, 95% CI: 1.197–3.390; P = 0.008; Figure 3H), N stage (N1–N2 vs. N0; HR: 3.839, 95% CI: 2.181–6.757; P < 0.001; Figure 3I), M stage (M1 vs. M0; HR: 3.788, 95% CI: 2.026–7.081; P < 0.001; Figure 3J), TNM stage (I–II vs. III-IV; HR: 4.581, 95% CI: 2.602–8.065; P < 0.001; Figure 3K), and FAM50A expression level (high vs. low; HR: 2.544, 95% CI: 1.504–4.303; P < 0.001; Figure 3L) were associated with poor OS. There was no statistically significant difference in other variables (Figure 3B-3G, all P > 0.05). The corresponding risk factors are presented in Table 2. Multivariate Cox regression analysis showed that advanced M stage, advanced TNM stage, and high FAM50A expression level were independent predictors of poor OS (all P < 0.05, Table 2). These findings indicated that FAM50A may be a prognostic biomarker in predicting the survival outcome of patients with CRC.

Next, forest plots were drawn to visually illustrate the prognostic impact of FAM50A. These revealed a significant association between high FAM50A expression levels and the poor survival of patients with CRC (Figure 4A and 4B, Table 2). Nomograms can employ multiple indicators to simultaneously diagnose or predict the onset or progression of disease . In the present study, a nomogram was constructed to display the prognostic ability of FAM50A levels in patients with CRC. Candidate prognostic variables were assessed using a proportional hazard plot. This revealed that T, N, M, and TNM stage, as well as FAM50A expression level, all fitted the inclusion criteria (Figure 5A–5E). A nomogram was subsequently constructed using these five variables and subsequently assessed using a calibration curve. A high level of FAM50A expression was associated with higher scores and with a decreased probability of patient survival (Figure 5F). The calibration curve also indicated that the nomogram was well-constructed, since the predicted values were in good agreement with the actual values (Figure 5G).

3.4. Construction of Cell Lines with Stable FAM50A Expression

Continuous proliferation signals and the evasion of growth inhibition are characteristic phenotypes of tumor cells. We, therefore, investigated the potential role of FAM50A in promoting tumor cell proliferation. The gene expression level of FAM50A was first analyzed in four CRC cell lines (SW480, HCT-8, HCT-116, and RKO) and in normal intestinal epithelial cells (FHC). Results from RT-qPCR and western blotting analyses revealed higher FAM50A expression levels in tumor cell lines than normal FHC cells, at both the mRNA (Figure 6A, upper panel) and protein levels (Figure 6A, lower panel). The SW480 and HCT-8 cell lines were used to construct knockdown cell lines as they revealed the highest expression levels. Furthermore, the expression level of FAM50A in RKO cells was relatively lower and therefore, these cells were used for overexpression. Cells were first transfected with lentivirus, followed by drug screening to obtain cells with stable knockdown or overexpression . The efficiencies of knockdown and overexpression were verified at the mRNA and protein level in both SW480 (Figure 6B), HCT-8 (Figure 6C), and RKO (Figure 7A) cells. These stably expressing cells were then used in experiments to study the cell phenotype, as described below.

3.5. FAM50A Promoted CRC Cell Proliferation

The efficiencies of knockdown were verified at the mRNA and protein level in SW480 (Figure 6B), and HCT-8 (Figure 6C) cells. These stably expressing cells were then used in experiments to study the cell phenotype, as described below. The proliferation of CRC cells was evaluated using the CCK-8 method. Knockdown of FAM50A significantly reduced the proliferation of SW480 (Figure 6D) and HCT-8 (Figure 6E) cells. EdU staining was used to evaluate the effect of FAM50A knockdown on cell proliferation ability. Knockdown of FAM50A significantly reduced the proportion of EdU-positive FAM50A-sh-1 and FAM50A-sh-2 cells (Figure 6F). In addition, the number of cell clones was greatly reduced in cells with low FAM50A expression levels (Figure 6G). Moreover, overexpression of FAM50A was verified (Figure 7A), and significantly enhanced the proliferation ability of the RKO cell line (Figure 7B). The results of the EdU staining and colony-formation experiments demonstrated that overexpression of FAM50A significantly increased the number of EdU-positive cells (Figure 7C) and the extent of colony formation (Figure 7D). Overall, these results confirmed the role of FAM50A in cell proliferation.

3.6. FAM50A Expression in CRC was Highly Correlated with the Levels of Numerous other Genes

To investigate the mechanism by which FAM50A is involved in tumor progression, TCGA database was utilized to identify molecules that could potentially interact with FAM50A. The correlation coefficients for gene expression levels were calculated individually, with the cut-off value being < -0.4 for a negative correlation coefficient, and > 0.6 for a positive correlation coefficient. Twenty-one candidate genes were identified, and a co-expression heatmap of these genes was generated based on the level of FAM50A expression. The co-expression heatmap revealed a negative correlation between the expression level of FAM50A and those of RO60, CDC42SE2, and NRIP1, whereas the expression levels of all remaining genes revealed positive correlations with the expression level of FAM50A (Figure 8A). The correlation heatmap was redrawn, revealing that three genes with a negative correlation and 18 genes with a positive correlation were closely linked to the regulation of FAM50A (Figure 8B). Subsequently, correlation graphs between the expression levels of 21 genes and FAM50A were individually generated and analyzed. The results for the three genes with negative correlations were as follows: RO60 (r = -0.431, P < 0.001, Figure 8C), CDC42SE2 (r = -0.408, P < 0.001, Figure 8D), and NRIP1 (r = -0.404, P < 0.001, Figure 8E). The results for the 18 genes with positive correlation coefficients are depicted in Figure 8F–W (P < 0.001). These findings suggested that FAM50A may be involved in the malignant progression of CRC through interactions with the abovementioned gene products.

4. Discussion

Ongoing research in CRC has demonstrated that for effective diagnosis and treatment, potential targets must exhibit significant differential expression at the mRNA or protein level between tumor and corresponding normal tissue. Additionally, the target should demonstrate specific functional effects in terms of promoting or suppressing tumor growth . We observed that the FAM50A expression level was significantly elevated in tumor tissues compared to adjacent normal tissues, suggesting that it may serve as a potential biomarker in CRC. The clinical significance of genes such as p53, C-myc, Ras, EGFR, PMS1, PMS2, COX-2, CD44, PD-1, PD-L1, and CTLA4 in cancer has been extensively investigated. This includes the study of inactivated tumor suppressor genes (p53), activated proto-oncogenes (C-myc, Ras, and EGFR), mismatch repair genes (PMS1 and PMS2), modifying genes (COX-2 and CD44), and immune-related genes (PD-1, PD-L1, and CTLA4). These genes play crucial roles in genetic diagnosis, prognostic assessment, immunotherapy, and targeted therapy [38,39,40,41,42]. The results of the current study revealed significant correlations between FAM50A expression levels and the clinical features of N, M, and TNM stage, which is consistent with the results of a related study . Together, they suggest that FAM50A may be a valuable prognostic marker for clinical use. Multivariate Cox regression analysis demonstrated that the FAM50A expression level was an independent prognostic factor in CRC.

There is limited research on FAM50A, and its function has yet to be fully elucidated. The expression level of FAM50A varies across different tissues and tumors in the human body, as revealed by pan-cancer analysis of public databases. This variability may be attributed to tissue specificity, the involvement of diverse signaling regulatory mechanisms, and distinct epigenetic modifications . The results for FAM50A expression in the TCGA database were validated by IHC, thus highlighting its potential oncogenic activity. To further investigate its potential role as an oncogene, we constructed HCT-8 and SW480 cell lines with reduced FAM50A expression levels and confirmed the efficacy of knockdown using qRT-PCR and western blotting methods. Given the importance of tumor cell proliferation , our initial focus was to examine the effect of FAM50A on the proliferation of malignant CRC cells. CCK-8, EdU, and colony-forming assays were performed to study the effects of FAM50A knockdown, with the results consistently demonstrating decreased cell proliferation.

The mechanism of action of FAM50A in tumor cell proliferation was subsequently explored using co-expression heat mapping and correlation analysis. A total of 21 genes with expression levels revealing high correlations with FAM50A expression levels were included in the analysis. RO60 is associated with inflammatory diseases , and thus FAM50A may play a role in the progression of CRC by influencing specific inflammatory factors. NRIP1 is associated with the cAMP signaling pathway and facilitates the progression of lung cancer . It may, therefore, have a similar impact on the progression of CRC. IRAK1 mediates the down-regulation of AKR1B10 in lung cancer, and hence FAM50A may also regulate AKR1B10 through IRAK1 in CRC. OTUD5 is involved in the progression of numerous tumor types. For examples, it promotes the growth of hepatocellular carcinoma by stabilizing SLC38A1 . Therefore, OTUD5 may play an important role in the development of CRC. NSDHL is involved in regulating the activity of EGFR, thereby enhancing tumor growth and can possibly interact with FAM50A in CRC. Elevated expression levels of LAGE3 in CRC have been associated with poor prognosis and immune infiltration . FAM50A may, therefore, regulate LAGE3 expression in CRC. A previous study reported that UBL4A inhibited autophagy and increased the proliferation of pancreatic ductal adenocarcinoma cells by targeting LAMP1 . Based on this earlier study and the present results, we hypothesize that FAM50A may play a role in CRC by targeting LAMP1. Other genes that were closely associated with FAM50A in the present study have been implicated in tumorigenesis in previous reports [53,54].

While our findings demonstrated the prognostic value of FAM50A in CRC and a possible role in promoting tumor cell proliferation, several issues warrant further exploration. First, the limited sample size and pathology information of our CRC patient cohort may have affected the reliability and generalizability of our results. Samples with relevant patient information will be collected for analysis in future research endeavors. Additionally, the cell lines used for in vitro studies may not replicate the intricate and diverse nature of human CRC because of the lack of appropriate immune and physiological environments. Therefore, we will consider incorporating further in vitro testing in more comprehensive investigations in future studies. Financial constraints also restricted the scale of the experiments, the number of replicates, and the use of more sophisticated techniques. If various factors, such as funding and time, permit, our team Further in vivo research will conduct more experiments in vivo, which could may provide more precise results and conclusions. ImportantlyNotably, the specific mechanism involving whereby FAM50A in contributed to CRC progression has not been fully elucidated and requires more in-depth study.

In summary, our results show revealed that FAM50A may be a novel biomarker for CRC prognosis; and moreover, FAM50A may participate be involved in regulating tumor progression by interacting with some certain highly correlated molecules to that drive CRC cell proliferation.

5. Data Accessibility

All data are available from the authors for reasonable purposes and means.

Author contributions: Conceptualization, Longhai Li, Jizhao Shuai, Guangyun Li, Hao Gu and Yan Sun; Data curation, Hao Gu; Investigation, Jizhao Shuai; Methodology, Longhai Li; Resources, Jizhao Shuai; Software, Yan Sun; Supervision, Yan Sun; Writing-original draft, Yan Sun; Writing-review & editing, Longhai Li.

Declaration of competing interest: All authors declared that there were no financial interests or personal relationships that could influence the work reported in this paper.