This article is mentioned in:

Introduction

Colorectal cancer (CRC) has high morbidity and mortality rates worldwide, at 10.2 and 9.2%, respectively (1). Despite the therapeutic advances and earlier detection, the 5-year survival rate of patients with CRC remains unsatisfactory (2). One of the main reasons for this is that the occurrence of CRC is a complex multi-stage process, and involves further investigation into the proliferation, differentiation, apoptosis and survival mechanism of intestinal epithelial cells (3). Therefore, biomarkers for early detection and targeted therapy are urgently required.

Biological rhythms are produced by conserved transcription and translation feedback loops of circadian clock genes within the cells (4). A circadian disruption has been recognized as a potential independent risk factor for cancer development (5). Circadian clock genes appear to have multifaceted functions during cancer development and can act to both suppress tumors and promote carcinogenesis (6). Research by the International Agency for Research on Cancer has also demonstrated that this disruption increases the risk of CRC (7). Several previous studies have also demonstrated that large variations in expression levels, both up- and downregulated, and the circadian clock genes are associated with tumor progression and mammalian tumorigenesis for several malignancies, such as breast cancer (8), liver cancer (9) and colorectal carcinoma (10). In addition, the association between single nucleotide polymorphisms (SNPs) in circadian clock genes and disease has also been analyzed (11,12). These studies indicate that mutations or deregulated expression of circadian clock genes are frequently detected in different tumors. Copy number variation (CNV) is a kind of structural variation at the submicroscopic level, which refers to the complex chromosomal structural variation forms derived from the deletion and/or duplication of DNA fragments longer than 1 kb (13). Increasing research has shown that CNV is closely associated with the risk of tumor occurrence (14,15). The mechanism of action of CNV and circadian clock genes in many cancer types has been intensively investigated, such as liver cancer (16) and lung cancer (14); however, the study of the mechanism of action of CNV-driven circadian clock genes in cancer (including CRC) has not yet been reported.

The aryl hydrocarbon receptor nuclear translocator-like 2 (ARNTL2) gene is also known as brain and muscle ARNT-like 2 (BMAL2), which is mapped to human chromosome 12p11.22–11.23 and shares 52% amino acid identity with zebrafish Bmal2 and 49% identity with human BMAL1 (17). Schoenhard et al (18) hypothesized that the different ARNTL2 spatiotemporal distributions allow intrinsic circadian clocks to modulate the amplitudes of their oscillators while maintaining circadian periodicity. Research on ARNTL2 in various complex diseases (19,20), particularly cancer, has gradually become accepted. Studies have shown that ARNTL2 is a potential biomarker for tumor invasion in colorectal cancer (21), and it is significantly associated with lung cancer risk (22). In addition, a previous study has analyzed SNPs associated with ARNTL2 expression in patients with breast cancer (23).

The aim of the present study, using somatic CNV, legacy clinical information and gene expression data from the Cancer Genome Atlas (TCGA; http://tcga-data.nci.nih.gov/tcga/), was to investigate the systematic association between somatic cell CNV and circadian clock gene expression in patients with CRC, and to identify ARNTL2 as a contributing gene in CRC development that may serve as a promising therapeutic strategy.

Materials and methods

Data source and preprocessing

The CNV data were downloaded from TCGA data portal on October 23, 2018. The data contained 979 files and 460 cases for CNV analysis by setting specific parameters: Data Type was Masked Copy Number Segment. In addition, mRNA expression profile data and the corresponding legacy clinical information of patients with CRC from TCGA were also downloaded and contained 480 CRC tumor specimens and 41 tumor-adjacent tissue specimens. Firstly, 18 samples without adequate clinical information were removed, which left 462 patients with CRC with complete survival information. Subsequently, low-abundance mRNA expression data were removed; mRNAs with expression value >1 in 90% samples were retained. For the duplication data in one sample, the average values of the mRNA expression were adopted. The 2,083 differentially expressed mRNAs were analyzed using R/Bioconductor package edgeR (version 3.26) (24), with the criteria of |log2fold-change (FC)| >1.5 and q-value <0.01. No patients were involved in clinical trials in this study.

Identification and functional analysis of CNV-driven circadian clock genes

Gene Ontology (GO) analysis was performed to explore the functional roles of the target genes using DAVID (http://www.david.abcc.ncifcrf.gov/) (25). Finally, the enriched GO terms with gene count >5 and P<0.05 were selected for further analyses. Cytoscape software (version 3.7.1) (26) (with ClueGO and CluePedia plugins) was used for the Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses, showing only pathways with P<0.05.

Circadian clocks exist endogenously in almost every organism (6). The Circadian Gene Database (CGDB; version 1.0; http://cgdb.biocuckoo.org/index.php) (27) was used to identify the circadian clock genes. Circadian genes were selected that had been identified experimentally. A literature search using PubMed database was performed to identify the latest candidate circadian clock genes.

To verify the expression profile of ARNTL2 in CRC tissues and their non-tumoral counterparts, a meta-analysis was performed using the Oncomine database (version 4.5; www.oncomine.org) by setting specific parameters: ‘ARNTL2’, ‘Cancer vs. Normal Analysis’, ‘Colorectal Cancer’ and ‘mRNA’.

The java software Gene Set Enrichment Analysis (GSEA; version 3.0) was employed to perform the statistical significance test between two phenotypes (http://software.broadinstitute.org/gsea/index.jsp), with gene expression data and phenotype data (high/low group of expression values of ARNTL2) to be prepared according to the GSEA guidelines (28). The parameters were set as follows: Using KEGG pathway as a reference, permutation type to be the phenotype, and at least 15 genes in a single pathway. The mean expression levels (905.75) of ARNTL2 in all cancer samples were obtained. In the GSEA analysis, the expression level higher than this value is considered to be high expression, and below this value is considered to be low expression.

Statistical analysis

The segment mean at (−0.2, 0.2) was generated by the error of the instrument measurement, so the copy number of such genes was confirmed as unchanged. A χ2 test was used to compare the number of CNVs in cancer tissue and paracancer tissue, with a criterion of false discovery rate (FDR) <0.01. A Kolmogorov-Smirnov test was used to identify genes with CNV and expression consistency, with the criterion of P<0.005. A Kolmogorov-Smirnov test is based on cumulative distribution functions to test whether a distribution conforms to a theoretical distribution or whether there is a significant difference between two empirical distributions. To assess the result set of genes, a hypergeometric test was used to verify whether known CRC-related genes were enriched on the set. To identify the associations between clinicopathological parameters and the presence of copy number loss or gain in the regions containing selected genes, a Pearson's χ2 test was performed. A Kaplan-Meier curve analysis was performed to analyze the association between the gene and survival time, and statistical significance was assessed using the R package ‘survival’ (29). P<0.05 (two-sided) was considered to indicate a statistically significant difference.

Results

Patient characteristics

The detailed clinical and pathological characteristics of the study population, including age, sex, pathological stage, pathological tumor (pathological T), pathological node (pathological N) and pathological metastasis (pathologic M), were summarized in Table I. All the 462 patients were pathologically diagnosed with colorectal cancer. The median age for all patients was 60 years (interquartile range, 31–90 years).

Table I. Clinicopathological features of the 462 patients with colorectal cancer.

Table I.

Clinicopathological features of the 462 patients with colorectal cancer.

Feature	Primary, n (%)	Metastatic, n (%)	NA, n (%)
Age, years
<60	81 (24.0)	42 (36.5)	4 (40.0)
>60	256 (76.0)	73 (63.5)	6 (60.0)
Sex
Male	177 (52.5)	53 (46.1)	6 (60.0)
Female	160 (47.5)	62 (53.9)	4 (40.0)
Pathological T
T1-T2	77 (22.8)	10 (8.7)	3 (30.0)
T3-T4	260 (77.2)	105 (91.3)	7 (70.0)
Pathological n stage
N0	231 (68.5)	35 (30.4)	5 (50.0)
N1-N2	106 (31.5)	80 (69.6)	5 (50.0)
Pathological stage
I–II	228 (67.7)	23 (20.0)	2 (20.0)
III–IV	106 (31.5)	87 (75.7)	5 (50.0)
NA	3 (0.9)	5 (4.3)	3 (30.0)
Vital status
Alive	288 (85.5)	73 (63.5)	7 (70.0)
Death	49 (14.5)	42 (36.5)	3 (30.0)

[i] T, tumor; N, node.

Screening of potential CRC-related gene CNVs

To identify potential candidate genes within the regions exhibiting CNVs in the TCGA dataset, the frequency of copy number loss and gain in the regions was obtained. First, the instrument measurement error was filtered, and the area where the CNV number was significantly different located, and finally the genes in these areas were identified. Finally, the χ2 test was conducted on CNV, and a total of 10,256 genes with significant differences in CNV expression were obtained. KEGG and GO enrichment analyses was then performed with a smaller set of genes (n=295). A detailed workflow chart of the methodology is illustrated in Fig. 1A. CNV occurred differently on each chromosome in patients with CRC. Large-scale losses of copy numbers occurred only on certain chromosomal regions, such as chromosomes 4, 11, 14, 15, 18, 21 and 22. However, on other chromosomes, such as chromosomes 7, 12 and 13, only gains occurred (Fig. 1B).

Figure 1. Analyzed workflow and CNV identification results. (A) Workflow for identification and analysis of CNV-driven circadian clock genes. (B) Circos plot of the human genome showing chromosome structure and CNV in CRC. The outer most layer is the chromosome model, and the inner layer illustrates the CNV. A gene with CNV is annotated in the inner layer, in which a blue dot near the center of the circle denotes the presence of copy number loss, and a black dot near the outer layer denotes the presence of copy number gain. CNV, copy number variation; CRC, colorectal cancer; CGBD, Circadian Genome Database; GSEA, Gene Set Enrichment Analysis; ARNTL2, aryl hydrocarbon receptor nuclear translocator-like 2.

Screening of differentially expressed mRNAs

Based on the threshold criteria of |log2FC| >1.5 and q-value <0.01, 2,083 mRNAs were identified as aberrantly expressed mRNAs in the CRC tissues compared with that in the adjacent non-tumorous tissues. It was found that a number of mRNAs were upregulated or downregulated >100-fold (Fig. 2A). To further investigate the mechanism of CNV in the development and progression of CRC, the intersection of genes involved in significant abnormal CNV and differentially expressed genes was obtained. Subsequently the association analysis of the expression profiles and copy number profiles for the aforementioned small gene set was performed, with a result that 295 mRNAs had statistically significant differential expression and a difference in CNV. Finally, GO enrichment analysis and a KEGG pathway analysis of these mRNAs were performed, suggesting that the mRNAs were primarily enriched in only one KEGG pathway (P<0.05; Fig. 2B) and eight GO terms (Benjamin P<0.01; Fig. 2C). Recent evidence suggests that the circadian system can influence the Wnt/β-catenin signaling pathway (30), which is a critical pathway for the development and progression of CRC (31). Known CRC-related genes were mapped to the set of 295 mRNAs, and the 73 CRC-related genes were significantly enriched in this gene set (hypergeometric test, P=4.725561×10−9).

Figure 2. Functional analysis of copy number variation-driven genes. (A) Volcano plots show the expression profiles of mRNAs. The vertical line represents 2.0-fold up- and downregulation between the colorectal cancer tissues and the adjacent non-tumorous tissues, while the horizontal line represents the q-value. The red dots in the figure represent the differentially expressed mRNAs, with the left side indicating downregulation of mRNA expression and the right side indicating upregulation of mRNA expression. (B) The Kyoto Encyclopedia of Genes and Genomes pathway analysis of the 295 mRNAs. (C) The Gene Ontology enrichment analysis of the 295 mRNAs.

PubMed and CGDB databases were searched, and 15 of the 295 mRNAs were found to be circadian clock genes (Table II). Among the 15 circadian clock genes, NR3C2 and P2RX1 were downregulated, and the remaining 13 genes were upregulated in patients with CRC. Gain was the predominant type of alteration for BIRC7, GNGT1, NFE2L3, PDX1 and UBE2C, while loss of APCDD1 and P2RX1 was found in >30% of cases. No significant changes in the expression levels of other important genes, such as PER and ARNTL1, in the circadian clock signaling pathway, were found.

Table II. Information of the 15 circadian clock genes.

Table II.

Information of the 15 circadian clock genes.

Gene	Location	Log2FC	FDR	Loss	Gain	Normal	FDR	P-valuea
APCDD1	Chr18: 10,454,628-10,489,948	2.483785	2.96×10−12	154	9	289	1.14×10−87	1.93×10−5
AQP9	Chr15: 58,138,169-58,185,911	2.470799	6.41×10−10	59	2	391	6.65×10−12	2.16×10−3
ARNTL2	Chr12: 27,332,854-27,425,289	2.484382	1.58×10−41	4	44	404	2.46×10−8	1.28×10−7
BIRC7	Chr20: 63,235,883-63,240,507	2.466525	1.09×10−11	0	272	180	3.08×10−87	1.37×10−6
GNGT1	Chr7: 93,591,573-93,911,265	3.164943	4.55×10−10	1	109	342	2.94×10−25	5.17×10−6
HSD11B2	Chr16: 67,430,652-67,437,553	−2.35877	5.09×10−65	2	27	423	1.08×10−4	1.44×10−3
IL23A	Chr12: 56,334,174-56,340,410	3.021143	4.60×10−23	0	35	417	6.42×10−5	4.20×10−3
KRT23	Chr17: 40,922,696-40,937,634	7.179667	2.02×10−34	9	38	405	8.50×10−9	7.79×10−7
NFE2L3	Chr7: 26,152,240-26,187,125	2.676112	1.01×10−85	0	161	291	1.53×10−41	3.89×10−17
NR3C2	Chr4: 148,078,762-148,444,698	−2.63761	4.00×10−84	29	2	421	1.74×10−3	4.86×10−5
ORM2	Chr9: 114,329,869-114,333,252	2.989211	1.39×10−11	9	19	424	1.42×10−3	3.11×10−3
P2RX1	Chr17: 3,896,592-3,916,500	−2.29178	2.75×10−55	138	3	311	1.42×10−36	5.26×10−7
PDX1	Chr13: 27,920,020-27,926,231	4.797965	1.54×10−55	0	200	252	1.37×10−58	8.04×10−18
RNFT2	Chr12: 116,738,178-116,853,631	2.006432	4.31×10−30	3	31	418	1.08×10−4	4.09×10−3
UBE2C	Chr20: 45,812,576-45,816,957	2.15803	5.49×10−43	1	276	175	5.72×10−89	6.51×10−36

a Kolmogorov-Smirnov test. FC, fold-change; FDR, false discovery rate.

Subsequently, a meta-analysis on the expression of the 15 clock genes in CRC using public microarray datasets from the Oncomine database was performed. As presented in Fig. 3, the expression patterns of the clock gene ARNTL2 in 10 independent microarray datasets and TCGA datasets were consistent with previous analyses (32,33). Overexpression was found in all CRC tissues compared with that in the tumor-adjacent tissue (gene median rank, 86.0; P=9.39×10−7).

Figure 3. Oncomine analysis of ARNTL2 mRNA expression levels in the 10 independent microarray datasets and The Cancer Genome Atlas datasets. Data are shown as the median rank of ARNTL2 through each dataset analysis. P-value for ARNTL2 is obtained using the median ranked analysis. ARNTL2, aryl hydrocarbon receptor nuclear translocator-like 2.

Function analysis of the clock gene ARNTL2 driven by CNV in CRC

The expression of the gene ARNTL2 was found in the 452 patients with CRC, among which a total of 48 CNVs occurred, with the presence of copy number gain in 44 patients and copy number loss in 4 patients. ARNTL2 was null in 10 samples, which were consequently removed from the study. The association of ARNTL2 mRNA expression levels with CNV type was identified. As shown in Fig. 4A, single gain and amplification of ARNTL2 were associated with increased mRNA expression, and deletion of ARNTL2 was associated with decreased mRNA expression. Therefore, ARNTL2 gene expression and CNV in CRC tissues show the same trend.

Figure 4. Functional analysis of copy number variation-driven clock gene ARNTL2 in CRC. (A) Box plot of ARNTL2 mRNA expression levels associated with the corresponding gene status. (B) Kaplan-Meier survival analysis of ARNTL2 gene status. (C) Kaplan-Meier survival analysis of ARNTL2 expression levels. (D) ARNTL2 expression associated with the development and progression of CRC. Y-axis represents the expression value following ARNTL2 correction. (E) Paired difference analysis of 41 pairs of paracancerous and cancer samples of ARNTL2. Y-axis indicates the uncorrected expression value of ARNTL2. CRC, colorectal cancer; ARNTL2, aryl hydrocarbon receptor nuclear translocator-like 2; CI, confidence interval; N, node.

A Kaplan-Meier curve analysis was performed to investigate the overall survival time for ARNTL2 in patients with CRC. Compared with that of the patients with normal copy number, the survival rate of the patients with abnormal copy number (gain or loss) of ARNTL2 was significantly decreased (Fig. 4B), whereas the overall survival of patients with CRC with ARNTL2 CNV was significantly decreased. The expression levels of ARNTL2 were also associated with the overall patient survival; higher expression levels indicated greater survival time (Fig. 4C).

To further investigate whether ARNTL2 is involved in the development and progression of CRC, the tumor tissue samples were divided into several subgroups based on pathological TNM (T3+T4 vs. T1+T2, N2+N3 vs. N0+N1, M1 vs. M0) and pathological stages (I–II vs. III–IV) (34). A comparative analysis of ARNTL2 expression profiles was performed. As a result, ARNTL2 expression demonstrated a statistically significant association with pathological stages (P<0.001) and pathological N (P<0.001) (Fig. 4D).

To gain a clearer understanding of the expression of ARNTL2 in patients with cancer and adjacent tissues, a paired difference analysis of ARNTL2 from 41 patients with cancer and adjacent tissues was performed. The expression of ARNTL2 in cancer tissues was significantly higher compared with that in adjacent tissues (P=1.058×108; Fig. 4E). This is consistent with the results of our previous analysis of the difference.

To investigate the biological characteristics shared by the different ARNTL2 expression levels, a GSEA was performed. The most significant pathways for the upregulated gene sets in the significance order (nominal P<0.05) are shown in Fig. 5. The six pathways, including ‘natural killer cell-mediated cytotoxicity’, ‘oocyte meiosis’, the ‘p53 signaling pathway’, ‘pancreatic cancer’, ‘prostate cancer’ and the ‘toll like receptor signaling pathway’ were significant in the ARNTL2 high expression phenotype. Among these pathways, some were directly linked to cancer pathogenesis, such as ‘pancreatic cancer’, the ‘p53 signaling pathway’ (35) and ‘prostate cancer’. There were no significant pathways for downregulated gene sets with nominal P<0.05.

Figure 5. Gene Set Enrichment Analysis results showing the association of aryl hydrocarbon receptor nuclear translocator-like 2 expression levels and colorectal cancer-related gene sets.

Discussion

CRC is the third most commonly occurring cancer worldwide and the fourth most frequent cause of death having an oncological origin (1); it is considered to be a complex disease resulting from a combination of environmental factors, genetic/epigenetic predisposing variants and specific molecular mechanisms. Chromosomal instability (CIN) has been defined as a major factor contributing to CRC carcinogenesis (36). CNV exists as a genetic polymorphism in the human genome that is a type of CIN (37). The form of CNV directly affecting the expression of a gene is mainly the deletion or amplification of a copy number of a gene, causing an increase or decrease in the amount of gene expression and increasing the occurrence of the disease (38). A previous study found that tumor necrosis factor receptor superfamily member 10C CNV is associated with metastatic colorectal cancer (39). In the present study, an integrated analysis of CNV data and gene expression profile for CRC with a large sample size (n=503, including 462 patient samples and 41 tumor-adjacent tissue samples) was performed. A total of 10,256 genes with significantly different CNV and 2,083 aberrantly expressed mRNAs were obtained, of which 295 genes showed a statistically significant association between the gene expression and CNV; therefore, these 295 genes were regarded as CRC-related CNV-driven genes. The present findings may provide a new theoretical basis for the pathogenesis of CRC and also contribute to the development of new therapeutic strategies.

In the present study, CNV-driven genes were only enriched in the Wnt signaling pathway. The Wnt pathway is involved in the regulation of important physiological processes such as normal embryo development, and cell proliferation and differentiation, and its abnormal activation plays an important role in the process of tumor development, metastasis and therapeutic resistance (40). A previous study showed that >90% of colorectal cancer cases have abnormal activation of the Wnt classical signaling pathway (41). Meanwhile, some studies suggested that the regulation of circadian clock genes, such as CRY1 (42) and Rev-erbα (43), was mediated by the classical Wnt/β-catenin signaling pathway. Further study of the Wnt signaling pathway will help to develop new strategies for CRC treatment. The present findings provide a new clue to study this signaling pathway.

Studies in circadian clock genes may expand the knowledge regarding the mechanism of occurrence and development of tumors, and may provide a new approach for tumor therapy (44). Indeed, multiple epidemiological studies have shown that impaired function of the circadian clock promotes development of cancer (45). For circadian clock genes, including Per1, Per2, and Per3, the expression levels of which are often found to be decreased in pancreatic cancer (46) and gastric cancer (47), as well as the disruption of autonomic rhythm. Additionally, in a previous study, CRC showed lower expression of NPAS2 compared with that in healthy tissues, and this was negatively associated with tumor size, stage and metastasis (48). A previous study has also shown that varying degrees of biorhythm destruction are found in 50% of metastatic cancer cases (49). In the present study, 15 CNV-driven circadian clock genes in CRC tissues were identified, indicating that these circadian clock genes may play a role in cancer. However, this requires further validation at the protein level. ARNTL2 has been described as a candidate biomarker in different cancer types, including kidney cancer (50), colorectal cancer (21) and hepatocellular carcinoma (51), and similar results for ARNTL2 were obtained in the present study. As research continues to deepen, numerous studies have found that genomic alterations involving circadian clock genes, such as point mutations or CNV, are frequently found in different human cancer types. The rs1801260 SNP, in the 3′ untranslated region of the clock circadian regulator gene, was found to be associated with the development of CRC (52). The CNV form of the BMAL1 gene has also been found in multiple cancers, such as breast and colorectal cancer (35 gains and 7 losses in CNV numbers) (53). Previous studies have observed a close association between ARNTL2 expression and various types of cancer (21,22); however, no studies have characterized the association between CNV in ARNTL2 and cancer. In the present study, upregulation of ARNTL2 in patients with CRC was found, and ARNTL2 CNV has three forms: Single loss, single gain, and amplification. Further analysis found that the expression of ARNTL2 has the same trend as CNV. Our study showed that the expression level of ARNTL2 was abnormal due to the presence of CNV, which promoted the occurrence and development of CRC.

Genetic drifts in ARNTL2 polymorphisms have been described in the human population leading to variation in the circadian rhythm regulation (54). The expression of ARNTL2 was significantly associated with survival time in patients with CRC. A previous study found that high ARNTL2 expression predicted poor survival in patients with lung adenocarcinoma (55). However, in the present study, low ARNTL2 expression predicted poor survival in patients with CRC. This may be due to the heterogeneity between different types of cancer. Moreover, a previous report indicated that ARNTL2 high levels significantly influence mammary tumor metastasis (23). ARNTL2 levels were also significantly associated with pathological stage and N stage in patients with CRC in the present study. ARNTL2 CNV was also significantly associated with survival time in patients with CRC. These data suggest that ARNTL2 can be used as a prognostic factor for CRC, which may bring more personalized treatment to patients with CRC. GSEA analysis showed that ARNTL2 is enriched for gene sets associated with CRC pathogenesis, such as the ‘p53 signaling pathway’. These findings suggest that the CNV-driven clock gene ARNTL2 plays a crucial role in the development and progression of CRC. However, this study has some limitations as it was an in silico study. Further in vivo investigations would be beneficial to fully understand the roles of ARNTL1 in CRC initiation and development.

In summary, to the best of our knowledge, the present study demonstrates for the first time that circadian clock genes play an important role in CRC in the form of CNV, and that 15 CNV-driven clock genes are associated with the etiology and pathogenesis of CRC. Finally, it was concluded that CNV in the circadian clock gene ARNTL2 may be a useful genetic biomarker for the treatment of individualized CRC patients and may identify patients who may benefit from more aggressive systemic treatment strategies.

Acknowledgements

Not applicable.

Funding

This study was partly supported by the National Natural Science Foundation of China (grant no. 61775139 and grant no. 61373057), the Shanghai Natural Science Foundation (grant no. 15ZR1420800) and the Key Projects of Humanities and Social Sciences Research in Colleges and Universities in Anhui Province (grant no. SK2018A0630).

Availability of data and materials

The datasets used and/or analyzed during the present study are available from the corresponding author upon reasonable request.

Authors' contributions

WLY, SHP and LHJ conceived the study. WLY collected the data and performed the bioinformatics analyses. LL, XBL, CW, DS, CXD, and SCL performed quality control of the raw data and performed data analyses. WLY and SHP wrote the manuscript. SHP and LHJ supervised the study and agreed to be responsible for ensuring that all aspects of the study are accurate and have been appropriately investigated. All authors read and approved the final manuscript.

Ethics approval and consent to participate

Not applicable.

Patient consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

References