BDP1 Variants I1264M and V1347M Significantly Associated with Clinical Outcomes of Pediatric Neuroblastoma Patients Imply a New Prognostic Biomarker: A 121-Patient Cancer Genome Study

Background: Neuroblastoma (N.B.) is the most common tumor in children. The gene BDP1 (B Double Prime 1) plays a role in cancers but is less known in N.B. Thus, we conducted this study to investigate the value of BDP1 mutations in N.B. prognosis. Methods: A dataset of 121 NB patients from the Cancer Genome Atlas database was used to analyze BDP1 gene mutations by RNA sequencing. Kaplan-Meier estimates were performed for overall survival (O.S.) analysis on BDP1 variants, and Cox’s proportional hazards regression model was used for multivariate analysis. Results: In 121 NB patients, we identified two variants of BDP1 associated with N.B., located at chr5:71511131 and chr5:71510884. The prevalence of these BDP1 variants, I1264M and V1347M, was 52.9% (64/121) and 45.5% (55/121), respectively. O.S. analysis showed a significant difference between subgroups with or without BDP1 variants (p < 0.05). Multivariate analysis further revealed that BDP1ariants were independent prognostic variables in N.B. (p < 0.05). Conclusion: Our results suggest BDP1 variants are associated with significantly improved clinical outcomes in N.B., thus providing clinicians with a new tool.


Introduction
As the most common extracranial pediatric solid tumor, ever since defined as an embryonal tumor of the autonomic nervous system, 90% of neuroblastoma (N.B.) cases are diagnosed before the patient is 5-year old, with a median age at diagnosis of around 18 months [1]. Such a delayed diagnosis results in patients with high-risk N.B. who come below 50% for 5-year survival rates [2], especially in older children, with poor outcomes despite aggressive multimodal therapy [2]. This situation sparks the International Neuroblastoma Risk Group (INRG) to standardize diagnostic systems using age, histology,

Dataset
A dataset representing 121 NB patients was extracted from The Cancer Genome Atlas (TCGA) database (https://www.genome.gov/Funded-Programs-Projects/Cancer-Genome-Atlas, assessed on 6 October 2020) -The Cancer Genome Atlas (TCGA) was a joint effort of the United States National Cancer Institute (NCI) and the US National Human Genome Research Institute (NHGRI). The Cancer Genome Atlas (TCGA), a landmark cancer genomics program, molecularly characterized over 20,000 primary cancer and matched normal samples spanning 33 cancer types (https://www.cancer.gov/about-nci/ organization/ccg/research/structural-genomics/tcga, accessed on 6 October 2020).These patients all completed a long-term clinical observation period from three months to ten years past diagnosis. We downloaded RNA-seq data on these tumor samples to conduct a risk analysis of the death of N.B. patients.

Variant Calling for Detecting BDP1 Variants
We followed the standard variant calling process of GATK ( [20] (GATK website www.broadinstitute.org/gatk, accessed on 6 October 2020, for each patient's RNA-seq data, followed by variant annotation and detection of unique genes/variants in N.B. First, according to the variant calling of GATK, we used Picard to make a duplicate and then used GATK to split reads into exon samples. After that, we ran base recalibration according to the common SNPs and indels. Then, we used the HaplotypeCaller tool in GATK to call variants and filter variants with criteria of D.P. > 10, Q.D. < 2, F.S. > 30, and no SNPcluster. The filters automatically applied by HaplotypeCaller built-in tools of evaluation metrics came with the abbreviations defined as the depth of informative coverage for each sample (D.P.), false-positive (F.P.), FisherStrand (F.S.), QualByDepth (Q.D.), true positive (T.P.). We used AnnoVar [21] for BDP1 variant annotations and only analyzed exonic variants. To study the BDP1 gene, we generated gene-sample and variant-sample matrices for the 121 samples. Each row represented a sample in the gene-sample matrix. Each column represented a gene from the set (the set of genes of filtered exonic variants). If BDP1 had exonic variants in a sample, the corresponding cell was assigned the value 1 (or otherwise 0). Each row represented a sample for the variant-sample matrix, and each column represented a gene variant. Similar to the gene-sample matrix, if a BDP1 variant occurred in a sample, the corresponding cell was assigned the value 1. Otherwise, it was assigned 0.

Statistical Analysis
Overall survival (O.S.) was defined as the time from diagnosis to disease-related death or the date of the last follow-up. For O.S. calculation, patients alive at the end of the study period died from other causes or lost to follow-up were excluded in the Kaplan-Meier analyses and compared using the log-rank test. Cox's proportional hazards regression model was used for multivariate analysis.

BDP1 Variants Found in 121 NB Patients
To tie in the clinical relevance, we investigated 121 NB patients to search for the prevalence rate of BDP1 variants, I1264M and V1347M, showing 52.9% (64/121) and 45.5% (55/121), respectively. Notably, the rate of non-BDP1 variants was 44.6% (54/121), while the percentage of patients harboring both was 43.0% (52/121). We followed up with characteristics for correlation of these variants with the clinical outcomes.

Clinical Outcomes of N.B. Patients Carrying BDP1 Variants
We compared the clinical outcome in subgroups with specific variant patterns, biological markers, and clinical characteristics, including MYCN-amplification, age, gender, stage, Children's Oncology Group (COG) risk score, and grade in this large, independent population (Table 1). Specifically, we found two variants in BDP1, I1264M and V1347M, significantly associated with improved clinical outcomes ( Figure 1). The median O.S. of the patients without I1264M was 35.8 months, and those without V1347M, considerably shorter than the median O.S. of those patients with variant I1264M (87.9 months, p = 0.0007) or patients with V1347M (87.9 months, p = 0.0019).    Table 1 suggests that I1264M alone has a dynamic prognostic value. Moreover, the median O.S. of the subset with non-BDP1 variants was 33.8 months, much shorter than the median O.S. of the subset harboring both BDP1 variants (87.9 months, p = 0.0011).

BDP1 Variants May Serve for Neuroblastoma Prognosis
We conducted a multivariate Cox hazard analysis of risk factors for N.B. prognosis ( Table 2). As expected, stage, ploidy, and age at diagnosis were independent variables significantly associated with O.S. time; BDP1 variants V1347M and I1264M reduced patients' risk of death. Noteworthy, BDP1 variants V1347M and I1264M reduced patients' risk of death, and the hazard ratios (H.R.) for these variants were 0.488 (95% confidence interval (CI): 0.289-0.826, p = 0.007) and 0.58 (95% CI: 0.349-0.964, p = 0.035), respectively.  Table 1 suggests that I1264M alone has a dynamic prognostic value. Moreover, the median O.S. of the subset with non-BDP1 variants was 33.8 months, much shorter than the median O.S. of the subset harboring both BDP1 variants (87.9 months, p = 0.0011).

BDP1 Variants May Serve for Neuroblastoma Prognosis
We conducted a multivariate Cox hazard analysis of risk factors for N.B. prognosis ( Table 2). As expected, stage, ploidy, and age at diagnosis were independent variables significantly associated with O.S. time; BDP1 variants V1347M and I1264M reduced patients' risk of death. Noteworthy, BDP1 variants V1347M and I1264M reduced patients' risk of

Discussion
We provide valuable clinical data on a pediatric brain N.B. Specifically, we found that BDP1 variants, I1264M, V1347M, are associated with clinical outcomes, as measured by Kaplan-Meier overall survival (O.S.) curves, extrapolated from 121 patients of the TCGA data. The high-level amplification of MYCN on chromosome 2p24 occurs in approximately 20% of N.B. patients, indicating its link to aggressive disease with a poor prognosis. Inhibitors that downregulate MYCN/MYC proteins can suppress N.B. tumor growth [22]. Mutations in the tyrosine kinase domain of Anaplastic lymphoma kinase (ALK) have also been identified in N.B. FDA approved a series of ALK tyrosine kinase inhibitors (TKIs) for use in ALK-driven cancers [23]. Furthermore, the paired-like homeobox 2B gene (PHOX2B) was discovered as the first N.B. predisposition gene, accounting for 6.4% of heritable cases [24]. Despite the development of molecular strategies, the diagnostic and prognostic value of many gene mutations in N.B. patients is not very clear. Here, we first report the association of BDP1 mutations with N.B. prognosis, showing that two BDP1 variants, I1264M and V1347M, are located at chr5:71510884 and chr5:71511131, are significantly correlated with better outcomes in an independent N.B. population. We further verified the prevalence of these BDP1 variants, I1264M and V1347M, were 52.9% (64/121) and 45.5% (55/121), respectively. Interestingly, the rate of patients harboring both was 43.0% (52/121), suggesting a high rate of overlap between them. Consistently, the minor allele frequency (MAF) of BDP1 polymorphism variants in an NCBI database were detected by DNA-seq ranging from 17% to 22%, which was similar to our results using RNA-seq.
The weakness of the current approach was to make a comparison by evaluating the survival on the ground of gene expression on a single gene, which might be uninformative, given the fact that the variables that affect the clinical course of these patients are many and deserve an evaluation on a case history with clinical data collected (retrospective or prospective). However, the lack of such case history with clinical data collected (retrospective or prospective) is evident in the N.B.-databases. The solution to this problem might fall on conducting the work on registries in the future, which may be useful on gene expression studies but not on the outcome of the study only. Artificial intelligence-based medicine such as machine learning [25] might integrate case history, clinical information, molecular pathology, whole-genome sequencing, RNA seq., lifestyle courses, treatment-driven pathophysiological, and biomarker changes into comprehensive lifetime management. In previous reports, mutations in genes such as ALK or ATRX are associated with poor prognosis in N.B. [5]. In contrast, others, such as TAM1 variants, were associated with improved clinical outcomes [22]. Gene polymorphisms also confer N.B. susceptibility resulting in different clinical outcomes such as METTL14 [7] and ERCC1/XPF [8]. In our study, age at diagnosis (≥18 months), diploidy, and stage IV were all associated with poor prognosis, as expected. Interestingly, both BDP1 variants significantly improved overall survival time and reduced patients' risk of death from N.B. Therefore, our findings support that the protective phenotype of BDP1 variants may exist in subsets of N.B. patients with the following implications.
First, establishing the association of BDP1 variants with clinical outcomes warrants in vitro or in vivo experimental verifications. Second, new bioinformatic methodologies may be used to explore the networks of these variants with accuracy. Third, both variants may affect BPD1 mRNA expression levels in subsets of N.B. patients [26]. Fourth, NB heterogeneity surfaces as single-cell characterization of malignant phenotypes and developmental trajectories of adrenal N.B. [27] for subclonal evolution [28]. Lastly, we must expand to explore the mechanisms by which these two variants affect the functions of NB-involved BDP1 (Figure 2). Thus, the value of BDP1 mutations in N.B. deserves further investigation. many and deserve an evaluation on a case history with clinical data collected (retrospective or prospective). However, the lack of such case history with clinical data collected (retrospective or prospective) is evident in the N.B.-databases. The solution to this problem might fall on conducting the work on registries in the future, which may be useful on gene expression studies but not on the outcome of the study only. Artificial intelligencebased medicine such as machine learning [25] might integrate case history, clinical information, molecular pathology, whole-genome sequencing, RNA seq., lifestyle courses, treatment-driven pathophysiological, and biomarker changes into comprehensive lifetime management. In previous reports, mutations in genes such as ALK or ATRX are associated with poor prognosis in N.B. [5]. In contrast, others, such as TAM1 variants, were associated with improved clinical outcomes [22]. Gene polymorphisms also confer N.B. susceptibility resulting in different clinical outcomes such as METTL14 [7] and ERCC1/XPF [8]. In our study, age at diagnosis (≥18 months), diploidy, and stage IV were all associated with poor prognosis, as expected. Interestingly, both BDP1 variants significantly improved overall survival time and reduced patients' risk of death from N.B. Therefore, our findings support that the protective phenotype of BDP1 variants may exist in subsets of N.B. patients with the following implications.
First, establishing the association of BDP1 variants with clinical outcomes warrants in vitro or in vivo experimental verifications. Second, new bioinformatic methodologies may be used to explore the networks of these variants with accuracy. Third, both variants may affect BPD1 mRNA expression levels in subsets of N.B. patients [26]. Fourth, NB heterogeneity surfaces as single-cell characterization of malignant phenotypes and developmental trajectories of adrenal N.B. [27] for subclonal evolution [28]. Lastly, we must expand to explore the mechanisms by which these two variants affect the functions of NBinvolved BDP1 (Figure 2). Thus, the value of BDP1 mutations in N.B. deserves further investigation. Figure 2. Role of BDP1 in the TFIIIB complex of transcription machinery in cell proliferation, cell transformation, and cancer progression. TFIIIB is a complex of the transcription machinery of RNA Pol III gene transcription [12]. TFIIIB consists of Brf1, BDP1, and TBP. Oncogenic proteins, such as c-Jun, c-Myc, Ras, and c-Fos activate TFIIIB to enhance RNA Pol III gene transcription, resulting in cancer tumorigenesis, while tumor suppressors, such as BRCA1 [13], PTEN, p53, and pRb, Figure 2. Role of BDP1 in the TFIIIB complex of transcription machinery in cell proliferation, cell transformation, and cancer progression. TFIIIB is a complex of the transcription machinery of RNA Pol III gene transcription [12]. TFIIIB consists of Brf1, BDP1, and TBP. Oncogenic proteins, such as c-Jun, c-Myc, Ras, and c-Fos activate TFIIIB to enhance RNA Pol III gene transcription, resulting in cancer tumorigenesis, while tumor suppressors, such as BRCA1 [13], PTEN, p53, and pRb, inactivate its activity to decrease the transcription of Pol III genes [14], leading to repression of tumor development [BDP1, a.k.a., TATA Box Binding Protein (TBP)-Associated Factor] [16]. TFIIIB). The concept of BDP1's cancer-involved is emerging with colorectal cancers [16,17], lung cancer [18], and breast cancer [19].
Author Contributions: X.J., Z.W. and J.C. conceived of the study. X.L., L.S. and S.C.L. drafted the manuscript. A.S., L.T. and X.C. performed the data analysis. All other authors participated in its experimentation and revision of the manuscript. All authors have read and agreed to the published version of the manuscript.
Funding: This study was supported by the Natural Science Foundation of Chongqing, China (cstc2020jcyj-msxmX1063), graduate tutor team construction project of Chongqing Municipal Education Commission Foundation, China (dstd201801), and National Natural Science Foundation of China (81703063). This work was also supported in part by the CHOC Children's-UC Irvine Child Health Research Awards #16004004; CHOC-UCI Child Health Research grant #16004003; and CHOC CSO grant #16986004.

Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable. Patient consent was waived for the NIH public databases (The Cancer Genome Atlas).

Data Availability Statement:
Publicly available datasets were analyzed in this study. This data can be found here: https://portal.gdc.cancer.gov.