An autophagy-related gene expression signature for survival prediction in multiple cohorts of hepatocellular carcinoma patients

Prognostic signatures have been proposed as clinical tools to estimate prognosis in hepatocellular carcinoma (HCC), which is the second most common contributor to cancer-related death at present globally. Autophagy-related genes play a dynamic and fundamental role in HCC, but knowledge of their utility as prognostic markers is limited. Here, we facilitated univariate and multivariate Cox proportional hazards regression analyses to reveal that 3 autophagy-related genes (BIRC5, FOXO1 and SQSTM1) were closely related to the survival of HCC. Then, we generated a prognosis index (PI) for predicting overall survival (OS) based on the three genes, which was an independent prognostic indicator for the OS of HCC (HR = 1.930, 95% CI: 1.200–3.104, P = 0.007). The PI showed moderate performance for predicting the survival of HCC patients and its efficacy was validated by data from three microarrays (GSE10143, GSE10186 and GSE17856). Furthermore, we deeply mined the integrated large-scale datasets from public microarrays and immunohistochemistry to validate the overexpression of BIRC5 and SQSTM1 while down-regulated FOXO1 expression in HCC. Bioinformatic analysis offered the hypothesis that proliferative signals in high-risk HCC patients were disturbing and thereby facilitated inferior clinical outcomes. Collectively, the prognostic signature we proposed is a promising biomarker for monitoring outcome of HCC. Nevertheless, prospective experimental studies are needed to validate the clinical utility.


INTRODUCTION
Hepatocellular carcinoma (HCC), the predominant primary tumor of the liver, is the second most common contributor to cancer-related death at present globally [1][2][3]. The considerably high mortality rate of liver tumor seriously threatens humanity, where 40,710 of the estimated new liver cancer cases in the United States in 2017 [4]. Although surgical resection technology or transplantation rapidly improved, the 5-year survival rate of HCC patients remains relatively low [5,6]. The dismal clinical outcome of HCC is in part driven by the delay in diagnosis and the still-rudimentary prognosis monitoring [7]. Hepatocarcinogenesis is heterogeneous and contains a multi-step process, including genetic and epigenetic factors that forms its unique molecular fingerprint.
Research Paper www.oncotarget.com Therefore, extensive analyses for identifying reliable biomarkers with prognostic significance that target major oncogenes in HCC is imperative.
Autophagy is a fundamental cell-physiologic regulator to ensure intracellular quality control, similar to apoptosis, energy production and waste removal [8,9]. Given the prominent functions of autophagy, it is closely associated with diverse human pathologies, such as immune disorders [10,11], neurodegenerative diseases [12,13] and cancers [14,15]. Oftentimes, compared with their normal counterparts, tumor cells are metabolically reprogrammed when several growth-promoting pathways are disturbed to obtain sufficient energy or additional stimuli to master of their own destinies. Interestingly, autophagy probably is present in both cancer and in cancer prevention, as well as potentially contributing to the growth of cancer. Additionally, its roles during the course of various tumors progression, including HCC [16], are variable. However, subtle mechanisms of autophagy in HCC are perplexing and have not been fully understood until now [17]. For the dynamic and fundamental role of autophagy in HCC, the prognostic signatures and therapeutic targets of autophagyrelated genes could show promise for understanding the genetic control mechanisms of HCC and offering promising targets. Recently, several studies have been proposed to identify autophagy-related prognostic signature in pancreatic ductal adenocarcinoma [18] and gliomas [19]. However, potential autophagy-related prognostic biomarkers for HCC are still urgently needed.
Due to the recent advances in RNA sequencing (RNA-Seq) technology and several public databases, such as The Cancer Genome Atlas (TCGA, https:// cancergenome.nih.gov/) and Gene Expression Omnibus (GEO, http://www.ncbi.nlm.nih.gov/geo/), it is feasible to identify several genes that can predict the clinical outcomes of HCC patients based on the expression profiles. In the present study, we first extracted autophagyrelated genes from The Human Autophagy Database (HADb, http://www.autophagy.lu/index.html), which provides an informative and up-to-date list of human genes involved in autophagy directly or indirectly. [20] Additionally, we calculated and obtained a series of differentially expressed autophagy-related genes in HCC and elucidated the cellular and molecular characteristics, dynamic role and pathways of these genes in HCC. The highlight was that we placed special emphasis on the prognostic value of autophagy-related genes and constructed a specific prognosis index (PI) that accurately predicted the clinical outcomes of HCC patients. These results were evidently confirmed by comprehensive data sources. Through these means, we proposed the rationale for clinically available prognosis monitoring and new insights into the molecular mechanisms for HCC patients, with particular emphasis on autophagy.

Differentially expressed autophagy-related genes in HCC
A total of 234 genes involved directly or indirectly in autophagy were downloaded via the online database HADb. After extracting expression data of these 234 autophagy-related genes from TCGA, 33 differentially expressed genes in HCC were identified (Figures 1  and 2), which were used for further investigation on their prognostic value. The independent sample t-test demonstrated that 21 genes showed remarkably higher expressions than non-cancerous tissues, including BIRC5, CDKN2A, PEA15, TP73, HSP90AB1, RAB24, ITGA6, CLN3, HGS, DAPK2, BAX, BAK1, NRG2, SQSTM1, TMEM74, IKBKE, SPHK1, ITGB4, ITGA3, IRGM and NKX2-3. Conversely, 9 genes (FOS, FOXO1, DIRAS3, DLC1, NAMPT, JUN, CCL2, MYC and NRG1) had notably lower expressions than the adjacent tissues. (Figures 3 and 4). Due to the discrepancy in calculation, no statistical significance was observed in HSPB8, TP63 and GRID1. Subsequently, ROC analysis was conducted to further investigate the ability of these genes to distinguish cancerous tissues from the non-cancerous ones. The results showed that the AUC value of 24 genes was above 0.7 ( Figures 5-6), which confirmed that these genes excelled at differentiating between cancerous and non-cancerous tissues.

Gene-enrichment and functional annotation for differentially expressed autophagy genes in HCC
According to the gene functional annotation summarized by clusterProfiler package of R software, we obtained Gene Ontology (GO), Disease Ontology (DO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) annotations of differently expressed autophagy-related genes in HCC. In the biological process (BP) term of GO analysis, not surprisingly, we found these genes were significantly related to autophagy ( Figure 7A; Figure  8A-8B). In cellular component (CC), "autophagosome", "integrin complex" and "protein complex involved in cell adhesion" were significantly enriched by these genes ( Figure 7C; Figure 8C-8D). For molecular function (MF), "sequence-specific binding", "growth factor binding" and "transcription factor binding" were items for these genes to play important roles ( Figure 7E; Figure 8E-8F). Interestingly, DO analysis indicated that these genes were also mainly related to various types of tumors, such as neuroblastoma, colorectal and stomach cancers ( Figure  8B; Figure 9A-9B). Do analysis provided us clues that these genes could act as general oncogenes. We also found these genes enriched in several risk KEGG pathways www.oncotarget.com which are relevant to tumor initiation and progression ( Figure 7D; Figure 9C-9D).

Autophagy-related prognosis index and clinicopathological parameters
After we removed the patients without sufficient survival data, a final cohort of 371 HCC patients in TCGA was used for prognosis evaluation. The univariate Cox analysis revealed that 6 genes including BIRC5, FOXO1, DLC1, SQSTM1, BAK1 and IKBKE had a prognostic value for HCC ( Figure 10). Subsequently, multivariate Cox analysis was conducted and BIRC5, FOXO1 and SQSTM1 were screened as independent prognostic indicators for overall survival (OS) of HCC (Figure 11A-11C; Supplementary Table 1).  Red dots indicated autophagy-related genes which were high expression in HCC and blue for low expression. This volcano plot was drawn by the ggplot2 package of R language. www.oncotarget.com Furthermore, t-test analyses indicated that the three genes were differentially expressed in various clinicopathological parameters. As shown in Table 1, differential BIRC5 expression was found at different tumor status, pathological stages, histological grades, T stages and vascular tumor cell types. Differential expression of FOXO1 was observed in different histological grades and vascular tumor cell types. For SQSTM1, it showed different expression in terms of age, gender and vascular tumor cell type.
Then, we generated a risk score model PI for predicting OS based on the three genes ( Figure 12) using the formula: PI = (1.242 * expression level of BIRC5) + (-1.619 * expression level of FOXO1) + (2.214 * expression level of SQSTM1). Then, the PI significantly stratified each HCC patients into highrisk (n = 186) or low-risk (n = 185) groups in terms of the OS (HR = 2.179, 95% CI: 1.537-3.089, P < 0.001; Figure 11D). Multivariate analysis demonstrated a consistent HR of 1.930 (95% CI: 1.200-3.104, P =  Table 2), which verified that the PI composed of the three genes remained as an independent prognostic indicator for the OS of HCC patients. When the HCC patients were stratified by pathological stage, the PI remained a high prognostic value for stage I/II (HR = = 2.280, 95CI%:1.461-3.558, P < 0.001; Figure 13A) and stage III/IV (HR = 1.720, 95CI%:0.953-3.105, P = 0.072; Figure 13B). When restricted to HCC patients with histological grades, a higher PI was still relevant with patients in G1/G2 (HR = 2.121, 95CI%:1.333-3.377, P = 0.002; Figure 13C) and G3/G4 (HR = 2.269, 95CI%:1.273-4.044, P = 0.006; Figure 13D). We also assessed the expression pattern of the three genes between the high-and low-risk groups. The results of data analysis showed that a remarkably higher expression for BIRC5 and SQSTM1 in the high-risk groups, while they had lower expression for FOXO1 in the high-risk groups ( Figure 14).

Validation of PI as an effective prognostic factor
The online bioinformatics tool SurvExpress provides the differences in mRNA expression levels to draw Kaplan-Meier curves and to classify risk groups. Groups with low and high expression are represented by green and red, respectively. SurvExpress contained two datasets with OS data (GSE10143, N = 162 and GSE10186, N = 112) and 1 dataset with recurrence free survival (RFS; GSE17856, N = 95) data. We first exhibited the expression pattern of these 3 genes ( Figure 15A, 16A, 17A) and clinical information of patients included in datasets ( Figure 15B, 16B, 17B). And then, similar to from the strategy of TCGA data dealing, patients in the three datasets from GEO were stratified into high-or lowrisk groups by PI signature (Figures 15C-15D, 16C-16D, 17C-17D). For GSE10143 cohort, patients in high-risk group had inferior OS than patients in low-risk group  Figure  15E-15F). For GSE10186 cohort, patients in high-risk group also had inferior OS than patients in low-risk group (HR = 2.07, 95% CI = 0.99-4.32, P = 0.053; Figure  16E-16F). In the GSE17856 cohort, PI could function as an indicator for the RFS of HCC patients. (HR = 2.12, 95% CI = 1.01-4.41, P = 0.046; Figure 17E-17F). These findings provided a valuable message that combination of these three genes could act as an indicator of prognosis monitoring with moderate degree.

Expression profiling of identified prognostic genes
After searching and selecting HCC-related datasets in the GEO, a total of 34 GEO datasets were used to identify the expression pattern of BIRC5, FOXO1 and SQSTM1. The characteristics of included datasets were shown in Table 3. After integrating the TCGA data and these datasets by using Standardized Mean Difference (SMD), we found that the BIRC5 expression level was significantly increased in the tumor group (SMD = 1.52, 95% CI: 1.23-1.83, P < 0.001; Figure 18) when using the random-effects model. And FOXO1 was validated to be lower in HCC (SMD = -1.04, 95% CI: -1.32-(-0.76), P < 0.001; Figure 19) when using the random-effects model. For SQSTM1, its expression level was significantly upregulated in tumor group (SMD = 1.19, 95% CI: 1.00-1.37, P < 0.001; Figure 20). Moreover, remarkable overexpression of BIRC5 and SQSTM1 and opposite expression of FOXO1 protein in HCC was confirmed by Proteinatlas ( Figure 21).

Identification of clinical significances by in-house immunohistochemistry
BIRC5, FOXO1 and SQSTM1 protein expression level were detected by immunohistochemistry (IHC) in 302 HCC tissues and in 41 non-tumor tissues. For BIRC5, 97 samples (32.12%) were assigned to the positive group and 205 cases (67.88%) were assigned to the negative group according to the immunoreactive score (IRS) of the tissue samples ( Figure 22C-22D). In non-tumor tissues, 2 samples (4.9%) were exhibited positive BIRC5 expression, and 39 (95.12%) were negative ( Figure 22A-22B). HCC and non-tumor tissues exerted statistical difference in BIRC5 protein expression (P < 0.001). Similarly, the IHC results of SQSTM1 also confirmed its higher expression in HCC tissues. A total of 152 (50.33%) HCC tissues were positive while 150 (49.67%) samples were negative ( Figure  22K-22L). The positive rate was significantly higher than non-tumor tissues (P = 0.001), which showed low positive rate (9/41, 21.95%) ( Figure 22I  Although there was no significant difference, which may be owed to the limited samples, FOXO1 still showed a high positive rate in non-tumor tissues. The correlations between these genes expression and clinicopathological features of HCC are presented in Table 4. We observed significant positive correlations between BIRC5 expression and age (P = 0.019) and the pathologic T stage (P = 0.017). Up-regulated SQSTM1 also showed notable relationship with advanced pathologic T stage (P = 0.029). And lower FOXO1 expression was significantly related to advanced TNM stage (P = 0.002).

Molecular pathways disturbed between highand low-risk groups
Gene Set Enrichment Analysis was computed to pick up the molecular pathways disturbed between the highand low-risk groups. According to the results of the GSEA  analysis, genes co-expressed in the high-risk group were significantly enriched in the biological pathways associated with the E2F targets, MYC targets and G2M checkpoint ( Figure 23).

DISCUSSION
HCC patients are at risk for frequently returning cancer and high mortality, even after complete surgical debulking. Reliable molecular biomarkers for HCC prognosis prediction are hence significant for selecting patients who might be sensitive to additional targeted therapy. Using RNA-Seq technology, several novel prognostic indicators based on gene expression profiles have been proposed. Nevertheless, most studies only focused on a signal biomarker, which may not yield sufficient results that are inherently unpredictable and risky. In the present study, we constructed a novel prognostic model based on the expression levels of autophagy-related genes. Furthermore, the prognostic signature was effectively validated against multiple distinct datasets.
In recent years, elaborate molecular research for autophagy in tumors has proposed profound changes in our understanding of tumor management. Katheder NS et al. [21] discovered that malignant tumor cells induced autophagy of the surrounding normal cells in the tumor microenvironment and then release amino acids that are absorbed to support tumor growth. The precise mechanism of autophagy in HCC is elusive and changes frequently [22]. Unbalanced autophagy could disturb various vital signaling pathways such as PI3K/AKT/mTOR      [23][24][25], ERK/MAPK [26,27] and the Wnt/β-catenin [28,29] and therefore act as an accomplice of malignant hepatocyte to facilitate the HCC progression. In our study, we first identified differentially expressed autophagyrelated genes in HCC. Subsequent several gene functional enrichment analyses were performed to assess biological process, disease or pathways influenced by these genes. GO functional enrichment analysis indicated that differentially expressed autophagy-related genes were enriched in autophagy and several proliferation-related biological processes in HCC. KEGG and DO analyzed revealed that dysregulated autophagy-related genes may also play a critical role in the various cancers, which is coincided with other previous studies. These findings suggested that autophagy could be a driver in the process of onset and proliferation of HCC. Through univariate and multivariate survival analysis, we identified that BIRC5, FOXO1 and SQSTM1 were significantly related with the HCC patients' overall survival. BIRC5, a member of the anti-apoptosis family, mainly exerted its influence on HCC cells in inhibiting apoptosis [30], promoting proliferation [31], enhancing  resistance to radiotherapy and chemotherapy [32] and inducing the tumor stromal angiogenesis [33]. It has been reported that BIRC5 is directly correlated with autophagosome formation and contributed to HCC cells survival [34]. FOXO1 is thought to be a tumor suppressor that is down-regulated in HCC to exert its positive role in reversing the epithelial-to-mesenchymal transition program [35]. However, we still know relatively little about the autophagy-related mechanisms of FOXO1 that control the course of tumorigenesis and tumor development. SQSTM1 is autophagy-adaptor gene and accumulates in autophagy-deficient tumors, including HCC. Its aberrant accumulation and phosphorylation have been a concern in HCC proliferation potency [36,37] and could represent a potential chemotherapeutic approach against HCC [38,39]. Our results also evidently corroborated that there are sharp upregulations of BIRC5 and SQSTM1 in HCC while FOXO1 was significantly downregulated by an integrated analytic approach. We then developed a PI based on the expression of the three genes. The prognostic signature related to autophagy may offer broad prospects for modifying clinical management strategies or at least to prolong the lifespan of HCC patients in the precision medicine era. Due to such evidence, PI may also provide a better program for patients treated with a combination of multiply targeted agents.
Our novel risk score model was constructed and displayed the ability to stratify patients in different pathological stages and histological grades into subgroups with two distinct survival outcomes. More importantly, the predictive model was confirmed by GEO datasets. The TCGA program, other large-scale data and IHC results provided a comprehensive and reliable source to assess the molecular features which were most interrelated with the clinical outcomes of HCC patients. In our previous articles, Wang ZH et al. [40] mined RNA-Seq data of HCC patients form the TCGA program and managed to construct a four-lncRNA signature, which was significantly associated with HCC prognosis. Similarly, Li B et al. [41] proposed a three-gene predictive signature BASED ON the expression data of three genes (UPB1, SOCS2 and RTN3). However, both these investigations rely on the TCGA dataset and were unable to validate the efficiency by USING other cohorts. Furthermore, we analyzed from the perspective  of autophagy, which could offer more precise information in clinical management. As such, our PI model provided potential directions for survival estimates of HCC and future clinical practice may come in the foreseeable future. Moreover, GSEA was performed to analyze the coordinate expression of genes between the high-and lowrisk group. The hallmarks of high-risk were distinguished from low-risk group, by the E2F targets, MYC targets and G2M checkpoint, which revealed that the prognostic signature might be inclined to sustain chronic proliferation. These hallmarks of the high-risk group fully indicated that the proliferative signals were relatively disturbing. Arguably, the ability to proliferate outside of cell cycle control is a fundamental trait of cancer cells [42]. Tumors are prone to autophagy because of the important role of autophagy in supplying amino acids and fatty acids to meet the needs the survival and proliferation needs of the cell [14].
Thus, in conclusion, we assessed the gene expression profiles data of autophagy-related genes based on the TCGA database and proposed a risk score model which had a moderate efficacy in predicting the OS of HCC patients both in univariate and multivariate survival analysis. These results shown that the autophagy-related genes signature is a promising prognostic indictor and provided a better understanding of autophagy in HCC. However, prospective studies are needed to further verify the clinical utility as well as the biological function of the signature using more experiments.

Differently expressed autophagy-related genes
We obtained autophagy-related genes via database HADb. RNA-seq data of autophagy-related genes from each individual and their clinical information was download from the available TCGA database, which containing 374 HCC and 50 adjacent non-tumor tissues. To test for differential expression between HCC and their non-tumor counterparts, the R package edgeR program was analyzed raw count data with the criteria of Padj < 0.05 and the log FC>1. Then, all the data were converted into [log2 (data+1)] for further analysis.

Gene functional annotation
To reveal the biological function and disturbing pathways of these differently expressed autophagy-related genes, we used clusterProfiler package of R software for GO, DO and KEGG pathways [43]. The visualize enrichment maps of annotation analysis results were drawn by R with the "ggplot2" and "GOplot" packages.

Manufacture of potential prognostic signature for HCC
To identify the prognostic genes, these patients were analyzed by univariate Cox proportional hazards regression methods based on expression data as well as the clinicopathological features obtained from TCGA. The procedure was analyzed by using Survival package in R. All statistically significant indicators were selected as candidates for multivariate Cox regression analysis and used to construct the PI score model. The coefficient of each prognostic genes in the risk score model was derived from the corresponding multivariate analysis of these identified indicators. All patients were divided into high-risk or low-risk group based on the median value of the PI, which was calculated according to the expression profiles of the autophagy-related genes and the estimated regression coefficient. PI = exp BIRC5 *β BIRC5 +exp FOXO1 *β FOXO1 +exp SQSTM1 *Β SQSTM1 .

Verification of the PI via GEO datasets
The online multi-gene biomarkers validation tool SurvExpress (http://bioinformatica.mty.itesm.mx:8080/ Biomatec/SurvivaX.jsp) is a gene expression database that includes datasets collected form TCGA and GEO with clinical information. SurvExpress is an online bioinformatics tool to validate the prognostic performance based on a set of genes in various human cancers [44]. Through SurvExpress, we facilitated a Cox regression analysis and separated HCC patients into high-and lowrisk groups based on three liver cancer-related datasets: GSE10143, GSE10186 and GSE17856.

Validation for the expression values of genes
HCC-related microarray expression datasets were searched and downloaded from GEO. A computer-aided microarray search was performed by using the following search terms: (malignan* OR cancer OR tumor OR tumour OR neoplas* OR carcinoma) AND (hepatocellular OR liver OR hepatic OR HCC). The eligible datasets were required to match the following criteria: (1) human Figure 18: Forest plots for meta-analysis of evaluating BIRC5 expression between HCC and non-tumor tissues. SMD and associated 95% confidence interval were calculated using the random effects model. tissues; (2) proven diagnosis of HCC; (3) measuring the expression of genes included in PI in HCC and their normal counterparts; (4) providing at least three cases of HCC samples. These datasets were used to determine the expression pattern of the genes included in the PI model between HCC and non-tumor tissues. All the relevant data were extracted and presented as the mean ± SD. To give a comprehensive and visual display of our results, integrated analytic approach was implemented in the form of a meta-analysis. SMD was applied as the effect quantity to evaluate the association between the gene expression levels and HCC using the software STATA 12.0 software. SMD > 0 and its 95% CI not crossing the integer 0 suggested that genes are significantly upregulated in HCC. When SMD < 0 and its 95% CI not crossing the integer 0, it indicated genes are significantly downregulated in HCC. We chose a fixed-effect model or random-effect model to pool the SMD across the GEO datasets depending on the Q statistical analysis. What's more, we also observed the protein expression pattern of the three genes via The Human Protein Atlas detected by immunohistochemistry.

Immunohistochemistry
Two tissue microarrays obtained from Pantomics, Inc. (Richmond, CA) containing a total of 302 HCC samples and 11 non-tumor samples were used for the analysis of the present study. We also collected 30 nontumor tissues from the First Affiliated Hospital of Guangxi Medical University, People's Republic of China from SMD and associated 95% confidence interval were calculated using the random effects model. January 2016 to November 2017. All the collected tissues are treated under the strict confidentiality and according to appropriate applicable laws that protect the confidentiality of personal information. IHC analysis was performed according to the procedure in the manufacturer's instruction. The polyclonal rabbit anti-SQSTM1 primary antibody (1:1000 dilution; Abcam, ab207305), monoclonal rabbit anti-SQSTM1 primary antibody (1:300 dilution; Abcam, ab52857) and Rabbit anti-Survivin Monoclonal Antibody (1:80 dilution; ZSGB-BIO, ZA-0530) were used in the IHC analysis. Two pathologists (Yi-wu Dang and Gang Chen) independently evaluated the expression level via the IRS. IRS = (staining intensity) * (percentage of marked tumor cells). For percentage of marked tumor cells: 1-10% positive cells scored 1; 11-50% scored 2; 51-80% scored 3; and >80% positive cells scored 4. Staining intensity was scored as: negative staining (0), weak staining (1), moderate staining (2) and strong staining (3). All of the HCC samples were then divided into positive group (IRS ≥ 6) or negative group (IRS < 6) based on the IRS.

Gene set enrichment analysis
To reveal the meaning of the biological states from the obtained gene expression data between high-and low- SMD and associated 95% confidence interval were calculated using the random effects model.   risk group, a GSEA (v3.0, Broad Institute, Cambridge, USA) was performed. The phenotype label was set to high-risk versus low-risk. One thousand permutations were conducted to calculate for each analysis. And enrichment map was used for the visualization of the GSEA results. The false discovery rate (FDR) value and normalized enrichment score (NES) were used to identify the hallmarks significantly enriched in each phenotype.

Statistical analysis
All the statistical analyses were conducted using SPSS 22.0 (Chicago, IL, USA) and STATA 12.0 software (StataCorp, College Station, TX, USA). We used R software, OriginPro 2017 software (Northampton, Massachusetts, USA) and GraphPad Prism 5 (San Diego, CA, USA) for diagram drawing. Independent sample t-tests were used to analyze the expression patterns of genes in HCC and the clinicopathological parameters. The survival curves were plotted by Kaplan-Meier (K-M) method, and differences in the survival rates were assessed using the log-rank test. Univariate and multivariate Cox regression analyses were also performed to analyze the independent prognostic value by Survival package of R. All the statistical results were considered to be significant when the p value was less than 0.05.