Clinical Significance and Potential Mechanisms of ATP Binding Cassette Subfamily C Genes in Hepatocellular Carcinoma

The purpose of this investigation was to assess the diagnostic and prognostic significance of ATP binding cassette subfamily C (ABCC) genes in hepatocellular carcinoma (HCC). The Student t-test was used to compare the expression level of ABCCs between HCC and paraneoplastic tissues. Receiver operating characteristic curve (ROC) analysis was applied for diagnostic efficiency assessment. The Kaplan–Meier method and Cox proportional hazards model were respectively applied for survival analysis. Genes with prognostic significance were subsequently used to construct prognostic models. From the perspective of genome-wide enrichment analysis, the mechanisms of prognosis-related ABCC genes were attempted to be elaborated by gene set enrichment analysis (GSEA). It was observed in the TCGA database that ABCC1, ABCC4, ABCC5, and ABCC10 were significantly upregulated in tumor tissues, while ABCC6 and ABCC7 were downregulated in HCC tissues. Receiver operating characteristic analysis revealed that ABCC7 might be a potential diagnostic biomarker in HCC. ABCC1, ABCC4, ABCC5, and ABCC6 were significantly related to the prognosis of HCC in the TCGA database. The prognostic significance of ABCC1, ABCC4, ABCC5, and ABCC6 was also observed in the Guangxi cohort. In the Guangxi cohort, both polymerase chain reaction and IHC (immunohistochemical) assays demonstrated higher expression of ABCC1, ABCC4, and ABCC5 in HCC compared to liver tissues, while the opposite was true for ABCC6. GSEA analysis indicated that ABCC1 was associated with tumor differentiation, nod-like receptor signal pathway, and so forth. It also revealed that ABCC4 might play a role in HCC by regulating epithelial-mesenchymal transition, cytidine analog pathway, met pathway, and so forth. ABCC5 might be associated with the fatty acid metabolism and KRT19 in HCC. ABCC6 might impact the cell cycle in HCC by regulating E2F1 and myc. The relationship between ABCC genes and immune infiltration was explored, and ABCC1,4,5 were found to be positively associated with infiltration of multiple immune cells, while ABCC6 was found to be the opposite. In conclusion, ABCC1, ABCC4, ABCC5, and ABCC6 might be prognostic biomarkers in HCC. The prognostic models constructed with ABCC1, ABCC4, ABCC5, and ABCC6 had satisfactory efficacy.

Keywords: HCC (hepatic cellular carcinoma), ABCC gene family, prognosis (carcinoma), nomogram, GSEA (gene set enrichment analysis) BACKGROUND Hepatocellular carcinoma (HCC) generally followed cirrhosis given rise by metabolic disorder (Yang et al., 2019), chronic ethanol intake (Llovet et al., 2016), and hepatitis virus infection (Fujiwara et al., 2018). The leading metabolic risk factor for HCC is non-alcoholic fatty liver disease (NAFLD) (Zhang, 2018), which is mainly related to obesity and type 2 diabetes. Currently, it is acknowledged that hepatitis B virus (HBV) and hepatitis C virus (HCV) were the main infectious etiologies for cirrhosis and HCC. It is worth mentioning that viral hepatitis could skip cirrhosis and induced HCC directly and independently (El-Serag, 2012;Levrero and Zucman-Rossi, 2016). Besides the aforementioned factors, the intake of Aflatoxin B1 (AFB1) was also demonstrated to be related to HCC (Rushing and Selim, 2019). People in specific regions entailing relatively high exposure to Aflatoxin B1 were accompanied by high incidence and mortality of HCC (Long et al., 2008;Wogan et al., 2012;. More than 8 million new cases of liver cancer occurred worldwide each year, which directly or indirectly gave rise to more than 4 million deaths worldwide each year (Torre et al., 2015;Bray et al., 2018). Asia has the highest incidence of liver cancer in the world, particularly in China, which accounts for almost half the global cases (Akinyemiju et al., 2017). At the same time, Asia is the high-incidence area of HBV and HCV (Gower et al., 2014;Polaris Observatory Collaborators, 2018). The main treatment methods of liver cancer mainly include surgical resection, transcatheter arterial chemoembolization (TACE), ablation, liver transplantation, radiotherapy, and so forth (Fattovich et al., 2004). Sorafenib, the multi-kinase inhibitor, is one of first-line drugs approved for the treatment of advanced HCC. Although it can improve survival, the long-term survival of HCC patients is limited due to the drug resistance. Hence, the discovery of new hub genes for developing HCC-targeted drugs and specific genes that improve and maintain drug susceptibility might be hopeful for advanced-stage HCC patients.
The ATP binding cassette subfamily C (ABCC) subfamily includes 13 members whose protein products take effect in transporters with different functional profiles, including ion transport, cell surface receptor, and toxin secretion activity (Childs and Ling, 1994;Robey et al., 2018;Yamada et al., 2018). The ATP-binding domain of the ABCC product possesses distinctive conserved motifs (Walker A and B motifs), which are separated by an uncertain sequence of around 100 amino acids . The distinctive interval and conserved motifs distinguish ABCC members from other ATP-binding proteins (Higgins et al., 1986). Genetic variations in these genes are substantiated in numerous research studies to be the cause or contributor to a variety of complex human diseases, including cystic fibrosis, neurological diseases, defects in cholesterol and bile transport, and drug responses. The ABCC subfamily plays an important role in the pharmacokinetics of endogenous and exogenous compounds. Studies have shown that the members of the ABCC family could transport drugs to the extracellular substances by virtue of ATP energy (Chen and Tiwari, 2011;Keppler, 2011;Leslie, 2012).

Data Acquisition and Specimen Collection
RNA-Seq data (FPKM) of 412 samples, 362 tumors, and 50 paraneoplastic tissues were acquired from the TCGA database (https://portal.gdc.cancer.gov/, accessed on 22 December 2019). The limma package was employed for normalization of this RNA-Seq data in R. Matched prognostic/clinicopathologic data of these 362 patients were acquired from UCSC Xena (http://xena.ucsc. edu/, accessed on 23 December 2019).
The HCC tissues and matched paracancer tissues of 102 patients hospitalized in the first affiliated hospital of Guangxi Medical University from September 2016 to December 2018 were collected after informed consent was obtained. Among them, excised tissues during surgery of 72 patients were well preserved in the Department of Pathology. Tissue slices of these patients were obtained from the Department of Pathology.

Expression Difference and Diagnostic Efficiency Analysis of ABCC Genes
The expression levels of ABCCs in HCC and paraneoplastic tissues were extracted from the RNA-Seq Chip matrix in the TCGA database. The normality test was assessed using the Kolmogorov-Smirnov normality test. Student's t-test was used to assess the statistical significance of ABCC genes' expression between HCC and paraneoplastic tissues. The area under the curve (AUC) of the receiver operating characteristic curve (ROC) was used to access the diagnostic efficiency of each ABCC gene in HCC. AUC > 0.8 with p < 0.05 was considered as satisfactory diagnostic performance (Hosmer et al., 2013).

Immunohistochemistry
Tissue sections were sequentially placed in xylene and graded concentrations of ethanol to achieve hydration. Antigens were repaired with a pH 6.0 citrate repair solution (ZSGB-BIO, Beijing, China). Subsequent antigen-antibody reactions and color development reactions were performed with the help of a universal two-step detection kit (Mouse/Rabbit Enhanced Polymer Detection System). Immunohistochemical scores were assessed by two experienced pathologists. Antibodies for ABCC1, ABCC4, ABCC5, and ABCC6 were diluted according to the recommended concentrations of the manufacturer (Proteintech, Wuhan, China).

Prognostic Significance Assessment of ABCC Genes
The patients in the TCGA database were divided into two groups in terms of the median value of each ABCC gene expression for survival analysis. The Kaplan-Meier method with a log-rank test was applied to assess the prognostic significance of each ABCC gene. The Cox proportional hazards model was applied to adjust the bias caused by prognosis-related clinicopathologic factors.
In terms of survival analysis results in the TCGA database, the prognostic significance of ABCC1, ABCC4, ABCC5, and ABCC6 was further validated in the Guangxi cohort.
For better predicting the prognosis and evaluating the combined effect of ABCCs, prognosis-related ABCCs (ABCC1, ABCC4, ABCC5, and ABCC6) were integrated in pairs into combined effect survival analysis. The patients were divided into four groups in terms of the expression level of ABCCs with details displayed in Table 2. The Kaplan-Meier method with the log-rank test and Cox proportional hazards model were applied to assess the prognostic significance.

Nomogram
Independent prognostic factors, including ABCCs and clinicopathologic features, were integrated to construct the nomogram in R with the rms package (Iasonos et al., 2008). In the nomogram, the risk degree of each variable in the nomogram was displayed by the integration line, and the total risk score is obtained by adding up the risk value of each variable (Zhang Z et al., 2017). The model was validated for calibration and discrimination using the bootstrap method (Wang et al., 2013).

Prognostic Signature Construction
The Cox proportional hazards model was used to assess the risk coefficient of ABCCs in overall survival. Then, the prognostic signature was constructed in terms of the expression of ABCCs and the corresponding risk coefficient. The formula of prognostic signature construction is as follows: Risk score N 1 (ExpVlueipβi) (Chen M et al., 2017). N is the number of prognostic genes. ExpVluei is the expression value of each ABCC gene. βi is the risk coefficient of the corresponding ABCC gene. A time-dependent ROC curve was constructed in R (version 3.6.2; www.r-project.org) with the survivalROC package to evaluate the availability of this prognostic signature (Chen M et al., 2017).

Biological Functional Exploration of ABCC Genes
The Gene Ontology (GO) database, the integrated database of calculable information about the functions of genes, was comprehensively used for identifying unique biological properties of high-throughput transcriptome or genome data (The Gene Ontology Consortium, 2017;Chen L et al., 2017). KEGG is a collection of databases dealing with genomes, diseases, biological pathways, drugs, and chemical materials (Kanehisa et al., 2017). DAVID (The Database for Annotation, Visualization, and Integrated Discovery, https://david.ncifcrf.gov/) is an online bioinformatics tool to access the GO database and the KEGG database (Long et al., 2008). DAVID was used to access the enrichment of biological functions and pathways of ABCC genes in this investigation. Then, the enrichment biological functions and pathways were visualized in R Studio (Version 1. 2.5033) with packages Goplot, Hmisc, and ggplot2 (Nolan et al., 2013;Ito and Murphy, 2013). The Biological Networks Gene Ontology tool (BiNGO) is an open-source online database, which was employed to determine the significantly overrepresented GO terms of ABCC genes (Maere et al., 2005). Functions and interactions of ABCC genes were performed in Genemania (http://genemania.org/. accessed on 11 August 2020) and STRING (https://string-db.org/. accessed on 11 August 2020), respectively (Szklarczyk et al., 2015;Luo et al., 2020).

Gene Set Enrichment Analysis (GSEA)
GSEA is software with additional resources for analyzing, annotating, and interpreting standardized chip matrices. In this investigation, GSEA enrichment was used to analyze the enriched biological pathways of ABCC1, ABCC4, ABCC5, and ABCC6 in the TCGA database. The Oncogenic Signatures c2.all.
v7.1.symbols.gmt data set was adopted as the reference data set. The biological pathways exported from GSEA with p < 0.05 and FDR < 0.25 were considered as significant results.

Correlation Analysis of Tumor-Infiltrating Immune Cells and ABCC Gene Expression
TIMER (http://timer.cistrome.org/) is a comprehensive resource for the systematical analysis of immune infiltrates across diverse cancer types, which provides immune infiltrates abundances estimated by multiple immune deconvolution methods. In this investigation, TIMER was accessed to explore the correlation between infiltrating immune cells and ABCC expression in HCC.

Expression and Diagnostic Efficiency of ABCC Genes in HCC
Several ABCC genes were discovered to be differentially expressed in HCC and paraneoplastic tissues based on the RNA-seq data of the TCGA database. ABCC1, ABCC4, ABCC5, and ABCC10 (Figures 1A,C,D,H) were significantly higher expressed in HCC tissues, but ABCC2, ABCC6, ABCC7, and ABCC9 (Figures

Prognostic Significance of ABCC Genes
Subsequently, the prognostic significance of ABCC genes was systematically discussed. The clinicopathologic characteristics of 362 HCC tissues in the TCGA database are displayed in Supplementary Table S1. The expression levels of ABCC1 (log-rank p = 0.002, adjusted p = 0.008, adjusted HR = 1.656), ABCC4 (log-rank p = 0.026, adjusted p = 0.038, adjusted HR = 1.479), ABCC5 (log-rank p = 0.002, adjusted p = 0.001, adjusted HR = 1.928), and ABCC6 (log-rank p < 0.001, adjusted p = 0.001, adjusted HR = 0.534) were significantly associated with the overall survival of HCC patients in univariate and multivariate survival analysis (  effect survival analysis. In the combined effect survival analysis, it was observed that there was an even bigger prognosis difference among groups in combined effect survival analysis ( Table 2).
Group C with high expression of ABCC1 and ABCC4 was significantly correlated to bad outcome (p = 0.001, Figure 4A), so were group 3 with high expression of ABCC1 and ABCC5 (p <  0.001, Figure 4B), group c with high expression of ABCC1 and low expression of ABCC4 (p < 0.001, Figure 4C), group Ⅲ with high expression of ABCC4 and ABCC5 (p = 0.001, Figure 4D), group ⅲ with high expression of ABCC4 and low expression of ABCC6 (p < 0.001, Figure 4E), and group γ with high expression of ABCC5 and low expression of ABCC6 (p < 0.001, Figure 4F).
Nomogram Based on ABCC1, 4, 5, and 6 and Tumor Stage In the survival analysis, we found that ABCC1.4.5.6 was strongly associated with the prognosis of HCC. In addition, the clinical factor tumor stage could also partially distinguish patients with good and bad prognoses. Thus, a nomogram integrating clinical elements and ABCC gene expression was constructed in terms of the COX proportional hazards model. In the nomogram, the contribution of ABCC1, ABCC4, ABCC5, ABCC6, and clinicopathologic features to the overall survival of HCC patients was displayed by virtue of the length of the scales ( Figure 5A). The calibration plot for 1-, 3-, and 5-year survival after surgery revealed a satisfactory overlap between calculation and reality ( Figures 5B-D).

Prognostic Signature Based on the TCGA Database
In terms of the expressions of ABCC1, ABCC4, ABCC5, and ABCC6, the prognostic signature for HCC patients was built in the TCGA database and Guangxi cohort. Each HCC patient was assigned with a risk score in terms of the expression of ABCC1, ABCC4, ABCC5, and ABCC6. In the prognostic signature built for the TCGA database, the risk score for each patient was displayed in the upper scatter plot, and the patients were divided into two groups based on the median value ( Figure 5E). The survival time and survival status of specific patients can be observed from the middle scatter plot, which showed that the dots representing patients in the highrisk group tended to cluster lower ( Figure 5F). The expression levels of ABCC1, ABCC4, ABCC5, and ABCC6 in patients were presented in the form of heat maps ( Figure 5G). A significant difference in overall survival was observed between the highrisk and low-risk groups ( Figure 5H, p = 0.003). The AUC value of the prognostic signature for 1-year, 3-year, and 5-year overall survival prediction was 0.689, 0.619, and 0.598, respectively ( Figure 5I).

Validation in the Guanxi HCC Cohort
A total of 102patients who were hospitalized in the first affiliated hospital of Guangxi Medical University from September 2016 to December 2018 were taken into the group for validation. The baseline information for these patients is presented in Table 3. The expressions of ABCC in HCC tissues and in paraneoplastic tissues were detected by immunohistochemical (IHC) and polymerase chain reaction (PCR) assays, respectively. In the IHC assay, the expressions of ABCC1, ABCC4, and ABCC5 in HCC tissues were significantly higher than that of paraneoplastic tissues, while ABCC6 was higher expressed in paraneoplastic tissues ( Figure 6A). The same expression trends of ABCC1, ABCC4, ABCC5, and ABCC6 were observed at the mRNA level ( Figures 6B-E). The prognostic significance of ABCC1, ABCC4, ABCC5, and ABCC6 was  Table 2.
In the prognostic signature built for the Guangxi cohort, patients were divided into two groups in terms of the risk score ( Figure 6J). The same as above-mentioned, the dots representing patients in the high-risk group also tended to cluster lower ( Figure 6K). The expression levels of ABCC1, ABCC4, ABCC5, and ABCC6 in patients were presented in the form of heat maps ( Figure 6L). The prognosis of the high-risk group was significantly worse than that of the lowexpression group ( Figure 6M, p = 0.001). The AUC value of the prognostic signature for 1-year, 2-year, 3-year, and 4year overall survival prediction was 0.787, 0.772, 0.838, and 0.788, respectively ( Figure 6N). In the nomogram, the contribution of ABCC1, ABCC4, ABCC5, ABCC6, and clinicopathologic features to overall survival was displayed by the length of the corresponding scales (Supplementary Figure S4A). The calibration plot for 1-, 2-and 3-year survival after the surgery revealed a satisfactory overlap between calculation and reality (Supplementary Figures  S4B-D).

Biological Functional Exploration of ABCCs
The enrichment analysis of the ABCC gene by setting Homo sapiens as the background was performed on the DAVID online database for obtaining enrichment information about GO terms. The corresponding relationship between ABCCs and GO terms is displayed in Figure 7A. The enrichment analysis of GO showed that ABCCs were mainly related to ATP binding, ATP activity, transmembrane, and other biological functions ( Figure 7B). The bubble color from red to green represents the biological function of the -log (p-value) from high to low. The network diagram of the relationship between enriched GO terms is shown in Figures 7C,D. Interactions of ABCCs which were analyzed from STRING and Genemania are respectively displayed in Figures 7E,F.

GSEA
The GSEA results revealed that the expression of ABCC1 was associated with tumor differentiation, nod-like receptor signal pathway, resistance to the bcl2 inhibitor up, and so on ( Figures  8A-F)   acid metabolism and the expression of kt19 and myc ( Figures  8M-O). The result of GSEA revealed that high expression of ABCC6 was accompanied with lower HCC late recurrence ( Figure 8P). It also illustrated that ABCC6 might impact HCC by regulating E2F1 and myc ( Figures 8Q,R).

Correlation Analysis of ABCC Gene Expression and Tumor-Infiltrating Immune Cells
The estimation of the abundance of immune cell infiltration showed that ABCC1, ABCC4, and ABCC5 were significantly positively associated with infiltration of immune cells, which include B cells, CD8 + T cells, CD4 + T cells, macrophages, neutrophils, and dendritic cells ( Figures 9A-C). However, ABCC6 was negatively associated with the infiltration of immune cells ( Figure 9D).

DISCUSSION
ABCC expressions were analyzed in two data sets, and consistent results were obtained. Compared with normal tissues, they revealed that ABCC1, ABCC4, ABCC5, and ABCC10 were significantly upregulated in HCC tissues, while ABCC6 and ABCC7 were significantly downregulated in HCC tissues. In the TCGA database, ABCC5, ABCC7, ABCC9, and ABCC10 were equipped with high diagnostic efficacy for HCC (AUC > 0.8). In GSE76427, the good diagnostic efficacy for HCC was only discovered in ABCC7. Combining the results of the two data sets, we consider ABCC7 as a potential diagnostic marker for HCC. In the TCGA database, ABCC1, ABCC4, ABCC5, and ABCC6 were found to be associated with the prognosis of HCC, while further verification in GSE14250 indicated that only ABCC6 was significantly correlated to the prognosis. The results of survival analysis in the two data sets were very similar, although not identical. We observed that the expression of ABCC1 and ABCC5 was associated with the prognosis of liver cancer in both data sets. The reason for the different conclusions may lie in the population difference and inconformity in the causes of neoplasm. HCC patients in GSE76427 were mainly in the Asian population, and the proportion of hepatitis B virus infection was high. However, the majority of HCC patients in the TCGA database were Caucasian and the proportion of hepatitis B virus infection was low. The sample size of both databases is relatively large, and the follow-up data were also of high quality. The results from both databases should be reliable but may apply to different populations. Both clinicopathologic features and biomarker expression were included in the nomogram as prognostic dependent variables, with the length of each variable clearly reflecting its contribution to the prognosis of liver cancer.
Based on the four prognostic biomarkers obtained from the survival analysis, we further performed combined effect survival analysis, nomogram, and prognostic signature based on biomarker expression. The combined survival analysis had obvious advantages, and the prognostic difference between groups was more remarkable. The length of each variable in the nomogram clearly reflects its contribution to the prognosis of liver cancer.
ABCC1 transports drugs to the extracellular substances, thereby reducing the drug concentration and generating drug resistance in cancer (Wlcek and Stieger, 2014). In the liver, ABCC1 undertakes excretion of the drugs into the bile (Zhou, 2008). The ontogeny, localization, expression, and function of ABCC1 in HCC were reported in several research studies, and the previous reports mainly focused on the role of ABCC1 in HCC drug resistance (Flens et al., 1996;Nies et al., 2001;Vander Borght et al., 2005). It was reported that ABCC1 was significantly upregulated in the tissues of oxaliplatin-resistant, 5fluorouracil-resistant, and sorafenib-resistant HCC patients (Ding et al., 2017;Huang et al., 2018;Ding et al., 2019). In HCC, increased ABCC1 expression was related to increasing dedifferentiation, tumor size, and microvascular invasion (Vander Borght et al., 2008;Zhou, 2008).
Located on the inner surface of the basal side of the liver cells, ABCC4 undertakes bile salt transport (Borst et al., 2007). Previous studies have shown that ABCC4 expression is extremely low in the normal adult liver and fetal liver (Sharma  et al., 2013), and ABCC4 expression is significantly increased in cholestatic hepatocyte cell membranes (Gradhand et al., 2008;Sharma et al., 2013). Studies have shown that ABCC4 is highly expressed in HCC tissues (Sekine et al., 2011;Borel et al., 2012;Luo et al., 2020). Recently, ABCC4 was found to play an important role in HCC oncogenesis and development promoted by decreasing the haploid of p53 (Luo et al., 2020). In addition, ABCC4 could specifically and independently distinguish the aggressive subtypes of HCC (Gradhand et al., 2008).
Here are a few reports on ABCC5 in HCC, with the existing relevant study indicating that ABCC5 is highly expressed in the liver cancer tissues. Our findings in this investigation also confirm this conclusion.
T lymphocytes are known as the main cells of the tumor immunity. Cytotoxic CD8+T cells play a particularly vital role in anticancer immune response (Vesely et al., 2011;Raskov et al., 2021). Once successfully activated, CD8+T cells secreted death-inducing granules to enhance the killing effect of target cells (Basu et al., 2016). Accumulating evidence indicates that TRM (tissue-resident CD8 + memory T cells) is essential for suppressing cancer growth. In a mouse model, whether generated during tumorigenesis or prior to tumor challenge, antitumor TRM cells revealed suppression in cancer growth (Park et al., 2019). Regulatory T cells inhibit anticancer immunity by preventing the protective immunosurveillance of neoplasia and hindering antitumor immune responses in tumor-bearing hosts, thereby promoting the tumor progression Wing and Sakaguchi, 2010;Togashi and Nishikawa, 2017).
B cells have a crucial part in the regulation of T cell response against tumors (Olkhanud et al., 2011;Tadmor et al., 2011). There is a crosstalk between the B and T lymphocytes in antitumor immunity (Blair et al., 2010;DiLillo et al., 2010). Natural killer cells (NK cells) in cancer are involved in priming a multilayered immune response to achieving long-lasting immunity against tumors, in which T cells are involved (Morandi et al., 2012;Ferlazzo and Moretta, 2014). Moreover, NK cells generate cytokines and chemokines that regulate immune responses. The function of non-NK ILCs (innate lymphoid cells, ILCs) in cancer remains unclear.
Combining this investigation and previous research studies, we could preliminarily conclude that ABCC1, ABCC4, and ABCC5 reduce drug sensitivity by influencing drug transport out of cells, thus resulting in a poor prognosis in these patients with HCC. In this study, we also found a significant positive correlation between ABCC1, ABCC4, and ABCC5 expression and immune cell infiltration.
There were no reports on ABCC6 in HCC before. The role of ABCC6 in HCC is completely opposite to that of ABCC1, ABCC4, and ABCC5. We found that ABCC6 expression was decreased in the liver cancer tissues, and the patients with low ABCC6 expression had a better prognosis. We speculate that ABCC6 may function through a completely different mechanism, and the specific findings need to be further studied.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusion of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
XZ and J-mH formulated the research direction and inclusion criteria for this study. G-zZ and XZ were mainly responsible for data extraction and the elimination of cases that did not meet the inclusion criteria. X-pY and TP analyzed the data, drew charts, and wrote the manuscript. All authors read and approved the final manuscript. All authors made a significant contribution to the work reported, whether this is in the conception, study design, execution, acquisition of data, analysis, and interpretation or in all these areas; took part in drafting, revising or critically reviewing the article; gave final approval of the version to be published; have agreed on the journal to which the article has been submitted; and agree to be accountable for all aspects of the work.